pkg/generic_advdiff/gad_implicit_r.F

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