pkg/generic_advdiff/gad_advection.F

C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_advection.F,v 1.24 2004/06/28 21:10:55 dimitri Exp $
C $Name:  $

#include "GAD_OPTIONS.h"

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C !ROUTINE: GAD_ADVECTION

C !INTERFACE: ==========================================================
      SUBROUTINE GAD_ADVECTION(
     I     implicitAdvection, advectionScheme, vertAdvecScheme,
     I     tracerIdentity,
     I     uVel, vVel, wVel, tracer,
     O     gTracer,
     I     bi,bj, myTime,myIter,myThid)

C !DESCRIPTION:
C Calculates the tendancy of a tracer due to advection.
C It uses the multi-dimensional method given in \ref{sect:multiDimAdvection}
C and can only be used for the non-linear advection schemes such as the
C direct-space-time method and flux-limiters. 
C
C The algorithm is as follows:
C \begin{itemize}
C \item{$\theta^{(n+1/3)} = \theta^{(n)}
C      - \Delta t \partial_x (u\theta^{(n)}) + \theta^{(n)} \partial_x u$}
C \item{$\theta^{(n+2/3)} = \theta^{(n+1/3)}
C      - \Delta t \partial_y (v\theta^{(n+1/3)}) + \theta^{(n)} \partial_y v$}
C \item{$\theta^{(n+3/3)} = \theta^{(n+2/3)}
C      - \Delta t \partial_r (w\theta^{(n+2/3)}) + \theta^{(n)} \partial_r w$}
C \item{$G_\theta = ( \theta^{(n+3/3)} - \theta^{(n)} )/\Delta t$}
C \end{itemize}
C
C The tendancy (output) is over-written by this routine.

C !USES: ===============================================================
      IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "GAD.h"
#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
# include "tamc_keys.h"
#endif
#ifdef ALLOW_EXCH2
#include "W2_EXCH2_TOPOLOGY.h"
#include "W2_EXCH2_PARAMS.h"
#endif /* ALLOW_EXCH2 */

C !INPUT PARAMETERS: ===================================================
C  implicitAdvection :: implicit vertical advection (later on)
C  advectionScheme   :: advection scheme to use (Horizontal plane)
C  vertAdvecScheme   :: advection scheme to use (vertical direction)
C  tracerIdentity    :: tracer identifier (required only for OBCS)
C  uVel              :: velocity, zonal component
C  vVel              :: velocity, meridional component
C  wVel              :: velocity, vertical component
C  tracer            :: tracer field
C  bi,bj             :: tile indices
C  myTime            :: current time
C  myIter            :: iteration number
C  myThid            :: thread number
      LOGICAL implicitAdvection
      INTEGER advectionScheme, vertAdvecScheme
      INTEGER tracerIdentity
      _RL uVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
      _RL vVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
      _RL wVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
      _RL tracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
      INTEGER bi,bj
      _RL myTime
      INTEGER myIter
      INTEGER myThid

C !OUTPUT PARAMETERS: ==================================================
C  gTracer           :: tendancy array
      _RL gTracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)

C !LOCAL VARIABLES: ====================================================
C  maskUp        :: 2-D array for mask at W points
C  iMin,iMax,    :: loop range for called routines
C  jMin,jMax     :: loop range for called routines
C  i,j,k         :: loop indices
C  kup           :: index into 2 1/2D array, toggles between 1 and 2
C  kdown         :: index into 2 1/2D array, toggles between 2 and 1
C  kp1           :: =k+1 for k<Nr, =Nr for k=Nr
C  xA,yA         :: areas of X and Y face of tracer cells
C  uTrans,vTrans :: 2-D arrays of volume transports at U,V points
C  rTrans        :: 2-D arrays of volume transports at W points
C  rTransKp1     :: vertical volume transport at interface k+1
C  af            :: 2-D array for horizontal advective flux
C  fVerT         :: 2 1/2D arrays for vertical advective flux
C  localTij      :: 2-D array, temporary local copy of tracer fld
C  localTijk     :: 3-D array, temporary local copy of tracer fld
C  kp1Msk        :: flag (0,1) for over-riding mask for W levels
C  calc_fluxes_X :: logical to indicate to calculate fluxes in X dir
C  calc_fluxes_Y :: logical to indicate to calculate fluxes in Y dir
C  nipass        :: number of passes in multi-dimensional method
C  ipass         :: number of the current pass being made
C  myTile        :: variables used to determine which cube face 
C  nCFace        :: owns a tile for cube grid runs using
C                :: multi-dim advection.
      _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      INTEGER iMin,iMax,jMin,jMax
      INTEGER i,j,k,kup,kDown
      _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL af      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
      _RL localTij(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL localTijk(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL kp1Msk
      LOGICAL calc_fluxes_X,calc_fluxes_Y
      INTEGER nipass,ipass
      INTEGER myTile, nCFace
      LOGICAL southWestCorner
      LOGICAL southEastCorner
      LOGICAL northWestCorner
      LOGICAL northEastCorner
CEOP

#ifdef ALLOW_AUTODIFF_TAMC
          act0 = tracerIdentity - 1
          max0 = maxpass
          act1 = bi - myBxLo(myThid)
          max1 = myBxHi(myThid) - myBxLo(myThid) + 1
          act2 = bj - myByLo(myThid)
          max2 = myByHi(myThid) - myByLo(myThid) + 1
          act3 = myThid - 1
          max3 = nTx*nTy
          act4 = ikey_dynamics - 1
          igadkey = (act0 + 1) 
     &                      + act1*max0
     &                      + act2*max0*max1
     &                      + act3*max0*max1*max2
     &                      + act4*max0*max1*max2*max3
          if (tracerIdentity.GT.maxpass) then
             print *, 'ph-pass gad_advection ', maxpass, tracerIdentity
             STOP 'maxpass seems smaller than tracerIdentity'
          endif
#endif /* ALLOW_AUTODIFF_TAMC */

C--   Set up work arrays with valid (i.e. not NaN) values
C     These inital values do not alter the numerical results. They
C     just ensure that all memory references are to valid floating
C     point numbers. This prevents spurious hardware signals due to
C     uninitialised but inert locations.
      DO j=1-OLy,sNy+OLy
       DO i=1-OLx,sNx+OLx
        xA(i,j)      = 0. _d 0
        yA(i,j)      = 0. _d 0
        uTrans(i,j)  = 0. _d 0
        vTrans(i,j)  = 0. _d 0
        rTrans(i,j)  = 0. _d 0
        fVerT(i,j,1) = 0. _d 0
        fVerT(i,j,2) = 0. _d 0
        rTransKp1(i,j)= 0. _d 0
       ENDDO
      ENDDO

      iMin = 1-OLx
      iMax = sNx+OLx
      jMin = 1-OLy
      jMax = sNy+OLy

C--   Start of k loop for horizontal fluxes
      DO k=1,Nr
#ifdef ALLOW_AUTODIFF_TAMC 
         kkey = (igadkey-1)*Nr + k
CADJ STORE tracer(:,:,k,bi,bj) = 
CADJ &     comlev1_bibj_k_gad, key=kkey, byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */

C--   Get temporary terms used by tendency routines
      CALL CALC_COMMON_FACTORS (
     I         bi,bj,iMin,iMax,jMin,jMax,k,
     O         xA,yA,uTrans,vTrans,rTrans,maskUp,
     I         myThid)

#ifdef ALLOW_GMREDI
C--   Residual transp = Bolus transp + Eulerian transp
       IF (useGMRedi)
     &   CALL GMREDI_CALC_UVFLOW(
     &            uTrans, vTrans, bi, bj, k, myThid)
#endif /* ALLOW_GMREDI */

C--   Make local copy of tracer array
      DO j=1-OLy,sNy+OLy
       DO i=1-OLx,sNx+OLx
        localTij(i,j)=tracer(i,j,k,bi,bj)
       ENDDO
      ENDDO

cph  The following block is needed for useCubedSphereExchange only,
cph  but needs to be set for all cases to avoid spurious 
cph  TAF dependencies
       southWestCorner = .TRUE.
       southEastCorner = .TRUE.
       northWestCorner = .TRUE.
       northEastCorner = .TRUE.
#ifdef ALLOW_EXCH2
       myTile = W2_myTileList(bi)
       nCFace = exch2_myFace(myTile)
       southWestCorner = exch2_isWedge(myTile).EQ.1
     &             .AND. exch2_isSedge(myTile).EQ.1
       southEastCorner = exch2_isEedge(myTile).EQ.1
     &             .AND. exch2_isSedge(myTile).EQ.1
       northEastCorner = exch2_isEedge(myTile).EQ.1
     &             .AND. exch2_isNedge(myTile).EQ.1
       northWestCorner = exch2_isWedge(myTile).EQ.1
     &             .AND. exch2_isNedge(myTile).EQ.1
#else
       nCFace = bi
#endif
      IF (useCubedSphereExchange) THEN

       nipass=3
#ifdef ALLOW_AUTODIFF_TAMC
       if ( nipass.GT.maxcube )
     &      STOP 'maxcube needs to be = 3'
#endif
      ELSE
       nipass=1
      ENDIF
cph       nipass=1

C--   Multiple passes for different directions on different tiles
C--   For cube need one pass for each of red, green and blue axes.
      DO ipass=1,nipass
#ifdef ALLOW_AUTODIFF_TAMC
         passkey = ipass + (k-1)      *maxcube
     &                   + (igadkey-1)*maxcube*Nr
         IF (nipass .GT. maxpass) THEN
          STOP 'GAD_ADVECTION: nipass > maxcube. check tamc.h'
         ENDIF
#endif /* ALLOW_AUTODIFF_TAMC */

      IF (nipass.EQ.3) THEN
       calc_fluxes_X=.FALSE.
       calc_fluxes_Y=.FALSE.
       IF (ipass.EQ.1 .AND. (nCFace.EQ.1 .OR. nCFace.EQ.2) ) THEN
        calc_fluxes_X=.TRUE.
       ELSEIF (ipass.EQ.1 .AND. (nCFace.EQ.4 .OR. nCFace.EQ.5) ) THEN
        calc_fluxes_Y=.TRUE.
       ELSEIF (ipass.EQ.2 .AND. (nCFace.EQ.1 .OR. nCFace.EQ.6) ) THEN
        calc_fluxes_Y=.TRUE.
       ELSEIF (ipass.EQ.2 .AND. (nCFace.EQ.3 .OR. nCFace.EQ.4) ) THEN
        calc_fluxes_X=.TRUE.
       ELSEIF (ipass.EQ.3 .AND. (nCFace.EQ.2 .OR. nCFace.EQ.3) ) THEN
        calc_fluxes_Y=.TRUE.
       ELSEIF (ipass.EQ.3 .AND. (nCFace.EQ.5 .OR. nCFace.EQ.6) ) THEN
        calc_fluxes_X=.TRUE.
       ENDIF
      ELSE
       calc_fluxes_X=.TRUE.
       calc_fluxes_Y=.TRUE.
      ENDIF
 
C--   X direction
      IF (calc_fluxes_X) THEN

C--   Internal exchange for calculations in X
      IF (useCubedSphereExchange) THEN
C--    For cube face corners we need to duplicate the
C--    i-1 and i+1 values into the null space as follows:
C
C
C      o NW corner: copy T(    0,sNy  ) into T(    0,sNy+1) e.g.
C                      |
C         x T(0,sNy+1) |
C        /\            |
C      --||------------|-----------
C        ||            |
C         x T(0,sNy)   |   x T(1,sNy)
C                      |
C
C      o SW corner: copy T(0,1) into T(0,0) e.g.
C                      |
C         x T(0,1)     |  x T(1,1)
C        ||            |
C      --||------------|-----------
C        \/            |
C         x T(0,0)     |
C                      |
C
C      o NE corner: copy T(sNx+1,sNy  ) into T(sNx+1,sNy+1) e.g.
C                      |
C                      |   x T(sNx+1,sNy+1)
C                      |  /\
C      ----------------|--||-------
C                      |  ||
C         x T(sNx,sNy) |   x T(sNx+1,sNy  )
C                      |
C      o SE corner: copy T(sNx+1,1    ) into T(sNx+1,0    ) e.g.
C                      |
C         x T(sNx,1)   |   x T(sNx+1,    1) 
C                      |  ||
C      ----------------|--||-------
C                      |  \/
C                      |   x T(sNx+1,    0)
       IF ( southWestCorner ) THEN
        localTij(0    ,0    )= localTij(0    ,1  )
       ENDIF
       IF ( southEastCorner ) THEN
        localTij(sNx+1,0    )= localTij(sNx+1,1  )
       ENDIF
       IF ( northWestCorner ) THEN
        localTij(0    ,sNy+1)= localTij(0    ,sNy)
       ENDIF
       IF ( northEastCorner ) THEN
        localTij(sNx+1,sNy+1)= localTij(sNx+1,sNy)
       ENDIF
      ENDIF

C-    Advective flux in X
      DO j=1-Oly,sNy+Oly
       DO i=1-Olx,sNx+Olx
        af(i,j) = 0.
       ENDDO
      ENDDO

#ifdef ALLOW_AUTODIFF_TAMC
#ifndef DISABLE_MULTIDIM_ADVECTION
CADJ STORE localTij(:,:)  = 
CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
#endif
#endif /* ALLOW_AUTODIFF_TAMC */

      IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
       CALL GAD_FLUXLIMIT_ADV_X(
     &      bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)
      ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
       CALL GAD_DST3_ADV_X(
     &       bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)
      ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
       CALL GAD_DST3FL_ADV_X(
     &       bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)
      ELSE
       STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'
      ENDIF

      DO j=1-Oly,sNy+Oly
       DO i=1-Olx,sNx+Olx-1
        localTij(i,j)=localTij(i,j)-deltaTtracer*
     &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
     &    *recip_rA(i,j,bi,bj)
     &    *( af(i+1,j)-af(i,j)
     &      -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
     &     )
       ENDDO
      ENDDO

#ifdef ALLOW_OBCS
C--   Apply open boundary conditions
      IF (useOBCS) THEN
       IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
        CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )
       ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
        CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
       END IF
      END IF
#endif /* ALLOW_OBCS */

C--   End of X direction
      ENDIF

C--   Y direction
      IF (calc_fluxes_Y) THEN

      IF (useCubedSphereExchange) THEN
C--   Internal exchange for calculations in Y
C--    For cube face corners we need to duplicate the
C--    j-1 and j+1 values into the null space as follows:
C
C      o SW corner: copy T(0,1) into T(0,0) e.g.
C                      |
C                      |  x T(1,1)
C                      |
C      ----------------|-----------
C                      |
C         x T(0,0)<====== x T(1,0)
C                      |
C
C      o NW corner: copy T(    0,sNy  ) into T(    0,sNy+1) e.g.
C                      |
C         x T(0,sNy+1)<=== x T(1,sNy+1)
C                      |
C      ----------------|-----------
C                      |
C                      |   x T(1,sNy)
C                      |
C
C      o NE corner: copy T(sNx+1,sNy  ) into T(sNx+1,sNy+1) e.g.
C                      |
C      x T(sNx,sNy+1)=====>x T(sNx+1,sNy+1)
C                      |     
C      ----------------|-----------
C                      |     
C      x T(sNx,sNy)    |                       
C                      |
C      o SE corner: copy T(sNx+1,1    ) into T(sNx+1,0    ) e.g.
C                      |
C         x T(sNx,1)   |                    
C                      |    
C      ----------------|-----------
C                      |     
C         x T(sNx,0) =====>x T(sNx+1,    0)
       IF ( southWestCorner ) THEN
         localTij(    0,0    ) = localTij(  1,0    )
       ENDIF
       IF ( southEastCorner ) THEN
         localTij(sNx+1,0    ) = localTij(sNx,0    )
       ENDIF
       IF ( northWestCorner ) THEN
         localTij(0    ,sNy+1) = localTij(  1,sNy+1)
       ENDIF
       IF ( northEastCorner ) THEN
         localTij(sNx+1,sNy+1) = localTij(sNx,sNy+1)
       ENDIF
      ENDIF

C-    Advective flux in Y
      DO j=1-Oly,sNy+Oly
       DO i=1-Olx,sNx+Olx
        af(i,j) = 0.
       ENDDO
      ENDDO

#ifdef ALLOW_AUTODIFF_TAMC 
#ifndef DISABLE_MULTIDIM_ADVECTION
CADJ STORE localTij(:,:)  = 
CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
#endif
#endif /* ALLOW_AUTODIFF_TAMC */

      IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
       CALL GAD_FLUXLIMIT_ADV_Y(
     &       bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid)
      ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
       CALL GAD_DST3_ADV_Y(
     &       bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid)
      ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
       CALL GAD_DST3FL_ADV_Y(
     &       bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid)
      ELSE
       STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
      ENDIF

      DO j=1-Oly,sNy+Oly-1
       DO i=1-Olx,sNx+Olx
        localTij(i,j)=localTij(i,j)-deltaTtracer*
     &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
     &    *recip_rA(i,j,bi,bj)
     &    *( af(i,j+1)-af(i,j)
     &      -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
     &     )
       ENDDO
      ENDDO

#ifdef ALLOW_OBCS
C--   Apply open boundary conditions
      IF (useOBCS) THEN
       IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
        CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )
       ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
        CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
       END IF
      END IF
#endif /* ALLOW_OBCS */

C--   End of Y direction
      ENDIF

C--   End of ipass loop
      ENDDO

      IF ( implicitAdvection ) THEN
C-    explicit advection is done ; store tendency in gTracer:
        DO j=1-Oly,sNy+Oly
         DO i=1-Olx,sNx+Olx
          gTracer(i,j,k,bi,bj)=
     &     (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer
         ENDDO
        ENDDO
      ELSE
C-    horizontal advection done; store intermediate result in 3D array:
       DO j=1-Oly,sNy+Oly
        DO i=1-Olx,sNx+Olx
         localTijk(i,j,k)=localTij(i,j)
        ENDDO
       ENDDO
      ENDIF

C--   End of K loop for horizontal fluxes
      ENDDO

      IF ( .NOT.implicitAdvection ) THEN
C--   Start of k loop for vertical flux
       DO k=Nr,1,-1
#ifdef ALLOW_AUTODIFF_TAMC 
         kkey = (igadkey-1)*Nr + k
#endif /* ALLOW_AUTODIFF_TAMC */
C--   kup    Cycles through 1,2 to point to w-layer above
C--   kDown  Cycles through 2,1 to point to w-layer below
        kup  = 1+MOD(k+1,2)
        kDown= 1+MOD(k,2)
c       kp1=min(Nr,k+1)
        kp1Msk=1.
        if (k.EQ.Nr) kp1Msk=0.

C-- Compute Vertical transport
#ifdef ALLOW_AIM
C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
        IF ( k.EQ.1 .OR.
     &     (useAIM .AND. tracerIdentity.EQ.GAD_SALINITY .AND. k.EQ.Nr)
     &              ) THEN
#else
        IF ( k.EQ.1 ) THEN
#endif

C- Surface interface :
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
           rTransKp1(i,j) = kp1Msk*rTrans(i,j)
           rTrans(i,j) = 0.
           fVerT(i,j,kUp) = 0.
           af(i,j) = 0.
          ENDDO
         ENDDO

        ELSE
C- Interior interface :

         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
           rTransKp1(i,j) = kp1Msk*rTrans(i,j)
           rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)
     &                 *maskC(i,j,k-1,bi,bj)
           af(i,j) = 0.
          ENDDO
         ENDDO

#ifdef ALLOW_GMREDI
C--   Residual transp = Bolus transp + Eulerian transp
         IF (useGMRedi) 
     &   CALL GMREDI_CALC_WFLOW(
     &                    rTrans, bi, bj, k, myThid)
#endif /* ALLOW_GMREDI */

#ifdef ALLOW_AUTODIFF_TAMC 
CADJ STORE localTijk(:,:,k)  
CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
CADJ STORE rTrans(:,:)  
CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */

C-    Compute vertical advective flux in the interior:
         IF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
          CALL GAD_FLUXLIMIT_ADV_R(
     &        bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
         ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
          CALL GAD_DST3_ADV_R(
     &        bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
         ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
          CALL GAD_DST3FL_ADV_R(
     &        bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
         ELSE
          STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
         ENDIF
C-    add the advective flux to fVerT
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
           fVerT(i,j,kUp) = af(i,j)
          ENDDO
         ENDDO

C- end Surface/Interior if bloc
        ENDIF

#ifdef ALLOW_AUTODIFF_TAMC 
CADJ STORE rTrans(:,:)  
CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
CADJ STORE rTranskp1(:,:)  
CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */

C--   Divergence of vertical fluxes
        DO j=1-Oly,sNy+Oly
         DO i=1-Olx,sNx+Olx
          localTij(i,j)=localTijk(i,j,k)-deltaTtracer*
     &     _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
     &     *recip_rA(i,j,bi,bj)
     &     *( fVerT(i,j,kUp)-fVerT(i,j,kDown)
     &       -tracer(i,j,k,bi,bj)*(rTrans(i,j)-rTransKp1(i,j))
     &      )*rkFac
          gTracer(i,j,k,bi,bj)=
     &     (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer
         ENDDO
        ENDDO
 
C--   End of K loop for vertical flux
       ENDDO
C--   end of if not.implicitAdvection block
      ENDIF 

      RETURN
      END
1	C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_advection.F,v 1.24 2004/06/28 21:10:55 dimitri Exp $
2	C $Name: $
3
4	#include "GAD_OPTIONS.h"
5
6	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
7	CBOP
8	C !ROUTINE: GAD_ADVECTION
9
10	C !INTERFACE: ==========================================================
11	SUBROUTINE GAD_ADVECTION(
12	I implicitAdvection, advectionScheme, vertAdvecScheme,
13	I tracerIdentity,
14	I uVel, vVel, wVel, tracer,
15	O gTracer,
16	I bi,bj, myTime,myIter,myThid)
17
18	C !DESCRIPTION:
19	C Calculates the tendancy of a tracer due to advection.
20	C It uses the multi-dimensional method given in \ref{sect:multiDimAdvection}
21	C and can only be used for the non-linear advection schemes such as the
22	C direct-space-time method and flux-limiters.
23	C
24	C The algorithm is as follows:
25	C \begin{itemize}
26	C \item{$\theta^{(n+1/3)} = \theta^{(n)}
27	C - \Delta t \partial_x (u\theta^{(n)}) + \theta^{(n)} \partial_x u$}
28	C \item{$\theta^{(n+2/3)} = \theta^{(n+1/3)}
29	C - \Delta t \partial_y (v\theta^{(n+1/3)}) + \theta^{(n)} \partial_y v$}
30	C \item{$\theta^{(n+3/3)} = \theta^{(n+2/3)}
31	C - \Delta t \partial_r (w\theta^{(n+2/3)}) + \theta^{(n)} \partial_r w$}
32	C \item{$G_\theta = ( \theta^{(n+3/3)} - \theta^{(n)} )/\Delta t$}
33	C \end{itemize}
34	C
35	C The tendancy (output) is over-written by this routine.
36
37	C !USES: ===============================================================
38	IMPLICIT NONE
39	#include "SIZE.h"
40	#include "EEPARAMS.h"
41	#include "PARAMS.h"
42	#include "GRID.h"
43	#include "GAD.h"
44	#ifdef ALLOW_AUTODIFF_TAMC
45	# include "tamc.h"
46	# include "tamc_keys.h"
47	#endif
48	#ifdef ALLOW_EXCH2
49	#include "W2_EXCH2_TOPOLOGY.h"
50	#include "W2_EXCH2_PARAMS.h"
51	#endif /* ALLOW_EXCH2 */
52
53	C !INPUT PARAMETERS: ===================================================
54	C implicitAdvection :: implicit vertical advection (later on)
55	C advectionScheme :: advection scheme to use (Horizontal plane)
56	C vertAdvecScheme :: advection scheme to use (vertical direction)
57	C tracerIdentity :: tracer identifier (required only for OBCS)
58	C uVel :: velocity, zonal component
59	C vVel :: velocity, meridional component
60	C wVel :: velocity, vertical component
61	C tracer :: tracer field
62	C bi,bj :: tile indices
63	C myTime :: current time
64	C myIter :: iteration number
65	C myThid :: thread number
66	LOGICAL implicitAdvection
67	INTEGER advectionScheme, vertAdvecScheme
68	INTEGER tracerIdentity
69	_RL uVel (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
70	_RL vVel (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
71	_RL wVel (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
72	_RL tracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
73	INTEGER bi,bj
74	_RL myTime
75	INTEGER myIter
76	INTEGER myThid
77
78	C !OUTPUT PARAMETERS: ==================================================
79	C gTracer :: tendancy array
80	_RL gTracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
81
82	C !LOCAL VARIABLES: ====================================================
83	C maskUp :: 2-D array for mask at W points
84	C iMin,iMax, :: loop range for called routines
85	C jMin,jMax :: loop range for called routines
86	C i,j,k :: loop indices
87	C kup :: index into 2 1/2D array, toggles between 1 and 2
88	C kdown :: index into 2 1/2D array, toggles between 2 and 1
89	C kp1 :: =k+1 for k<Nr, =Nr for k=Nr
90	C xA,yA :: areas of X and Y face of tracer cells
91	C uTrans,vTrans :: 2-D arrays of volume transports at U,V points
92	C rTrans :: 2-D arrays of volume transports at W points
93	C rTransKp1 :: vertical volume transport at interface k+1
94	C af :: 2-D array for horizontal advective flux
95	C fVerT :: 2 1/2D arrays for vertical advective flux
96	C localTij :: 2-D array, temporary local copy of tracer fld
97	C localTijk :: 3-D array, temporary local copy of tracer fld
98	C kp1Msk :: flag (0,1) for over-riding mask for W levels
99	C calc_fluxes_X :: logical to indicate to calculate fluxes in X dir
100	C calc_fluxes_Y :: logical to indicate to calculate fluxes in Y dir
101	C nipass :: number of passes in multi-dimensional method
102	C ipass :: number of the current pass being made
103	C myTile :: variables used to determine which cube face
104	C nCFace :: owns a tile for cube grid runs using
105	C :: multi-dim advection.
106	_RS maskUp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
107	INTEGER iMin,iMax,jMin,jMax
108	INTEGER i,j,k,kup,kDown
109	_RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
110	_RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
111	_RL uTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
112	_RL vTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
113	_RL rTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
114	_RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
115	_RL af (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116	_RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
117	_RL localTij(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
118	_RL localTijk(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
119	_RL kp1Msk
120	LOGICAL calc_fluxes_X,calc_fluxes_Y
121	INTEGER nipass,ipass
122	INTEGER myTile, nCFace
123	LOGICAL southWestCorner
124	LOGICAL southEastCorner
125	LOGICAL northWestCorner
126	LOGICAL northEastCorner
127	CEOP
128
129	#ifdef ALLOW_AUTODIFF_TAMC
130	act0 = tracerIdentity - 1
131	max0 = maxpass
132	act1 = bi - myBxLo(myThid)
133	max1 = myBxHi(myThid) - myBxLo(myThid) + 1
134	act2 = bj - myByLo(myThid)
135	max2 = myByHi(myThid) - myByLo(myThid) + 1
136	act3 = myThid - 1
137	max3 = nTx*nTy
138	act4 = ikey_dynamics - 1
139	igadkey = (act0 + 1)
140	& + act1*max0
141	& + act2max0max1
142	& + act3max0max1*max2
143	& + act4max0max1max2max3
144	if (tracerIdentity.GT.maxpass) then
145	print *, 'ph-pass gad_advection ', maxpass, tracerIdentity
146	STOP 'maxpass seems smaller than tracerIdentity'
147	endif
148	#endif /* ALLOW_AUTODIFF_TAMC */
149
150	C-- Set up work arrays with valid (i.e. not NaN) values
151	C These inital values do not alter the numerical results. They
152	C just ensure that all memory references are to valid floating
153	C point numbers. This prevents spurious hardware signals due to
154	C uninitialised but inert locations.
155	DO j=1-OLy,sNy+OLy
156	DO i=1-OLx,sNx+OLx
157	xA(i,j) = 0. _d 0
158	yA(i,j) = 0. _d 0
159	uTrans(i,j) = 0. _d 0
160	vTrans(i,j) = 0. _d 0
161	rTrans(i,j) = 0. _d 0
162	fVerT(i,j,1) = 0. _d 0
163	fVerT(i,j,2) = 0. _d 0
164	rTransKp1(i,j)= 0. _d 0
165	ENDDO
166	ENDDO
167
168	iMin = 1-OLx
169	iMax = sNx+OLx
170	jMin = 1-OLy
171	jMax = sNy+OLy
172
173	C-- Start of k loop for horizontal fluxes
174	DO k=1,Nr
175	#ifdef ALLOW_AUTODIFF_TAMC
176	kkey = (igadkey-1)*Nr + k
177	CADJ STORE tracer(:,:,k,bi,bj) =
178	CADJ & comlev1_bibj_k_gad, key=kkey, byte=isbyte
179	#endif /* ALLOW_AUTODIFF_TAMC */
180
181	C-- Get temporary terms used by tendency routines
182	CALL CALC_COMMON_FACTORS (
183	I bi,bj,iMin,iMax,jMin,jMax,k,
184	O xA,yA,uTrans,vTrans,rTrans,maskUp,
185	I myThid)
186
187	#ifdef ALLOW_GMREDI
188	C-- Residual transp = Bolus transp + Eulerian transp
189	IF (useGMRedi)
190	& CALL GMREDI_CALC_UVFLOW(
191	& uTrans, vTrans, bi, bj, k, myThid)
192	#endif /* ALLOW_GMREDI */
193
194	C-- Make local copy of tracer array
195	DO j=1-OLy,sNy+OLy
196	DO i=1-OLx,sNx+OLx
197	localTij(i,j)=tracer(i,j,k,bi,bj)
198	ENDDO
199	ENDDO
200
201	cph The following block is needed for useCubedSphereExchange only,
202	cph but needs to be set for all cases to avoid spurious
203	cph TAF dependencies
204	southWestCorner = .TRUE.
205	southEastCorner = .TRUE.
206	northWestCorner = .TRUE.
207	northEastCorner = .TRUE.
208	#ifdef ALLOW_EXCH2
209	myTile = W2_myTileList(bi)
210	nCFace = exch2_myFace(myTile)
211	southWestCorner = exch2_isWedge(myTile).EQ.1
212	& .AND. exch2_isSedge(myTile).EQ.1
213	southEastCorner = exch2_isEedge(myTile).EQ.1
214	& .AND. exch2_isSedge(myTile).EQ.1
215	northEastCorner = exch2_isEedge(myTile).EQ.1
216	& .AND. exch2_isNedge(myTile).EQ.1
217	northWestCorner = exch2_isWedge(myTile).EQ.1
218	& .AND. exch2_isNedge(myTile).EQ.1
219	#else
220	nCFace = bi
221	#endif
222	IF (useCubedSphereExchange) THEN
223
224	nipass=3
225	#ifdef ALLOW_AUTODIFF_TAMC
226	if ( nipass.GT.maxcube )
227	& STOP 'maxcube needs to be = 3'
228	#endif
229	ELSE
230	nipass=1
231	ENDIF
232	cph nipass=1
233
234	C-- Multiple passes for different directions on different tiles
235	C-- For cube need one pass for each of red, green and blue axes.
236	DO ipass=1,nipass
237	#ifdef ALLOW_AUTODIFF_TAMC
238	passkey = ipass + (k-1) *maxcube
239	& + (igadkey-1)maxcubeNr
240	IF (nipass .GT. maxpass) THEN
241	STOP 'GAD_ADVECTION: nipass > maxcube. check tamc.h'
242	ENDIF
243	#endif /* ALLOW_AUTODIFF_TAMC */
244
245	IF (nipass.EQ.3) THEN
246	calc_fluxes_X=.FALSE.
247	calc_fluxes_Y=.FALSE.
248	IF (ipass.EQ.1 .AND. (nCFace.EQ.1 .OR. nCFace.EQ.2) ) THEN
249	calc_fluxes_X=.TRUE.
250	ELSEIF (ipass.EQ.1 .AND. (nCFace.EQ.4 .OR. nCFace.EQ.5) ) THEN
251	calc_fluxes_Y=.TRUE.
252	ELSEIF (ipass.EQ.2 .AND. (nCFace.EQ.1 .OR. nCFace.EQ.6) ) THEN
253	calc_fluxes_Y=.TRUE.
254	ELSEIF (ipass.EQ.2 .AND. (nCFace.EQ.3 .OR. nCFace.EQ.4) ) THEN
255	calc_fluxes_X=.TRUE.
256	ELSEIF (ipass.EQ.3 .AND. (nCFace.EQ.2 .OR. nCFace.EQ.3) ) THEN
257	calc_fluxes_Y=.TRUE.
258	ELSEIF (ipass.EQ.3 .AND. (nCFace.EQ.5 .OR. nCFace.EQ.6) ) THEN
259	calc_fluxes_X=.TRUE.
260	ENDIF
261	ELSE
262	calc_fluxes_X=.TRUE.
263	calc_fluxes_Y=.TRUE.
264	ENDIF
265
266	C-- X direction
267	IF (calc_fluxes_X) THEN
268
269	C-- Internal exchange for calculations in X
270	IF (useCubedSphereExchange) THEN
271	C-- For cube face corners we need to duplicate the
272	C-- i-1 and i+1 values into the null space as follows:
273	C
274	C
275	C o NW corner: copy T( 0,sNy ) into T( 0,sNy+1) e.g.
276	C \|
277	C x T(0,sNy+1) \|
278	C /\ \|
279	C --\|\|------------\|-----------
280	C \|\| \|
281	C x T(0,sNy) \| x T(1,sNy)
282	C \|
283	C
284	C o SW corner: copy T(0,1) into T(0,0) e.g.
285	C \|
286	C x T(0,1) \| x T(1,1)
287	C \|\| \|
288	C --\|\|------------\|-----------
289	C \/ \|
290	C x T(0,0) \|
291	C \|
292	C
293	C o NE corner: copy T(sNx+1,sNy ) into T(sNx+1,sNy+1) e.g.
294	C \|
295	C \| x T(sNx+1,sNy+1)
296	C \| /\
297	C ----------------\|--\|\|-------
298	C \| \|\|
299	C x T(sNx,sNy) \| x T(sNx+1,sNy )
300	C \|
301	C o SE corner: copy T(sNx+1,1 ) into T(sNx+1,0 ) e.g.
302	C \|
303	C x T(sNx,1) \| x T(sNx+1, 1)
304	C \| \|\|
305	C ----------------\|--\|\|-------
306	C \| \/
307	C \| x T(sNx+1, 0)
308	IF ( southWestCorner ) THEN
309	localTij(0 ,0 )= localTij(0 ,1 )
310	ENDIF
311	IF ( southEastCorner ) THEN
312	localTij(sNx+1,0 )= localTij(sNx+1,1 )
313	ENDIF
314	IF ( northWestCorner ) THEN
315	localTij(0 ,sNy+1)= localTij(0 ,sNy)
316	ENDIF
317	IF ( northEastCorner ) THEN
318	localTij(sNx+1,sNy+1)= localTij(sNx+1,sNy)
319	ENDIF
320	ENDIF
321
322	C- Advective flux in X
323	DO j=1-Oly,sNy+Oly
324	DO i=1-Olx,sNx+Olx
325	af(i,j) = 0.
326	ENDDO
327	ENDDO
328
329	#ifdef ALLOW_AUTODIFF_TAMC
330	#ifndef DISABLE_MULTIDIM_ADVECTION
331	CADJ STORE localTij(:,:) =
332	CADJ & comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
333	#endif
334	#endif /* ALLOW_AUTODIFF_TAMC */
335
336	IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
337	CALL GAD_FLUXLIMIT_ADV_X(
338	& bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)
339	ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
340	CALL GAD_DST3_ADV_X(
341	& bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)
342	ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
343	CALL GAD_DST3FL_ADV_X(
344	& bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)
345	ELSE
346	STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'
347	ENDIF
348
349	DO j=1-Oly,sNy+Oly
350	DO i=1-Olx,sNx+Olx-1
351	localTij(i,j)=localTij(i,j)-deltaTtracer*
352	& _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
353	& *recip_rA(i,j,bi,bj)
354	& *( af(i+1,j)-af(i,j)
355	& -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
356	& )
357	ENDDO
358	ENDDO
359
360	#ifdef ALLOW_OBCS
361	C-- Apply open boundary conditions
362	IF (useOBCS) THEN
363	IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
364	CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )
365	ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
366	CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
367	END IF
368	END IF
369	#endif /* ALLOW_OBCS */
370
371	C-- End of X direction
372	ENDIF
373
374	C-- Y direction
375	IF (calc_fluxes_Y) THEN
376
377	IF (useCubedSphereExchange) THEN
378	C-- Internal exchange for calculations in Y
379	C-- For cube face corners we need to duplicate the
380	C-- j-1 and j+1 values into the null space as follows:
381	C
382	C o SW corner: copy T(0,1) into T(0,0) e.g.
383	C \|
384	C \| x T(1,1)
385	C \|
386	C ----------------\|-----------
387	C \|
388	C x T(0,0)<====== x T(1,0)
389	C \|
390	C
391	C o NW corner: copy T( 0,sNy ) into T( 0,sNy+1) e.g.
392	C \|
393	C x T(0,sNy+1)<=== x T(1,sNy+1)
394	C \|
395	C ----------------\|-----------
396	C \|
397	C \| x T(1,sNy)
398	C \|
399	C
400	C o NE corner: copy T(sNx+1,sNy ) into T(sNx+1,sNy+1) e.g.
401	C \|
402	C x T(sNx,sNy+1)=====>x T(sNx+1,sNy+1)
403	C \|
404	C ----------------\|-----------
405	C \|
406	C x T(sNx,sNy) \|
407	C \|
408	C o SE corner: copy T(sNx+1,1 ) into T(sNx+1,0 ) e.g.
409	C \|
410	C x T(sNx,1) \|
411	C \|
412	C ----------------\|-----------
413	C \|
414	C x T(sNx,0) =====>x T(sNx+1, 0)
415	IF ( southWestCorner ) THEN
416	localTij( 0,0 ) = localTij( 1,0 )
417	ENDIF
418	IF ( southEastCorner ) THEN
419	localTij(sNx+1,0 ) = localTij(sNx,0 )
420	ENDIF
421	IF ( northWestCorner ) THEN
422	localTij(0 ,sNy+1) = localTij( 1,sNy+1)
423	ENDIF
424	IF ( northEastCorner ) THEN
425	localTij(sNx+1,sNy+1) = localTij(sNx,sNy+1)
426	ENDIF
427	ENDIF
428
429	C- Advective flux in Y
430	DO j=1-Oly,sNy+Oly
431	DO i=1-Olx,sNx+Olx
432	af(i,j) = 0.
433	ENDDO
434	ENDDO
435
436	#ifdef ALLOW_AUTODIFF_TAMC
437	#ifndef DISABLE_MULTIDIM_ADVECTION
438	CADJ STORE localTij(:,:) =
439	CADJ & comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
440	#endif
441	#endif /* ALLOW_AUTODIFF_TAMC */
442
443	IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
444	CALL GAD_FLUXLIMIT_ADV_Y(
445	& bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid)
446	ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
447	CALL GAD_DST3_ADV_Y(
448	& bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid)
449	ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
450	CALL GAD_DST3FL_ADV_Y(
451	& bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid)
452	ELSE
453	STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
454	ENDIF
455
456	DO j=1-Oly,sNy+Oly-1
457	DO i=1-Olx,sNx+Olx
458	localTij(i,j)=localTij(i,j)-deltaTtracer*
459	& _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
460	& *recip_rA(i,j,bi,bj)
461	& *( af(i,j+1)-af(i,j)
462	& -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
463	& )
464	ENDDO
465	ENDDO
466
467	#ifdef ALLOW_OBCS
468	C-- Apply open boundary conditions
469	IF (useOBCS) THEN
470	IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
471	CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )
472	ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
473	CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
474	END IF
475	END IF
476	#endif /* ALLOW_OBCS */
477
478	C-- End of Y direction
479	ENDIF
480
481	C-- End of ipass loop
482	ENDDO
483
484	IF ( implicitAdvection ) THEN
485	C- explicit advection is done ; store tendency in gTracer:
486	DO j=1-Oly,sNy+Oly
487	DO i=1-Olx,sNx+Olx
488	gTracer(i,j,k,bi,bj)=
489	& (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer
490	ENDDO
491	ENDDO
492	ELSE
493	C- horizontal advection done; store intermediate result in 3D array:
494	DO j=1-Oly,sNy+Oly
495	DO i=1-Olx,sNx+Olx
496	localTijk(i,j,k)=localTij(i,j)
497	ENDDO
498	ENDDO
499	ENDIF
500
501	C-- End of K loop for horizontal fluxes
502	ENDDO
503
504	IF ( .NOT.implicitAdvection ) THEN
505	C-- Start of k loop for vertical flux
506	DO k=Nr,1,-1
507	#ifdef ALLOW_AUTODIFF_TAMC
508	kkey = (igadkey-1)*Nr + k
509	#endif /* ALLOW_AUTODIFF_TAMC */
510	C-- kup Cycles through 1,2 to point to w-layer above
511	C-- kDown Cycles through 2,1 to point to w-layer below
512	kup = 1+MOD(k+1,2)
513	kDown= 1+MOD(k,2)
514	c kp1=min(Nr,k+1)
515	kp1Msk=1.
516	if (k.EQ.Nr) kp1Msk=0.
517
518	C-- Compute Vertical transport
519	#ifdef ALLOW_AIM
520	C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
521	IF ( k.EQ.1 .OR.
522	& (useAIM .AND. tracerIdentity.EQ.GAD_SALINITY .AND. k.EQ.Nr)
523	& ) THEN
524	#else
525	IF ( k.EQ.1 ) THEN
526	#endif
527
528	C- Surface interface :
529	DO j=1-Oly,sNy+Oly
530	DO i=1-Olx,sNx+Olx
531	rTransKp1(i,j) = kp1Msk*rTrans(i,j)
532	rTrans(i,j) = 0.
533	fVerT(i,j,kUp) = 0.
534	af(i,j) = 0.
535	ENDDO
536	ENDDO
537
538	ELSE
539	C- Interior interface :
540
541	DO j=1-Oly,sNy+Oly
542	DO i=1-Olx,sNx+Olx
543	rTransKp1(i,j) = kp1Msk*rTrans(i,j)
544	rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)
545	& *maskC(i,j,k-1,bi,bj)
546	af(i,j) = 0.
547	ENDDO
548	ENDDO
549
550	#ifdef ALLOW_GMREDI
551	C-- Residual transp = Bolus transp + Eulerian transp
552	IF (useGMRedi)
553	& CALL GMREDI_CALC_WFLOW(
554	& rTrans, bi, bj, k, myThid)
555	#endif /* ALLOW_GMREDI */
556
557	#ifdef ALLOW_AUTODIFF_TAMC
558	CADJ STORE localTijk(:,:,k)
559	CADJ & = comlev1_bibj_k_gad, key=kkey, byte=isbyte
560	CADJ STORE rTrans(:,:)
561	CADJ & = comlev1_bibj_k_gad, key=kkey, byte=isbyte
562	#endif /* ALLOW_AUTODIFF_TAMC */
563
564	C- Compute vertical advective flux in the interior:
565	IF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
566	CALL GAD_FLUXLIMIT_ADV_R(
567	& bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
568	ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
569	CALL GAD_DST3_ADV_R(
570	& bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
571	ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
572	CALL GAD_DST3FL_ADV_R(
573	& bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
574	ELSE
575	STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
576	ENDIF
577	C- add the advective flux to fVerT
578	DO j=1-Oly,sNy+Oly
579	DO i=1-Olx,sNx+Olx
580	fVerT(i,j,kUp) = af(i,j)
581	ENDDO
582	ENDDO
583
584	C- end Surface/Interior if bloc
585	ENDIF
586
587	#ifdef ALLOW_AUTODIFF_TAMC
588	CADJ STORE rTrans(:,:)
589	CADJ & = comlev1_bibj_k_gad, key=kkey, byte=isbyte
590	CADJ STORE rTranskp1(:,:)
591	CADJ & = comlev1_bibj_k_gad, key=kkey, byte=isbyte
592	#endif /* ALLOW_AUTODIFF_TAMC */
593
594	C-- Divergence of vertical fluxes
595	DO j=1-Oly,sNy+Oly
596	DO i=1-Olx,sNx+Olx
597	localTij(i,j)=localTijk(i,j,k)-deltaTtracer*
598	& _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
599	& *recip_rA(i,j,bi,bj)
600	& *( fVerT(i,j,kUp)-fVerT(i,j,kDown)
601	& -tracer(i,j,k,bi,bj)*(rTrans(i,j)-rTransKp1(i,j))
602	& )*rkFac
603	gTracer(i,j,k,bi,bj)=
604	& (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer
605	ENDDO
606	ENDDO
607
608	C-- End of K loop for vertical flux
609	ENDDO
610	C-- end of if not.implicitAdvection block
611	ENDIF
612
613	RETURN
614	END