pkg/generic_advdiff/gad_advection.F

C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_advection.F,v 1.23 2004/06/26 02:38:54 jmc 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

      IF (useCubedSphereExchange) THEN
       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

       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.23 2004/06/26 02:38:54 jmc 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	IF (useCubedSphereExchange) THEN
202	southWestCorner = .TRUE.
203	southEastCorner = .TRUE.
204	northWestCorner = .TRUE.
205	northEastCorner = .TRUE.
206	#ifdef ALLOW_EXCH2
207	myTile = W2_myTileList(bi)
208	nCFace = exch2_myFace(myTile)
209	southWestCorner = exch2_isWedge(myTile).EQ.1
210	& .AND. exch2_isSedge(myTile).EQ.1
211	southEastCorner = exch2_isEedge(myTile).EQ.1
212	& .AND. exch2_isSedge(myTile).EQ.1
213	northEastCorner = exch2_isEedge(myTile).EQ.1
214	& .AND. exch2_isNedge(myTile).EQ.1
215	northWestCorner = exch2_isWedge(myTile).EQ.1
216	& .AND. exch2_isNedge(myTile).EQ.1
217	#else
218	nCFace = bi
219	#endif
220
221	nipass=3
222	#ifdef ALLOW_AUTODIFF_TAMC
223	if ( nipass.GT.maxcube )
224	& STOP 'maxcube needs to be = 3'
225	#endif
226	ELSE
227	nipass=1
228	ENDIF
229	cph nipass=1
230
231	C-- Multiple passes for different directions on different tiles
232	C-- For cube need one pass for each of red, green and blue axes.
233	DO ipass=1,nipass
234	#ifdef ALLOW_AUTODIFF_TAMC
235	passkey = ipass + (k-1) *maxcube
236	& + (igadkey-1)maxcubeNr
237	IF (nipass .GT. maxpass) THEN
238	STOP 'GAD_ADVECTION: nipass > maxcube. check tamc.h'
239	ENDIF
240	#endif /* ALLOW_AUTODIFF_TAMC */
241
242	IF (nipass.EQ.3) THEN
243	calc_fluxes_X=.FALSE.
244	calc_fluxes_Y=.FALSE.
245	IF (ipass.EQ.1 .AND. (nCFace.EQ.1 .OR. nCFace.EQ.2) ) THEN
246	calc_fluxes_X=.TRUE.
247	ELSEIF (ipass.EQ.1 .AND. (nCFace.EQ.4 .OR. nCFace.EQ.5) ) THEN
248	calc_fluxes_Y=.TRUE.
249	ELSEIF (ipass.EQ.2 .AND. (nCFace.EQ.1 .OR. nCFace.EQ.6) ) THEN
250	calc_fluxes_Y=.TRUE.
251	ELSEIF (ipass.EQ.2 .AND. (nCFace.EQ.3 .OR. nCFace.EQ.4) ) THEN
252	calc_fluxes_X=.TRUE.
253	ELSEIF (ipass.EQ.3 .AND. (nCFace.EQ.2 .OR. nCFace.EQ.3) ) THEN
254	calc_fluxes_Y=.TRUE.
255	ELSEIF (ipass.EQ.3 .AND. (nCFace.EQ.5 .OR. nCFace.EQ.6) ) THEN
256	calc_fluxes_X=.TRUE.
257	ENDIF
258	ELSE
259	calc_fluxes_X=.TRUE.
260	calc_fluxes_Y=.TRUE.
261	ENDIF
262
263	C-- X direction
264	IF (calc_fluxes_X) THEN
265
266	C-- Internal exchange for calculations in X
267	IF (useCubedSphereExchange) THEN
268	C-- For cube face corners we need to duplicate the
269	C-- i-1 and i+1 values into the null space as follows:
270	C
271	C
272	C o NW corner: copy T( 0,sNy ) into T( 0,sNy+1) e.g.
273	C \|
274	C x T(0,sNy+1) \|
275	C /\ \|
276	C --\|\|------------\|-----------
277	C \|\| \|
278	C x T(0,sNy) \| x T(1,sNy)
279	C \|
280	C
281	C o SW corner: copy T(0,1) into T(0,0) e.g.
282	C \|
283	C x T(0,1) \| x T(1,1)
284	C \|\| \|
285	C --\|\|------------\|-----------
286	C \/ \|
287	C x T(0,0) \|
288	C \|
289	C
290	C o NE corner: copy T(sNx+1,sNy ) into T(sNx+1,sNy+1) e.g.
291	C \|
292	C \| x T(sNx+1,sNy+1)
293	C \| /\
294	C ----------------\|--\|\|-------
295	C \| \|\|
296	C x T(sNx,sNy) \| x T(sNx+1,sNy )
297	C \|
298	C o SE corner: copy T(sNx+1,1 ) into T(sNx+1,0 ) e.g.
299	C \|
300	C x T(sNx,1) \| x T(sNx+1, 1)
301	C \| \|\|
302	C ----------------\|--\|\|-------
303	C \| \/
304	C \| x T(sNx+1, 0)
305	IF ( southWestCorner ) THEN
306	localTij(0 ,0 )= localTij(0 ,1 )
307	ENDIF
308	IF ( southEastCorner ) THEN
309	localTij(sNx+1,0 )= localTij(sNx+1,1 )
310	ENDIF
311	IF ( northWestCorner ) THEN
312	localTij(0 ,sNy+1)= localTij(0 ,sNy)
313	ENDIF
314	IF ( northEastCorner ) THEN
315	localTij(sNx+1,sNy+1)= localTij(sNx+1,sNy)
316	ENDIF
317	ENDIF
318
319	C- Advective flux in X
320	DO j=1-Oly,sNy+Oly
321	DO i=1-Olx,sNx+Olx
322	af(i,j) = 0.
323	ENDDO
324	ENDDO
325
326	#ifdef ALLOW_AUTODIFF_TAMC
327	#ifndef DISABLE_MULTIDIM_ADVECTION
328	CADJ STORE localTij(:,:) =
329	CADJ & comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
330	#endif
331	#endif /* ALLOW_AUTODIFF_TAMC */
332
333	IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
334	CALL GAD_FLUXLIMIT_ADV_X(
335	& bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)
336	ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
337	CALL GAD_DST3_ADV_X(
338	& bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)
339	ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
340	CALL GAD_DST3FL_ADV_X(
341	& bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)
342	ELSE
343	STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'
344	ENDIF
345
346	DO j=1-Oly,sNy+Oly
347	DO i=1-Olx,sNx+Olx-1
348	localTij(i,j)=localTij(i,j)-deltaTtracer*
349	& _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
350	& *recip_rA(i,j,bi,bj)
351	& *( af(i+1,j)-af(i,j)
352	& -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
353	& )
354	ENDDO
355	ENDDO
356
357	#ifdef ALLOW_OBCS
358	C-- Apply open boundary conditions
359	IF (useOBCS) THEN
360	IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
361	CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )
362	ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
363	CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
364	END IF
365	END IF
366	#endif /* ALLOW_OBCS */
367
368	C-- End of X direction
369	ENDIF
370
371	C-- Y direction
372	IF (calc_fluxes_Y) THEN
373
374	IF (useCubedSphereExchange) THEN
375	C-- Internal exchange for calculations in Y
376	C-- For cube face corners we need to duplicate the
377	C-- j-1 and j+1 values into the null space as follows:
378	C
379	C o SW corner: copy T(0,1) into T(0,0) e.g.
380	C \|
381	C \| x T(1,1)
382	C \|
383	C ----------------\|-----------
384	C \|
385	C x T(0,0)<====== x T(1,0)
386	C \|
387	C
388	C o NW corner: copy T( 0,sNy ) into T( 0,sNy+1) e.g.
389	C \|
390	C x T(0,sNy+1)<=== x T(1,sNy+1)
391	C \|
392	C ----------------\|-----------
393	C \|
394	C \| x T(1,sNy)
395	C \|
396	C
397	C o NE corner: copy T(sNx+1,sNy ) into T(sNx+1,sNy+1) e.g.
398	C \|
399	C x T(sNx,sNy+1)=====>x T(sNx+1,sNy+1)
400	C \|
401	C ----------------\|-----------
402	C \|
403	C x T(sNx,sNy) \|
404	C \|
405	C o SE corner: copy T(sNx+1,1 ) into T(sNx+1,0 ) e.g.
406	C \|
407	C x T(sNx,1) \|
408	C \|
409	C ----------------\|-----------
410	C \|
411	C x T(sNx,0) =====>x T(sNx+1, 0)
412	IF ( southWestCorner ) THEN
413	localTij( 0,0 ) = localTij( 1,0 )
414	ENDIF
415	IF ( southEastCorner ) THEN
416	localTij(sNx+1,0 ) = localTij(sNx,0 )
417	ENDIF
418	IF ( northWestCorner ) THEN
419	localTij(0 ,sNy+1) = localTij( 1,sNy+1)
420	ENDIF
421	IF ( northEastCorner ) THEN
422	localTij(sNx+1,sNy+1) = localTij(sNx,sNy+1)
423	ENDIF
424	ENDIF
425
426	C- Advective flux in Y
427	DO j=1-Oly,sNy+Oly
428	DO i=1-Olx,sNx+Olx
429	af(i,j) = 0.
430	ENDDO
431	ENDDO
432
433	#ifdef ALLOW_AUTODIFF_TAMC
434	#ifndef DISABLE_MULTIDIM_ADVECTION
435	CADJ STORE localTij(:,:) =
436	CADJ & comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
437	#endif
438	#endif /* ALLOW_AUTODIFF_TAMC */
439
440	IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
441	CALL GAD_FLUXLIMIT_ADV_Y(
442	& bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid)
443	ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
444	CALL GAD_DST3_ADV_Y(
445	& bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid)
446	ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
447	CALL GAD_DST3FL_ADV_Y(
448	& bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid)
449	ELSE
450	STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
451	ENDIF
452
453	DO j=1-Oly,sNy+Oly-1
454	DO i=1-Olx,sNx+Olx
455	localTij(i,j)=localTij(i,j)-deltaTtracer*
456	& _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
457	& *recip_rA(i,j,bi,bj)
458	& *( af(i,j+1)-af(i,j)
459	& -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
460	& )
461	ENDDO
462	ENDDO
463
464	#ifdef ALLOW_OBCS
465	C-- Apply open boundary conditions
466	IF (useOBCS) THEN
467	IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
468	CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )
469	ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
470	CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
471	END IF
472	END IF
473	#endif /* ALLOW_OBCS */
474
475	C-- End of Y direction
476	ENDIF
477
478	C-- End of ipass loop
479	ENDDO
480
481	IF ( implicitAdvection ) THEN
482	C- explicit advection is done ; store tendency in gTracer:
483	DO j=1-Oly,sNy+Oly
484	DO i=1-Olx,sNx+Olx
485	gTracer(i,j,k,bi,bj)=
486	& (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer
487	ENDDO
488	ENDDO
489	ELSE
490	C- horizontal advection done; store intermediate result in 3D array:
491	DO j=1-Oly,sNy+Oly
492	DO i=1-Olx,sNx+Olx
493	localTijk(i,j,k)=localTij(i,j)
494	ENDDO
495	ENDDO
496	ENDIF
497
498	C-- End of K loop for horizontal fluxes
499	ENDDO
500
501	IF ( .NOT.implicitAdvection ) THEN
502	C-- Start of k loop for vertical flux
503	DO k=Nr,1,-1
504	#ifdef ALLOW_AUTODIFF_TAMC
505	kkey = (igadkey-1)*Nr + k
506	#endif /* ALLOW_AUTODIFF_TAMC */
507	C-- kup Cycles through 1,2 to point to w-layer above
508	C-- kDown Cycles through 2,1 to point to w-layer below
509	kup = 1+MOD(k+1,2)
510	kDown= 1+MOD(k,2)
511	c kp1=min(Nr,k+1)
512	kp1Msk=1.
513	if (k.EQ.Nr) kp1Msk=0.
514
515	C-- Compute Vertical transport
516	#ifdef ALLOW_AIM
517	C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
518	IF ( k.EQ.1 .OR.
519	& (useAIM .AND. tracerIdentity.EQ.GAD_SALINITY .AND. k.EQ.Nr)
520	& ) THEN
521	#else
522	IF ( k.EQ.1 ) THEN
523	#endif
524
525	C- Surface interface :
526	DO j=1-Oly,sNy+Oly
527	DO i=1-Olx,sNx+Olx
528	rTransKp1(i,j) = kp1Msk*rTrans(i,j)
529	rTrans(i,j) = 0.
530	fVerT(i,j,kUp) = 0.
531	af(i,j) = 0.
532	ENDDO
533	ENDDO
534
535	ELSE
536	C- Interior interface :
537
538	DO j=1-Oly,sNy+Oly
539	DO i=1-Olx,sNx+Olx
540	rTransKp1(i,j) = kp1Msk*rTrans(i,j)
541	rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)
542	& *maskC(i,j,k-1,bi,bj)
543	af(i,j) = 0.
544	ENDDO
545	ENDDO
546
547	#ifdef ALLOW_GMREDI
548	C-- Residual transp = Bolus transp + Eulerian transp
549	IF (useGMRedi)
550	& CALL GMREDI_CALC_WFLOW(
551	& rTrans, bi, bj, k, myThid)
552	#endif /* ALLOW_GMREDI */
553
554	#ifdef ALLOW_AUTODIFF_TAMC
555	CADJ STORE localTijk(:,:,k)
556	CADJ & = comlev1_bibj_k_gad, key=kkey, byte=isbyte
557	CADJ STORE rTrans(:,:)
558	CADJ & = comlev1_bibj_k_gad, key=kkey, byte=isbyte
559	#endif /* ALLOW_AUTODIFF_TAMC */
560
561	C- Compute vertical advective flux in the interior:
562	IF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
563	CALL GAD_FLUXLIMIT_ADV_R(
564	& bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
565	ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
566	CALL GAD_DST3_ADV_R(
567	& bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
568	ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
569	CALL GAD_DST3FL_ADV_R(
570	& bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
571	ELSE
572	STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
573	ENDIF
574	C- add the advective flux to fVerT
575	DO j=1-Oly,sNy+Oly
576	DO i=1-Olx,sNx+Olx
577	fVerT(i,j,kUp) = af(i,j)
578	ENDDO
579	ENDDO
580
581	C- end Surface/Interior if bloc
582	ENDIF
583
584	#ifdef ALLOW_AUTODIFF_TAMC
585	CADJ STORE rTrans(:,:)
586	CADJ & = comlev1_bibj_k_gad, key=kkey, byte=isbyte
587	CADJ STORE rTranskp1(:,:)
588	CADJ & = comlev1_bibj_k_gad, key=kkey, byte=isbyte
589	#endif /* ALLOW_AUTODIFF_TAMC */
590
591	C-- Divergence of vertical fluxes
592	DO j=1-Oly,sNy+Oly
593	DO i=1-Olx,sNx+Olx
594	localTij(i,j)=localTijk(i,j,k)-deltaTtracer*
595	& _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
596	& *recip_rA(i,j,bi,bj)
597	& *( fVerT(i,j,kUp)-fVerT(i,j,kDown)
598	& -tracer(i,j,k,bi,bj)*(rTrans(i,j)-rTransKp1(i,j))
599	& )*rkFac
600	gTracer(i,j,k,bi,bj)=
601	& (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer
602	ENDDO
603	ENDDO
604
605	C-- End of K loop for vertical flux
606	ENDDO
607	C-- end of if not.implicitAdvection block
608	ENDIF
609
610	RETURN
611	END