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