Back to home page

MITgcm
 
 

    

File indexing completed on 2026-05-23 05:08:40 UTC

view on githubraw file Latest commit 9b89fcf6 on 2026-05-22 13:35:26 UTC
b83c015289 Jean*0001 #include "GAD_OPTIONS.h"
                0002 
                0003 CBOP
                0004 C     !ROUTINE: GAD_IMPLICIT_R
                0005 C     !INTERFACE:
0845f1a203 Jean*0006       SUBROUTINE GAD_IMPLICIT_R(
b83c015289 Jean*0007      I      implicitAdvection, advectionScheme, tracerIdentity,
4e66ab0b67 Oliv*0008      I      deltaTLev,
09e49e8561 Jean*0009      I      kappaRX, recip_hFac, wFld, tracer,
b83c015289 Jean*0010      U      gTracer,
                0011      I      bi, bj, myTime, myIter, myThid )
1b5fb69d21 Ed H*0012 C     !DESCRIPTION:
                0013 C     Solve implicitly vertical advection and diffusion for one tracer.
b83c015289 Jean*0014 
                0015 C     !USES:
                0016       IMPLICIT NONE
                0017 C     == Global data ==
                0018 #include "SIZE.h"
                0019 #include "EEPARAMS.h"
                0020 #include "PARAMS.h"
                0021 #include "GRID.h"
9b9094601f Jean*0022 #include "SURFACE.h"
b83c015289 Jean*0023 #include "GAD.h"
9b89fcf692 antn*0024 #ifdef ALLOW_LAYERS
                0025 # include "LAYERS_P2SHARE.h"
                0026 #endif
b83c015289 Jean*0027 
1b5fb69d21 Ed H*0028 C !INPUT/OUTPUT PARAMETERS:
                0029 C == Routine Arguments ==
                0030 C implicitAdvection :: if True, treat vertical advection implicitly
                0031 C advectionScheme   :: advection scheme to use
                0032 C tracerIdentity    :: Identifier for the tracer
                0033 C kappaRX           :: 3-D array for vertical diffusion coefficient
9b9094601f Jean*0034 C recip_hFac        :: inverse of cell open-depth factor
09e49e8561 Jean*0035 C wFld              :: Advection velocity field, vertical component
1b5fb69d21 Ed H*0036 C tracer            :: tracer field at current time step
                0037 C gTracer           :: future tracer field
                0038 C bi,bj             :: tile indices
                0039 C myTime            :: current time
                0040 C myIter            :: current iteration number
                0041 C myThid            :: thread number
b83c015289 Jean*0042       LOGICAL implicitAdvection
                0043       INTEGER advectionScheme
                0044       INTEGER tracerIdentity
4e66ab0b67 Oliv*0045       _RL     deltaTLev(Nr)
935a76deec Jean*0046       _RL kappaRX   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0047       _RS recip_hFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0048       _RL wFld      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0049       _RL tracer    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0050       _RL gTracer   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
b83c015289 Jean*0051       INTEGER bi, bj
                0052       _RL myTime
                0053       INTEGER myIter, myThid
                0054 
9b9094601f Jean*0055 #ifdef ALLOW_DIAGNOSTICS
                0056 C !FUNCTIONS:
                0057       CHARACTER*4 GAD_DIAG_SUFX
                0058       EXTERNAL    GAD_DIAG_SUFX
                0059       LOGICAL     DIAGNOSTICS_IS_ON
                0060       EXTERNAL    DIAGNOSTICS_IS_ON
                0061 #endif
                0062 
1b5fb69d21 Ed H*0063 C !LOCAL VARIABLES:
                0064 C == Local variables ==
9b9094601f Jean*0065 C iMin,iMax,jMin,jMax :: computational domain
                0066 C i,j,k     :: loop indices
                0067 C a5d       :: 2nd  lower diagonal of the pentadiagonal matrix
                0068 C b5d       :: 1rst lower diagonal of the pentadiagonal matrix
                0069 C c5d       :: main diagonal       of the pentadiagonal matrix
                0070 C d5d       :: 1rst upper diagonal of the pentadiagonal matrix
                0071 C e5d       :: 2nd  upper diagonal of the pentadiagonal matrix
                0072 C rTrans    :: vertical volume transport at interface k
dd23c36d3b Jean*0073 C localTr   :: local copy of tracer (for Non-Lin Adv.Scheme)
1b5fb69d21 Ed H*0074 C diagonalNumber :: number of non-zero diagonals in the matrix
                0075 C errCode   :: > 0 if singular matrix
1a714a01de Jean*0076 C msgBuf    :: Informational/error message buffer
b83c015289 Jean*0077       INTEGER iMin,iMax,jMin,jMax
9b9094601f Jean*0078       PARAMETER( iMin = 1, iMax = sNx )
                0079       PARAMETER( jMin = 1, jMax = sNy )
b83c015289 Jean*0080       INTEGER i,j,k
                0081       INTEGER diagonalNumber, errCode
935a76deec Jean*0082       _RL a5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0083       _RL b5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0084       _RL c5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0085       _RL d5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0086       _RL e5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
dd23c36d3b Jean*0087       _RL rTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0088       _RL localTr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
81c8d7b9aa Jean*0089 #ifdef ALLOW_DIAGNOSTICS
                0090       CHARACTER*8 diagName
9b9094601f Jean*0091       CHARACTER*4 diagSufx
                0092       LOGICAL     diagDif, diagAdv
                0093       INTEGER km1, km2, kp1, kp2
935a76deec Jean*0094       _RL df (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0095       _RL af (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0096       _RL div(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0097       _RL flx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1a714a01de Jean*0098 # ifdef SOLVE_DIAGONAL_LOWMEMORY
                0099       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0100 # endif /* SOLVE_DIAGONAL_LOWMEMORY */
                0101 #endif /* ALLOW_DIAGNOSTICS */
b83c015289 Jean*0102 CEOP
                0103 
9de7a55d87 Jean*0104 C--   no need to solve anything with only 1 level:
                0105       IF (Nr.GT.1) THEN
b83c015289 Jean*0106 
                0107 C--   Initialise
                0108       DO k=1,Nr
935a76deec Jean*0109        DO j=1-OLy,sNy+OLy
                0110         DO i=1-OLx,sNx+OLx
b83c015289 Jean*0111          a5d(i,j,k) = 0. _d 0
                0112          b5d(i,j,k) = 0. _d 0
                0113          c5d(i,j,k) = 1. _d 0
                0114          d5d(i,j,k) = 0. _d 0
                0115          e5d(i,j,k) = 0. _d 0
                0116         ENDDO
                0117        ENDDO
                0118       ENDDO
                0119       diagonalNumber = 1
                0120 
                0121       IF (implicitDiffusion) THEN
                0122 C--   set the tri-diagonal matrix to solve the implicit diffusion problem
                0123        diagonalNumber = 3
                0124 C-     1rst lower diagonal :
                0125        DO k=2,Nr
                0126         DO j=jMin,jMax
                0127          DO i=iMin,iMax
4e66ab0b67 Oliv*0128            b5d(i,j,k) = -deltaTLev(k)*maskC(i,j,k-1,bi,bj)
9b9094601f Jean*0129      &                  *recip_hFac(i,j,k)*recip_drF(k)
dd23c36d3b Jean*0130      &                  *recip_deepFac2C(k)*recip_rhoFacC(k)
b83c015289 Jean*0131      &                  *kappaRX(i,j, k )*recip_drC( k )
dd23c36d3b Jean*0132      &                  *deepFac2F( k )*rhoFacF( k )
b83c015289 Jean*0133          ENDDO
                0134         ENDDO
                0135        ENDDO
                0136 C-     1rst upper diagonal :
                0137        DO k=1,Nr-1
                0138         DO j=jMin,jMax
                0139          DO i=iMin,iMax
4e66ab0b67 Oliv*0140            d5d(i,j,k) = -deltaTLev(k)*maskC(i,j,k+1,bi,bj)
dd23c36d3b Jean*0141      &                  *recip_hFac(i,j,k)*recip_drF(k)
                0142      &                  *recip_deepFac2C(k)*recip_rhoFacC(k)
                0143      &                  *kappaRX(i,j,k+1)*recip_drC(k+1)
                0144      &                  *deepFac2F(k+1)*rhoFacF(k+1)
b83c015289 Jean*0145          ENDDO
                0146         ENDDO
                0147        ENDDO
                0148 C-     Main diagonal :
                0149        DO k=1,Nr
                0150         DO j=jMin,jMax
                0151          DO i=iMin,iMax
93a010d2da Jean*0152            c5d(i,j,k) = 1. _d 0 - ( b5d(i,j,k) + d5d(i,j,k) )
                0153 C      to recover older (prior to 2016-10-05) results:
                0154 c          c5d(i,j,k) = 1. _d 0 - b5d(i,j,k) - d5d(i,j,k)
b83c015289 Jean*0155          ENDDO
                0156         ENDDO
                0157        ENDDO
                0158 
                0159 C--   end if implicitDiffusion
                0160       ENDIF
                0161 
                0162       IF (implicitAdvection) THEN
                0163 
dd23c36d3b Jean*0164 C--   Non-Linear Advection scheme: keep a local copy of tracer field
                0165        IF ( advectionScheme.EQ.ENUM_FLUX_LIMIT .OR.
                0166      &      advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
                0167         IF ( multiDimAdvection ) THEN
                0168          DO k=1,Nr
                0169           DO j=1-OLy,sNy+OLy
                0170            DO i=1-OLx,sNx+OLx
                0171             localTr(i,j,k) = gTracer(i,j,k)
                0172            ENDDO
                0173           ENDDO
                0174          ENDDO
                0175         ELSE
                0176          DO k=1,Nr
                0177           DO j=1-OLy,sNy+OLy
                0178            DO i=1-OLx,sNx+OLx
                0179             localTr(i,j,k) = tracer(i,j,k)
                0180            ENDDO
                0181           ENDDO
                0182          ENDDO
                0183         ENDIF
                0184        ENDIF
                0185 
75abb052f1 Jean*0186        DO k=Nr,1,-1
b83c015289 Jean*0187 
75abb052f1 Jean*0188 C--    Compute transport
                0189         IF (k.EQ.1) THEN
935a76deec Jean*0190          DO j=1-OLy,sNy+OLy
                0191           DO i=1-OLx,sNx+OLx
0845f1a203 Jean*0192             rTrans(i,j) = 0. _d 0
75abb052f1 Jean*0193           ENDDO
                0194          ENDDO
                0195         ELSE
935a76deec Jean*0196          DO j=1-OLy,sNy+OLy
                0197           DO i=1-OLx,sNx+OLx
09e49e8561 Jean*0198             rTrans(i,j) = wFld(i,j,k)*rA(i,j,bi,bj)
dd23c36d3b Jean*0199      &                               *deepFac2F(k)*rhoFacF(k)
09e49e8561 Jean*0200      &                               *maskC(i,j,k-1,bi,bj)
b83c015289 Jean*0201           ENDDO
75abb052f1 Jean*0202          ENDDO
                0203         ENDIF
b83c015289 Jean*0204 
b7f518a1e0 Jean*0205 #ifdef ALLOW_AIM
                0206 C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
9b9094601f Jean*0207         IF ( k.GE.2 .AND.
                0208      &     (.NOT.useAIM .OR.tracerIdentity.NE.GAD_SALINITY .OR.k.LT.Nr)
b7f518a1e0 Jean*0209      &              ) THEN
                0210 #else
9b9094601f Jean*0211         IF ( k.GE.2 ) THEN
b7f518a1e0 Jean*0212 #endif
75abb052f1 Jean*0213 
                0214          IF ( advectionScheme.EQ.ENUM_CENTERED_2ND ) THEN
                0215           diagonalNumber = 3
                0216           CALL GAD_C2_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
9b9094601f Jean*0217      I                        deltaTLev, rTrans, recip_hFac,
75abb052f1 Jean*0218      U                        b5d, c5d, d5d,
381eb13d88 Jean*0219      I                        myThid )
                0220          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
                0221      &       .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
                0222           diagonalNumber = 3
                0223           CALL GAD_DST2U1_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
9b9094601f Jean*0224      I                        advectionScheme, deltaTLev,
                0225      I                        rTrans, recip_hFac,
381eb13d88 Jean*0226      U                        b5d, c5d, d5d,
                0227      I                        myThid )
                0228          ELSEIF ( advectionScheme.EQ.ENUM_FLUX_LIMIT ) THEN
75abb052f1 Jean*0229           diagonalNumber = 3
                0230           CALL GAD_FLUXLIMIT_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
dd23c36d3b Jean*0231      I                        deltaTLev, rTrans, recip_hFac, localTr,
75abb052f1 Jean*0232      U                        b5d, c5d, d5d,
381eb13d88 Jean*0233      I                        myThid )
                0234          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_3RD
                0235      &       .OR. advectionScheme.EQ.ENUM_CENTERED_4TH
                0236      &       .OR. advectionScheme.EQ.ENUM_DST3 ) THEN
75abb052f1 Jean*0237           diagonalNumber = 5
                0238           CALL GAD_U3C4_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
9b9094601f Jean*0239      I                        advectionScheme, deltaTLev,
                0240      I                        rTrans, recip_hFac,
75abb052f1 Jean*0241      U                        a5d, b5d, c5d, d5d, e5d,
381eb13d88 Jean*0242      I                        myThid )
                0243          ELSEIF ( advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
                0244           diagonalNumber = 5
                0245           CALL GAD_DST3FL_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
dd23c36d3b Jean*0246      I                        deltaTLev, rTrans, recip_hFac, localTr,
381eb13d88 Jean*0247      U                        a5d, b5d, c5d, d5d, e5d,
                0248      I                        myThid )
75abb052f1 Jean*0249          ELSE
                0250           STOP 'GAD_IMPLICIT_R: Adv.Scheme in Impl form not yet coded'
                0251          ENDIF
                0252 
                0253         ENDIF
                0254 
                0255 C--     end k loop
                0256        ENDDO
b83c015289 Jean*0257 
                0258 C--   end if implicitAdvection
                0259       ENDIF
                0260 
                0261       IF ( diagonalNumber .EQ. 3 ) THEN
                0262 C--   Solve tri-diagonal system :
9c3fc51fd7 Jean*0263         errCode = -1
b83c015289 Jean*0264         CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
                0265      I                          b5d, c5d, d5d,
                0266      U                          gTracer,
                0267      O                          errCode,
                0268      I                          bi, bj, myThid )
                0269         IF (errCode.GE.1) THEN
                0270           STOP 'GAD_IMPLICIT_R: error when solving 3-Diag problem'
                0271         ENDIF
75abb052f1 Jean*0272       ELSEIF ( diagonalNumber .EQ. 5 ) THEN
                0273 C--   Solve penta-diagonal system :
9c3fc51fd7 Jean*0274         errCode = -1
75abb052f1 Jean*0275         CALL SOLVE_PENTADIAGONAL( iMin,iMax, jMin,jMax,
                0276      I                            a5d, b5d, c5d, d5d, e5d,
                0277      U                            gTracer,
                0278      O                            errCode,
                0279      I                            bi, bj, myThid )
                0280         IF (errCode.GE.1) THEN
                0281           STOP 'GAD_IMPLICIT_R: error when solving 5-Diag problem'
                0282         ENDIF
b83c015289 Jean*0283       ELSEIF ( diagonalNumber .NE. 1 ) THEN
                0284         STOP 'GAD_IMPLICIT_R: no solver available'
                0285       ENDIF
                0286 
81c8d7b9aa Jean*0287 #ifdef ALLOW_DIAGNOSTICS
                0288 C--   Set diagnostic suffix for the current tracer
9b9094601f Jean*0289       IF ( useDiagnostics ) THEN
81c8d7b9aa Jean*0290         diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )
                0291         diagName = 'DFrI'//diagSufx
9b9094601f Jean*0292         diagDif = implicitDiffusion
                0293         IF ( diagDif ) diagDif = DIAGNOSTICS_IS_ON(diagName,myThid)
9b89fcf692 antn*0294 #ifdef ALLOW_LAYERS
                0295      &                      .OR. layers_useThermo
                0296 #endif
9b9094601f Jean*0297         diagName = 'ADVr'//diagSufx
                0298         diagAdv = implicitAdvection
                0299         IF ( diagAdv ) diagAdv = DIAGNOSTICS_IS_ON(diagName,myThid)
9b89fcf692 antn*0300 #ifdef ALLOW_LAYERS
                0301      &                      .OR. layers_useThermo
                0302 #endif
9b9094601f Jean*0303 
                0304         IF ( diagDif .OR. diagAdv ) THEN
                0305          DO j=1-OLy,sNy+OLy
                0306           DO i=1-OLx,sNx+OLx
919b9da988 Gael*0307             flx(i,j) = 0. _d 0
9b9094601f Jean*0308           ENDDO
                0309          ENDDO
1a714a01de Jean*0310 C--      start diagnostics k loop
9b9094601f Jean*0311          DO k= Nr,1,-1
1a714a01de Jean*0312 
9b9094601f Jean*0313           IF ( implicitDiffusion .AND. k.GE.2 ) THEN
                0314             DO j=jMin,jMax
                0315              DO i=iMin,iMax
                0316                df(i,j) =
6e2ef715a4 Patr*0317 cc#ifdef ALLOW_AUTODIFF_OPENAD
                0318 cc     &             -rA(i,j,bi,bj)%v
                0319 cc#else
dd23c36d3b Jean*0320      &             -rA(i,j,bi,bj)*deepFac2F(k)*rhoFacF(k)
6e2ef715a4 Patr*0321 cc#endif
dd23c36d3b Jean*0322      &            * kappaRX(i,j,k)*recip_drC(k)*rkSign
935a76deec Jean*0323      &            * (gTracer(i,j,k) - gTracer(i,j,k-1))
9b9094601f Jean*0324      &            * maskC(i,j,k,bi,bj)
                0325      &            * maskC(i,j,k-1,bi,bj)
                0326              ENDDO
                0327             ENDDO
                0328           ELSE
81c8d7b9aa Jean*0329             DO j=1-OLy,sNy+OLy
                0330              DO i=1-OLx,sNx+OLx
                0331                df(i,j) = 0. _d 0
                0332              ENDDO
                0333             ENDDO
9b9094601f Jean*0334           ENDIF
1a714a01de Jean*0335 
9b9094601f Jean*0336 C-  Note: Needs to explicitly increment counter (call DIAGNOSTICS_COUNT)
                0337 C         since skipping k=1 DIAGNOSTICS_FILL call.
                0338           IF ( diagDif .AND. k.GE.2 ) THEN
                0339            diagName = 'DFrI'//diagSufx
                0340            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
                0341            IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)
cf336ab6c5 Ryan*0342 #ifdef ALLOW_LAYERS
9b89fcf692 antn*0343            IF ( layers_useThermo ) THEN
50d8304171 Ryan*0344              CALL LAYERS_FILL( df, tracerIdentity, 'DFR',
cf336ab6c5 Ryan*0345      &                           k, 1, 2,bi,bj, myThid )
935a76deec Jean*0346            ENDIF
cf336ab6c5 Ryan*0347 #endif /* ALLOW_LAYERS */
9b9094601f Jean*0348           ENDIF
1a714a01de Jean*0349 
9b9094601f Jean*0350           IF ( diagAdv ) THEN
1a714a01de Jean*0351 #ifdef SOLVE_DIAGONAL_LOWMEMORY
                0352            diagName = 'ADVr'//diagSufx
                0353            WRITE(msgBuf,'(4A)') 'GAD_IMPLICIT_R: ',
                0354      &      'unable to compute Diagnostic "', diagName, '" with'
                0355            CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
                0356      &                         SQUEEZE_RIGHT, myThid )
                0357            WRITE(msgBuf,'(4A)') 'GAD_IMPLICIT_R: ',
                0358      &      '#define SOLVE_DIAGONAL_LOWMEMORY (in CPP_OPTIONS.h)'
                0359            CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
                0360      &                         SQUEEZE_RIGHT, myThid )
                0361            STOP 'ABNORMAL END: S/R GAD_IMPLICIT_R'
                0362 #endif /* SOLVE_DIAGONAL_LOWMEMORY */
9b9094601f Jean*0363            km1=MAX(1,k-1)
                0364            km2=MAX(1,k-2)
                0365            kp1=MIN(Nr,k+1)
                0366            kp2=MIN(Nr,k+2)
                0367 C--   Flux_divergence*deltaT = Tr^n - Tr^n+1 = [A-I](Tr^n+1)
                0368 C                            = deltaT*rkSign*[ Flx_k+1 - Flx_k ]/dz
                0369            DO j=jMin,jMax
                0370             DO i=iMin,iMax
935a76deec Jean*0371               div(i,j) = gTracer(i,j,k)*( c5d(i,j,k) - 1. _d 0 )
                0372      &                 + gTracer(i,j,km1)*b5d(i,j,k)
                0373      &                 + gTracer(i,j,kp1)*d5d(i,j,k)
9b9094601f Jean*0374             ENDDO
                0375            ENDDO
                0376            IF ( diagonalNumber .EQ. 5 ) THEN
                0377             DO j=jMin,jMax
                0378              DO i=iMin,iMax
                0379               div(i,j) = div(i,j)
935a76deec Jean*0380      &                 + gTracer(i,j,km2)*a5d(i,j,k)
                0381      &                 + gTracer(i,j,kp2)*e5d(i,j,k)
81c8d7b9aa Jean*0382              ENDDO
                0383             ENDDO
9b9094601f Jean*0384            ENDIF
                0385 #ifdef NONLIN_FRSURF
                0386            IF ( nonlinFreeSurf.GT.0 ) THEN
                0387 C--    use future hFac to stay consistent with solver matrix
                0388             IF ( select_rStar.GT.0 ) THEN
                0389              DO j=jMin,jMax
                0390               DO i=iMin,iMax
                0391                div(i,j) = div(i,j)*h0FacC(i,j,k,bi,bj)*drF(k)
                0392      &                            *rStarFacC(i,j,bi,bj)
                0393               ENDDO
                0394              ENDDO
                0395             ELSEIF ( selectSigmaCoord.NE.0 ) THEN
                0396              DO j=jMin,jMax
                0397               DO i=iMin,iMax
                0398                div(i,j) = div(i,j)*(
                0399      &               _hFacC(i,j,k,bi,bj)*drF(k)
935a76deec Jean*0400      &              + dBHybSigF(k)*dEtaHdt(i,j,bi,bj)*deltaTFreeSurf
9b9094601f Jean*0401      &                             )
                0402               ENDDO
                0403              ENDDO
                0404             ELSE
                0405              DO j=jMin,jMax
                0406               DO i=iMin,iMax
                0407                IF ( k.EQ.kSurfC(i,j,bi,bj) ) THEN
                0408                 div(i,j) = div(i,j)*hFac_surfC(i,j,bi,bj)*drF(k)
                0409                ELSE
                0410                 div(i,j) = div(i,j)*_hFacC(i,j,k,bi,bj)*drF(k)
                0411                ENDIF
                0412               ENDDO
                0413              ENDDO
                0414             ENDIF
                0415            ELSE
                0416 #else /* NONLIN_FRSURF */
                0417            IF ( .TRUE. ) THEN
                0418 #endif /* NONLIN_FRSURF */
                0419 C--    use current hFac (consistent with solver matrix)
                0420              DO j=jMin,jMax
                0421               DO i=iMin,iMax
                0422                div(i,j) = div(i,j)*_hFacC(i,j,k,bi,bj)*drF(k)
                0423               ENDDO
                0424              ENDDO
                0425            ENDIF
                0426            DO j=jMin,jMax
                0427             DO i=iMin,iMax
919b9da988 Gael*0428               flx(i,j) = flx(i,j)
6e2ef715a4 Patr*0429 cc#ifdef ALLOW_AUTODIFF_OPENAD
                0430 cc     &                - rkSign*div(i,j)*rA(i,j,bi,bj)%v/deltaTLev(k)
                0431 cc#else
dd23c36d3b Jean*0432      &                - rkSign*div(i,j)*rA(i,j,bi,bj)
                0433      &                        *deepFac2C(k)*rhoFacC(k)/deltaTLev(k)
6e2ef715a4 Patr*0434 cc#endif
919b9da988 Gael*0435               af(i,j) = flx(i,j) - df(i,j)
9b9094601f Jean*0436             ENDDO
                0437            ENDDO
                0438            diagName = 'ADVr'//diagSufx
                0439            CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
50d8304171 Ryan*0440 #ifdef ALLOW_LAYERS
9b89fcf692 antn*0441            IF ( layers_useThermo ) THEN
50d8304171 Ryan*0442              CALL LAYERS_FILL(af,tracerIdentity,'AFR',
                0443      &                            k,1,2,bi,bj,myThid)
                0444            ENDIF
                0445 #endif /* ALLOW_LAYERS */
81c8d7b9aa Jean*0446           ENDIF
1a714a01de Jean*0447 
                0448 C--      end diagnostics k loop
81c8d7b9aa Jean*0449          ENDDO
                0450         ENDIF
                0451       ENDIF
                0452 #endif /* ALLOW_DIAGNOSTICS */
                0453 
9de7a55d87 Jean*0454 C--   end if Nr > 1
                0455       ENDIF
                0456 
b83c015289 Jean*0457       RETURN
                0458       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!