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