pkg/generic_advdiff/gad_som_advect.F

d7ce0d34f8 Jean*0001 #include "GAD_OPTIONS.h"
7448700841 Mart*0002 #ifdef ALLOW_PTRACERS
                0003 # include "PTRACERS_OPTIONS.h"
                0004 #endif
1574069d50 Mart*0005 #ifdef ALLOW_AUTODIFF
                0006 # include "AUTODIFF_OPTIONS.h"
                0007 #endif
d7ce0d34f8 Jean*0008
                0009 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0010 CBOP
                0011 C !ROUTINE: GAD_SOM_ADVECT
                0012
                0013 C !INTERFACE: ==========================================================
                0014       SUBROUTINE GAD_SOM_ADVECT(
                0015      I     implicitAdvection, advectionScheme, vertAdvecScheme,
4e66ab0b67 Oliv*0016      I     tracerIdentity, deltaTLev,
09e49e8561 Jean*0017      I     uFld, vFld, wFld, tracer,
d7ce0d34f8 Jean*0018      U     smTr,
                0019      O     gTracer,
                0020      I     bi,bj, myTime,myIter,myThid)
                0021
                0022 C !DESCRIPTION:
                0023 C Calculates the tendency of a tracer due to advection.
                0024 C It uses the 2nd-Order moment advection scheme with multi-dimensional method
                0025 C  see Prather, 1986, JGR, v.91, D-6, pp.6671-6681.
                0026 C
                0027 C The tendency (output) is over-written by this routine.
                0028
                0029 C !USES: ===============================================================
                0030       IMPLICIT NONE
                0031 #include "SIZE.h"
                0032 #include "EEPARAMS.h"
                0033 #include "PARAMS.h"
                0034 #include "GRID.h"
                0035 #include "GAD.h"
1574069d50 Mart*0036 #ifdef ALLOW_AUTODIFF_TAMC
                0037 # include "tamc.h"
                0038 # ifdef ALLOW_PTRACERS
                0039 #  include "PTRACERS_SIZE.h"
                0040 # endif
                0041 #endif
b79a2b44f2 Jean*0042 #ifdef ALLOW_EXCH2
f9f661930b Jean*0043 #include "W2_EXCH2_SIZE.h"
b79a2b44f2 Jean*0044 #include "W2_EXCH2_TOPOLOGY.h"
                0045 #endif /* ALLOW_EXCH2 */
d7ce0d34f8 Jean*0046
                0047 C !INPUT PARAMETERS: ===================================================
                0048 C  implicitAdvection :: implicit vertical advection (later on)
                0049 C  advectionScheme   :: advection scheme to use (Horizontal plane)
                0050 C  vertAdvecScheme   :: advection scheme to use (vertical direction)
                0051 C  tracerIdentity    :: tracer identifier (required only for OBCS)
09e49e8561 Jean*0052 C  uFld              :: Advection velocity field, zonal component
                0053 C  vFld              :: Advection velocity field, meridional component
                0054 C  wFld              :: Advection velocity field, vertical component
d7ce0d34f8 Jean*0055 C  tracer            :: tracer field
                0056 C  bi,bj             :: tile indices
                0057 C  myTime            :: current time
                0058 C  myIter            :: iteration number
                0059 C  myThid            :: thread number
                0060       LOGICAL implicitAdvection
                0061       INTEGER advectionScheme, vertAdvecScheme
                0062       INTEGER tracerIdentity
4e66ab0b67 Oliv*0063       _RL deltaTLev(Nr)
09e49e8561 Jean*0064       _RL uFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0065       _RL vFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0066       _RL wFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
d7ce0d34f8 Jean*0067       _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0068       INTEGER bi,bj
                0069       _RL myTime
                0070       INTEGER myIter
                0071       INTEGER myThid
                0072
                0073 C !OUTPUT PARAMETERS: ==================================================
                0074 C  smTr              :: tracer 1rst & 2nd Order moments
                0075 C  gTracer           :: tendency array
                0076       _RL smTr   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy,nSOM)
935a76deec Jean*0077       _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
d7ce0d34f8 Jean*0078
7448700841 Mart*0079 #if ( defined GAD_ALLOW_TS_SOM_ADV || defined PTRACERS_ALLOW_DYN_STATE )
d7ce0d34f8 Jean*0080 C !LOCAL VARIABLES: ====================================================
                0081 C  maskUp        :: 2-D array mask for W points
                0082 C  i,j,k         :: loop indices
                0083 C  kUp           :: index into 2 1/2D array, toggles between 1 and 2
                0084 C  kDown         :: index into 2 1/2D array, toggles between 2 and 1
                0085 C  xA,yA         :: areas of X and Y face of tracer cells
                0086 C  uTrans,vTrans :: 2-D arrays of volume transports at U,V points
                0087 C  rTrans        :: 2-D arrays of volume transports at W points
                0088 C  afx           :: 2-D array for horizontal advective flux, x direction
                0089 C  afy           :: 2-D array for horizontal advective flux, y direction
                0090 C  afr           :: 2-D array for vertical advective flux
                0091 C  fVerT         :: 2 1/2D arrays for vertical advective flux
                0092 C  calc_fluxes_X :: logical to indicate to calculate fluxes in X dir
                0093 C  calc_fluxes_Y :: logical to indicate to calculate fluxes in Y dir
                0094 C  interiorOnly  :: only update the interior of myTile, but not the edges
                0095 C  overlapOnly   :: only update the edges of myTile, but not the interior
                0096 C  npass         :: number of passes in multi-dimensional method
                0097 C  ipass         :: number of the current pass being made
                0098 C  myTile        :: variables used to determine which cube face
                0099 C  nCFace        :: owns a tile for cube grid runs using
                0100 C                :: multi-dim advection.
                0101 C [N,S,E,W]_edge :: true if N,S,E,W edge of myTile is an Edge of the cube
56c999258f Jean*0102 C  msgBuf        :: Informational/error message buffer
d7ce0d34f8 Jean*0103       _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0104       INTEGER i,j,k,km1,kUp,kDown
                0105       _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0106       _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0107       _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0108       _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0109       _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0110       _RL afx     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0111       _RL afy     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0112       _RL afr     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0113       _RL  smVol  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0114       _RL  smTr0  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0115       _RL  alp    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0116       _RL  aln    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0117       _RL  fp_v   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0118       _RL  fn_v   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0119       _RL  fp_o   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0120       _RL  fn_o   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0121       _RL  fp_x   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0122       _RL  fn_x   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0123       _RL  fp_y   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0124       _RL  fn_y   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0125       _RL  fp_z   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0126       _RL  fn_z   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0127       _RL  fp_xx  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0128       _RL  fn_xx  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0129       _RL  fp_yy  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0130       _RL  fn_yy  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0131       _RL  fp_zz  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0132       _RL  fn_zz  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0133       _RL  fp_xy  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0134       _RL  fn_xy  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0135       _RL  fp_xz  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0136       _RL  fn_xz  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0137       _RL  fp_yz  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0138       _RL  fn_yz  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
b79a2b44f2 Jean*0139       _RL  smCorners(OLx,OLy,4,-1:nSOM)
f08c4b8486 Jean*0140 c     _RL  localTr
d7ce0d34f8 Jean*0141       LOGICAL calc_fluxes_X, calc_fluxes_Y
b79a2b44f2 Jean*0142       LOGICAL interiorOnly, overlapOnly
d7ce0d34f8 Jean*0143       INTEGER limiter
                0144       INTEGER npass, ipass
b79a2b44f2 Jean*0145       INTEGER nCFace, n
                0146       LOGICAL N_edge, S_edge, E_edge, W_edge
56c999258f Jean*0147       LOGICAL noFlowAcrossSurf
d7ce0d34f8 Jean*0148       CHARACTER*(MAX_LEN_MBUF) msgBuf
1574069d50 Mart*0149 #ifdef ALLOW_AUTODIFF_TAMC
                0150       INTEGER ipasskey
                0151 #endif
b79a2b44f2 Jean*0152 #ifdef ALLOW_EXCH2
                0153       INTEGER myTile
                0154 #endif
d7ce0d34f8 Jean*0155 #ifdef ALLOW_DIAGNOSTICS
                0156       CHARACTER*8 diagName
b79a2b44f2 Jean*0157       CHARACTER*4 diagSufx
                0158       LOGICAL     doDiagAdvX, doDiagAdvY, doDiagAdvR
                0159 C-    Functions:
                0160       CHARACTER*4 GAD_DIAG_SUFX
d7ce0d34f8 Jean*0161       EXTERNAL    GAD_DIAG_SUFX
b79a2b44f2 Jean*0162       LOGICAL  DIAGNOSTICS_IS_ON
                0163       EXTERNAL DIAGNOSTICS_IS_ON
d7ce0d34f8 Jean*0164 #endif
                0165 CEOP
                0166
                0167 #ifdef ALLOW_DIAGNOSTICS
b79a2b44f2 Jean*0168 C--   Set diagnostics flags and suffix for the current tracer
                0169       doDiagAdvX = .FALSE.
                0170       doDiagAdvY = .FALSE.
                0171       doDiagAdvR = .FALSE.
d7ce0d34f8 Jean*0172       IF ( useDiagnostics ) THEN
                0173         diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )
b79a2b44f2 Jean*0174         diagName = 'ADVx'//diagSufx
                0175         doDiagAdvX = DIAGNOSTICS_IS_ON( diagName, myThid )
                0176         diagName = 'ADVy'//diagSufx
                0177         doDiagAdvY = DIAGNOSTICS_IS_ON( diagName, myThid )
                0178         diagName = 'ADVr'//diagSufx
                0179         doDiagAdvR = DIAGNOSTICS_IS_ON( diagName, myThid )
d7ce0d34f8 Jean*0180       ENDIF
                0181 #endif
                0182
                0183 C--   Set up work arrays with valid (i.e. not NaN) values
                0184 C     These inital values do not alter the numerical results. They
                0185 C     just ensure that all memory references are to valid floating
                0186 C     point numbers. This prevents spurious hardware signals due to
                0187 C     uninitialised but inert locations.
b79a2b44f2 Jean*0188       DO j=1-OLy,sNy+OLy
                0189        DO i=1-OLx,sNx+OLx
                0190          afx(i,j) = 0.
                0191          afy(i,j) = 0.
09e49e8561 Jean*0192 C-    xA,yA,uTrans,vTrans are set over the full domain
                0193 C      => no need for extra initialisation
d7ce0d34f8 Jean*0194 c       xA(i,j)      = 0. _d 0
                0195 c       yA(i,j)      = 0. _d 0
                0196 c       uTrans(i,j)  = 0. _d 0
                0197 c       vTrans(i,j)  = 0. _d 0
                0198 C-    rTrans is set over the full domain: no need for extra initialisation
                0199 c       rTrans(i,j)  = 0. _d 0
b79a2b44f2 Jean*0200        ENDDO
                0201       ENDDO
                0202       DO n=-1,nSOM
                0203        DO k=1,4
                0204         DO j=1,OLy
                0205          DO i=1,OLx
                0206           smCorners(i,j,k,n) = 0.
                0207          ENDDO
                0208         ENDDO
                0209        ENDDO
                0210       ENDDO
d7ce0d34f8 Jean*0211
                0212       IF ( implicitAdvection ) THEN
                0213         WRITE(msgBuf,'(2A)') 'S/R GAD_SOM_ADVECT: ',
                0214      &     'not coded for implicit-vertical Advection'
                0215         CALL PRINT_ERROR( msgBuf, myThid )
                0216         STOP 'ABNORMAL END: S/R GAD_SOM_ADVECT'
                0217       ENDIF
                0218       IF ( vertAdvecScheme .NE. advectionScheme ) THEN
                0219         WRITE(msgBuf,'(2A)') 'S/R GAD_SOM_ADVECT: ',
                0220      &     'not coded for different vertAdvecScheme'
                0221         CALL PRINT_ERROR( msgBuf, myThid )
                0222         STOP 'ABNORMAL END: S/R GAD_SOM_ADVECT'
                0223       ENDIF
                0224
b79a2b44f2 Jean*0225 C--   Set tile-specific parameters for horizontal fluxes
                0226       IF (useCubedSphereExchange) THEN
                0227        npass  = 3
                0228 #ifdef ALLOW_EXCH2
c424ee7cc7 Jean*0229        myTile = W2_myTileList(bi,bj)
b79a2b44f2 Jean*0230        nCFace = exch2_myFace(myTile)
                0231        N_edge = exch2_isNedge(myTile).EQ.1
                0232        S_edge = exch2_isSedge(myTile).EQ.1
                0233        E_edge = exch2_isEedge(myTile).EQ.1
                0234        W_edge = exch2_isWedge(myTile).EQ.1
                0235 #else
                0236        nCFace = bi
                0237        N_edge = .TRUE.
                0238        S_edge = .TRUE.
                0239        E_edge = .TRUE.
                0240        W_edge = .TRUE.
                0241 #endif
                0242       ELSE
                0243        npass  = 2
                0244        nCFace = 0
                0245        N_edge = .FALSE.
                0246        S_edge = .FALSE.
                0247        E_edge = .FALSE.
                0248        W_edge = .FALSE.
                0249       ENDIF
                0250
d7ce0d34f8 Jean*0251       limiter = MOD(advectionScheme, 10)
                0252
1574069d50 Mart*0253 #ifdef ALLOW_AUTODIFF_TAMC
                0254       IF ( npass.GT.maxcube ) THEN
                0255         WRITE(msgBuf,'(A,2(I3,A))') 'S/R GAD_SOM_ADVECT: npass =',
                0256      &     npass, ' >', maxcube, ' = maxcube, ==> check "tamc.h"'
                0257         CALL PRINT_ERROR( msgBuf, myThid )
                0258         STOP 'ABNORMAL END: S/R GAD_SOM_ADVECT'
                0259       ENDIF
                0260 CADJ INIT somtape = COMMON, 1
                0261 CADJ INIT somtape_k = COMMON, Nr
                0262 CADJ INIT somtape_k_pass = COMMON, maxcube*Nr
                0263 #endif
d7ce0d34f8 Jean*0264 C--   Start of k loop for horizontal fluxes
                0265       DO k=1,Nr
                0266
                0267 C--   Get temporary terms used by tendency routines
09e49e8561 Jean*0268        DO j=1-OLy,sNy+OLy
                0269         DO i=1-OLx,sNx+OLx
                0270           xA(i,j) = _dyG(i,j,bi,bj)*deepFacC(k)
                0271      &            *drF(k)*_hFacW(i,j,k,bi,bj)
                0272           yA(i,j) = _dxG(i,j,bi,bj)*deepFacC(k)
                0273      &            *drF(k)*_hFacS(i,j,k,bi,bj)
                0274         ENDDO
                0275        ENDDO
                0276 C--   Calculate "volume transports" through tracer cell faces.
                0277 C     anelastic: scaled by rhoFacC (~ mass transport)
                0278        DO j=1-OLy,sNy+OLy
                0279         DO i=1-OLx,sNx+OLx
                0280           uTrans(i,j) = uFld(i,j,k)*xA(i,j)*rhoFacC(k)
                0281           vTrans(i,j) = vFld(i,j,k)*yA(i,j)*rhoFacC(k)
                0282         ENDDO
                0283        ENDDO
d7ce0d34f8 Jean*0284
                0285 C--   grid-box volume and tracer content (zero order moment)
b79a2b44f2 Jean*0286        DO j=1-OLy,sNy+OLy
                0287         DO i=1-OLx,sNx+OLx
d7ce0d34f8 Jean*0288          smVol(i,j,k) = rA(i,j,bi,bj)*deepFac2C(k)
                0289      &                *drF(k)*hFacC(i,j,k,bi,bj)
                0290      &                *rhoFacC(k)
                0291          smTr0(i,j,k) = tracer(i,j,k,bi,bj)*smVol(i,j,k)
                0292 C-    fill empty grid-box:
                0293          smVol(i,j,k) = smVol(i,j,k)
                0294      &                + (1. _d 0 - maskC(i,j,k,bi,bj))
b79a2b44f2 Jean*0295         ENDDO
d7ce0d34f8 Jean*0296        ENDDO
                0297
1574069d50 Mart*0298 #ifdef ALLOW_AUTODIFF_TAMC
                0299 CADJ STORE smtr(:,:,k,bi,bj,:) = somtape_k, key=k, kind=isbyte
                0300 CADJ STORE smtr0(:,:,k)        = somtape_k, key=k, kind=isbyte
                0301 CADJ STORE smVol(:,:,k)        = somtape_k, key=k, kind=isbyte
                0302 CADJ STORE smcorners           = somtape_k, key=k, kind=isbyte
                0303 #endif
d7ce0d34f8 Jean*0304 C--   Multiple passes for different directions on different tiles
                0305 C--   For cube need one pass for each of red, green and blue axes.
b79a2b44f2 Jean*0306        DO ipass=1,npass
1574069d50 Mart*0307 #ifdef ALLOW_AUTODIFF_TAMC
                0308         ipasskey = ipass + (k-1) * maxcube
                0309 #endif /* ALLOW_AUTODIFF_TAMC */
b79a2b44f2 Jean*0310
                0311         interiorOnly = .FALSE.
                0312         overlapOnly  = .FALSE.
                0313         IF (useCubedSphereExchange) THEN
                0314 C-    CubedSphere : pass 3 times, with partial update of local tracer field
                0315          IF (ipass.EQ.1) THEN
                0316           overlapOnly   = MOD(nCFace,3).EQ.0
                0317           interiorOnly  = MOD(nCFace,3).NE.0
                0318           calc_fluxes_X = nCFace.EQ.6 .OR. nCFace.EQ.1 .OR. nCFace.EQ.2
                0319           calc_fluxes_Y = nCFace.EQ.3 .OR. nCFace.EQ.4 .OR. nCFace.EQ.5
                0320          ELSEIF (ipass.EQ.2) THEN
                0321           overlapOnly   = MOD(nCFace,3).EQ.2
                0322           interiorOnly  = MOD(nCFace,3).EQ.1
                0323           calc_fluxes_X = nCFace.EQ.2 .OR. nCFace.EQ.3 .OR. nCFace.EQ.4
                0324           calc_fluxes_Y = nCFace.EQ.5 .OR. nCFace.EQ.6 .OR. nCFace.EQ.1
                0325          ELSE
                0326           interiorOnly  = .TRUE.
                0327           calc_fluxes_X = nCFace.EQ.5 .OR. nCFace.EQ.6
                0328           calc_fluxes_Y = nCFace.EQ.2 .OR. nCFace.EQ.3
                0329          ENDIF
                0330         ELSE
d7ce0d34f8 Jean*0331 C-    not CubedSphere
b79a2b44f2 Jean*0332           calc_fluxes_X = MOD(ipass,2).EQ.1
                0333           calc_fluxes_Y = .NOT.calc_fluxes_X
                0334         ENDIF
1574069d50 Mart*0335 #ifdef ALLOW_AUTODIFF_TAMC
                0336 CADJ STORE smtr(:,:,k,bi,bj,:) =somtape_k_pass,key=ipasskey,kind=isbyte
                0337 CADJ STORE smtr0(:,:,k)        =somtape_k_pass,key=ipasskey,kind=isbyte
                0338 CADJ STORE smVol(:,:,k)        =somtape_k_pass,key=ipasskey,kind=isbyte
                0339 CADJ STORE smcorners           =somtape_k_pass,key=ipasskey,kind=isbyte
                0340 #endif
d7ce0d34f8 Jean*0341
                0342 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0343
b79a2b44f2 Jean*0344 C--   X direction
                0345 C-     Do not compute fluxes if
                0346 C       a) needed in overlap only
                0347 C   and b) the overlap of myTile are not cube-face Edges
                0348         IF ( calc_fluxes_X .AND.
                0349      &      (.NOT.overlapOnly .OR. N_edge .OR. S_edge)
                0350      &     ) THEN
                0351
                0352 C-     Internal exchange for calculations in X
                0353           IF ( useCubedSphereExchange .AND. .NOT.interiorOnly ) THEN
                0354            CALL GAD_SOM_PREP_CS_CORNER(
                0355      U                     smVol, smTr0, smTr, smCorners,
                0356      I                     .TRUE., overlapOnly, interiorOnly,
                0357      I                     N_edge, S_edge, E_edge, W_edge,
                0358      I                     ipass, k, Nr, bi, bj, myThid )
                0359           ENDIF
                0360
                0361 C-     Solve advection in X and update moments
                0362           IF ( advectionScheme.EQ.ENUM_SOM_PRATHER
                0363      &       .OR. advectionScheme.EQ.ENUM_SOM_LIMITER ) THEN
                0364            CALL GAD_SOM_ADV_X(
d7ce0d34f8 Jean*0365      I                     bi,bj,k, limiter,
b79a2b44f2 Jean*0366      I                     overlapOnly, interiorOnly,
                0367      I                     N_edge, S_edge, E_edge, W_edge,
4e66ab0b67 Oliv*0368      I                     deltaTLev(k), uTrans,
72de869c1b Jean*0369      I                     maskInC(1-OLx,1-OLy,bi,bj),
d7ce0d34f8 Jean*0370      U                     smVol(1-OLx,1-OLy,k),
                0371      U                     smTr0(1-OLx,1-OLy,k),
                0372      U                     smTr(1-OLx,1-OLy,k,bi,bj,1),
                0373      U                     smTr(1-OLx,1-OLy,k,bi,bj,2),
                0374      U                     smTr(1-OLx,1-OLy,k,bi,bj,3),
                0375      U                     smTr(1-OLx,1-OLy,k,bi,bj,4),
                0376      U                     smTr(1-OLx,1-OLy,k,bi,bj,5),
                0377      U                     smTr(1-OLx,1-OLy,k,bi,bj,6),
                0378      U                     smTr(1-OLx,1-OLy,k,bi,bj,7),
                0379      U                     smTr(1-OLx,1-OLy,k,bi,bj,8),
                0380      U                     smTr(1-OLx,1-OLy,k,bi,bj,9),
                0381      O                     afx, myThid )
b79a2b44f2 Jean*0382           ELSE
                0383            STOP 'GAD_SOM_ADVECT: adv. scheme incompatibale with SOM'
                0384           ENDIF
d7ce0d34f8 Jean*0385
                0386 C--   End of X direction
b79a2b44f2 Jean*0387         ENDIF
d7ce0d34f8 Jean*0388
                0389 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0390
b79a2b44f2 Jean*0391 C--   Y direction
d7ce0d34f8 Jean*0392 C-     Do not compute fluxes if
                0393 C       a) needed in overlap only
                0394 C   and b) the overlap of myTile are not cube-face edges
b79a2b44f2 Jean*0395         IF ( calc_fluxes_Y .AND.
                0396      &      (.NOT.overlapOnly .OR. E_edge .OR. W_edge)
                0397      &     ) THEN
                0398
                0399 C-     Internal exchange for calculations in Y
                0400           IF ( useCubedSphereExchange .AND. .NOT.interiorOnly ) THEN
                0401            CALL GAD_SOM_PREP_CS_CORNER(
                0402      U                     smVol, smTr0, smTr, smCorners,
                0403      I                     .FALSE., overlapOnly, interiorOnly,
                0404      I                     N_edge, S_edge, E_edge, W_edge,
                0405      I                     iPass, k, Nr, bi, bj, myThid )
                0406           ENDIF
                0407
                0408 C-     Solve advection in Y and update moments
                0409           IF ( advectionScheme.EQ.ENUM_SOM_PRATHER
                0410      &       .OR. advectionScheme.EQ.ENUM_SOM_LIMITER ) THEN
                0411            CALL GAD_SOM_ADV_Y(
d7ce0d34f8 Jean*0412      I                     bi,bj,k, limiter,
b79a2b44f2 Jean*0413      I                     overlapOnly, interiorOnly,
                0414      I                     N_edge, S_edge, E_edge, W_edge,
4e66ab0b67 Oliv*0415      I                     deltaTLev(k), vTrans,
72de869c1b Jean*0416      I                     maskInC(1-OLx,1-OLy,bi,bj),
d7ce0d34f8 Jean*0417      U                     smVol(1-OLx,1-OLy,k),
                0418      U                     smTr0(1-OLx,1-OLy,k),
                0419      U                     smTr(1-OLx,1-OLy,k,bi,bj,1),
                0420      U                     smTr(1-OLx,1-OLy,k,bi,bj,2),
                0421      U                     smTr(1-OLx,1-OLy,k,bi,bj,3),
                0422      U                     smTr(1-OLx,1-OLy,k,bi,bj,4),
                0423      U                     smTr(1-OLx,1-OLy,k,bi,bj,5),
                0424      U                     smTr(1-OLx,1-OLy,k,bi,bj,6),
                0425      U                     smTr(1-OLx,1-OLy,k,bi,bj,7),
                0426      U                     smTr(1-OLx,1-OLy,k,bi,bj,8),
                0427      U                     smTr(1-OLx,1-OLy,k,bi,bj,9),
                0428      O                     afy, myThid )
b79a2b44f2 Jean*0429           ELSE
                0430            STOP 'GAD_SOM_ADVECT: adv. scheme incompatibale with SOM'
                0431           ENDIF
d7ce0d34f8 Jean*0432
                0433 C--   End of Y direction
b79a2b44f2 Jean*0434         ENDIF
d7ce0d34f8 Jean*0435 C--   End of ipass loop
b79a2b44f2 Jean*0436        ENDDO
d7ce0d34f8 Jean*0437
b79a2b44f2 Jean*0438        IF ( implicitAdvection ) THEN
d7ce0d34f8 Jean*0439 C-    explicit advection is done ; store tendency in gTracer:
                0440         DO j=1-OLy,sNy+OLy
                0441          DO i=1-OLx,sNx+OLx
f08c4b8486 Jean*0442 C--  without rescaling of tendencies:
                0443 c         localTr = smTr0(i,j,k)/smVol(i,j,k)
935a76deec Jean*0444 c         gTracer(i,j,k) = ( localTr - tracer(i,j,k,bi,bj) )
                0445 c    &                   / deltaTLev(k)
f08c4b8486 Jean*0446 C--  consistent with rescaling of tendencies (in FREESURF_RESCALE_G):
935a76deec Jean*0447           gTracer(i,j,k) =
f08c4b8486 Jean*0448      &          ( smTr0(i,j,k) - tracer(i,j,k,bi,bj)*smVol(i,j,k) )
                0449      &            *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
                0450      &            *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
                0451      &            *recip_rhoFacC(k)
4e66ab0b67 Oliv*0452      &            /deltaTLev(k)
d7ce0d34f8 Jean*0453          ENDDO
                0454         ENDDO
b79a2b44f2 Jean*0455        ENDIF
d7ce0d34f8 Jean*0456
                0457 #ifdef ALLOW_DIAGNOSTICS
b79a2b44f2 Jean*0458        IF ( doDiagAdvX ) THEN
                0459          diagName = 'ADVx'//diagSufx
                0460          CALL DIAGNOSTICS_FILL(afx,diagName, k,1, 2,bi,bj, myThid )
                0461        ENDIF
                0462        IF ( doDiagAdvY ) THEN
                0463          diagName = 'ADVy'//diagSufx
                0464          CALL DIAGNOSTICS_FILL(afy,diagName, k,1, 2,bi,bj, myThid )
                0465        ENDIF
50d8304171 Ryan*0466 #ifdef ALLOW_LAYERS
                0467        IF ( useLayers ) THEN
                0468          CALL LAYERS_FILL(afx,tracerIdentity,'AFX',k,1,2,bi,bj,myThid)
                0469          CALL LAYERS_FILL(afy,tracerIdentity,'AFY',k,1,2,bi,bj,myThid)
                0470        ENDIF
                0471 #endif /* ALLOW_LAYERS */
d7ce0d34f8 Jean*0472 #endif
                0473
                0474 #ifdef ALLOW_DEBUG
188f784a21 Jean*0475        IF ( debugLevel .GE. debLevC
d7ce0d34f8 Jean*0476      &   .AND. tracerIdentity.EQ.GAD_TEMPERATURE
                0477      &   .AND. k.LE.3 .AND. myIter.EQ.1+nIter0
                0478      &   .AND. nPx.EQ.1 .AND. nPy.EQ.1
                0479      &   .AND. useCubedSphereExchange ) THEN
                0480         CALL DEBUG_CS_CORNER_UV( ' afx,afy from GAD_SOM_ADVECT',
                0481      &             afx,afy, k, standardMessageUnit,bi,bj,myThid )
b79a2b44f2 Jean*0482        ENDIF
d7ce0d34f8 Jean*0483 #endif /* ALLOW_DEBUG */
                0484
                0485 C--   End of K loop for horizontal fluxes
                0486       ENDDO
                0487
                0488 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0489
56c999258f Jean*0490       noFlowAcrossSurf = rigidLid .OR. nonlinFreeSurf.GE.1
                0491      &                            .OR. select_rStar.NE.0
                0492
d7ce0d34f8 Jean*0493       IF ( .NOT.implicitAdvection ) THEN
                0494 C--   Apply limiter (if any):
1574069d50 Mart*0495 #ifdef ALLOW_AUTODIFF_TAMC
                0496 CADJ STORE smtr(:,:,:,bi,bj,:) = somtape, key=1, kind=isbyte
                0497 CADJ STORE smtr0               = somtape, key=1, kind=isbyte
                0498 #endif
d7ce0d34f8 Jean*0499        CALL GAD_SOM_LIM_R( bi,bj, limiter,
                0500      U                     smVol,
                0501      U                     smTr0,
                0502      U                     smTr(1-OLx,1-OLy,1,bi,bj,1),
                0503      U                     smTr(1-OLx,1-OLy,1,bi,bj,2),
                0504      U                     smTr(1-OLx,1-OLy,1,bi,bj,3),
                0505      U                     smTr(1-OLx,1-OLy,1,bi,bj,4),
                0506      U                     smTr(1-OLx,1-OLy,1,bi,bj,5),
                0507      U                     smTr(1-OLx,1-OLy,1,bi,bj,6),
                0508      U                     smTr(1-OLx,1-OLy,1,bi,bj,7),
                0509      U                     smTr(1-OLx,1-OLy,1,bi,bj,8),
                0510      U                     smTr(1-OLx,1-OLy,1,bi,bj,9),
                0511      I                     myThid )
                0512
                0513 C--   Start of k loop for vertical flux
                0514        DO k=Nr,1,-1
1574069d50 Mart*0515 #ifdef ALLOW_AUTODIFF_TAMC
                0516 CADJ STORE smtr(:,:,k,bi,bj,:) = somtape_k, key=k, kind=isbyte
                0517 CADJ STORE smtr0(:,:,k)        = somtape_k, key=k, kind=isbyte
                0518 CADJ STORE smVol(:,:,k)        = somtape_k, key=k, kind=isbyte
                0519 CADJ STORE alp, aln   = somtape_k, key=k, kind=isbyte
                0520 CADJ STORE fp_v,fn_v,fp_o,fn_o,fp_x,fn_x,fp_y,fn_y,fp_z,fn_z
                0521 CADJ &                         = somtape_k, key=k, kind=isbyte
                0522 CADJ STORE fp_xx, fn_xx, fp_yy, fn_yy, fp_zz, fn_zz
                0523 CADJ &                         = somtape_k, key=k, kind=isbyte
                0524 CADJ STORE fp_xy, fn_xy, fp_xz, fn_xz, fp_yz, fn_yz
                0525 CADJ &                         = somtape_k, key=k, kind=isbyte
                0526 #endif
d7ce0d34f8 Jean*0527 C--   kUp    Cycles through 1,2 to point to w-layer above
                0528 C--   kDown  Cycles through 2,1 to point to w-layer below
                0529         kUp  = 1+MOD(Nr-k,2)
                0530         kDown= 1+MOD(Nr-k+1,2)
                0531         IF (k.EQ.Nr) THEN
                0532 C--   Set advective fluxes at the very bottom:
                0533          DO j=1-OLy,sNy+OLy
                0534           DO i=1-OLx,sNx+OLx
                0535            alp  (i,j,kDown) = 0. _d 0
                0536            aln  (i,j,kDown) = 0. _d 0
                0537            fp_v (i,j,kDown) = 0. _d 0
                0538            fn_v (i,j,kDown) = 0. _d 0
                0539            fp_o (i,j,kDown) = 0. _d 0
                0540            fn_o (i,j,kDown) = 0. _d 0
                0541            fp_x (i,j,kDown) = 0. _d 0
                0542            fn_x (i,j,kDown) = 0. _d 0
                0543            fp_y (i,j,kDown) = 0. _d 0
                0544            fn_y (i,j,kDown) = 0. _d 0
                0545            fp_z (i,j,kDown) = 0. _d 0
                0546            fn_z (i,j,kDown) = 0. _d 0
                0547            fp_xx(i,j,kDown) = 0. _d 0
                0548            fn_xx(i,j,kDown) = 0. _d 0
                0549            fp_yy(i,j,kDown) = 0. _d 0
                0550            fn_yy(i,j,kDown) = 0. _d 0
                0551            fp_zz(i,j,kDown) = 0. _d 0
                0552            fn_zz(i,j,kDown) = 0. _d 0
                0553            fp_xy(i,j,kDown) = 0. _d 0
                0554            fn_xy(i,j,kDown) = 0. _d 0
                0555            fp_xz(i,j,kDown) = 0. _d 0
                0556            fn_xz(i,j,kDown) = 0. _d 0
                0557            fp_yz(i,j,kDown) = 0. _d 0
                0558            fn_yz(i,j,kDown) = 0. _d 0
                0559           ENDDO
                0560          ENDDO
                0561         ENDIF
                0562
                0563 C-- Compute Vertical transport
                0564 #ifdef ALLOW_AIM
                0565 C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
                0566 c       IF ( k.EQ.1 .OR.
                0567 c    &     (useAIM .AND. tracerIdentity.EQ.GAD_SALINITY .AND. k.EQ.Nr)
                0568 c    &              ) THEN
                0569 #else
                0570 c       IF ( k.EQ.1 ) THEN
                0571 #endif
56c999258f Jean*0572         IF ( noFlowAcrossSurf .AND. k.EQ.1 ) THEN
d7ce0d34f8 Jean*0573 C- Surface interface :
                0574          DO j=1-OLy,sNy+OLy
                0575           DO i=1-OLx,sNx+OLx
                0576            rTrans(i,j) = 0.
                0577            maskUp(i,j) = 0.
                0578           ENDDO
                0579          ENDDO
                0580
56c999258f Jean*0581         ELSEIF ( noFlowAcrossSurf ) THEN
d7ce0d34f8 Jean*0582 C- Interior interface :
                0583          DO j=1-OLy,sNy+OLy
                0584           DO i=1-OLx,sNx+OLx
09e49e8561 Jean*0585            rTrans(i,j) = wFld(i,j,k)*rA(i,j,bi,bj)
d7ce0d34f8 Jean*0586      &                 *deepFac2F(k)*rhoFacF(k)
                0587      &                 *maskC(i,j,k-1,bi,bj)
                0588            maskUp(i,j) = 1.
                0589           ENDDO
                0590          ENDDO
                0591
                0592         ELSE
                0593 C- Linear Free-Surface: do not mask rTrans :
                0594          km1= MAX(k-1,1)
                0595          DO j=1-OLy,sNy+OLy
                0596           DO i=1-OLx,sNx+OLx
09e49e8561 Jean*0597            rTrans(i,j) = wFld(i,j,k)*rA(i,j,bi,bj)
d7ce0d34f8 Jean*0598      &                 *deepFac2F(k)*rhoFacF(k)
                0599            maskUp(i,j) = maskC(i,j,km1,bi,bj)*maskC(i,j,k,bi,bj)
                0600           ENDDO
                0601          ENDDO
                0602
                0603 C- end Surface/Interior if bloc
                0604         ENDIF
                0605
                0606 C-    Compute vertical advective flux in the interior:
                0607         IF ( vertAdvecScheme.EQ.ENUM_SOM_PRATHER
b79a2b44f2 Jean*0608      &     .OR. vertAdvecScheme.EQ.ENUM_SOM_LIMITER ) THEN
                0609            CALL GAD_SOM_ADV_R(
d7ce0d34f8 Jean*0610      I                     bi,bj,k, kUp, kDown,
4e66ab0b67 Oliv*0611      I                     deltaTLev(k), rTrans, maskUp,
72de869c1b Jean*0612      I                     maskInC(1-OLx,1-OLy,bi,bj),
d7ce0d34f8 Jean*0613      U                     smVol,
                0614      U                     smTr0,
                0615      U                     smTr(1-OLx,1-OLy,1,bi,bj,1),
                0616      U                     smTr(1-OLx,1-OLy,1,bi,bj,2),
                0617      U                     smTr(1-OLx,1-OLy,1,bi,bj,3),
                0618      U                     smTr(1-OLx,1-OLy,1,bi,bj,4),
                0619      U                     smTr(1-OLx,1-OLy,1,bi,bj,5),
                0620      U                     smTr(1-OLx,1-OLy,1,bi,bj,6),
                0621      U                     smTr(1-OLx,1-OLy,1,bi,bj,7),
                0622      U                     smTr(1-OLx,1-OLy,1,bi,bj,8),
                0623      U                     smTr(1-OLx,1-OLy,1,bi,bj,9),
                0624      U                     alp,   aln,   fp_v,  fn_v,  fp_o,  fn_o,
                0625      U                     fp_x,  fn_x,  fp_y,  fn_y,  fp_z,  fn_z,
                0626      U                     fp_xx, fn_xx, fp_yy, fn_yy, fp_zz, fn_zz,
                0627      U                     fp_xy, fn_xy, fp_xz, fn_xz, fp_yz, fn_yz,
                0628      O                     afr, myThid )
                0629         ELSE
b79a2b44f2 Jean*0630            STOP 'GAD_SOM_ADVECT: adv. scheme incompatibale with SOM'
d7ce0d34f8 Jean*0631         ENDIF
                0632
b79a2b44f2 Jean*0633 C--   Compute new tracer value and store tracer tendency
d7ce0d34f8 Jean*0634         DO j=1-OLy,sNy+OLy
                0635          DO i=1-OLx,sNx+OLx
f08c4b8486 Jean*0636 C--  without rescaling of tendencies:
                0637 c         localTr = smTr0(i,j,k)/smVol(i,j,k)
935a76deec Jean*0638 c         gTracer(i,j,k) = ( localTr - tracer(i,j,k,bi,bj) )
                0639 c    &                   / deltaTLev(k)
56c999258f Jean*0640 C--  Non-Lin Free-Surf: consistent with rescaling of tendencies
                0641 C     (in FREESURF_RESCALE_G) and RealFreshFlux/addMass.
                0642 C    Also valid for linear Free-Surf (r & r* coords) except that surf tracer
                0643 C    loss/gain is computed (in GAD_SOM_ADV_R) from partially updated tracer
                0644 C     (instead of from Tr^n as fresh-water dilution effect) resulting in
                0645 C    inaccurate linFSConserveTr and "surfExpan_" monitor.
935a76deec Jean*0646           gTracer(i,j,k) =
f08c4b8486 Jean*0647      &          ( smTr0(i,j,k) - tracer(i,j,k,bi,bj)*smVol(i,j,k) )
d7ce0d34f8 Jean*0648      &            *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
                0649      &            *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
                0650      &            *recip_rhoFacC(k)
4e66ab0b67 Oliv*0651      &            /deltaTLev(k)
d7ce0d34f8 Jean*0652          ENDDO
                0653         ENDDO
                0654
                0655 #ifdef ALLOW_DIAGNOSTICS
b79a2b44f2 Jean*0656         IF ( doDiagAdvR ) THEN
                0657          diagName = 'ADVr'//diagSufx
                0658          CALL DIAGNOSTICS_FILL( afr,
                0659      &                          diagName, k,1, 2,bi,bj, myThid )
d7ce0d34f8 Jean*0660         ENDIF
50d8304171 Ryan*0661 #ifdef ALLOW_LAYERS
                0662         IF ( useLayers ) THEN
                0663           CALL LAYERS_FILL(afr,tracerIdentity,'AFR',
                0664      &                     k,1,2,bi,bj,myThid)
                0665         ENDIF
                0666 #endif /* ALLOW_LAYERS */
d7ce0d34f8 Jean*0667 #endif
                0668
b79a2b44f2 Jean*0669 C--   End of k loop for vertical flux
d7ce0d34f8 Jean*0670        ENDDO
                0671 C--   end of if not.implicitAdvection block
                0672       ENDIF
7448700841 Mart*0673 #endif /* GAD_ALLOW_TS_SOM_ADV or PTRACERS_ALLOW_DYN_STATE */
d7ce0d34f8 Jean*0674
                0675       RETURN
                0676       END