Back to home page

MITgcm
 
 

    

File indexing completed on 2019-12-22 06:11:14 UTC

view on githubraw file Latest commit 5511e94c on 2019-12-11 04:48:28 UTC
dcaf8922f9 Jean*0001 #include "MOM_COMMON_OPTIONS.h"
b83c015289 Jean*0002 
                0003 CBOP
                0004 C     !ROUTINE: MOM_V_IMPLICIT_R
                0005 C     !INTERFACE:
dcaf8922f9 Jean*0006       SUBROUTINE MOM_V_IMPLICIT_R(
b83c015289 Jean*0007      I                 kappaRV,
                0008      I                 bi, bj, myTime, myIter, myThid )
                0009 C     !DESCRIPTION: \bv
                0010 C     *==========================================================*
                0011 C     | S/R MOM_V_IMPLICIT_R
                0012 C     | o Solve implicitly vertical advection & diffusion
                0013 C     |   of momentum, meridional component
                0014 C     *==========================================================*
                0015 C     *==========================================================*
                0016 C     \ev
                0017 
                0018 C     !USES:
                0019       IMPLICIT NONE
                0020 C     == Global data ==
                0021 #include "SIZE.h"
                0022 #include "EEPARAMS.h"
                0023 #include "PARAMS.h"
                0024 #include "GRID.h"
                0025 #include "DYNVARS.h"
79074ef66b Jean*0026 #include "FFIELDS.h"
b83c015289 Jean*0027 
                0028 C     !INPUT/OUTPUT PARAMETERS:
                0029 C     == Routine Arguments ==
c44d420e3b Jean*0030       _RL kappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
b83c015289 Jean*0031       INTEGER bi, bj
                0032       _RL myTime
                0033       INTEGER myIter, myThid
                0034 
                0035 C     !LOCAL VARIABLES:
                0036 C     == Local variables ==
e60362cf6f Jean*0037 C iMin,iMax,jMin,jMax :: computational domain
                0038 C i,j,k     :: loop indices
                0039 C a5d       :: 2nd  lower diagonal of the pentadiagonal matrix
                0040 C b5d       :: 1rst lower diagonal of the pentadiagonal matrix
                0041 C c5d       :: main diagonal       of the pentadiagonal matrix
                0042 C d5d       :: 1rst upper diagonal of the pentadiagonal matrix
                0043 C e5d       :: 2nd  upper diagonal of the pentadiagonal matrix
                0044 C rTrans    :: vertical volume transport at interface k
79074ef66b Jean*0045 C cDrag     :: bottom drag coefficient
                0046 C KE2d      :: horizontal Kinetic Energy
e60362cf6f Jean*0047 C diagonalNumber :: number of non-zero diagonals in the matrix
                0048 C errCode   :: > 0 if singular matrix
b83c015289 Jean*0049       INTEGER iMin,iMax,jMin,jMax
e60362cf6f Jean*0050       PARAMETER( iMin = 1, iMax = sNx )
                0051       PARAMETER( jMin = 1, jMax = sNy+1 )
b83c015289 Jean*0052       INTEGER i,j,k
                0053       INTEGER diagonalNumber, errCode
8a58850ca8 Jean*0054 c     _RL a5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0055       _RL b5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0056       _RL c5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0057       _RL d5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0058 c     _RL e5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0059       _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79074ef66b Jean*0060       _RL cDrag (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0061       _RL KE2d  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
6d4c91318f Jean*0062       _RL rCenter, rUpwind, upwindFac
1c99f96b44 Jean*0063 #ifdef ALLOW_DIAGNOSTICS
                0064       CHARACTER*8 diagName
                0065       LOGICAL     DIAGNOSTICS_IS_ON
                0066       EXTERNAL    DIAGNOSTICS_IS_ON
79074ef66b Jean*0067       _RL recip_dT
                0068       _RL delVdt(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0069       _RL vf    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1c99f96b44 Jean*0070 #endif
b83c015289 Jean*0071 CEOP
                0072 
                0073 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0074 C     Solve for V-component :
                0075 C----------------------------
                0076 
                0077 C--   Initialise
79074ef66b Jean*0078       diagonalNumber = 1
b83c015289 Jean*0079       DO k=1,Nr
8a58850ca8 Jean*0080        DO j=1-OLy,sNy+OLy
                0081         DO i=1-OLx,sNx+OLx
b83c015289 Jean*0082 c        a5d(i,j,k) = 0. _d 0
                0083          b5d(i,j,k) = 0. _d 0
                0084          c5d(i,j,k) = 1. _d 0
                0085          d5d(i,j,k) = 0. _d 0
                0086 c        e5d(i,j,k) = 0. _d 0
                0087         ENDDO
                0088        ENDDO
                0089       ENDDO
79074ef66b Jean*0090       DO j=1-OLy,sNy+OLy
                0091        DO i=1-OLx,sNx+OLx
                0092          KE2d (i,j) = 0. _d 0
                0093          cDrag(i,j) = 0. _d 0
                0094        ENDDO
                0095       ENDDO
                0096 #ifdef ALLOW_DIAGNOSTICS
                0097       IF ( useDiagnostics .AND. ( implicitViscosity .OR.
                0098      &                            selectImplicitDrag.GE.1 ) ) THEN
                0099        DO k=1,Nr
                0100         DO j=1-OLy,sNy+OLy
                0101          DO i=1-OLx,sNx+OLx
                0102            delVdt(i,j,k) = gV(i,j,k,bi,bj)
                0103          ENDDO
                0104         ENDDO
                0105        ENDDO
                0106       ENDIF
                0107 #endif /* ALLOW_DIAGNOSTICS */
b83c015289 Jean*0108 
362dc18408 Jean*0109       IF ( implicitViscosity .AND. Nr.GT.1 ) THEN
                0110 
b83c015289 Jean*0111 C--   set the tri-diagonal matrix to solve the implicit viscosity
                0112        diagonalNumber = 3
                0113 C-     1rst lower diagonal :
                0114        DO k=2,Nr
                0115         DO j=jMin,jMax
                0116          DO i=iMin,iMax
e60362cf6f Jean*0117           IF (maskS(i,j,k-1,bi,bj).EQ.oneRS)
8a58850ca8 Jean*0118      &     b5d(i,j,k) = -deltaTMom
616600b8d2 Patr*0119      &                  *_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
e60362cf6f Jean*0120      &                  *recip_deepFac2C(k)*recip_rhoFacC(k)
b83c015289 Jean*0121      &                  *kappaRV(i,j, k )*recip_drC( k )
e60362cf6f Jean*0122      &                  *deepFac2F( k )*rhoFacF( k )
b83c015289 Jean*0123          ENDDO
                0124         ENDDO
                0125        ENDDO
                0126 C-     1rst upper diagonal :
                0127        DO k=1,Nr-1
                0128         DO j=jMin,jMax
                0129          DO i=iMin,iMax
e60362cf6f Jean*0130           IF (maskS(i,j,k+1,bi,bj).EQ.oneRS)
8a58850ca8 Jean*0131      &     d5d(i,j,k) = -deltaTMom
e60362cf6f Jean*0132      &                  *_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
                0133      &                  *recip_deepFac2C(k)*recip_rhoFacC(k)
                0134      &                  *kappaRV(i,j,k+1)*recip_drC(k+1)
                0135      &                  *deepFac2F(k+1)*rhoFacF(k+1)
b83c015289 Jean*0136          ENDDO
                0137         ENDDO
                0138        ENDDO
                0139 C-     Main diagonal :
                0140        DO k=1,Nr
                0141         DO j=jMin,jMax
                0142          DO i=iMin,iMax
e60362cf6f Jean*0143            c5d(i,j,k) = 1. _d 0 - ( b5d(i,j,k) + d5d(i,j,k) )
b83c015289 Jean*0144          ENDDO
                0145         ENDDO
                0146        ENDDO
                0147 
                0148 C--   end if implicitDiffusion
                0149       ENDIF
                0150 
b76b90dc9b Jean*0151       IF ( selectImplicitDrag.GE.1 ) THEN
5511e94c2d Jean*0152 C--   Add bottom and top drag coeff to main diagonal "c5d":
                0153 
b76b90dc9b Jean*0154        IF ( no_slip_bottom .OR. selectBotDragQuadr.GE.0
                0155      &                     .OR. bottomDragLinear.NE.0. ) THEN
79074ef66b Jean*0156         DO k=1,Nr
a125ab7be7 Jean*0157           CALL MOM_V_BOTDRAG_COEFF( bi, bj, k, .FALSE.,
79074ef66b Jean*0158      I                     uVel(1-OLx,1-OLy,k,bi,bj),
                0159      I                     vVel(1-OLx,1-OLy,k,bi,bj),
b76b90dc9b Jean*0160      I                     kappaRV,
79074ef66b Jean*0161      U                     KE2d,
                0162      O                     cDrag,
                0163      I                     myIter, myThid )
                0164           DO j=jMin,jMax
                0165            DO i=iMin,iMax
                0166              c5d(i,j,k) = c5d(i,j,k)
                0167      &                  + cDrag(i,j)*deltaTMom
                0168      &                  *_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
                0169            ENDDO
                0170           ENDDO
                0171 C--   Save (in 2-D array) drag-coeff for diagnostics of bottom-stress
                0172           IF ( useDiagnostics ) THEN
                0173            DO j=1-OLy,sNy+OLy
                0174             DO i=1-OLx+1,sNx+OLx
                0175               botDragV(i,j,bi,bj) = botDragV(i,j,bi,bj)
                0176      &                            + cDrag(i,j)*rUnit2mass
                0177             ENDDO
                0178            ENDDO
                0179           ENDIF
                0180 C-   End k loop
                0181         ENDDO
b76b90dc9b Jean*0182        ENDIF
5511e94c2d Jean*0183 
                0184 #ifdef ALLOW_SHELFICE
                0185        IF ( useShelfIce ) THEN
                0186 C--   Add shelfice drag coeff to main diagonal "c5d":
                0187         DO k=1,Nr
                0188           CALL SHELFICE_V_DRAG_COEFF( bi, bj, k, .FALSE.,
                0189      I                     uVel(1-OLx,1-OLy,k,bi,bj),
                0190      I                     vVel(1-OLx,1-OLy,k,bi,bj),
                0191      I                     kappaRV,
                0192      U                     KE2d,
                0193      O                     cDrag,
                0194      I                     myIter, myThid )
                0195           DO j=jMin,jMax
                0196            DO i=iMin,iMax
                0197              c5d(i,j,k) = c5d(i,j,k)
                0198      &                  + cDrag(i,j)*deltaTMom
                0199      &                  *_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
                0200            ENDDO
                0201           ENDDO
                0202 C-   End k loop
                0203         ENDDO
                0204        ENDIF
                0205 #endif /* ALLOW_SHELFICE */
                0206 
                0207 C--   end if selectImplicitDrag > 0
b76b90dc9b Jean*0208       ENDIF
                0209 
362dc18408 Jean*0210       IF ( momImplVertAdv .AND. Nr.GT.1 ) THEN
b83c015289 Jean*0211 
                0212         diagonalNumber = 3
                0213         DO k=2,Nr
                0214 
                0215           DO j=jMin,jMax
                0216            DO i=iMin,iMax
                0217             rTrans(i,j) = 0.5 _d 0 * (
                0218      &                wVel(i, j ,k,bi,bj)*rA(i, j ,bi,bj)
                0219      &                                 *maskC(i, j ,k-1,bi,bj)
                0220      &              + wVel(i,j-1,k,bi,bj)*rA(i,j-1,bi,bj)
                0221      &                                 *maskC(i,j-1,k-1,bi,bj)
e60362cf6f Jean*0222      &                               )*deepFac2F(k)*rhoFacF(k)
b83c015289 Jean*0223            ENDDO
                0224           ENDDO
                0225 
                0226           IF ( vectorInvariantMomentum ) THEN
6d4c91318f Jean*0227 C-          space Centered/Upwind advection scheme, Advective form:
                0228             IF ( upwindShear ) THEN
                0229               upwindFac = 1. _d 0
                0230             ELSE
                0231               upwindFac = 0. _d 0
                0232             ENDIF
b83c015289 Jean*0233             DO j=jMin,jMax
                0234              DO i=iMin,iMax
8a58850ca8 Jean*0235               rCenter = 0.5 _d 0 *deltaTMom*rTrans(i,j)
6d4c91318f Jean*0236      &                           *recip_rAs(i,j,bi,bj)*rkSign
                0237               rUpwind = ABS(rCenter)*upwindFac
b83c015289 Jean*0238               b5d(i,j,k) = b5d(i,j,k)
6d4c91318f Jean*0239      &                   - (rCenter+rUpwind)
616600b8d2 Patr*0240      &                     *_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
e60362cf6f Jean*0241      &                     *recip_deepFac2C(k)*recip_rhoFacC(k)
b83c015289 Jean*0242               c5d(i,j,k) = c5d(i,j,k)
6d4c91318f Jean*0243      &                   + (rCenter+rUpwind)
616600b8d2 Patr*0244      &                     *_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
e60362cf6f Jean*0245      &                     *recip_deepFac2C(k)*recip_rhoFacC(k)
b83c015289 Jean*0246               c5d(i,j,k-1) = c5d(i,j,k-1)
6d4c91318f Jean*0247      &                   - (rCenter-rUpwind)
616600b8d2 Patr*0248      &                     *_recip_hFacS(i,j,k-1,bi,bj)*recip_drF(k-1)
e60362cf6f Jean*0249      &                     *recip_deepFac2C(k-1)*recip_rhoFacC(k-1)
b83c015289 Jean*0250               d5d(i,j,k-1) = d5d(i,j,k-1)
6d4c91318f Jean*0251      &                   + (rCenter-rUpwind)
616600b8d2 Patr*0252      &                     *_recip_hFacS(i,j,k-1,bi,bj)*recip_drF(k-1)
e60362cf6f Jean*0253      &                     *recip_deepFac2C(k-1)*recip_rhoFacC(k-1)
b83c015289 Jean*0254              ENDDO
                0255             ENDDO
                0256           ELSE
                0257 C-          space Centered advection scheme, Flux form:
                0258             DO j=jMin,jMax
                0259              DO i=iMin,iMax
8a58850ca8 Jean*0260               rCenter = 0.5 _d 0 *deltaTMom*rTrans(i,j)
6d4c91318f Jean*0261      &                           *recip_rAs(i,j,bi,bj)*rkSign
b83c015289 Jean*0262               b5d(i,j,k) = b5d(i,j,k)
616600b8d2 Patr*0263      &            - rCenter*_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
e60362cf6f Jean*0264      &                     *recip_deepFac2C(k)*recip_rhoFacC(k)
b83c015289 Jean*0265               c5d(i,j,k) = c5d(i,j,k)
616600b8d2 Patr*0266      &            - rCenter*_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
e60362cf6f Jean*0267      &                     *recip_deepFac2C(k)*recip_rhoFacC(k)
b83c015289 Jean*0268               c5d(i,j,k-1) = c5d(i,j,k-1)
616600b8d2 Patr*0269      &            + rCenter*_recip_hFacS(i,j,k-1,bi,bj)*recip_drF(k-1)
e60362cf6f Jean*0270      &                     *recip_deepFac2C(k-1)*recip_rhoFacC(k-1)
b83c015289 Jean*0271               d5d(i,j,k-1) = d5d(i,j,k-1)
616600b8d2 Patr*0272      &            + rCenter*_recip_hFacS(i,j,k-1,bi,bj)*recip_drF(k-1)
e60362cf6f Jean*0273      &                     *recip_deepFac2C(k-1)*recip_rhoFacC(k-1)
b83c015289 Jean*0274              ENDDO
                0275             ENDDO
dcaf8922f9 Jean*0276             STOP 'MOM_IMPLICIT_R: Flux Form not yet finished.'
b83c015289 Jean*0277           ENDIF
                0278 
                0279 C--     end k loop
                0280         ENDDO
                0281 
                0282 C--   end if momImplVertAdv
                0283       ENDIF
                0284 
b76b90dc9b Jean*0285       IF ( diagonalNumber .EQ. 1 ) THEN
                0286         errCode = 0
                0287         DO k=1,Nr
                0288          DO j=jMin,jMax
                0289           DO i=iMin,iMax
                0290             IF ( c5d(i,j,k).NE.zeroRL ) THEN
                0291              c5d(i,j,k) = 1. _d 0 / c5d(i,j,k)
                0292             ELSE
                0293              c5d(i,j,k) = 0. _d 0
                0294              errCode = 1
                0295             ENDIF
                0296           ENDDO
                0297          ENDDO
                0298          DO j=jMin,jMax
                0299           DO i=iMin,iMax
                0300             gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)*c5d(i,j,k)
                0301           ENDDO
                0302          ENDDO
                0303         ENDDO
                0304         IF (errCode.GE.1) THEN
                0305           STOP 'MOM_IMPLICIT_R: error when solving 1-Diag problem.'
                0306         ENDIF
                0307       ELSEIF ( diagonalNumber .EQ. 3 ) THEN
b83c015289 Jean*0308 C--   Solve tri-diagonal system :
fa8d87e2db Jean*0309         errCode = -1
b83c015289 Jean*0310         CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
                0311      I                          b5d, c5d, d5d,
8a58850ca8 Jean*0312      U                          gV(1-OLx,1-OLy,1,bi,bj),
b83c015289 Jean*0313      O                          errCode,
                0314      I                          bi, bj, myThid )
                0315         IF (errCode.GE.1) THEN
dcaf8922f9 Jean*0316           STOP 'MOM_IMPLICIT_R: error when solving 3-Diag problem.'
b83c015289 Jean*0317         ENDIF
b76b90dc9b Jean*0318       ELSE
dcaf8922f9 Jean*0319         STOP 'MOM_IMPLICIT_R: no solver available.'
b83c015289 Jean*0320       ENDIF
                0321 
b76b90dc9b Jean*0322 #ifdef ALLOW_SOLVE4_PS_AND_DRAG
                0323       IF ( selectImplicitDrag.EQ.2 ) THEN
                0324        DO k=1,Nr
                0325         DO j=jMin,jMax
                0326          DO i=iMin,iMax
                0327            dV_psFacY(i,j,k,bi,bj) = maskS(i,j,k,bi,bj)
                0328      &                             *recip_deepFacC(k)*recip_rhoFacC(k)
                0329          ENDDO
                0330         ENDDO
                0331        ENDDO
                0332        IF ( diagonalNumber .EQ. 1 ) THEN
                0333         DO k=1,Nr
                0334          DO j=jMin,jMax
                0335           DO i=iMin,iMax
                0336            dV_psFacY(i,j,k,bi,bj) = dV_psFacY(i,j,k,bi,bj)*c5d(i,j,k)
                0337           ENDDO
                0338          ENDDO
                0339         ENDDO
                0340        ELSEIF ( diagonalNumber .EQ. 3 ) THEN
                0341         CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
                0342      I                          b5d, c5d, d5d,
                0343      U                          dV_psFacY(1-OLx,1-OLy,1,bi,bj),
                0344      O                          errCode,
                0345      I                          bi, bj, myThid )
                0346        ENDIF
                0347       ENDIF
                0348 #endif /* ALLOW_SOLVE4_PS_AND_DRAG */
                0349 
1c99f96b44 Jean*0350 #ifdef ALLOW_DIAGNOSTICS
79074ef66b Jean*0351 C--   Diagnostics of momentum dissipation/viscous tendency, implicit part:
                0352       IF ( useDiagnostics .AND. ( implicitViscosity .OR.
                0353      &                            selectImplicitDrag.GE.1 ) ) THEN
                0354        diagName = 'Vm_ImplD'
                0355        IF ( DIAGNOSTICS_IS_ON(diagName,myThid) ) THEN
                0356          recip_dT = 0. _d 0
                0357          IF ( deltaTMom.GT.zeroRL ) recip_dT = oneRL / deltaTMom
                0358          DO k=1,Nr
                0359           DO j=1-OLy,sNy+OLy
                0360            DO i=1-OLx,sNx+OLx
                0361             delVdt(i,j,k) = ( gV(i,j,k,bi,bj)-delVdt(i,j,k) )*recip_dT
                0362            ENDDO
                0363           ENDDO
                0364          ENDDO
                0365          CALL DIAGNOSTICS_FILL(delVdt,diagName, 0,Nr, 2,bi,bj, myThid)
                0366        ENDIF
                0367       ENDIF
                0368 
1c99f96b44 Jean*0369 C--   Diagnostics of vertical viscous flux:
                0370       IF ( useDiagnostics .AND. implicitViscosity ) THEN
79074ef66b Jean*0371        diagName = 'VISrI_Vm'
                0372        IF ( DIAGNOSTICS_IS_ON(diagName,myThid) ) THEN
1c99f96b44 Jean*0373          DO k= 1,Nr
                0374           IF ( k.EQ.1 ) THEN
                0375 C-  Note: Needs to call DIAGNOSTICS_FILL at level k=1 even if array == 0
                0376 C         otherwise counter is not incremented !!
                0377             DO j=1-OLy,sNy+OLy
                0378              DO i=1-OLx,sNx+OLx
                0379                vf(i,j) = 0. _d 0
                0380              ENDDO
                0381             ENDDO
                0382           ELSE
                0383             DO j=jMin,jMax
                0384              DO i=iMin,iMax
e60362cf6f Jean*0385                vf(i,j) = -rAs(i,j,bi,bj)*deepFac2F(k)*rhoFacF(k)
                0386      &            * kappaRV(i,j,k)*recip_drC(k)*rkSign
                0387      &            * (gV(i,j,k,bi,bj) - gV(i,j,k-1,bi,bj))
1c99f96b44 Jean*0388      &            *_maskS(i,j,k,bi,bj)
                0389      &            *_maskS(i,j,k-1,bi,bj)
                0390              ENDDO
                0391             ENDDO
                0392           ENDIF
                0393           CALL DIAGNOSTICS_FILL(vf,diagName, k,1, 2,bi,bj, myThid)
                0394          ENDDO
79074ef66b Jean*0395        ENDIF
1c99f96b44 Jean*0396       ENDIF
                0397 #endif /* ALLOW_DIAGNOSTICS */
                0398 
b83c015289 Jean*0399       RETURN
                0400       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!