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