pkg/generic_advdiff/gad_advection.F

24016b3859 Alis*0001 #include "GAD_OPTIONS.h"
d8efe0aa23 Jean*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
24016b3859 Alis*0005
2b4c849245 Ed H*0006 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
28d97917ae Alis*0007 CBOP
                0008 C !ROUTINE: GAD_ADVECTION
                0009
                0010 C !INTERFACE: ==========================================================
6fd9490fbd Jean*0011       SUBROUTINE GAD_ADVECTION(
0d75a51072 Mart*0012      I     implicitAdvection, advectionSchArg, vertAdvecSchArg,
cb0d108b91 Jean*0013      I     trIdentity, deltaTLev,
09e49e8561 Jean*0014      I     uFld, vFld, wFld, tracer,
1b5fb69d21 Ed H*0015      O     gTracer,
                0016      I     bi,bj, myTime,myIter,myThid)
24016b3859 Alis*0017
28d97917ae Alis*0018 C !DESCRIPTION:
29fd21a3ae Jean*0019 C Calculates the tendency of a tracer due to advection.
28d97917ae Alis*0020 C It uses the multi-dimensional method given in \ref{sect:multiDimAdvection}
                0021 C and can only be used for the non-linear advection schemes such as the
0845f1a203 Jean*0022 C direct-space-time method and flux-limiters.
28d97917ae Alis*0023 C
                0024 C The algorithm is as follows:
                0025 C \begin{itemize}
                0026 C \item{$\theta^{(n+1/3)} = \theta^{(n)}
5ebd7636db Ed D*0027 C      - \Delta t (\partial_x (u\theta^{(n)}) - \theta^{(n)} \partial_x u$)}
28d97917ae Alis*0028 C \item{$\theta^{(n+2/3)} = \theta^{(n+1/3)}
5ebd7636db Ed D*0029 C      - \Delta t (\partial_y (v\theta^{(n+1/3)}) - \theta^{(n)} \partial_y v$)}
28d97917ae Alis*0030 C \item{$\theta^{(n+3/3)} = \theta^{(n+2/3)}
5ebd7636db Ed D*0031 C      - \Delta t (\partial_r (w\theta^{(n+2/3)}) - \theta^{(n)} \partial_r w$)}
28d97917ae Alis*0032 C \item{$G_\theta = ( \theta^{(n+3/3)} - \theta^{(n)} )/\Delta t$}
                0033 C \end{itemize}
                0034 C
29fd21a3ae Jean*0035 C The tendency (output) is over-written by this routine.
28d97917ae Alis*0036
                0037 C !USES: ===============================================================
                0038       IMPLICIT NONE
24016b3859 Alis*0039 #include "SIZE.h"
                0040 #include "EEPARAMS.h"
                0041 #include "PARAMS.h"
                0042 #include "GRID.h"
                0043 #include "GAD.h"
d8efe0aa23 Jean*0044 #ifdef ALLOW_AUTODIFF
0d75a51072 Mart*0045 # include "AUTODIFF_PARAMS.h"
7c50f07931 Mart*0046 # ifdef ALLOW_AUTODIFF_TAMC
                0047 #  include "tamc.h"
                0048 # endif
27cc6013c1 Patr*0049 # ifdef ALLOW_PTRACERS
                0050 #  include "PTRACERS_SIZE.h"
                0051 # endif
67a1e439d8 Patr*0052 #endif
9b89fcf692 antn*0053 #ifdef ALLOW_LAYERS
                0054 # include "LAYERS_P2SHARE.h"
                0055 #endif
e73e2206ff Dimi*0056 #ifdef ALLOW_EXCH2
9b89fcf692 antn*0057 # include "W2_EXCH2_SIZE.h"
                0058 # include "W2_EXCH2_TOPOLOGY.h"
e73e2206ff Dimi*0059 #endif /* ALLOW_EXCH2 */
24016b3859 Alis*0060
28d97917ae Alis*0061 C !INPUT PARAMETERS: ===================================================
1b5fb69d21 Ed H*0062 C  implicitAdvection :: implicit vertical advection (later on)
0d75a51072 Mart*0063 C  advectionSchArg   :: advection scheme to use (Horizontal plane)
                0064 C  vertAdvecSchArg   :: advection scheme to use (vertical direction)
cb0d108b91 Jean*0065 C  trIdentity        :: tracer identifier
09e49e8561 Jean*0066 C  uFld              :: Advection velocity field, zonal component
                0067 C  vFld              :: Advection velocity field, meridional component
                0068 C  wFld              :: Advection velocity field, vertical component
1b5fb69d21 Ed H*0069 C  tracer            :: tracer field
                0070 C  bi,bj             :: tile indices
                0071 C  myTime            :: current time
                0072 C  myIter            :: iteration number
                0073 C  myThid            :: thread number
6fd9490fbd Jean*0074       LOGICAL implicitAdvection
0d75a51072 Mart*0075       INTEGER advectionSchArg, vertAdvecSchArg
cb0d108b91 Jean*0076       INTEGER trIdentity
4e66ab0b67 Oliv*0077       _RL deltaTLev(Nr)
09e49e8561 Jean*0078       _RL uFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0079       _RL vFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0080       _RL wFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
cb0d108b91 Jean*0081       _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
6fd9490fbd Jean*0082       INTEGER bi,bj
24016b3859 Alis*0083       _RL myTime
                0084       INTEGER myIter
                0085       INTEGER myThid
                0086
28d97917ae Alis*0087 C !OUTPUT PARAMETERS: ==================================================
29fd21a3ae Jean*0088 C  gTracer           :: tendency array
935a76deec Jean*0089       _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
28d97917ae Alis*0090
20a156cdce Mart*0091 #ifdef DISABLE_MULTIDIM_ADVECTION
                0092       STOP 'GAD_ADVECTION is empty with DISABLE_MULTIDIM_ADVECTION'
                0093 #else /*  DISABLE_MULTIDIM_ADVECTION */
d077f4a702 Jean*0094 C !FUNCTIONS: ==========================================================
                0095 #ifdef ALLOW_DIAGNOSTICS
                0096       CHARACTER*4 GAD_DIAG_SUFX
                0097       EXTERNAL    GAD_DIAG_SUFX
                0098       LOGICAL  DIAGNOSTICS_IS_ON
                0099       EXTERNAL DIAGNOSTICS_IS_ON
                0100 #endif
                0101
28d97917ae Alis*0102 C !LOCAL VARIABLES: ====================================================
1b5fb69d21 Ed H*0103 C  maskUp        :: 2-D array for mask at W points
616d5d70aa Jean*0104 C  maskLocW      :: 2-D array for mask at West points
                0105 C  maskLocS      :: 2-D array for mask at South points
77cabae5d9 Jean*0106 C [iMin,iMax]Upd :: loop range to update tracer field
                0107 C [jMin,jMax]Upd :: loop range to update tracer field
1b5fb69d21 Ed H*0108 C  i,j,k         :: loop indices
0845f1a203 Jean*0109 C  kUp           :: index into 2 1/2D array, toggles between 1 and 2
                0110 C  kDown         :: index into 2 1/2D array, toggles between 2 and 1
0d75a51072 Mart*0111 C advectionScheme:: local copy of routine argument advectionSchArg
                0112 C vertAdvecScheme:: local copy of routine argument vertAdvecSchArg
1b5fb69d21 Ed H*0113 C  xA,yA         :: areas of X and Y face of tracer cells
                0114 C  uTrans,vTrans :: 2-D arrays of volume transports at U,V points
                0115 C  rTrans        :: 2-D arrays of volume transports at W points
d077f4a702 Jean*0116 C  rTransKp      :: vertical volume transport at interface k+1
8e4c181d69 Jean*0117 C  rTran3d       :: 3-D array of volume transport at W points
                0118 C  afr           :: 3-D array of vertical advective flux
77cabae5d9 Jean*0119 C  af            :: 2-D array for horizontal advective flux
616d5d70aa Jean*0120 C  afx           :: 2-D array for horizontal advective flux, x direction
                0121 C  afy           :: 2-D array for horizontal advective flux, y direction
1b5fb69d21 Ed H*0122 C  fVerT         :: 2 1/2D arrays for vertical advective flux
d077f4a702 Jean*0123 C  localTij      :: 2-D array, temporary local copy of tracer field
                0124 C  localT3d      :: 3-D array, temporary local copy of tracer field
1b5fb69d21 Ed H*0125 C  kp1Msk        :: flag (0,1) for over-riding mask for W levels
                0126 C  calc_fluxes_X :: logical to indicate to calculate fluxes in X dir
                0127 C  calc_fluxes_Y :: logical to indicate to calculate fluxes in Y dir
77cabae5d9 Jean*0128 C  interiorOnly  :: only update the interior of myTile, but not the edges
                0129 C  overlapOnly   :: only update the edges of myTile, but not the interior
da54d4c5af Jean*0130 C  npass         :: number of passes in multi-dimensional method
1b5fb69d21 Ed H*0131 C  ipass         :: number of the current pass being made
0845f1a203 Jean*0132 C  myTile        :: variables used to determine which cube face
e73e2206ff Dimi*0133 C  nCFace        :: owns a tile for cube grid runs using
                0134 C                :: multi-dim advection.
77cabae5d9 Jean*0135 C [N,S,E,W]_edge :: true if N,S,E,W edge of myTile is an Edge of the cube
0845f1a203 Jean*0136 c     _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
616d5d70aa Jean*0137       _RS maskLocW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0138       _RS maskLocS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
77cabae5d9 Jean*0139       INTEGER iMinUpd,iMaxUpd,jMinUpd,jMaxUpd
0845f1a203 Jean*0140       INTEGER i,j,k,kUp,kDown
0d75a51072 Mart*0141       INTEGER advectionScheme, vertAdvecScheme
24016b3859 Alis*0142       _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0143       _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0144       _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0145       _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0146       _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
d077f4a702 Jean*0147       _RL rTransKp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
20a156cdce Mart*0148 #ifndef ALLOW_AUTODIFF
                0149        LOGICAL usePPMvertAdv, usePQMvertAdv
8e4c181d69 Jean*0150       _RL rTran3d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0151       _RL afr     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
20a156cdce Mart*0152 #endif
77cabae5d9 Jean*0153       _RL af      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
616d5d70aa Jean*0154       _RL afx     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0155       _RL afy     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
24016b3859 Alis*0156       _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0157       _RL localTij(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
d077f4a702 Jean*0158       _RL localT3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0159 #ifdef GAD_MULTIDIM_COMPRESSIBLE
                0160       _RL tmpTrac
                0161       _RL localVol(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0162       _RL locVol3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0163 #endif
24016b3859 Alis*0164       _RL kp1Msk
616d5d70aa Jean*0165       LOGICAL calc_fluxes_X, calc_fluxes_Y, withSigns
77cabae5d9 Jean*0166       LOGICAL interiorOnly, overlapOnly
da54d4c5af Jean*0167       INTEGER npass, ipass
8996cf5a3c Jean*0168       INTEGER nCFace
77cabae5d9 Jean*0169       LOGICAL N_edge, S_edge, E_edge, W_edge
cd3a97f874 Jean*0170 #ifdef ALLOW_AUTODIFF_TAMC
0f757ba827 Patr*0171 C     msgBuf     :: Informational/error message buffer
edb6656069 Mart*0172 C     dkey       :: tape key (direction dependent)
0f757ba827 Patr*0173       CHARACTER*(MAX_LEN_MBUF) msgBuf
edb6656069 Mart*0174       INTEGER dkey
20a156cdce Mart*0175 # ifdef GAD_MULTIDIM_COMPRESSIBLE
edb6656069 Mart*0176 C     ijkey      :: tape key (depends on indices i,j, and direction)
20a156cdce Mart*0177       INTEGER ijkey
                0178 # endif
cd3a97f874 Jean*0179 #endif
8996cf5a3c Jean*0180 #ifdef ALLOW_EXCH2
                0181       INTEGER myTile
                0182 #endif
81c8d7b9aa Jean*0183 #ifdef ALLOW_DIAGNOSTICS
                0184       CHARACTER*8 diagName
78c3e29ef3 Jean*0185       CHARACTER*4 diagSufx
                0186       LOGICAL     doDiagAdvX, doDiagAdvY, doDiagAdvR
d077f4a702 Jean*0187 #endif /* ALLOW_DIAGNOSTICS */
28d97917ae Alis*0188 CEOP
24016b3859 Alis*0189
0d75a51072 Mart*0190 C     make local copies to be tampered with if necessary
                0191       advectionScheme = advectionSchArg
                0192       vertAdvecScheme = vertAdvecSchArg
67a1e439d8 Patr*0193 #ifdef ALLOW_AUTODIFF_TAMC
0d75a51072 Mart*0194       IF ( inAdMode .AND. useApproxAdvectionInAdMode ) THEN
                0195 C     In AD-mode, we change non-linear, potentially unstable AD advection
                0196 C     schemes to linear schemes with more stability. So far only DST3 with
                0197 C     flux limiting is replaced by DST3 without flux limiting, but any
                0198 C     combination is possible.
                0199        IF ( advectionSchArg.EQ.ENUM_DST3_FLUX_LIMIT )
                0200      &      advectionScheme = ENUM_DST3
                0201        IF ( vertAdvecSchArg.EQ.ENUM_DST3_FLUX_LIMIT )
                0202      &      vertAdvecScheme = ENUM_DST3
                0203 C     here is room for more advection schemes as this becomes necessary
                0204       ENDIF
67a1e439d8 Patr*0205 #endif /* ALLOW_AUTODIFF_TAMC */
                0206
81c8d7b9aa Jean*0207 #ifdef ALLOW_DIAGNOSTICS
78c3e29ef3 Jean*0208 C--   Set diagnostics flags and suffix for the current tracer
                0209       doDiagAdvX = .FALSE.
                0210       doDiagAdvY = .FALSE.
                0211       doDiagAdvR = .FALSE.
81c8d7b9aa Jean*0212       IF ( useDiagnostics ) THEN
20a156cdce Mart*0213        diagSufx = GAD_DIAG_SUFX( trIdentity, myThid )
                0214        diagName = 'ADVx'//diagSufx
                0215        doDiagAdvX = DIAGNOSTICS_IS_ON( diagName, myThid )
                0216        diagName = 'ADVy'//diagSufx
                0217        doDiagAdvY = DIAGNOSTICS_IS_ON( diagName, myThid )
                0218        diagName = 'ADVr'//diagSufx
                0219        doDiagAdvR = DIAGNOSTICS_IS_ON( diagName, myThid )
81c8d7b9aa Jean*0220       ENDIF
d077f4a702 Jean*0221 #endif /* ALLOW_DIAGNOSTICS */
81c8d7b9aa Jean*0222
24016b3859 Alis*0223 C--   Set up work arrays with valid (i.e. not NaN) values
                0224 C     These inital values do not alter the numerical results. They
                0225 C     just ensure that all memory references are to valid floating
                0226 C     point numbers. This prevents spurious hardware signals due to
                0227 C     uninitialised but inert locations.
                0228       DO j=1-OLy,sNy+OLy
                0229        DO i=1-OLx,sNx+OLx
09e49e8561 Jean*0230 C-    xA,yA,vFld,uTrans,vTrans are set over the full domain
                0231 C      => no need for extra initialisation
                0232 c       xA(i,j)      = 0. _d 0
                0233 c       yA(i,j)      = 0. _d 0
                0234 c       uTrans(i,j)  = 0. _d 0
                0235 c       vTrans(i,j)  = 0. _d 0
                0236 C-    rTransKp is set over the full domain: no need for extra initialisation
                0237 c       rTransKp(i,j)= 0. _d 0
                0238 C-    rTrans and fVerT need to be initialised to zero:
24016b3859 Alis*0239         rTrans(i,j)  = 0. _d 0
                0240         fVerT(i,j,1) = 0. _d 0
                0241         fVerT(i,j,2) = 0. _d 0
d8efe0aa23 Jean*0242 #ifdef ALLOW_AUTODIFF
d077f4a702 Jean*0243 # ifdef GAD_MULTIDIM_COMPRESSIBLE
                0244         localVol(i,j) = 0. _d 0
                0245 # endif
cdc9f269ae Patr*0246         localTij(i,j) = 0. _d 0
d8efe0aa23 Jean*0247 #endif /* ALLOW_AUTODIFF */
24016b3859 Alis*0248        ENDDO
                0249       ENDDO
                0250
77cabae5d9 Jean*0251 C--   Set tile-specific parameters for horizontal fluxes
                0252       IF (useCubedSphereExchange) THEN
da54d4c5af Jean*0253        npass  = 3
77cabae5d9 Jean*0254 #ifdef ALLOW_EXCH2
c424ee7cc7 Jean*0255        myTile = W2_myTileList(bi,bj)
77cabae5d9 Jean*0256        nCFace = exch2_myFace(myTile)
                0257        N_edge = exch2_isNedge(myTile).EQ.1
                0258        S_edge = exch2_isSedge(myTile).EQ.1
                0259        E_edge = exch2_isEedge(myTile).EQ.1
                0260        W_edge = exch2_isWedge(myTile).EQ.1
                0261 #else
                0262        nCFace = bi
                0263        N_edge = .TRUE.
                0264        S_edge = .TRUE.
                0265        E_edge = .TRUE.
                0266        W_edge = .TRUE.
                0267 #endif
                0268       ELSE
da54d4c5af Jean*0269        npass  = 2
                0270        nCFace = 0
77cabae5d9 Jean*0271        N_edge = .FALSE.
                0272        S_edge = .FALSE.
                0273        E_edge = .FALSE.
                0274        W_edge = .FALSE.
                0275       ENDIF
                0276
1574069d50 Mart*0277 #ifdef ALLOW_AUTODIFF_TAMC
                0278       IF ( npass.GT.maxcube ) THEN
                0279         WRITE(msgBuf,'(A,2(I3,A))') 'S/R GAD_ADVECTION: npass =',
                0280      &     npass, ' >', maxcube, ' = maxcube, ==> check "tamc.h"'
                0281         CALL PRINT_ERROR( msgBuf, myThid )
                0282         STOP 'ABNORMAL END: S/R GAD_ADVECTION'
                0283       ENDIF
20a156cdce Mart*0284 C     Define local tapes. Using local tapes makes it necessary to recompute
                0285 C     the forward routine (once) in AD-mode, but saves a lot of memory.
                0286 CADJ INIT loctape_gad_adv = COMMON, 1
                0287 CADJ INIT loctape_gad_adv_k = COMMON, Nr
                0288 CADJ INIT loctape_gad_adv_k_pass = COMMON, maxcube*Nr
                0289 # ifdef GAD_MULTIDIM_COMPRESSIBLE
                0290 CADJ INIT loctape_gad_adv_ijk_pass
                0291 CADJ &  = COMMON, maxcube*Nr*(sNx+2*OLx)*nSx*(sNy+2*OLy)*nSy
                0292 # endif
1574069d50 Mart*0293 #endif
24016b3859 Alis*0294 C--   Start of k loop for horizontal fluxes
                0295       DO k=1,Nr
                0296
                0297 C--   Get temporary terms used by tendency routines
20a156cdce Mart*0298        DO j=1-OLy,sNy+OLy
                0299         DO i=1-OLx,sNx+OLx
09e49e8561 Jean*0300          xA(i,j) = _dyG(i,j,bi,bj)*deepFacC(k)
                0301      &           *drF(k)*_hFacW(i,j,k,bi,bj)
                0302          yA(i,j) = _dxG(i,j,bi,bj)*deepFacC(k)
                0303      &           *drF(k)*_hFacS(i,j,k,bi,bj)
20a156cdce Mart*0304         ENDDO
09e49e8561 Jean*0305        ENDDO
                0306 C--   Calculate "volume transports" through tracer cell faces.
                0307 C     anelastic: scaled by rhoFacC (~ mass transport)
20a156cdce Mart*0308        DO j=1-OLy,sNy+OLy
                0309         DO i=1-OLx,sNx+OLx
09e49e8561 Jean*0310          uTrans(i,j) = uFld(i,j,k)*xA(i,j)*rhoFacC(k)
                0311          vTrans(i,j) = vFld(i,j,k)*yA(i,j)*rhoFacC(k)
20a156cdce Mart*0312         ENDDO
09e49e8561 Jean*0313        ENDDO
b4413f3c59 Jean*0314
616d5d70aa Jean*0315 C--   Make local copy of tracer array and mask West & South
20a156cdce Mart*0316        DO j=1-OLy,sNy+OLy
                0317         DO i=1-OLx,sNx+OLx
cb0d108b91 Jean*0318          localTij(i,j) = tracer(i,j,k,bi,bj)
d077f4a702 Jean*0319 #ifdef GAD_MULTIDIM_COMPRESSIBLE
                0320          localVol(i,j) = rA(i,j,bi,bj)*deepFac2C(k)
                0321      &                  *rhoFacC(k)*drF(k)*hFacC(i,j,k,bi,bj)
                0322      &                 + ( oneRS - maskC(i,j,k,bi,bj) )
                0323 #endif
b98672ba3a Jean*0324 #ifdef ALLOW_OBCS
                0325          maskLocW(i,j) = _maskW(i,j,k,bi,bj)*maskInW(i,j,bi,bj)
                0326          maskLocS(i,j) = _maskS(i,j,k,bi,bj)*maskInS(i,j,bi,bj)
                0327 #else /* ALLOW_OBCS */
                0328          maskLocW(i,j) = _maskW(i,j,k,bi,bj)
                0329          maskLocS(i,j) = _maskS(i,j,k,bi,bj)
                0330 #endif /* ALLOW_OBCS */
20a156cdce Mart*0331         ENDDO
24016b3859 Alis*0332        ENDDO
                0333
20a156cdce Mart*0334        IF (useCubedSphereExchange) THEN
616d5d70aa Jean*0335         withSigns = .FALSE.
0845f1a203 Jean*0336         CALL FILL_CS_CORNER_UV_RS(
616d5d70aa Jean*0337      &            withSigns, maskLocW,maskLocS, bi,bj, myThid )
20a156cdce Mart*0338        ENDIF
7711d7d079 Alis*0339
                0340 C--   Multiple passes for different directions on different tiles
e73e2206ff Dimi*0341 C--   For cube need one pass for each of red, green and blue axes.
20a156cdce Mart*0342        DO ipass=1,npass
67a1e439d8 Patr*0343 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0344         dkey = ipass + (k-1) * maxcube
20a156cdce Mart*0345 #endif
7711d7d079 Alis*0346
20a156cdce Mart*0347         interiorOnly = .FALSE.
                0348         overlapOnly  = .FALSE.
                0349         IF (useCubedSphereExchange) THEN
77cabae5d9 Jean*0350 C-    CubedSphere : pass 3 times, with partial update of local tracer field
20a156cdce Mart*0351          IF (ipass.EQ.1) THEN
                0352           overlapOnly  = MOD(nCFace,3).EQ.0
                0353           interiorOnly = MOD(nCFace,3).NE.0
                0354           calc_fluxes_X = nCFace.EQ.6 .OR. nCFace.EQ.1 .OR. nCFace.EQ.2
                0355           calc_fluxes_Y = nCFace.EQ.3 .OR. nCFace.EQ.4 .OR. nCFace.EQ.5
                0356          ELSEIF (ipass.EQ.2) THEN
                0357           overlapOnly  = MOD(nCFace,3).EQ.2
                0358           interiorOnly = MOD(nCFace,3).EQ.1
                0359           calc_fluxes_X = nCFace.EQ.2 .OR. nCFace.EQ.3 .OR. nCFace.EQ.4
                0360           calc_fluxes_Y = nCFace.EQ.5 .OR. nCFace.EQ.6 .OR. nCFace.EQ.1
                0361          ELSE
                0362           interiorOnly = .TRUE.
                0363           calc_fluxes_X = nCFace.EQ.5 .OR. nCFace.EQ.6
                0364           calc_fluxes_Y = nCFace.EQ.2 .OR. nCFace.EQ.3
                0365          ENDIF
                0366         ELSE
77cabae5d9 Jean*0367 C-    not CubedSphere
20a156cdce Mart*0368          calc_fluxes_X = MOD(ipass,2).EQ.1
                0369          calc_fluxes_Y = .NOT.calc_fluxes_X
                0370         ENDIF
0845f1a203 Jean*0371
616d5d70aa Jean*0372 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
7711d7d079 Alis*0373 C--   X direction
d077f4a702 Jean*0374
d8efe0aa23 Jean*0375 #ifdef ALLOW_AUTODIFF
d077f4a702 Jean*0376 C-     Always reset advective flux in X
cb0d108b91 Jean*0377         DO j=1-OLy,sNy+OLy
                0378          DO i=1-OLx,sNx+OLx
cdc9f269ae Patr*0379           af(i,j) = 0.
                0380          ENDDO
                0381         ENDDO
d8efe0aa23 Jean*0382 #endif /* ALLOW_AUTODIFF */
d077f4a702 Jean*0383
20a156cdce Mart*0384         IF (calc_fluxes_X) THEN
7711d7d079 Alis*0385
77cabae5d9 Jean*0386 C-     Do not compute fluxes if
0845f1a203 Jean*0387 C       a) needed in overlap only
77cabae5d9 Jean*0388 C   and b) the overlap of myTile are not cube-face Edges
20a156cdce Mart*0389          IF ( .NOT.overlapOnly .OR. N_edge .OR. S_edge ) THEN
77cabae5d9 Jean*0390
                0391 C-     Internal exchange for calculations in X
20a156cdce Mart*0392           IF ( overlapOnly ) THEN
                0393            CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
                0394      &                                localTij, bi,bj, myThid )
                0395           ENDIF
7711d7d079 Alis*0396
d077f4a702 Jean*0397 C-     Advective flux in X
d8efe0aa23 Jean*0398 #ifndef ALLOW_AUTODIFF
20a156cdce Mart*0399           DO j=1-OLy,sNy+OLy
                0400            DO i=1-OLx,sNx+OLx
                0401             af(i,j) = 0.
                0402            ENDDO
                0403           ENDDO
d8efe0aa23 Jean*0404 #else /* ALLOW_AUTODIFF */
20a156cdce Mart*0405 # ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0406 CADJ STORE localTij = loctape_gad_adv_k_pass, key = dkey, kind = isbyte
cdc9f269ae Patr*0407 # endif
d8efe0aa23 Jean*0408 #endif /* ALLOW_AUTODIFF */
67a1e439d8 Patr*0409
20a156cdce Mart*0410           IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
                0411      &         .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
                0412            CALL GAD_DST2U1_ADV_X(    bi,bj,k, advectionScheme, .TRUE.,
                0413      I              deltaTLev(k),uTrans,uFld(1-OLx,1-OLy,k), localTij,
                0414      O              af, myThid )
                0415           ELSEIF ( advectionScheme.EQ.ENUM_FLUX_LIMIT ) THEN
                0416            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
                0417      I              uTrans, uFld(1-OLx,1-OLy,k), maskLocW, localTij,
                0418      O              af, myThid )
                0419           ELSEIF ( advectionScheme.EQ.ENUM_DST3 ) THEN
                0420            CALL GAD_DST3_ADV_X(      bi,bj,k, .TRUE., deltaTLev(k),
                0421      I              uTrans, uFld(1-OLx,1-OLy,k), maskLocW, localTij,
                0422      O              af, myThid )
                0423           ELSEIF ( advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
                0424            CALL GAD_DST3FL_ADV_X(    bi,bj,k, .TRUE., deltaTLev(k),
                0425      I              uTrans, uFld(1-OLx,1-OLy,k), maskLocW, localTij,
                0426      O              af, myThid )
                0427           ELSEIF ( advectionScheme.EQ.ENUM_OS7MP ) THEN
                0428            CALL GAD_OS7MP_ADV_X(     bi,bj,k, .TRUE., deltaTLev(k),
                0429      I              uTrans, uFld(1-OLx,1-OLy,k), maskLocW, localTij,
                0430      O              af, myThid )
598aebfcee Mart*0431 #ifndef ALLOW_AUTODIFF
20a156cdce Mart*0432           ELSEIF ( advectionScheme.EQ.ENUM_PPM_NULL_LIMIT  .OR.
                0433      &             advectionScheme.EQ.ENUM_PPM_MONO_LIMIT  .OR.
                0434      &             advectionScheme.EQ.ENUM_PPM_WENO_LIMIT ) THEN
                0435            CALL GAD_PPM_ADV_X( advectionScheme, bi, bj, k , .TRUE.,
                0436      I              deltaTLev(k), uFld(1-OLx,1-OLy,k), uTrans, localTij,
                0437      O              af, myThid )
                0438           ELSEIF ( advectionScheme.EQ.ENUM_PQM_NULL_LIMIT  .OR.
                0439      &             advectionScheme.EQ.ENUM_PQM_MONO_LIMIT  .OR.
                0440      &             advectionScheme.EQ.ENUM_PQM_WENO_LIMIT ) THEN
                0441            CALL GAD_PQM_ADV_X( advectionScheme, bi, bj, k , .TRUE.,
                0442      I              deltaTLev(k), uFld(1-OLx,1-OLy,k), uTrans, localTij,
                0443      O              af, myThid )
8e4c181d69 Jean*0444 #endif /* ndef ALLOW_AUTODIFF */
20a156cdce Mart*0445           ELSE
                0446            STOP 'GAD_ADVECTION: adv. scheme incompatible with multi-dim'
                0447           ENDIF
77cabae5d9 Jean*0448
cb0d108b91 Jean*0449 #ifdef ALLOW_OBCS
20a156cdce Mart*0450           IF ( useOBCS ) THEN
cb0d108b91 Jean*0451 C-      replace advective flux with 1st order upwind scheme estimate
20a156cdce Mart*0452            CALL OBCS_U1_ADV_TRACER( .TRUE., trIdentity, bi, bj, k,
                0453      I                              maskLocW, uTrans, localTij,
                0454      U                              af, myThid )
                0455           ENDIF
cb0d108b91 Jean*0456 #endif /* ALLOW_OBCS */
                0457
77cabae5d9 Jean*0458 C-     Internal exchange for next calculations in Y
20a156cdce Mart*0459           IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
                0460            CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
                0461      &                                localTij, bi,bj, myThid )
                0462           ENDIF
da54d4c5af Jean*0463
                0464 C-     Advective flux in X : done
20a156cdce Mart*0465          ENDIF
77cabae5d9 Jean*0466
7448700841 Mart*0467 #if ( defined ALLOW_AUTODIFF_TAMC && defined NONLIN_FRSURF )
edb6656069 Mart*0468 CADJ STORE af = loctape_gad_adv_k_pass, key = dkey, kind = isbyte
20a156cdce Mart*0469 #endif
77cabae5d9 Jean*0470 C-     Update the local tracer field where needed:
bce2a21b32 Jean*0471 C      use "maksInC" to prevent updating tracer field in OB regions
77cabae5d9 Jean*0472
                0473 C      update in overlap-Only
20a156cdce Mart*0474          IF ( overlapOnly ) THEN
                0475           iMinUpd = 1-OLx+1
                0476           iMaxUpd = sNx+OLx-1
0845f1a203 Jean*0477 C- notes: these 2 lines below have no real effect (because recip_hFac=0
77cabae5d9 Jean*0478 C         in corner region) but safer to keep them.
20a156cdce Mart*0479           IF ( W_edge ) iMinUpd = 1
                0480           IF ( E_edge ) iMaxUpd = sNx
77cabae5d9 Jean*0481
20a156cdce Mart*0482           IF ( S_edge ) THEN
                0483 #if defined ALLOW_AUTODIFF_TAMC && defined GAD_MULTIDIM_COMPRESSIBLE
edb6656069 Mart*0484 CADJ STORE localTij = loctape_gad_adv_k_pass, key = dkey, kind = isbyte
20a156cdce Mart*0485 #endif
                0486            DO j=1-OLy,0
                0487             DO i=iMinUpd,iMaxUpd
d077f4a702 Jean*0488 #ifdef GAD_MULTIDIM_COMPRESSIBLE
20a156cdce Mart*0489 # ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0490              ijkey = i + ((j-1) + (dkey-1)*(sNy+2*OLy))*(sNx+2*OLx)
20a156cdce Mart*0491 CADJ STORE localVol(i,j) = loctape_gad_adv_ijk_pass, key=ijkey
                0492 CADJ &     , kind = isbyte
                0493 # endif
                0494              tmpTrac = localTij(i,j)*localVol(i,j)
                0495      &            -deltaTLev(k)*( af(i+1,j) - af(i,j) )
                0496      &                         *maskInC(i,j,bi,bj)
                0497              localVol(i,j) = localVol(i,j)
                0498      &            -deltaTLev(k)*( uTrans(i+1,j) - uTrans(i,j) )
                0499      &                         *maskInC(i,j,bi,bj)
                0500              localTij(i,j) = tmpTrac/localVol(i,j)
d077f4a702 Jean*0501 #else /* GAD_MULTIDIM_COMPRESSIBLE */
20a156cdce Mart*0502              localTij(i,j) = localTij(i,j)
                0503      &           -deltaTLev(k)*recip_rhoFacC(k)
                0504      &            *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
                0505      &            *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
                0506      &            *( af(i+1,j)-af(i,j)
                0507      &              -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
                0508      &             )*maskInC(i,j,bi,bj)
d077f4a702 Jean*0509 #endif /* GAD_MULTIDIM_COMPRESSIBLE */
20a156cdce Mart*0510             ENDDO
                0511            ENDDO
                0512           ENDIF
                0513           IF ( N_edge ) THEN
                0514 #if defined ALLOW_AUTODIFF_TAMC && defined GAD_MULTIDIM_COMPRESSIBLE
edb6656069 Mart*0515 CADJ STORE localTij = loctape_gad_adv_k_pass, key = dkey, kind = isbyte
20a156cdce Mart*0516 #endif
                0517            DO j=sNy+1,sNy+OLy
                0518             DO i=iMinUpd,iMaxUpd
d077f4a702 Jean*0519 #ifdef GAD_MULTIDIM_COMPRESSIBLE
20a156cdce Mart*0520 # ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0521              ijkey = i + ((j-1) + (dkey-1)*(sNy+2*OLy))*(sNx+2*OLx)
20a156cdce Mart*0522 CADJ STORE localVol(i,j) = loctape_gad_adv_ijk_pass, key=ijkey
                0523 CADJ &     , kind = isbyte
                0524 # endif
                0525              tmpTrac = localTij(i,j)*localVol(i,j)
                0526      &            -deltaTLev(k)*( af(i+1,j) - af(i,j) )
                0527      &                         *maskInC(i,j,bi,bj)
                0528              localVol(i,j) = localVol(i,j)
                0529      &            -deltaTLev(k)*( uTrans(i+1,j) - uTrans(i,j) )
                0530      &                         *maskInC(i,j,bi,bj)
                0531              localTij(i,j) = tmpTrac/localVol(i,j)
d077f4a702 Jean*0532 #else /* GAD_MULTIDIM_COMPRESSIBLE */
20a156cdce Mart*0533              localTij(i,j) = localTij(i,j)
                0534      &           -deltaTLev(k)*recip_rhoFacC(k)
                0535      &            *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
                0536      &            *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
                0537      &            *( af(i+1,j)-af(i,j)
                0538      &              -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
                0539      &             )*maskInC(i,j,bi,bj)
d077f4a702 Jean*0540 #endif /* GAD_MULTIDIM_COMPRESSIBLE */
20a156cdce Mart*0541             ENDDO
                0542            ENDDO
                0543           ENDIF
77cabae5d9 Jean*0544
20a156cdce Mart*0545          ELSE
                0546 #if defined ALLOW_AUTODIFF_TAMC && defined GAD_MULTIDIM_COMPRESSIBLE
edb6656069 Mart*0547 CADJ STORE localTij = loctape_gad_adv_k_pass, key = dkey, kind = isbyte
20a156cdce Mart*0548 #endif
77cabae5d9 Jean*0549 C      do not only update the overlap
20a156cdce Mart*0550           jMinUpd = 1-OLy
                0551           jMaxUpd = sNy+OLy
                0552           IF ( interiorOnly .AND. S_edge ) jMinUpd = 1
                0553           IF ( interiorOnly .AND. N_edge ) jMaxUpd = sNy
                0554           DO j=jMinUpd,jMaxUpd
                0555            DO i=1-OLx+1,sNx+OLx-1
d077f4a702 Jean*0556 #ifdef GAD_MULTIDIM_COMPRESSIBLE
20a156cdce Mart*0557 # ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0558              ijkey = i + ((j-1) + (dkey-1)*(sNy+2*OLy))*(sNx+2*OLx)
20a156cdce Mart*0559 CADJ STORE localVol(i,j) = loctape_gad_adv_ijk_pass, key=ijkey
                0560 CADJ &     , kind = isbyte
                0561 # endif
                0562             tmpTrac = localTij(i,j)*localVol(i,j)
                0563      &           -deltaTLev(k)*( af(i+1,j) - af(i,j) )
                0564      &                        *maskInC(i,j,bi,bj)
                0565             localVol(i,j) = localVol(i,j)
                0566      &           -deltaTLev(k)*( uTrans(i+1,j) - uTrans(i,j) )
                0567      &                        *maskInC(i,j,bi,bj)
                0568             localTij(i,j) = tmpTrac/localVol(i,j)
d077f4a702 Jean*0569 #else /* GAD_MULTIDIM_COMPRESSIBLE */
20a156cdce Mart*0570             localTij(i,j) = localTij(i,j)
                0571      &          -deltaTLev(k)*recip_rhoFacC(k)
                0572      &           *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
                0573      &           *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
                0574      &           *( af(i+1,j)-af(i,j)
                0575      &             -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
                0576      &            )*maskInC(i,j,bi,bj)
d077f4a702 Jean*0577 #endif /* GAD_MULTIDIM_COMPRESSIBLE */
20a156cdce Mart*0578            ENDDO
                0579           ENDDO
77cabae5d9 Jean*0580 C-      keep advective flux (for diagnostics)
20a156cdce Mart*0581           DO j=1-OLy,sNy+OLy
                0582            DO i=1-OLx,sNx+OLx
                0583             afx(i,j) = af(i,j)
                0584            ENDDO
                0585           ENDDO
24016b3859 Alis*0586
77cabae5d9 Jean*0587 C-     end if/else update overlap-Only
20a156cdce Mart*0588          ENDIF
0845f1a203 Jean*0589
7711d7d079 Alis*0590 C--   End of X direction
20a156cdce Mart*0591         ENDIF
7711d7d079 Alis*0592
616d5d70aa Jean*0593 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
7711d7d079 Alis*0594 C--   Y direction
d077f4a702 Jean*0595
d8efe0aa23 Jean*0596 #ifdef ALLOW_AUTODIFF
d077f4a702 Jean*0597 C-     Always reset advective flux in Y
cb0d108b91 Jean*0598         DO j=1-OLy,sNy+OLy
                0599          DO i=1-OLx,sNx+OLx
cdc9f269ae Patr*0600           af(i,j) = 0.
                0601          ENDDO
                0602         ENDDO
d8efe0aa23 Jean*0603 #endif /* ALLOW_AUTODIFF */
d077f4a702 Jean*0604
20a156cdce Mart*0605         IF (calc_fluxes_Y) THEN
7711d7d079 Alis*0606
77cabae5d9 Jean*0607 C-     Do not compute fluxes if
                0608 C       a) needed in overlap only
                0609 C   and b) the overlap of myTile are not cube-face edges
20a156cdce Mart*0610          IF ( .NOT.overlapOnly .OR. E_edge .OR. W_edge ) THEN
77cabae5d9 Jean*0611
                0612 C-     Internal exchange for calculations in Y
20a156cdce Mart*0613           IF ( overlapOnly ) THEN
                0614            CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
                0615      &                                localTij, bi,bj, myThid )
                0616           ENDIF
7711d7d079 Alis*0617
77cabae5d9 Jean*0618 C-     Advective flux in Y
d8efe0aa23 Jean*0619 #ifndef ALLOW_AUTODIFF
20a156cdce Mart*0620           DO j=1-OLy,sNy+OLy
                0621            DO i=1-OLx,sNx+OLx
                0622             af(i,j) = 0.
                0623            ENDDO
                0624           ENDDO
d8efe0aa23 Jean*0625 #else /* ALLOW_AUTODIFF */
20a156cdce Mart*0626 # ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0627 CADJ STORE localTij = loctape_gad_adv_k_pass, key = dkey, kind = isbyte
20a156cdce Mart*0628 # endif
d8efe0aa23 Jean*0629 #endif /* ALLOW_AUTODIFF */
67a1e439d8 Patr*0630
20a156cdce Mart*0631           IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
                0632      &         .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
                0633            CALL GAD_DST2U1_ADV_Y(    bi,bj,k, advectionScheme, .TRUE.,
                0634      I              deltaTLev(k),vTrans,vFld(1-OLx,1-OLy,k), localTij,
                0635      O              af, myThid )
                0636           ELSEIF ( advectionScheme.EQ.ENUM_FLUX_LIMIT ) THEN
                0637            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
                0638      I              vTrans, vFld(1-OLx,1-OLy,k), maskLocS, localTij,
                0639      O              af, myThid )
                0640           ELSEIF ( advectionScheme.EQ.ENUM_DST3 ) THEN
                0641            CALL GAD_DST3_ADV_Y(      bi,bj,k, .TRUE., deltaTLev(k),
                0642      I              vTrans, vFld(1-OLx,1-OLy,k), maskLocS, localTij,
                0643      O              af, myThid )
                0644           ELSEIF ( advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
                0645            CALL GAD_DST3FL_ADV_Y(    bi,bj,k, .TRUE., deltaTLev(k),
                0646      I              vTrans, vFld(1-OLx,1-OLy,k), maskLocS, localTij,
                0647      O              af, myThid )
                0648           ELSEIF ( advectionScheme.EQ.ENUM_OS7MP ) THEN
                0649            CALL GAD_OS7MP_ADV_Y(     bi,bj,k, .TRUE., deltaTLev(k),
                0650      I              vTrans, vFld(1-OLx,1-OLy,k), maskLocS, localTij,
                0651      O              af, myThid )
598aebfcee Mart*0652 #ifndef ALLOW_AUTODIFF
20a156cdce Mart*0653           ELSEIF ( advectionScheme.EQ.ENUM_PPM_NULL_LIMIT  .OR.
                0654      &             advectionScheme.EQ.ENUM_PPM_MONO_LIMIT  .OR.
                0655      &             advectionScheme.EQ.ENUM_PPM_WENO_LIMIT ) THEN
                0656            CALL GAD_PPM_ADV_Y( advectionScheme, bi, bj, k , .TRUE.,
                0657      I              deltaTLev(k), vFld(1-OLX,1-OLy,k), vTrans, localTij,
                0658      O              af, myThid )
                0659           ELSEIF ( advectionScheme.EQ.ENUM_PQM_NULL_LIMIT  .OR.
                0660      &             advectionScheme.EQ.ENUM_PQM_MONO_LIMIT  .OR.
                0661      &             advectionScheme.EQ.ENUM_PQM_WENO_LIMIT ) THEN
                0662            CALL GAD_PQM_ADV_Y( advectionScheme, bi, bj, k , .TRUE.,
                0663      I              deltaTLev(k), vFld(1-OLX,1-OLy,k), vTrans, localTij,
                0664      O              af, myThid )
8e4c181d69 Jean*0665 #endif /* ndef ALLOW_AUTODIFF */
20a156cdce Mart*0666           ELSE
                0667            STOP 'GAD_ADVECTION: adv. scheme incompatible with mutli-dim'
                0668           ENDIF
77cabae5d9 Jean*0669
cb0d108b91 Jean*0670 #ifdef ALLOW_OBCS
20a156cdce Mart*0671           IF ( useOBCS ) THEN
cb0d108b91 Jean*0672 C-      replace advective flux with 1st order upwind scheme estimate
20a156cdce Mart*0673            CALL OBCS_U1_ADV_TRACER( .FALSE., trIdentity, bi, bj, k,
                0674      I                              maskLocS, vTrans, localTij,
                0675      U                              af, myThid )
                0676           ENDIF
cb0d108b91 Jean*0677 #endif /* ALLOW_OBCS */
                0678
77cabae5d9 Jean*0679 C-     Internal exchange for next calculations in X
20a156cdce Mart*0680           IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
                0681            CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
                0682      &                                localTij, bi,bj, myThid )
                0683           ENDIF
da54d4c5af Jean*0684
                0685 C-     Advective flux in Y : done
20a156cdce Mart*0686          ENDIF
77cabae5d9 Jean*0687
7448700841 Mart*0688 #if ( defined ALLOW_AUTODIFF_TAMC && defined NONLIN_FRSURF )
edb6656069 Mart*0689 CADJ STORE af = loctape_gad_adv_k_pass, key = dkey, kind = isbyte
20a156cdce Mart*0690 #endif
77cabae5d9 Jean*0691 C-     Update the local tracer field where needed:
bce2a21b32 Jean*0692 C      use "maksInC" to prevent updating tracer field in OB regions
77cabae5d9 Jean*0693
                0694 C      update in overlap-Only
20a156cdce Mart*0695          IF ( overlapOnly ) THEN
                0696           jMinUpd = 1-OLy+1
                0697           jMaxUpd = sNy+OLy-1
0845f1a203 Jean*0698 C- notes: these 2 lines below have no real effect (because recip_hFac=0
77cabae5d9 Jean*0699 C         in corner region) but safer to keep them.
20a156cdce Mart*0700           IF ( S_edge ) jMinUpd = 1
                0701           IF ( N_edge ) jMaxUpd = sNy
77cabae5d9 Jean*0702
20a156cdce Mart*0703           IF ( W_edge ) THEN
                0704 #if defined ALLOW_AUTODIFF_TAMC && defined GAD_MULTIDIM_COMPRESSIBLE
edb6656069 Mart*0705 CADJ STORE localTij = loctape_gad_adv_k_pass, key = dkey, kind = isbyte
20a156cdce Mart*0706 #endif
                0707            DO j=jMinUpd,jMaxUpd
                0708             DO i=1-OLx,0
d077f4a702 Jean*0709 #ifdef GAD_MULTIDIM_COMPRESSIBLE
20a156cdce Mart*0710 # ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0711              ijkey = i + ((j-1) + (dkey-1)*(sNy+2*OLy))*(sNx+2*OLx)
20a156cdce Mart*0712 CADJ STORE localVol(i,j) = loctape_gad_adv_ijk_pass, key = ijkey
                0713 CADJ &     , kind = isbyte
                0714 # endif
                0715              tmpTrac = localTij(i,j)*localVol(i,j)
                0716      &            -deltaTLev(k)*( af(i,j+1) - af(i,j) )
                0717      &                         *maskInC(i,j,bi,bj)
                0718              localVol(i,j) = localVol(i,j)
                0719      &            -deltaTLev(k)*( vTrans(i,j+1) - vTrans(i,j) )
                0720      &                         *maskInC(i,j,bi,bj)
                0721              localTij(i,j) = tmpTrac/localVol(i,j)
d077f4a702 Jean*0722 #else /* GAD_MULTIDIM_COMPRESSIBLE */
20a156cdce Mart*0723              localTij(i,j) = localTij(i,j)
                0724      &           -deltaTLev(k)*recip_rhoFacC(k)
                0725      &            *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
                0726      &            *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
                0727      &            *( af(i,j+1)-af(i,j)
                0728      &              -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
                0729      &             )*maskInC(i,j,bi,bj)
d077f4a702 Jean*0730 #endif /* GAD_MULTIDIM_COMPRESSIBLE */
20a156cdce Mart*0731             ENDDO
                0732            ENDDO
                0733           ENDIF
                0734           IF ( E_edge ) THEN
                0735 #if defined ALLOW_AUTODIFF_TAMC && defined GAD_MULTIDIM_COMPRESSIBLE
edb6656069 Mart*0736 CADJ STORE localTij = loctape_gad_adv_k_pass, key = dkey
20a156cdce Mart*0737 CADJ &     , kind = isbyte
                0738 #endif
                0739            DO j=jMinUpd,jMaxUpd
                0740             DO i=sNx+1,sNx+OLx
d077f4a702 Jean*0741 #ifdef GAD_MULTIDIM_COMPRESSIBLE
20a156cdce Mart*0742 # ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0743              ijkey = i + ((j-1) + (dkey-1)*(sNy+2*OLy))*(sNx+2*OLx)
20a156cdce Mart*0744 CADJ STORE localVol(i,j) = loctape_gad_adv_ijk_pass, key=ijkey
                0745 CADJ &     , kind = isbyte
                0746 # endif
                0747              tmpTrac = localTij(i,j)*localVol(i,j)
                0748      &            -deltaTLev(k)*( af(i,j+1) - af(i,j) )
                0749      &                         *maskInC(i,j,bi,bj)
                0750              localVol(i,j) = localVol(i,j)
                0751      &            -deltaTLev(k)*( vTrans(i,j+1) - vTrans(i,j) )
                0752      &                         *maskInC(i,j,bi,bj)
                0753              localTij(i,j) = tmpTrac/localVol(i,j)
d077f4a702 Jean*0754 #else /* GAD_MULTIDIM_COMPRESSIBLE */
20a156cdce Mart*0755              localTij(i,j) = localTij(i,j)
                0756      &           -deltaTLev(k)*recip_rhoFacC(k)
                0757      &            *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
                0758      &            *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
                0759      &            *( af(i,j+1)-af(i,j)
                0760      &              -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
                0761      &             )*maskInC(i,j,bi,bj)
d077f4a702 Jean*0762 #endif /* GAD_MULTIDIM_COMPRESSIBLE */
20a156cdce Mart*0763             ENDDO
                0764            ENDDO
                0765           ENDIF
77cabae5d9 Jean*0766
20a156cdce Mart*0767          ELSE
                0768 #if defined ALLOW_AUTODIFF_TAMC && defined GAD_MULTIDIM_COMPRESSIBLE
edb6656069 Mart*0769 CADJ STORE localTij = loctape_gad_adv_k_pass, key = dkey
20a156cdce Mart*0770 CADJ &     , kind = isbyte
                0771 #endif
77cabae5d9 Jean*0772 C      do not only update the overlap
20a156cdce Mart*0773           iMinUpd = 1-OLx
                0774           iMaxUpd = sNx+OLx
                0775           IF ( interiorOnly .AND. W_edge ) iMinUpd = 1
                0776           IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx
                0777           DO j=1-OLy+1,sNy+OLy-1
                0778            DO i=iMinUpd,iMaxUpd
d077f4a702 Jean*0779 #ifdef GAD_MULTIDIM_COMPRESSIBLE
20a156cdce Mart*0780 # ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0781              ijkey = i + ((j-1) + (dkey-1)*(sNy+2*OLy))*(sNx+2*OLx)
20a156cdce Mart*0782 CADJ STORE localVol(i,j) = loctape_gad_adv_ijk_pass, key=ijkey
                0783 CADJ &     , kind = isbyte
                0784 # endif
                0785             tmpTrac = localTij(i,j)*localVol(i,j)
                0786      &           -deltaTLev(k)*( af(i,j+1) - af(i,j) )
                0787      &                        *maskInC(i,j,bi,bj)
                0788             localVol(i,j) = localVol(i,j)
                0789      &           -deltaTLev(k)*( vTrans(i,j+1) - vTrans(i,j) )
                0790      &                        *maskInC(i,j,bi,bj)
                0791             localTij(i,j) = tmpTrac/localVol(i,j)
d077f4a702 Jean*0792 #else /* GAD_MULTIDIM_COMPRESSIBLE */
20a156cdce Mart*0793             localTij(i,j) = localTij(i,j)
                0794      &          -deltaTLev(k)*recip_rhoFacC(k)
                0795      &           *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
                0796      &           *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
                0797      &           *( af(i,j+1)-af(i,j)
                0798      &             -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
                0799      &            )*maskInC(i,j,bi,bj)
d077f4a702 Jean*0800 #endif /* GAD_MULTIDIM_COMPRESSIBLE */
20a156cdce Mart*0801            ENDDO
                0802           ENDDO
77cabae5d9 Jean*0803 C-      keep advective flux (for diagnostics)
20a156cdce Mart*0804           DO j=1-OLy,sNy+OLy
                0805            DO i=1-OLx,sNx+OLx
                0806             afy(i,j) = af(i,j)
                0807            ENDDO
                0808           ENDDO
7711d7d079 Alis*0809
77cabae5d9 Jean*0810 C      end if/else update overlap-Only
20a156cdce Mart*0811          ENDIF
77cabae5d9 Jean*0812
7711d7d079 Alis*0813 C--   End of Y direction
20a156cdce Mart*0814         ENDIF
7711d7d079 Alis*0815
                0816 C--   End of ipass loop
20a156cdce Mart*0817        ENDDO
24016b3859 Alis*0818
20a156cdce Mart*0819        IF ( implicitAdvection ) THEN
75abb052f1 Jean*0820 C-    explicit advection is done ; store tendency in gTracer:
d077f4a702 Jean*0821 #ifdef GAD_MULTIDIM_COMPRESSIBLE
                0822         STOP 'GAD_ADVECTION: missing code for implicitAdvection'
                0823 #endif /* GAD_MULTIDIM_COMPRESSIBLE */
cb0d108b91 Jean*0824         DO j=1-OLy,sNy+OLy
                0825          DO i=1-OLx,sNx+OLx
935a76deec Jean*0826           gTracer(i,j,k) =
                0827      &     ( localTij(i,j) - tracer(i,j,k,bi,bj) )/deltaTLev(k)
75abb052f1 Jean*0828          ENDDO
                0829         ENDDO
20a156cdce Mart*0830        ELSE
75abb052f1 Jean*0831 C-    horizontal advection done; store intermediate result in 3D array:
cb0d108b91 Jean*0832         DO j=1-OLy,sNy+OLy
                0833          DO i=1-OLx,sNx+OLx
d077f4a702 Jean*0834 #ifdef GAD_MULTIDIM_COMPRESSIBLE
                0835           locVol3d(i,j,k) = localVol(i,j)
                0836 #endif /* GAD_MULTIDIM_COMPRESSIBLE */
                0837           localT3d(i,j,k) = localTij(i,j)
eaba2fd266 Jean*0838          ENDDO
75abb052f1 Jean*0839         ENDDO
20a156cdce Mart*0840        ENDIF
75abb052f1 Jean*0841
81c8d7b9aa Jean*0842 #ifdef ALLOW_DIAGNOSTICS
20a156cdce Mart*0843        IF ( doDiagAdvX ) THEN
                0844         diagName = 'ADVx'//diagSufx
                0845         CALL DIAGNOSTICS_FILL( afx, diagName, k,1, 2,bi,bj, myThid )
                0846        ENDIF
                0847        IF ( doDiagAdvY ) THEN
                0848         diagName = 'ADVy'//diagSufx
                0849         CALL DIAGNOSTICS_FILL( afy, diagName, k,1, 2,bi,bj, myThid )
                0850        ENDIF
50d8304171 Ryan*0851 #ifdef ALLOW_LAYERS
9b89fcf692 antn*0852        IF ( layers_useThermo ) THEN
20a156cdce Mart*0853         CALL LAYERS_FILL(afx,trIdentity,'AFX',k,1,2,bi,bj,myThid)
                0854         CALL LAYERS_FILL(afy,trIdentity,'AFY',k,1,2,bi,bj,myThid)
                0855        ENDIF
50d8304171 Ryan*0856 #endif /* ALLOW_LAYERS */
d077f4a702 Jean*0857 #endif /* ALLOW_DIAGNOSTICS */
81c8d7b9aa Jean*0858
616d5d70aa Jean*0859 #ifdef ALLOW_DEBUG
20a156cdce Mart*0860        IF ( debugLevel .GE. debLevC
cb0d108b91 Jean*0861      &   .AND. trIdentity.EQ.GAD_TEMPERATURE
77cabae5d9 Jean*0862      &   .AND. k.LE.3 .AND. myIter.EQ.1+nIter0
616d5d70aa Jean*0863      &   .AND. nPx.EQ.1 .AND. nPy.EQ.1
                0864      &   .AND. useCubedSphereExchange ) THEN
                0865         CALL DEBUG_CS_CORNER_UV( ' afx,afy from GAD_ADVECTION',
                0866      &             afx,afy, k, standardMessageUnit,bi,bj,myThid )
20a156cdce Mart*0867        ENDIF
616d5d70aa Jean*0868 #endif /* ALLOW_DEBUG */
                0869
24016b3859 Alis*0870 C--   End of K loop for horizontal fluxes
                0871       ENDDO
                0872
616d5d70aa Jean*0873 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0874
75abb052f1 Jean*0875       IF ( .NOT.implicitAdvection ) THEN
8e4c181d69 Jean*0876
                0877 #ifndef ALLOW_AUTODIFF
20a156cdce Mart*0878        usePPMvertAdv = vertAdvecScheme.EQ.ENUM_PPM_NULL_LIMIT
                0879      &            .OR. vertAdvecScheme.EQ.ENUM_PPM_MONO_LIMIT
                0880      &            .OR. vertAdvecScheme.EQ.ENUM_PPM_WENO_LIMIT
                0881        usePQMvertAdv = vertAdvecScheme.EQ.ENUM_PQM_NULL_LIMIT
                0882      &            .OR. vertAdvecScheme.EQ.ENUM_PQM_MONO_LIMIT
                0883      &            .OR. vertAdvecScheme.EQ.ENUM_PQM_WENO_LIMIT
                0884        IF ( usePPMvertAdv .OR. usePQMvertAdv ) THEN
                0885 C-- PPM-style or PQM-style vertical advection
8e4c181d69 Jean*0886         DO k=1,Nr
                0887          IF (k.EQ.1) THEN
                0888 C-- vertical transport:  surface
                0889           DO j=1-OLy,sNy+OLy
                0890            DO i=1-OLx,sNx+OLx
                0891             rTran3d(i,j,k) = 0. _d 0
                0892            ENDDO
                0893           ENDDO
                0894          ELSE
                0895 C-- vertical transport: interior
                0896           DO j=1-OLy,sNy+OLy
                0897            DO i=1-OLx,sNx+OLx
                0898             rTran3d(i,j,k) = wFld(i,j,k)*rA(i,j,bi,bj)
                0899      &                      *deepFac2F(k)*rhoFacF(k)
                0900      &                      *maskC(i,j,k-1,bi,bj)
                0901            ENDDO
                0902           ENDDO
                0903          ENDIF
                0904         ENDDO
20a156cdce Mart*0905         IF ( usePPMvertAdv ) THEN
8e4c181d69 Jean*0906 C-- calc. PPM vertical flux data
20a156cdce Mart*0907          CALL GAD_PPM_ADV_R( vertAdvecScheme, bi, bj,
                0908      I                deltaTLev, wFld, rTran3d, localT3d,
                0909      O                afr, myThid )
8e4c181d69 Jean*0910         ENDIF
20a156cdce Mart*0911         IF ( usePQMvertAdv ) THEN
8e4c181d69 Jean*0912 C-- calc. PQM vertical flux data
20a156cdce Mart*0913          CALL GAD_PQM_ADV_R( vertAdvecScheme, bi, bj,
                0914      I                deltaTLev, wFld, rTran3d, localT3d,
                0915      O                afr, myThid )
                0916         ENDIF
                0917 C-- end PPM / PQM style vertical advection
8e4c181d69 Jean*0918        ENDIF
                0919 #endif /* ndef ALLOW_AUTODIFF */
                0920
20a156cdce Mart*0921 #ifdef ALLOW_AUTODIFF_TAMC
                0922 CADJ STORE localT3d = loctape_gad_adv, key=1, kind = isbyte
                0923
                0924 # ifdef GAD_MULTIDIM_COMPRESSIBLE
                0925 CADJ STORE locVol3d = loctape_gad_adv, key=1, kind = isbyte
                0926 # endif /* GAD_MULTIDIM_COMPRESSIBLE */
                0927 #endif
                0928 #ifdef ALLOW_AUTODIFF
                0929        DO j=1-OLy,sNy+OLy
                0930         DO i=1-OLx,sNx+OLx
                0931          fVerT(i,j,1) = 0.
                0932          fVerT(i,j,2) = 0.
                0933         ENDDO
                0934        ENDDO
                0935 #endif
                0936
24016b3859 Alis*0937 C--   Start of k loop for vertical flux
75abb052f1 Jean*0938        DO k=Nr,1,-1
0845f1a203 Jean*0939 C--   kUp    Cycles through 1,2 to point to w-layer above
24016b3859 Alis*0940 C--   kDown  Cycles through 2,1 to point to w-layer below
0845f1a203 Jean*0941         kUp  = 1+MOD(k+1,2)
75abb052f1 Jean*0942         kDown= 1+MOD(k,2)
                0943         kp1Msk=1.
09e49e8561 Jean*0944         IF (k.EQ.Nr) kp1Msk=0.
24016b3859 Alis*0945
da54d4c5af Jean*0946 #ifdef ALLOW_AUTODIFF_TAMC
20a156cdce Mart*0947 C     This directive suppresses a warning and hence makes a store
                0948 C     directive unnecessary. Before inserting this directive, the store
                0949 C     directive did not resolve a warning about a potential conflict due
                0950 C     to the incomplete variable rtrans.
                0951 CADJ INCOMPLETE rtrans
ba8396a821 Patr*0952 #endif
                0953
b4413f3c59 Jean*0954 C-- Compute Vertical transport
ca239d4b54 Jean*0955 #ifdef ALLOW_AIM
                0956 C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
                0957         IF ( k.EQ.1 .OR.
cb0d108b91 Jean*0958      &     (useAIM .AND. trIdentity.EQ.GAD_SALINITY .AND. k.EQ.Nr)
ca239d4b54 Jean*0959      &              ) THEN
                0960 #else
                0961         IF ( k.EQ.1 ) THEN
                0962 #endif
b4413f3c59 Jean*0963
                0964 C- Surface interface :
cb0d108b91 Jean*0965          DO j=1-OLy,sNy+OLy
                0966           DO i=1-OLx,sNx+OLx
d077f4a702 Jean*0967            rTransKp(i,j) = kp1Msk*rTrans(i,j)
75abb052f1 Jean*0968            rTrans(i,j) = 0.
                0969            fVerT(i,j,kUp) = 0.
                0970           ENDDO
                0971          ENDDO
b4413f3c59 Jean*0972
75abb052f1 Jean*0973         ELSE
                0974
7db058c202 Patr*0975 C- Interior interface :
cb0d108b91 Jean*0976          DO j=1-OLy,sNy+OLy
                0977           DO i=1-OLx,sNx+OLx
d077f4a702 Jean*0978            rTransKp(i,j) = kp1Msk*rTrans(i,j)
09e49e8561 Jean*0979            rTrans(i,j) = wFld(i,j,k)*rA(i,j,bi,bj)
eaba2fd266 Jean*0980      &                 *deepFac2F(k)*rhoFacF(k)
75abb052f1 Jean*0981      &                 *maskC(i,j,k-1,bi,bj)
616d5d70aa Jean*0982            fVerT(i,j,kUp) = 0.
75abb052f1 Jean*0983           ENDDO
                0984          ENDDO
b4413f3c59 Jean*0985
24016b3859 Alis*0986 C-    Compute vertical advective flux in the interior:
d299295b47 Jean*0987          IF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
                0988      &      .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
20a156cdce Mart*0989           CALL GAD_DST2U1_ADV_R(    bi,bj,k, advectionScheme,
09e49e8561 Jean*0990      I              deltaTLev(k),rTrans,wFld(1-OLx,1-OLy,k),localT3d,
                0991      O              fVerT(1-OLx,1-OLy,kUp), myThid )
0d75a51072 Mart*0992          ELSEIF ( vertAdvecScheme.EQ.ENUM_FLUX_LIMIT ) THEN
20a156cdce Mart*0993           CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, deltaTLev(k),
09e49e8561 Jean*0994      I              rTrans, wFld(1-OLx,1-OLy,k), localT3d,
                0995      O              fVerT(1-OLx,1-OLy,kUp), myThid )
0d75a51072 Mart*0996          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3) THEN
20a156cdce Mart*0997           CALL GAD_DST3_ADV_R(      bi,bj,k, deltaTLev(k),
09e49e8561 Jean*0998      I              rTrans, wFld(1-OLx,1-OLy,k), localT3d,
                0999      O              fVerT(1-OLx,1-OLy,kUp), myThid )
0d75a51072 Mart*1000          ELSEIF ( vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
20a156cdce Mart*1001           CALL GAD_DST3FL_ADV_R(    bi,bj,k, deltaTLev(k),
09e49e8561 Jean*1002      I              rTrans, wFld(1-OLx,1-OLy,k), localT3d,
                1003      O              fVerT(1-OLx,1-OLy,kUp), myThid )
0d75a51072 Mart*1004          ELSEIF ( vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
20a156cdce Mart*1005           CALL GAD_OS7MP_ADV_R(     bi,bj,k, deltaTLev(k),
09e49e8561 Jean*1006      I              rTrans, wFld(1-OLx,1-OLy,k), localT3d,
                1007      O              fVerT(1-OLx,1-OLy,kUp), myThid )
598aebfcee Mart*1008 #ifndef ALLOW_AUTODIFF
20a156cdce Mart*1009          ELSEIF ( usePPMvertAdv .OR. usePQMvertAdv ) THEN
8e4c181d69 Jean*1010 C- copy level from 3d flux data
20a156cdce Mart*1011           DO j = 1-OLy,sNy+OLy
                1012            DO i = 1-OLx,sNx+OLx
                1013             fVerT(i,j,kUp) = afr(i,j,k)
8e4c181d69 Jean*1014            ENDDO
20a156cdce Mart*1015           ENDDO
8e4c181d69 Jean*1016 #endif /* ndef ALLOW_AUTODIFF */
75abb052f1 Jean*1017          ELSE
20a156cdce Mart*1018           STOP 'GAD_ADVECTION: adv. scheme incompatible with mutli-dim'
75abb052f1 Jean*1019          ENDIF
b4413f3c59 Jean*1020
                1021 C- end Surface/Interior if bloc
75abb052f1 Jean*1022         ENDIF
24016b3859 Alis*1023
da54d4c5af Jean*1024 #ifdef ALLOW_AUTODIFF_TAMC
20a156cdce Mart*1025 CADJ STORE rtrans   = loctape_gad_adv_k, key = k, kind = isbyte
                1026 CADJ STORE rtransKp = loctape_gad_adv_k, key = k, kind = isbyte
7448700841 Mart*1027 # ifdef NONLIN_FRSURF
cc249757ef Patr*1028 cph --- following storing of fVerT is critical for correct
                1029 cph --- gradient with multiDimAdvection
                1030 cph --- Without it, kDown component is not properly recomputed
                1031 cph --- This is a TAF bug (and no warning available)
20a156cdce Mart*1032 CADJ STORE fVerT(:,:,:) = loctape_gad_adv_k, key = k, kind = isbyte
7448700841 Mart*1033 # endif
20a156cdce Mart*1034 #endif
105822914f Patr*1035
75abb052f1 Jean*1036 C--   Divergence of vertical fluxes
d077f4a702 Jean*1037 #ifdef GAD_MULTIDIM_COMPRESSIBLE
cb0d108b91 Jean*1038         DO j=1-OLy,sNy+OLy
                1039          DO i=1-OLx,sNx+OLx
d077f4a702 Jean*1040           tmpTrac = localT3d(i,j,k)*locVol3d(i,j,k)
                1041      &      -deltaTLev(k)*( fVerT(i,j,kDown)-fVerT(i,j,kUp) )
                1042      &                   *rkSign*maskInC(i,j,bi,bj)
                1043           localVol(i,j) = locVol3d(i,j,k)
                1044      &      -deltaTLev(k)*( rTransKp(i,j) - rTrans(i,j) )
                1045      &                   *rkSign*maskInC(i,j,bi,bj)
                1046 C- localTij only needed for Variance Bugget: can be move there
                1047           localTij(i,j) = tmpTrac/localVol(i,j)
                1048 C--  without rescaling of tendencies:
935a76deec Jean*1049 c         gTracer(i,j,k) =
                1050 c    &     ( localTij(i,j) - tracer(i,j,k,bi,bj) )/deltaTLev(k)
d077f4a702 Jean*1051 C--  Non-Lin Free-Surf: consistent with rescaling of tendencies
                1052 C     (in FREESURF_RESCALE_G) and RealFreshFlux/addMass.
                1053 C    Also valid for linear Free-Surf (r & r* coords) w/wout RealFreshFlux
                1054 C     and consistent with linFSConserveTr and "surfExpan_" monitor.
935a76deec Jean*1055           gTracer(i,j,k) =
d077f4a702 Jean*1056      &          ( tmpTrac - tracer(i,j,k,bi,bj)*localVol(i,j) )
                1057      &            *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
                1058      &            *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
                1059      &            *recip_rhoFacC(k)
                1060      &            /deltaTLev(k)
                1061          ENDDO
                1062         ENDDO
                1063 #else /* GAD_MULTIDIM_COMPRESSIBLE */
                1064         DO j=1-OLy,sNy+OLy
                1065          DO i=1-OLx,sNx+OLx
                1066           localTij(i,j) = localT3d(i,j,k)
4e66ab0b67 Oliv*1067      &      -deltaTLev(k)*recip_rhoFacC(k)
eaba2fd266 Jean*1068      &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
                1069      &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
                1070      &       *( fVerT(i,j,kDown)-fVerT(i,j,kUp)
d077f4a702 Jean*1071      &         -tracer(i,j,k,bi,bj)*(rTransKp(i,j)-rTrans(i,j))
4d24da31b2 Mart*1072      &        )*rkSign*maskInC(i,j,bi,bj)
935a76deec Jean*1073           gTracer(i,j,k) =
                1074      &     ( localTij(i,j) - tracer(i,j,k,bi,bj) )/deltaTLev(k)
75abb052f1 Jean*1075          ENDDO
                1076         ENDDO
d077f4a702 Jean*1077 #endif /* GAD_MULTIDIM_COMPRESSIBLE */
0845f1a203 Jean*1078
81c8d7b9aa Jean*1079 #ifdef ALLOW_DIAGNOSTICS
78c3e29ef3 Jean*1080         IF ( doDiagAdvR ) THEN
20a156cdce Mart*1081          diagName = 'ADVr'//diagSufx
                1082          CALL DIAGNOSTICS_FILL( fVerT(1-OLx,1-OLy,kUp),
                1083      &                          diagName, k,1, 2,bi,bj, myThid )
81c8d7b9aa Jean*1084         ENDIF
50d8304171 Ryan*1085 #ifdef ALLOW_LAYERS
9b89fcf692 antn*1086         IF ( layers_useThermo ) THEN
20a156cdce Mart*1087          CALL LAYERS_FILL( fVerT(1-OLx,1-OLy,kUp), trIdentity,
                1088      &                     'AFR', k, 1, 2,bi,bj, myThid)
50d8304171 Ryan*1089         ENDIF
                1090 #endif /* ALLOW_LAYERS */
d077f4a702 Jean*1091 #endif /* ALLOW_DIAGNOSTICS */
81c8d7b9aa Jean*1092
24016b3859 Alis*1093 C--   End of K loop for vertical flux
75abb052f1 Jean*1094        ENDDO
                1095 C--   end of if not.implicitAdvection block
0845f1a203 Jean*1096       ENDIF
20a156cdce Mart*1097 #endif /* DISABLE_MULTIDIM_ADVECTION */
24016b3859 Alis*1098
                1099       RETURN
                1100       END