pkg/generic_advdiff/gad_calc_rhs.F

31566b6684 Alis*0001 #include "GAD_OPTIONS.h"
1574069d50 Mart*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
31566b6684 Alis*0005
28d97917ae Alis*0006 CBOP
                0007 C !ROUTINE: GAD_CALC_RHS
                0008
                0009 C !INTERFACE: ==========================================================
0845f1a203 Jean*0010       SUBROUTINE GAD_CALC_RHS(
31566b6684 Alis*0011      I           bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
0845f1a203 Jean*0012      I           xA, yA, maskUp, uFld, vFld, wFld,
                0013      I           uTrans, vTrans, rTrans, rTransKp1,
6b7e107e49 Jean*0014      I           diffKh, diffK4, KappaR, diffKr4, TracerN, TracAB,
cb0d108b91 Jean*0015      I           deltaTLev, trIdentity,
0d75a51072 Mart*0016      I           advectionSchArg, vertAdvecSchArg,
77c3815cb8 Jean*0017      I           calcAdvection, implicitAdvection, applyAB_onTracer,
46918f1b26 Jean*0018      I           trUseDiffKr4, trUseGMRedi, trUseKPP, trUseSmolHack,
a7a3adb269 Jean*0019      O           fZon, fMer,
31566b6684 Alis*0020      U           fVerT, gTracer,
616d5d70aa Jean*0021      I           myTime, myIter, myThid )
31566b6684 Alis*0022
28d97917ae Alis*0023 C !DESCRIPTION:
1bb133c00c Jean*0024 C Calculates the tendency of a tracer due to advection and diffusion.
28d97917ae Alis*0025 C It calculates the fluxes in each direction indepentently and then
1bb133c00c Jean*0026 C sets the tendency to the divergence of these fluxes. The advective
28d97917ae Alis*0027 C fluxes are only calculated here when using the linear advection schemes
                0028 C otherwise only the diffusive and parameterized fluxes are calculated.
                0029 C
                0030 C Contributions to the flux are calculated and added:
                0031 C \begin{equation*}
                0032 C {\bf F} = {\bf F}_{adv} + {\bf F}_{diff} +{\bf F}_{GM} + {\bf F}_{KPP}
                0033 C \end{equation*}
                0034 C
1bb133c00c Jean*0035 C The tendency is the divergence of the fluxes:
28d97917ae Alis*0036 C \begin{equation*}
                0037 C G_\theta = G_\theta + \nabla \cdot {\bf F}
                0038 C \end{equation*}
                0039 C
1bb133c00c Jean*0040 C The tendency is assumed to contain data on entry.
28d97917ae Alis*0041
                0042 C !USES: ===============================================================
                0043       IMPLICIT NONE
31566b6684 Alis*0044 #include "SIZE.h"
                0045 #include "EEPARAMS.h"
                0046 #include "PARAMS.h"
                0047 #include "GRID.h"
e0c3eb6fa1 Jean*0048 #include "SURFACE.h"
31566b6684 Alis*0049 #include "GAD.h"
5127d1d91b Jean*0050 #ifdef ALLOW_AUTODIFF
                0051 # include "AUTODIFF_PARAMS.h"
                0052 #endif /* ALLOW_AUTODIFF */
09c4028516 Patr*0053
28d97917ae Alis*0054 C !INPUT PARAMETERS: ===================================================
b92db9c400 Jean*0055 C bi, bj           :: tile indices
                0056 C iMin, iMax       :: for called routines, to get valid output "gTracer"
                0057 C jMin, jMax       ::                      over this range of indices
0845f1a203 Jean*0058 C k                :: vertical index
                0059 C kM1              :: =k-1 for k>1, =1 for k=1
                0060 C kUp              :: index into 2 1/2D array, toggles between 1|2
                0061 C kDown            :: index into 2 1/2D array, toggles between 2|1
b92db9c400 Jean*0062 C xA, yA           :: areas of X and Y face of tracer cells
0845f1a203 Jean*0063 C maskUp           :: 2-D array for mask at W points
b92db9c400 Jean*0064 C uFld, vFld, wFld :: Local copy of velocity field (3 components)
                0065 C uTrans, vTrans   :: 2-D arrays of volume transports at U,V points
1b5fb69d21 Ed H*0066 C rTrans           :: 2-D arrays of volume transports at W points
                0067 C rTransKp1        :: 2-D array of volume trans at W pts, interf k+1
                0068 C diffKh           :: horizontal diffusion coefficient
6b7e107e49 Jean*0069 C diffK4           :: horizontal bi-harmonic diffusion coefficient
5a39f5fdb9 Jean*0070 C KappaR           :: 2-D array for vertical diffusion coefficient, interf k
6b7e107e49 Jean*0071 C diffKr4          :: 1-D array for vertical bi-harmonic diffusion coefficient
77c3815cb8 Jean*0072 C TracerN          :: tracer field @ time-step n (Note: only used
                0073 C                     if applying AB on tracer field rather than on tendency gTr)
                0074 C TracAB           :: current tracer field (@ time-step n if applying AB on gTr
a8c74226d9 Jean*0075 C                     or extrapolated fwd in time to n+1/2 if applying AB on Tr)
cb0d108b91 Jean*0076 C trIdentity       :: tracer identifier (required for KPP,GM)
0d75a51072 Mart*0077 C advectionSchArg  :: advection scheme to use (Horizontal plane)
                0078 C vertAdvecSchArg  :: advection scheme to use (Vertical direction)
1b5fb69d21 Ed H*0079 C calcAdvection    :: =False if Advec computed with multiDim scheme
                0080 C implicitAdvection:: =True if vertical Advec computed implicitly
77c3815cb8 Jean*0081 C applyAB_onTracer :: apply Adams-Bashforth on Tracer (rather than on gTr)
6b7e107e49 Jean*0082 C trUseDiffKr4     :: true if this tracer uses vertical bi-harmonic diffusion
6ea562f7f0 Jean*0083 C trUseGMRedi      :: true if this tracer uses GM-Redi
                0084 C trUseKPP         :: true if this tracer uses KPP
46918f1b26 Jean*0085 C trUseSmolHack    :: true if this tracer uses Smolarkiewicz-Hack to remain > 0
616d5d70aa Jean*0086 C myTime           :: current time
                0087 C myIter           :: iteration number
1b5fb69d21 Ed H*0088 C myThid           :: thread number
31566b6684 Alis*0089       INTEGER bi,bj,iMin,iMax,jMin,jMax
28d97917ae Alis*0090       INTEGER k,kUp,kDown,kM1
31566b6684 Alis*0091       _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0092       _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
0845f1a203 Jean*0093       _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0094       _RL uFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0095       _RL vFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0096       _RL wFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
31566b6684 Alis*0097       _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0098       _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0099       _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
75abb052f1 Jean*0100       _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
31566b6684 Alis*0101       _RL diffKh, diffK4
5a39f5fdb9 Jean*0102       _RL KappaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
6b7e107e49 Jean*0103       _RL diffKr4(Nr)
e9de1d7682 Jean*0104       _RL TracerN(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0105       _RL TracAB (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
4e66ab0b67 Oliv*0106       _RL deltaTLev(Nr)
cb0d108b91 Jean*0107       INTEGER trIdentity
0d75a51072 Mart*0108       INTEGER advectionSchArg, vertAdvecSchArg
6ea562f7f0 Jean*0109       LOGICAL calcAdvection
77c3815cb8 Jean*0110       LOGICAL implicitAdvection, applyAB_onTracer
46918f1b26 Jean*0111       LOGICAL trUseDiffKr4, trUseGMRedi, trUseKPP, trUseSmolHack
616d5d70aa Jean*0112       _RL     myTime
                0113       INTEGER myIter, myThid
31566b6684 Alis*0114
28d97917ae Alis*0115 C !OUTPUT PARAMETERS: ==================================================
1bb133c00c Jean*0116 C gTracer          :: tendency array
a7a3adb269 Jean*0117 C fZon             :: zonal flux
                0118 C fMer             :: meridional flux
1b5fb69d21 Ed H*0119 C fVerT            :: 2 1/2D arrays for vertical advective flux
935a76deec Jean*0120       _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
a7a3adb269 Jean*0121       _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0122       _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
28d97917ae Alis*0123       _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0124
cb0d108b91 Jean*0125 C !FUNCTIONS:       ====================================================
                0126 #ifdef ALLOW_DIAGNOSTICS
                0127       CHARACTER*4 GAD_DIAG_SUFX
                0128       EXTERNAL    GAD_DIAG_SUFX
                0129 #endif /* ALLOW_DIAGNOSTICS */
                0130
28d97917ae Alis*0131 C !LOCAL VARIABLES: ====================================================
1b5fb69d21 Ed H*0132 C i,j              :: loop indices
0d75a51072 Mart*0133 C advectionScheme  :: local copy of routine argument advectionSchArg
                0134 C vertAdvecScheme  :: local copy of routine argument vertAdvecSchArg
1b5fb69d21 Ed H*0135 C df4              :: used for storing del^2 T for bi-harmonic term
                0136 C af               :: advective flux
                0137 C df               :: diffusive flux
                0138 C localT           :: local copy of tracer field
1bb133c00c Jean*0139 C locABT           :: local copy of (AB-extrapolated) tracer field
31566b6684 Alis*0140       INTEGER i,j
0d75a51072 Mart*0141       INTEGER advectionScheme, vertAdvecScheme
2a09aadbbc Jean*0142       _RS maskLocW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0143       _RS maskLocS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
31566b6684 Alis*0144       _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0145       _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0146       _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0147       _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1bb133c00c Jean*0148       _RL locABT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
75abb052f1 Jean*0149       _RL advFac, rAdvFac
998cbd211e Oliv*0150 #ifdef GAD_SMOLARKIEWICZ_HACK
622cf2fb44 Jean*0151       _RL outFlux, trac, gTrFac
998cbd211e Oliv*0152 #endif
cb0d108b91 Jean*0153 #ifdef ALLOW_DIAGNOSTICS
                0154       CHARACTER*8 diagName
                0155       CHARACTER*4 diagSufx
                0156 #endif
28d97917ae Alis*0157 CEOP
31566b6684 Alis*0158
0d75a51072 Mart*0159 C     make local copies to be tampered with if necessary
                0160       advectionScheme = advectionSchArg
                0161       vertAdvecScheme = vertAdvecSchArg
64442af1fe Jean*0162 #ifdef ALLOW_AUTODIFF
1574069d50 Mart*0163 #ifdef ALLOW_AUTODIFF_TAMC
                0164 CADJ INIT gad_local_tape = COMMON, 1
                0165 C     This store directive just suppresses a recomputation warning.
                0166 C     TAF generates an extra field with or without this directive.
                0167 CADJ STORE fvert = gad_local_tape
                0168 #endif
31566b6684 Alis*0169 C--   only the kUp part of fverT is set in this subroutine
                0170 C--   the kDown is still required
                0171       fVerT(1,1,kDown) = fVerT(1,1,kDown)
0d75a51072 Mart*0172 C
                0173       IF ( inAdMode .AND. useApproxAdvectionInAdMode ) THEN
                0174 C     In AD-mode, we change non-linear, potentially unstable AD advection
                0175 C     schemes to linear schemes with more stability. So far only DST3 with
                0176 C     flux limiting is replaced by DST3 without flux limiting, but any
                0177 C     combination is possible.
                0178        IF ( advectionSchArg.EQ.ENUM_DST3_FLUX_LIMIT )
                0179      &      advectionScheme = ENUM_DST3
                0180        IF ( vertAdvecSchArg.EQ.ENUM_DST3_FLUX_LIMIT )
                0181      &      vertAdvecScheme = ENUM_DST3
                0182 C     here is room for more advection schemes as this becomes necessary
                0183       ENDIF
                0184 #endif /* ALLOW_AUTODIFF */
09c4028516 Patr*0185
81c8d7b9aa Jean*0186 #ifdef ALLOW_DIAGNOSTICS
                0187 C--   Set diagnostic suffix for the current tracer
                0188       IF ( useDiagnostics ) THEN
cb0d108b91 Jean*0189         diagSufx = GAD_DIAG_SUFX( trIdentity, myThid )
81c8d7b9aa Jean*0190       ENDIF
                0191 #endif
                0192
75abb052f1 Jean*0193       advFac  = 0. _d 0
                0194       IF (calcAdvection) advFac = 1. _d 0
bb6c554092 Jean*0195       rAdvFac = rkSign*advFac
95826db035 Jean*0196       IF (implicitAdvection) rAdvFac = rkSign
75abb052f1 Jean*0197
31566b6684 Alis*0198       DO j=1-OLy,sNy+OLy
                0199        DO i=1-OLx,sNx+OLx
67a1e439d8 Patr*0200         fZon(i,j)      = 0. _d 0
                0201         fMer(i,j)      = 0. _d 0
                0202         fVerT(i,j,kUp) = 0. _d 0
09c4028516 Patr*0203         df(i,j)        = 0. _d 0
                0204         df4(i,j)       = 0. _d 0
31566b6684 Alis*0205        ENDDO
                0206       ENDDO
                0207
                0208 C--   Make local copy of tracer array
77c3815cb8 Jean*0209       IF ( applyAB_onTracer ) THEN
                0210         DO j=1-OLy,sNy+OLy
                0211          DO i=1-OLx,sNx+OLx
e9de1d7682 Jean*0212           localT(i,j)=TracerN(i,j,k)
                0213           locABT(i,j)= TracAB(i,j,k)
77c3815cb8 Jean*0214          ENDDO
                0215         ENDDO
                0216       ELSE
                0217         DO j=1-OLy,sNy+OLy
                0218          DO i=1-OLx,sNx+OLx
e9de1d7682 Jean*0219           localT(i,j)=TracerN(i,j,k)
                0220           locABT(i,j)=TracerN(i,j,k)
77c3815cb8 Jean*0221          ENDDO
                0222         ENDDO
                0223       ENDIF
31566b6684 Alis*0224
                0225 C--   Pre-calculate del^2 T if bi-harmonic coefficient is non-zero
                0226       IF (diffK4 .NE. 0.) THEN
1574069d50 Mart*0227 #ifdef ALLOW_AUTODIFF_TAMC
                0228 C     These store directives suppress two recomputation warnings and
                0229 C     avoid the corresponding recomputations of fZon, fMer at the cost
                0230 C     of only four extra local 2D-fields.
                0231 CADJ STORE localT = gad_local_tape
                0232 #endif
31566b6684 Alis*0233        CALL GAD_GRAD_X(bi,bj,k,xA,localT,fZon,myThid)
1574069d50 Mart*0234 #ifdef ALLOW_AUTODIFF_TAMC
                0235 CADJ STORE localT = gad_local_tape
                0236 #endif
31566b6684 Alis*0237        CALL GAD_GRAD_Y(bi,bj,k,yA,localT,fMer,myThid)
1574069d50 Mart*0238 #ifdef ALLOW_AUTODIFF_TAMC
                0239 CADJ STORE fZon, fMer = gad_local_tape
                0240 #endif
31566b6684 Alis*0241        CALL GAD_DEL2(bi,bj,k,fZon,fMer,df4,myThid)
                0242       ENDIF
                0243
                0244 C--   Initialize net flux in X direction
cb0d108b91 Jean*0245       DO j=1-OLy,sNy+OLy
                0246        DO i=1-OLx,sNx+OLx
67a1e439d8 Patr*0247         fZon(i,j) = 0. _d 0
31566b6684 Alis*0248        ENDDO
                0249       ENDDO
                0250
                0251 C-    Advective flux in X
c17bf9e7ce Jean*0252       IF (calcAdvection) THEN
0d75a51072 Mart*0253         IF ( advectionScheme.EQ.ENUM_CENTERED_2ND ) THEN
2a09aadbbc Jean*0254            CALL GAD_C2_ADV_X( bi,bj,k, uTrans, locABT, af, myThid )
0845f1a203 Jean*0255         ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
a90759ed9a Jean*0256      &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
2a09aadbbc Jean*0257            CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme, .TRUE.,
                0258      I              deltaTLev(k), uTrans, uFld, locABT,
                0259      O              af, myThid )
                0260         ELSE
                0261          DO j=1-OLy,sNy+OLy
                0262           DO i=1-OLx,sNx+OLx
                0263 #ifdef ALLOW_OBCS
                0264            maskLocW(i,j) = _maskW(i,j,k,bi,bj)*maskInW(i,j,bi,bj)
                0265 #else /* ALLOW_OBCS */
                0266            maskLocW(i,j) = _maskW(i,j,k,bi,bj)
                0267 #endif /* ALLOW_OBCS */
                0268           ENDDO
                0269          ENDDO
0d75a51072 Mart*0270          IF ( advectionScheme.EQ.ENUM_FLUX_LIMIT ) THEN
2a09aadbbc Jean*0271            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
                0272      I              uTrans, uFld, maskLocW, locABT,
                0273      O              af, myThid )
0d75a51072 Mart*0274          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
188f784a21 Jean*0275            CALL GAD_U3_ADV_X( bi,bj,k, uTrans, maskLocW, locABT,
2a09aadbbc Jean*0276      O              af, myThid )
0d75a51072 Mart*0277          ELSEIF ( advectionScheme.EQ.ENUM_CENTERED_4TH ) THEN
2a09aadbbc Jean*0278            CALL GAD_C4_ADV_X( bi,bj,k, uTrans, maskLocW, locABT,
                0279      O              af, myThid )
0d75a51072 Mart*0280          ELSEIF ( advectionScheme.EQ.ENUM_DST3 ) THEN
2a09aadbbc Jean*0281            CALL GAD_DST3_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
                0282      I              uTrans, uFld, maskLocW, locABT,
                0283      O              af, myThid )
0d75a51072 Mart*0284          ELSEIF ( advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
2a09aadbbc Jean*0285            CALL GAD_DST3FL_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
                0286      I              uTrans, uFld, maskLocW, locABT,
                0287      O              af, myThid )
0d75a51072 Mart*0288          ELSEIF ( advectionScheme.EQ.ENUM_OS7MP ) THEN
2a09aadbbc Jean*0289            CALL GAD_OS7MP_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
                0290      I              uTrans, uFld, maskLocW, locABT,
                0291      O              af, myThid )
                0292          ELSE
                0293           STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'
                0294          ENDIF
81c8d7b9aa Jean*0295         ENDIF
cb0d108b91 Jean*0296 #ifdef ALLOW_OBCS
                0297         IF ( useOBCS ) THEN
                0298 C-      replace advective flux with 1st order upwind scheme estimate
                0299           CALL OBCS_U1_ADV_TRACER( .TRUE., trIdentity, bi, bj, k,
                0300      I                             maskW(1-OLx,1-OLy,k,bi,bj),
                0301      I                             uTrans, locABT,
                0302      U                             af, myThid )
                0303         ENDIF
                0304 #endif /* ALLOW_OBCS */
                0305         DO j=1-OLy,sNy+OLy
                0306          DO i=1-OLx,sNx+OLx
81c8d7b9aa Jean*0307           fZon(i,j) = fZon(i,j) + af(i,j)
                0308          ENDDO
                0309         ENDDO
                0310 #ifdef ALLOW_DIAGNOSTICS
                0311         IF ( useDiagnostics ) THEN
                0312           diagName = 'ADVx'//diagSufx
cb0d108b91 Jean*0313           CALL DIAGNOSTICS_FILL( af, diagName, k,1, 2,bi,bj, myThid )
81c8d7b9aa Jean*0314         ENDIF
50d8304171 Ryan*0315 #ifdef ALLOW_LAYERS
                0316         IF ( useLayers ) THEN
9897ae9d92 Jean*0317           CALL LAYERS_FILL( af, trIdentity, 'AFX',
50d8304171 Ryan*0318      &                      k, 1, 2,bi,bj, myThid )
                0319         ENDIF
                0320 #endif /* ENDIF */
81c8d7b9aa Jean*0321 #endif
24016b3859 Alis*0322       ENDIF
31566b6684 Alis*0323
                0324 C-    Diffusive flux in X
                0325       IF (diffKh.NE.0.) THEN
                0326        CALL GAD_DIFF_X(bi,bj,k,xA,diffKh,localT,df,myThid)
                0327       ELSE
cb0d108b91 Jean*0328        DO j=1-OLy,sNy+OLy
                0329         DO i=1-OLx,sNx+OLx
67a1e439d8 Patr*0330          df(i,j) = 0. _d 0
31566b6684 Alis*0331         ENDDO
                0332        ENDDO
                0333       ENDIF
                0334
81c8d7b9aa Jean*0335 C-    Add bi-harmonic diffusive flux in X
                0336       IF (diffK4 .NE. 0.) THEN
                0337        CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)
                0338       ENDIF
                0339
31566b6684 Alis*0340 #ifdef ALLOW_GMREDI
                0341 C-    GM/Redi flux in X
8e4e737494 Jean*0342       IF ( trUseGMRedi ) THEN
972c0130ec Jean*0343         CALL GMREDI_XTRANSPORT(
b92db9c400 Jean*0344      I         trIdentity, bi, bj, k, iMin, iMax+1, jMin, jMax,
8233d0ceb9 Jean*0345      I         xA, maskUp, TracerN,
77c3815cb8 Jean*0346      U         df,
e9de1d7682 Jean*0347      I         myThid )
31566b6684 Alis*0348       ENDIF
                0349 #endif
eaba2fd266 Jean*0350 C     anelastic: advect.fluxes are scaled by rhoFac but hor.diff. flx are not
cb0d108b91 Jean*0351       DO j=1-OLy,sNy+OLy
                0352        DO i=1-OLx,sNx+OLx
eaba2fd266 Jean*0353         fZon(i,j) = fZon(i,j) + df(i,j)*rhoFacC(k)
31566b6684 Alis*0354        ENDDO
                0355       ENDDO
                0356
81c8d7b9aa Jean*0357 #ifdef ALLOW_DIAGNOSTICS
                0358 C-    Diagnostics of Tracer flux in X dir (mainly Diffusive term),
                0359 C       excluding advective terms:
                0360       IF ( useDiagnostics .AND.
8e4e737494 Jean*0361      &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. trUseGMRedi) ) THEN
f68e77b60d Jean*0362           diagName = 'DFxE'//diagSufx
cb0d108b91 Jean*0363           CALL DIAGNOSTICS_FILL( df, diagName, k,1, 2,bi,bj, myThid )
cf336ab6c5 Ryan*0364 #ifdef ALLOW_LAYERS
ee16a2cae4 Ryan*0365           IF ( useLayers ) THEN
50d8304171 Ryan*0366            CALL LAYERS_FILL( df, trIdentity, 'DFX',
                0367      &                       k, 1, 2,bi,bj, myThid )
ee16a2cae4 Ryan*0368           ENDIF
cf336ab6c5 Ryan*0369 #endif /* ALLOW_LAYERS */
31566b6684 Alis*0370       ENDIF
81c8d7b9aa Jean*0371 #endif
31566b6684 Alis*0372
                0373 C--   Initialize net flux in Y direction
cb0d108b91 Jean*0374       DO j=1-OLy,sNy+OLy
                0375        DO i=1-OLx,sNx+OLx
67a1e439d8 Patr*0376         fMer(i,j) = 0. _d 0
31566b6684 Alis*0377        ENDDO
                0378       ENDDO
                0379
                0380 C-    Advective flux in Y
c17bf9e7ce Jean*0381       IF (calcAdvection) THEN
0d75a51072 Mart*0382         IF ( advectionScheme.EQ.ENUM_CENTERED_2ND ) THEN
2a09aadbbc Jean*0383            CALL GAD_C2_ADV_Y( bi,bj,k, vTrans, locABT, af, myThid )
0845f1a203 Jean*0384         ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
a90759ed9a Jean*0385      &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
2a09aadbbc Jean*0386            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .TRUE.,
                0387      I              deltaTLev(k), vTrans, vFld, locABT,
                0388      O              af, myThid )
                0389         ELSE
                0390          DO j=1-OLy,sNy+OLy
                0391           DO i=1-OLx,sNx+OLx
                0392 #ifdef ALLOW_OBCS
                0393            maskLocS(i,j) = _maskS(i,j,k,bi,bj)*maskInS(i,j,bi,bj)
                0394 #else /* ALLOW_OBCS */
                0395            maskLocS(i,j) = _maskS(i,j,k,bi,bj)
                0396 #endif /* ALLOW_OBCS */
                0397           ENDDO
                0398          ENDDO
0d75a51072 Mart*0399          IF ( advectionScheme.EQ.ENUM_FLUX_LIMIT ) THEN
2a09aadbbc Jean*0400            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
                0401      I              vTrans, vFld, maskLocS, locABT,
                0402      O              af, myThid )
0d75a51072 Mart*0403          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
2a09aadbbc Jean*0404            CALL GAD_U3_ADV_Y( bi,bj,k, vTrans, maskLocS, locABT,
                0405      O              af, myThid )
0d75a51072 Mart*0406          ELSEIF ( advectionScheme.EQ.ENUM_CENTERED_4TH ) THEN
2a09aadbbc Jean*0407            CALL GAD_C4_ADV_Y( bi,bj,k, vTrans, maskLocS, locABT,
                0408      O              af, myThid )
0d75a51072 Mart*0409          ELSEIF ( advectionScheme.EQ.ENUM_DST3 ) THEN
2a09aadbbc Jean*0410            CALL GAD_DST3_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
                0411      I              vTrans, vFld, maskLocS, locABT,
                0412      O              af, myThid )
0d75a51072 Mart*0413          ELSEIF ( advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
2a09aadbbc Jean*0414            CALL GAD_DST3FL_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
                0415      I              vTrans, vFld, maskLocS, locABT,
                0416      O              af, myThid )
0d75a51072 Mart*0417          ELSEIF ( advectionScheme.EQ.ENUM_OS7MP ) THEN
2a09aadbbc Jean*0418            CALL GAD_OS7MP_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
                0419      I              vTrans, vFld, maskLocS, locABT,
                0420      O              af, myThid )
                0421          ELSE
                0422            STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'
                0423          ENDIF
81c8d7b9aa Jean*0424         ENDIF
cb0d108b91 Jean*0425 #ifdef ALLOW_OBCS
                0426         IF ( useOBCS ) THEN
                0427 C-      replace advective flux with 1st order upwind scheme estimate
                0428           CALL OBCS_U1_ADV_TRACER( .FALSE., trIdentity, bi, bj, k,
                0429      I                             maskS(1-OLx,1-OLy,k,bi,bj),
                0430      I                             vTrans, locABT,
                0431      U                             af, myThid )
                0432         ENDIF
                0433 #endif /* ALLOW_OBCS */
                0434         DO j=1-OLy,sNy+OLy
                0435          DO i=1-OLx,sNx+OLx
81c8d7b9aa Jean*0436           fMer(i,j) = fMer(i,j) + af(i,j)
                0437          ENDDO
                0438         ENDDO
                0439 #ifdef ALLOW_DIAGNOSTICS
                0440         IF ( useDiagnostics ) THEN
                0441           diagName = 'ADVy'//diagSufx
cb0d108b91 Jean*0442           CALL DIAGNOSTICS_FILL( af, diagName, k,1, 2,bi,bj, myThid )
81c8d7b9aa Jean*0443         ENDIF
50d8304171 Ryan*0444 #ifdef ALLOW_LAYERS
                0445         IF ( useLayers ) THEN
                0446           CALL LAYERS_FILL( af, trIdentity, 'AFY',
                0447      &                          k, 1, 2,bi,bj, myThid )
                0448         ENDIF
                0449 #endif /* ALLOW_LAYES */
81c8d7b9aa Jean*0450 #endif
24016b3859 Alis*0451       ENDIF
31566b6684 Alis*0452
                0453 C-    Diffusive flux in Y
                0454       IF (diffKh.NE.0.) THEN
                0455        CALL GAD_DIFF_Y(bi,bj,k,yA,diffKh,localT,df,myThid)
                0456       ELSE
cb0d108b91 Jean*0457        DO j=1-OLy,sNy+OLy
                0458         DO i=1-OLx,sNx+OLx
67a1e439d8 Patr*0459          df(i,j) = 0. _d 0
31566b6684 Alis*0460         ENDDO
                0461        ENDDO
                0462       ENDIF
                0463
81c8d7b9aa Jean*0464 C-    Add bi-harmonic flux in Y
                0465       IF (diffK4 .NE. 0.) THEN
                0466        CALL GAD_BIHARM_Y(bi,bj,k,yA,df4,diffK4,df,myThid)
                0467       ENDIF
                0468
31566b6684 Alis*0469 #ifdef ALLOW_GMREDI
                0470 C-    GM/Redi flux in Y
8e4e737494 Jean*0471       IF ( trUseGMRedi ) THEN
972c0130ec Jean*0472         CALL GMREDI_YTRANSPORT(
b92db9c400 Jean*0473      I         trIdentity, bi, bj, k, iMin, iMax, jMin, jMax+1,
8233d0ceb9 Jean*0474      I         yA, maskUp, TracerN,
77c3815cb8 Jean*0475      U         df,
e9de1d7682 Jean*0476      I         myThid )
31566b6684 Alis*0477       ENDIF
                0478 #endif
eaba2fd266 Jean*0479 C     anelastic: advect.fluxes are scaled by rhoFac but hor.diff. flx are not
cb0d108b91 Jean*0480       DO j=1-OLy,sNy+OLy
                0481        DO i=1-OLx,sNx+OLx
eaba2fd266 Jean*0482         fMer(i,j) = fMer(i,j) + df(i,j)*rhoFacC(k)
31566b6684 Alis*0483        ENDDO
                0484       ENDDO
                0485
81c8d7b9aa Jean*0486 #ifdef ALLOW_DIAGNOSTICS
                0487 C-    Diagnostics of Tracer flux in Y dir (mainly Diffusive terms),
                0488 C       excluding advective terms:
                0489       IF ( useDiagnostics .AND.
8e4e737494 Jean*0490      &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. trUseGMRedi) ) THEN
f68e77b60d Jean*0491           diagName = 'DFyE'//diagSufx
cb0d108b91 Jean*0492           CALL DIAGNOSTICS_FILL( df, diagName, k,1, 2,bi,bj, myThid )
cf336ab6c5 Ryan*0493 #ifdef ALLOW_LAYERS
ee16a2cae4 Ryan*0494           IF ( useLayers ) THEN
50d8304171 Ryan*0495            CALL LAYERS_FILL( df, trIdentity, 'DFY',
                0496      &                           k, 1, 2,bi,bj, myThid )
ee16a2cae4 Ryan*0497           ENDIF
cf336ab6c5 Ryan*0498 #endif /* ALLOW_LAYERS */
31566b6684 Alis*0499       ENDIF
81c8d7b9aa Jean*0500 #endif
31566b6684 Alis*0501
e0c3eb6fa1 Jean*0502 C--   Compute vertical flux fVerT(kUp) at interface k (between k-1 & k):
31566b6684 Alis*0503 C-    Advective flux in R
ca239d4b54 Jean*0504 #ifdef ALLOW_AIM
                0505 C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
77c3815cb8 Jean*0506       IF (calcAdvection .AND. .NOT.implicitAdvection .AND. k.GE.2 .AND.
cb0d108b91 Jean*0507      &     (.NOT.useAIM .OR. trIdentity.NE.GAD_SALINITY .OR. k.LT.Nr)
ca239d4b54 Jean*0508      &   ) THEN
                0509 #else
77c3815cb8 Jean*0510       IF (calcAdvection .AND. .NOT.implicitAdvection .AND. k.GE.2) THEN
ca239d4b54 Jean*0511 #endif
972c0130ec Jean*0512        IF ( applyAB_onTracer ) THEN
                0513 C-    Compute vertical advective flux in the interior using TracAB:
0d75a51072 Mart*0514         IF ( vertAdvecScheme.EQ.ENUM_CENTERED_2ND ) THEN
cb0d108b91 Jean*0515            CALL GAD_C2_ADV_R( bi,bj,k, rTrans, TracAB, af, myThid )
0845f1a203 Jean*0516         ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
a90759ed9a Jean*0517      &          .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
cb0d108b91 Jean*0518            CALL GAD_DST2U1_ADV_R( bi,bj,k,vertAdvecScheme,deltaTLev(k),
e9de1d7682 Jean*0519      I              rTrans, wFld, TracAB,
cb0d108b91 Jean*0520      O              af, myThid )
0d75a51072 Mart*0521         ELSEIF ( vertAdvecScheme.EQ.ENUM_FLUX_LIMIT ) THEN
cb0d108b91 Jean*0522            CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, deltaTLev(k),
e9de1d7682 Jean*0523      I              rTrans, wFld, TracAB,
cb0d108b91 Jean*0524      O              af, myThid )
0d75a51072 Mart*0525         ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN
cb0d108b91 Jean*0526            CALL GAD_U3_ADV_R( bi,bj,k, rTrans, TracAB, af, myThid )
81c8d7b9aa Jean*0527         ELSEIF (vertAdvecScheme.EQ.ENUM_CENTERED_4TH) THEN
cb0d108b91 Jean*0528            CALL GAD_C4_ADV_R( bi,bj,k, rTrans, TracAB, af, myThid )
0d75a51072 Mart*0529         ELSEIF ( vertAdvecScheme.EQ.ENUM_DST3 ) THEN
cb0d108b91 Jean*0530            CALL GAD_DST3_ADV_R( bi,bj,k, deltaTLev(k),
e9de1d7682 Jean*0531      I              rTrans, wFld, TracAB,
cb0d108b91 Jean*0532      O              af, myThid )
0d75a51072 Mart*0533         ELSEIF ( vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
cb0d108b91 Jean*0534            CALL GAD_DST3FL_ADV_R( bi,bj,k, deltaTLev(k),
e9de1d7682 Jean*0535      I              rTrans, wFld, TracAB,
cb0d108b91 Jean*0536      O              af, myThid )
0d75a51072 Mart*0537         ELSEIF ( vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
cb0d108b91 Jean*0538            CALL GAD_OS7MP_ADV_R( bi,bj,k, deltaTLev(k),
e9de1d7682 Jean*0539      I              rTrans, wFld, TracAB,
cb0d108b91 Jean*0540      O              af, myThid )
81c8d7b9aa Jean*0541         ELSE
                0542           STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'
                0543         ENDIF
972c0130ec Jean*0544        ELSE
                0545 C-    Compute vertical advective flux in the interior using TracerN:
0d75a51072 Mart*0546         IF ( vertAdvecScheme.EQ.ENUM_CENTERED_2ND ) THEN
972c0130ec Jean*0547            CALL GAD_C2_ADV_R( bi,bj,k, rTrans, TracerN, af, myThid )
                0548         ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
                0549      &          .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
                0550            CALL GAD_DST2U1_ADV_R( bi,bj,k,vertAdvecScheme,deltaTLev(k),
e9de1d7682 Jean*0551      I              rTrans, wFld, TracerN,
972c0130ec Jean*0552      O              af, myThid )
0d75a51072 Mart*0553         ELSEIF ( vertAdvecScheme.EQ.ENUM_FLUX_LIMIT ) THEN
972c0130ec Jean*0554            CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, deltaTLev(k),
e9de1d7682 Jean*0555      I              rTrans, wFld, TracerN,
972c0130ec Jean*0556      O              af, myThid )
0d75a51072 Mart*0557         ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN
972c0130ec Jean*0558            CALL GAD_U3_ADV_R( bi,bj,k, rTrans, TracerN, af, myThid )
0d75a51072 Mart*0559         ELSEIF ( vertAdvecScheme.EQ.ENUM_CENTERED_4TH ) THEN
972c0130ec Jean*0560            CALL GAD_C4_ADV_R( bi,bj,k, rTrans, TracerN, af, myThid )
                0561         ELSEIF( vertAdvecScheme.EQ.ENUM_DST3 ) THEN
                0562            CALL GAD_DST3_ADV_R( bi,bj,k, deltaTLev(k),
e9de1d7682 Jean*0563      I              rTrans, wFld, TracerN,
972c0130ec Jean*0564      O              af, myThid )
0d75a51072 Mart*0565         ELSEIF ( vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
972c0130ec Jean*0566            CALL GAD_DST3FL_ADV_R( bi,bj,k, deltaTLev(k),
e9de1d7682 Jean*0567      I              rTrans, wFld, TracerN,
972c0130ec Jean*0568      O              af, myThid )
0d75a51072 Mart*0569         ELSEIF ( vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
972c0130ec Jean*0570            CALL GAD_OS7MP_ADV_R( bi,bj,k, deltaTLev(k),
e9de1d7682 Jean*0571      I              rTrans, wFld, TracerN,
972c0130ec Jean*0572      O              af, myThid )
                0573         ELSE
                0574           STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'
                0575         ENDIF
                0576        ENDIF
75abb052f1 Jean*0577 C-     add the advective flux to fVerT
972c0130ec Jean*0578        DO j=1-OLy,sNy+OLy
                0579         DO i=1-OLx,sNx+OLx
4d24da31b2 Mart*0580           fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)*maskInC(i,j,bi,bj)
fb78c12fcf Jean*0581         ENDDO
972c0130ec Jean*0582        ENDDO
81c8d7b9aa Jean*0583 #ifdef ALLOW_DIAGNOSTICS
972c0130ec Jean*0584        IF ( useDiagnostics ) THEN
81c8d7b9aa Jean*0585           diagName = 'ADVr'//diagSufx
cb0d108b91 Jean*0586           CALL DIAGNOSTICS_FILL( af, diagName, k,1, 2,bi,bj, myThid )
ed8ddd9828 Jean*0587 C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
                0588 C        does it only if k=1 (never the case here)
                0589           IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)
50d8304171 Ryan*0590 #ifdef ALLOW_LAYERS
                0591           IF ( useLayers ) THEN
                0592             CALL LAYERS_FILL(af,trIdentity,'AFR',k,1,2,bi,bj,myThid)
                0593           ENDIF
                0594 #endif /* ALLOW_LAYERS */
972c0130ec Jean*0595        ENDIF
81c8d7b9aa Jean*0596 #endif
24016b3859 Alis*0597       ENDIF
31566b6684 Alis*0598
                0599 C-    Diffusive flux in R
                0600 C     Note: For K=1 then KM1=1 and this gives a dT/dr = 0 upper
                0601 C           boundary condition.
                0602       IF (implicitDiffusion) THEN
cb0d108b91 Jean*0603        DO j=1-OLy,sNy+OLy
                0604         DO i=1-OLx,sNx+OLx
67a1e439d8 Patr*0605          df(i,j) = 0. _d 0
31566b6684 Alis*0606         ENDDO
                0607        ENDDO
                0608       ELSE
0e9223926c Jean*0609         CALL GAD_DIFF_R(
                0610      I                    bi,bj,k, maskup, KappaR, TracerN,
                0611      O                    df, myThid )
31566b6684 Alis*0612       ENDIF
                0613
6b7e107e49 Jean*0614       IF ( trUseDiffKr4 ) THEN
0e9223926c Jean*0615         CALL GAD_BIHARM_R(
                0616      I                    bi,bj,k, maskUp, diffKr4, TracerN,
                0617      U                    df, myThid )
6b7e107e49 Jean*0618       ENDIF
                0619
31566b6684 Alis*0620 #ifdef ALLOW_GMREDI
                0621 C-    GM/Redi flux in R
8e4e737494 Jean*0622       IF ( trUseGMRedi ) THEN
972c0130ec Jean*0623         CALL GMREDI_RTRANSPORT(
e9de1d7682 Jean*0624      I         trIdentity, bi, bj, k, iMin, iMax, jMin, jMax,
0e9223926c Jean*0625      I         maskUp, TracerN,
77c3815cb8 Jean*0626      U         df,
e9de1d7682 Jean*0627      I         myThid )
31566b6684 Alis*0628       ENDIF
                0629 #endif
                0630
cb0d108b91 Jean*0631       DO j=1-OLy,sNy+OLy
                0632        DO i=1-OLx,sNx+OLx
0e9223926c Jean*0633         fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)
31566b6684 Alis*0634        ENDDO
                0635       ENDDO
                0636
81c8d7b9aa Jean*0637 #ifdef ALLOW_DIAGNOSTICS
0845f1a203 Jean*0638 C-    Diagnostics of Tracer flux in R dir (mainly Diffusive terms),
81c8d7b9aa Jean*0639 C       Explicit terms only & excluding advective terms:
                0640       IF ( useDiagnostics .AND.
9897ae9d92 Jean*0641      &    (.NOT.implicitDiffusion .OR. trUseDiffKr4 .OR. trUseGMRedi)
                0642      &   ) THEN
81c8d7b9aa Jean*0643           diagName = 'DFrE'//diagSufx
cb0d108b91 Jean*0644           CALL DIAGNOSTICS_FILL( df, diagName, k,1, 2,bi,bj, myThid )
cf336ab6c5 Ryan*0645 #ifdef ALLOW_LAYERS
ee16a2cae4 Ryan*0646           IF ( useLayers ) THEN
50d8304171 Ryan*0647            CALL LAYERS_FILL(df,trIdentity,'DFR',k,1,2,bi,bj,myThid)
ee16a2cae4 Ryan*0648           ENDIF
cf336ab6c5 Ryan*0649 #endif /* ALLOW_LAYERS */
81c8d7b9aa Jean*0650       ENDIF
                0651 #endif
                0652
31566b6684 Alis*0653 #ifdef ALLOW_KPP
60fa91c126 Jean*0654 C-    Set non local KPP transport term (ghat):
8e4e737494 Jean*0655       IF ( trUseKPP .AND. k.GE.2 ) THEN
cb0d108b91 Jean*0656        DO j=1-OLy,sNy+OLy
                0657         DO i=1-OLx,sNx+OLx
67a1e439d8 Patr*0658          df(i,j) = 0. _d 0
31566b6684 Alis*0659         ENDDO
                0660        ENDDO
cb0d108b91 Jean*0661        IF (trIdentity.EQ.GAD_TEMPERATURE) THEN
31566b6684 Alis*0662         CALL KPP_TRANSPORT_T(
eb76443793 Jean*0663      I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
                0664      O           df,
                0665      I           myTime, myIter, myThid )
cb0d108b91 Jean*0666        ELSEIF (trIdentity.EQ.GAD_SALINITY) THEN
31566b6684 Alis*0667         CALL KPP_TRANSPORT_S(
eb76443793 Jean*0668      I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
                0669      O           df,
                0670      I           myTime, myIter, myThid )
5174b865b3 Mart*0671 #ifdef ALLOW_PTRACERS
cb0d108b91 Jean*0672        ELSEIF (trIdentity .GE. GAD_TR1) THEN
5174b865b3 Mart*0673         CALL KPP_TRANSPORT_PTR(
eb76443793 Jean*0674      I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
cb0d108b91 Jean*0675      I           trIdentity-GAD_TR1+1,
eb76443793 Jean*0676      O           df,
                0677      I           myTime, myIter, myThid )
5174b865b3 Mart*0678 #endif
31566b6684 Alis*0679        ELSE
cfa7d1c982 Jean*0680         WRITE(errorMessageUnit,*)
cb0d108b91 Jean*0681      &    'tracer identity =', trIdentity, ' is not valid => STOP'
cfa7d1c982 Jean*0682         STOP 'ABNORMAL END: S/R GAD_CALC_RHS: invalid tracer identity'
31566b6684 Alis*0683        ENDIF
cb0d108b91 Jean*0684        DO j=1-OLy,sNy+OLy
                0685         DO i=1-OLx,sNx+OLx
eaba2fd266 Jean*0686          fVerT(i,j,kUp) = fVerT(i,j,kUp)
                0687      &                  + df(i,j)*maskUp(i,j)*rhoFacF(k)
31566b6684 Alis*0688         ENDDO
                0689        ENDDO
cfa7d1c982 Jean*0690 #ifdef ALLOW_DIAGNOSTICS
                0691 C-    Diagnostics of Non-Local Tracer (vertical) flux
                0692        IF ( useDiagnostics ) THEN
                0693          diagName = 'KPPg'//diagSufx
                0694          CALL DIAGNOSTICS_FILL( df, diagName, k,1, 2,bi,bj, myThid )
                0695 C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
                0696 C        does it only if k=1 (never the case here)
                0697          IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)
cf336ab6c5 Ryan*0698 #ifdef ALLOW_LAYERS
ee16a2cae4 Ryan*0699          IF ( useLayers ) THEN
50d8304171 Ryan*0700           CALL LAYERS_FILL(df,trIdentity,'DFR',k,1,2,bi,bj,myThid)
ee16a2cae4 Ryan*0701          ENDIF
cf336ab6c5 Ryan*0702 #endif /* ALLOW_LAYERS */
cfa7d1c982 Jean*0703        ENDIF
31566b6684 Alis*0704 #endif
cfa7d1c982 Jean*0705       ENDIF
                0706 #endif /* ALLOW_KPP */
31566b6684 Alis*0707
998cbd211e Oliv*0708 #ifdef GAD_SMOLARKIEWICZ_HACK
2bd077d2cd Oliv*0709 coj   Hack to make redi (and everything else in this s/r) positive
                0710 coj   (see Smolarkiewicz MWR 1989 and Bott MWR 1989).
                0711 coj   Only works if 'down' is k+1 and k loop in thermodynamics is k=Nr,1,-1
46918f1b26 Jean*0712       IF ( trUseSmolHack ) THEN
cb0d108b91 Jean*0713        DO j=1-OLy,sNy+OLy-1
                0714         DO i=1-OLx,sNx+OLx-1
2bd077d2cd Oliv*0715 coj   Add outgoing fluxes
4e66ab0b67 Oliv*0716          outFlux=deltaTLev(k)*
998cbd211e Oliv*0717      &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
                0718      &   *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
                0719      &    *( MAX(0. _d 0,fZon(i+1,j)) + MAX(0. _d 0,-fZon(i,j))
                0720      &      +MAX(0. _d 0,fMer(i,j+1)) + MAX(0. _d 0,-fMer(i,j))
                0721      &      +MAX(0. _d 0,fVerT(i,j,kDown)*rkSign)
                0722      &      +MAX(0. _d 0,-fVerT(i,j,kUp)*rkSign)
                0723      &     )
3ad42fd705 Jean*0724          trac = localT(i,j)
2bd077d2cd Oliv*0725 coj   If they would reduce tracer by a fraction of more than
                0726 coj   SmolarkiewiczMaxFrac, scale them down
998cbd211e Oliv*0727          IF (outFlux.GT.0. _d 0 .AND.
                0728      &       outFlux.GT.SmolarkiewiczMaxFrac*trac) THEN
2bd077d2cd Oliv*0729 coj   If tracer is already negative, scale flux to zero
622cf2fb44 Jean*0730            gTrFac = MAX(0. _d 0,SmolarkiewiczMaxFrac*trac/outFlux)
2bd077d2cd Oliv*0731
622cf2fb44 Jean*0732            IF (fZon(i+1,j).GT.0. _d 0) fZon(i+1,j)=gTrFac*fZon(i+1,j)
                0733            IF (-fZon(i,j) .GT.0. _d 0) fZon(i,j)  =gTrFac*fZon(i,j)
                0734            IF (fMer(i,j+1).GT.0. _d 0) fMer(i,j+1)=gTrFac*fMer(i,j+1)
                0735            IF (-fMer(i,j) .GT.0. _d 0) fMer(i,j)  =gTrFac*fMer(i,j)
998cbd211e Oliv*0736            IF (-fVerT(i,j,kUp)*rkSign .GT.0. _d 0)
622cf2fb44 Jean*0737      &       fVerT(i,j,kUp)=gTrFac*fVerT(i,j,kUp)
2bd077d2cd Oliv*0738
                0739            IF (k.LT.Nr .AND. fVerT(i,j,kDown)*rkSign.GT.0. _d 0) THEN
                0740 coj   Down flux is special: it has already been applied in lower layer,
                0741 coj   so we have to readjust this.
20a3ab10b2 Oliv*0742 coj   undo down flux, ...
935a76deec Jean*0743              gTracer(i,j,k+1) = gTracer(i,j,k+1)
20a3ab10b2 Oliv*0744      &        +_recip_hFacC(i,j,k+1,bi,bj)*recip_drF(k+1)
2bd077d2cd Oliv*0745      &         *recip_rA(i,j,bi,bj)*recip_deepFac2C(k+1)
                0746      &         *recip_rhoFacC(k+1)
                0747      &         *( -fVerT(i,j,kDown)*rkSign )
                0748 coj   ... scale ...
622cf2fb44 Jean*0749              fVerT(i,j,kDown)=gTrFac*fVerT(i,j,kDown)
2bd077d2cd Oliv*0750 coj   ... and reapply
935a76deec Jean*0751              gTracer(i,j,k+1) = gTracer(i,j,k+1)
20a3ab10b2 Oliv*0752      &        +_recip_hFacC(i,j,k+1,bi,bj)*recip_drF(k+1)
2bd077d2cd Oliv*0753      &         *recip_rA(i,j,bi,bj)*recip_deepFac2C(k+1)
                0754      &         *recip_rhoFacC(k+1)
                0755      &         *( fVerT(i,j,kDown)*rkSign )
998cbd211e Oliv*0756            ENDIF
2bd077d2cd Oliv*0757
998cbd211e Oliv*0758          ENDIF
                0759         ENDDO
                0760        ENDDO
                0761       ENDIF
46918f1b26 Jean*0762 #endif /* GAD_SMOLARKIEWICZ_HACK */
998cbd211e Oliv*0763
31566b6684 Alis*0764 C--   Divergence of fluxes
eaba2fd266 Jean*0765 C     Anelastic: scale vertical fluxes by rhoFac and leave Horizontal fluxes unchanged
4d24da31b2 Mart*0766 C     for Stevens OBC: keep only vertical diffusive contribution on boundaries
cb0d108b91 Jean*0767       DO j=1-OLy,sNy+OLy-1
                0768        DO i=1-OLx,sNx+OLx-1
935a76deec Jean*0769         gTracer(i,j,k) = gTracer(i,j,k)
eaba2fd266 Jean*0770      &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
                0771      &   *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
4d24da31b2 Mart*0772      &   *( (fZon(i+1,j)-fZon(i,j))*maskInC(i,j,bi,bj)
                0773      &     +(fMer(i,j+1)-fMer(i,j))*maskInC(i,j,bi,bj)
bb6c554092 Jean*0774      &     +(fVerT(i,j,kDown)-fVerT(i,j,kUp))*rkSign
95826db035 Jean*0775      &     -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))*advFac
                0776      &                   +(vTrans(i,j+1)-vTrans(i,j))*advFac
bb6c554092 Jean*0777      &                   +(rTransKp1(i,j)-rTrans(i,j))*rAdvFac
4d24da31b2 Mart*0778      &                  )*maskInC(i,j,bi,bj)
31566b6684 Alis*0779      &    )
                0780        ENDDO
                0781       ENDDO
                0782
616d5d70aa Jean*0783 #ifdef ALLOW_DEBUG
188f784a21 Jean*0784       IF ( debugLevel .GE. debLevC
cb0d108b91 Jean*0785      &   .AND. trIdentity.EQ.GAD_TEMPERATURE
616d5d70aa Jean*0786      &   .AND. k.EQ.2 .AND. myIter.EQ.1+nIter0
                0787      &   .AND. nPx.EQ.1 .AND. nPy.EQ.1
                0788      &   .AND. useCubedSphereExchange ) THEN
                0789         CALL DEBUG_CS_CORNER_UV( ' fZon,fMer from GAD_CALC_RHS',
                0790      &             fZon,fMer, k, standardMessageUnit,bi,bj,myThid )
                0791       ENDIF
                0792 #endif /* ALLOW_DEBUG */
0845f1a203 Jean*0793
31566b6684 Alis*0794       RETURN
                0795       END