pkg/generic_advdiff/gad_advection.F

C $Header: /u/u0/gcmpack/MITgcm/pkg/generic_advdiff/gad_advection.F,v 1.25 2004/06/30 23:45:35 heimbach 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
        DO J=1,OLy
         DO I=1,OLx
          localTij(1-I, 1-J   )= localTij(1-J  ,1  )
         ENDDO
        ENDDO
       ENDIF
       IF ( southEastCorner ) THEN
        DO J=1,OLy
         DO I=1,OLx
          localTij(sNx+I, 1-J )=localTij(sNx+J, I  )
         ENDDO
        ENDDO
       ENDIF
       IF ( northWestCorner ) THEN
        DO J=1,OLy
         DO I=1,OLx
          localTij( 1-I ,sNy+J)=localTij( 1-J , sNy+1-I )
         ENDDO
        ENDDO
       ENDIF
       IF ( northEastCorner ) THEN
        DO J=1,OLy
         DO I=1,OLx
          localTij(sNx+I,sNy+J)=localTij(sNx+J, sNy+1-I )
         ENDDO
        ENDDO
       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
        DO J=1,Oly
         DO I=1,Olx
          localTij( 1-i , 1-j ) = localTij(j   , 1-i )
         ENDDO
        ENDDO
       ENDIF
       IF ( southEastCorner ) THEN
        DO J=1,Oly
         DO I=1,Olx
          localTij(sNx+i, 1-j ) = localTij(sNx+1-j, 1-i )
         ENDDO
        ENDDO
       ENDIF
       IF ( northWestCorner ) THEN
        DO J=1,Oly
         DO I=1,Olx
          localTij( 1-i ,sNy+j) = localTij(j   ,sNy+i)
         ENDDO
        ENDDO
       ENDIF
       IF ( northEastCorner ) THEN
        DO J=1,Oly
         DO I=1,Olx
          localTij(sNx+i,sNy+j) = localTij(sNx+1-j,sNy+i)
         ENDDO
        ENDDO
       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	cnh	1.26	C $Header: /u/u0/gcmpack/MITgcm/pkg/generic_advdiff/gad_advection.F,v 1.25 2004/06/30 23:45:35 heimbach Exp $
2	adcroft	1.2	C $Name: $
3	adcroft	1.4
4	adcroft	1.1	#include "GAD_OPTIONS.h"
5
6	edhill	1.19	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
7	adcroft	1.4	CBOP
8			C !ROUTINE: GAD_ADVECTION
9
10			C !INTERFACE: ==========================================================
11	jmc	1.17	SUBROUTINE GAD_ADVECTION(
12	jmc	1.23	I implicitAdvection, advectionScheme, vertAdvecScheme,
13			I tracerIdentity,
14	edhill	1.21	I uVel, vVel, wVel, tracer,
15			O gTracer,
16			I bi,bj, myTime,myIter,myThid)
17	adcroft	1.4
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	adcroft	1.5	C - \Delta t \partial_x (u\theta^{(n)}) + \theta^{(n)} \partial_x u$}
28	adcroft	1.4	C \item{$\theta^{(n+2/3)} = \theta^{(n+1/3)}
29	adcroft	1.5	C - \Delta t \partial_y (v\theta^{(n+1/3)}) + \theta^{(n)} \partial_y v$}
30	adcroft	1.4	C \item{$\theta^{(n+3/3)} = \theta^{(n+2/3)}
31	adcroft	1.5	C - \Delta t \partial_r (w\theta^{(n+2/3)}) + \theta^{(n)} \partial_r w$}
32	adcroft	1.4	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	adcroft	1.1	IMPLICIT NONE
39			#include "SIZE.h"
40			#include "EEPARAMS.h"
41			#include "PARAMS.h"
42			#include "GRID.h"
43			#include "GAD.h"
44	heimbach	1.6	#ifdef ALLOW_AUTODIFF_TAMC
45			# include "tamc.h"
46			# include "tamc_keys.h"
47			#endif
48	dimitri	1.24	#ifdef ALLOW_EXCH2
49			#include "W2_EXCH2_TOPOLOGY.h"
50			#include "W2_EXCH2_PARAMS.h"
51			#endif /* ALLOW_EXCH2 */
52	adcroft	1.1
53	adcroft	1.4	C !INPUT PARAMETERS: ===================================================
54	edhill	1.21	C implicitAdvection :: implicit vertical advection (later on)
55	jmc	1.23	C advectionScheme :: advection scheme to use (Horizontal plane)
56			C vertAdvecScheme :: advection scheme to use (vertical direction)
57	edhill	1.21	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	jmc	1.17	LOGICAL implicitAdvection
67	jmc	1.23	INTEGER advectionScheme, vertAdvecScheme
68	adcroft	1.1	INTEGER tracerIdentity
69	jmc	1.17	_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	adcroft	1.1	_RL myTime
75			INTEGER myIter
76			INTEGER myThid
77
78	adcroft	1.4	C !OUTPUT PARAMETERS: ==================================================
79	edhill	1.21	C gTracer :: tendancy array
80	adcroft	1.9	_RL gTracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
81	adcroft	1.4
82			C !LOCAL VARIABLES: ====================================================
83	edhill	1.21	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	dimitri	1.24	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	adcroft	1.1	_RS maskUp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
107			INTEGER iMin,iMax,jMin,jMax
108	jmc	1.11	INTEGER i,j,k,kup,kDown
109	adcroft	1.1	_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	jmc	1.11	_RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
115	adcroft	1.1	_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	adcroft	1.3	LOGICAL calc_fluxes_X,calc_fluxes_Y
121			INTEGER nipass,ipass
122	dimitri	1.24	INTEGER myTile, nCFace
123			LOGICAL southWestCorner
124			LOGICAL southEastCorner
125			LOGICAL northWestCorner
126			LOGICAL northEastCorner
127	adcroft	1.4	CEOP
128	adcroft	1.1
129	heimbach	1.6	#ifdef ALLOW_AUTODIFF_TAMC
130	heimbach	1.14	act0 = tracerIdentity - 1
131			max0 = maxpass
132	heimbach	1.6	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	heimbach	1.14	igadkey = (act0 + 1)
140			& + act1*max0
141			& + act2max0max1
142			& + act3max0max1*max2
143			& + act4max0max1max2max3
144	heimbach	1.15	if (tracerIdentity.GT.maxpass) then
145			print *, 'ph-pass gad_advection ', maxpass, tracerIdentity
146			STOP 'maxpass seems smaller than tracerIdentity'
147			endif
148	heimbach	1.6	#endif /* ALLOW_AUTODIFF_TAMC */
149
150	adcroft	1.1	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	jmc	1.11	rTransKp1(i,j)= 0. _d 0
165	adcroft	1.1	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	heimbach	1.6	#ifdef ALLOW_AUTODIFF_TAMC
176	heimbach	1.14	kkey = (igadkey-1)*Nr + k
177			CADJ STORE tracer(:,:,k,bi,bj) =
178			CADJ & comlev1_bibj_k_gad, key=kkey, byte=isbyte
179	heimbach	1.6	#endif /* ALLOW_AUTODIFF_TAMC */
180	adcroft	1.1
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	jmc	1.11	#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	adcroft	1.1	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	heimbach	1.25	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	dimitri	1.24	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	heimbach	1.25	IF (useCubedSphereExchange) THEN
223	dimitri	1.24
224	adcroft	1.3	nipass=3
225	heimbach	1.14	#ifdef ALLOW_AUTODIFF_TAMC
226			if ( nipass.GT.maxcube )
227			& STOP 'maxcube needs to be = 3'
228			#endif
229	adcroft	1.3	ELSE
230			nipass=1
231			ENDIF
232	heimbach	1.6	cph nipass=1
233	adcroft	1.3
234			C-- Multiple passes for different directions on different tiles
235	dimitri	1.24	C-- For cube need one pass for each of red, green and blue axes.
236	adcroft	1.3	DO ipass=1,nipass
237	heimbach	1.6	#ifdef ALLOW_AUTODIFF_TAMC
238	heimbach	1.14	passkey = ipass + (k-1) *maxcube
239			& + (igadkey-1)maxcubeNr
240	heimbach	1.6	IF (nipass .GT. maxpass) THEN
241	heimbach	1.14	STOP 'GAD_ADVECTION: nipass > maxcube. check tamc.h'
242	heimbach	1.6	ENDIF
243			#endif /* ALLOW_AUTODIFF_TAMC */
244	adcroft	1.3
245			IF (nipass.EQ.3) THEN
246			calc_fluxes_X=.FALSE.
247			calc_fluxes_Y=.FALSE.
248	dimitri	1.24	IF (ipass.EQ.1 .AND. (nCFace.EQ.1 .OR. nCFace.EQ.2) ) THEN
249	adcroft	1.3	calc_fluxes_X=.TRUE.
250	dimitri	1.24	ELSEIF (ipass.EQ.1 .AND. (nCFace.EQ.4 .OR. nCFace.EQ.5) ) THEN
251	adcroft	1.3	calc_fluxes_Y=.TRUE.
252	dimitri	1.24	ELSEIF (ipass.EQ.2 .AND. (nCFace.EQ.1 .OR. nCFace.EQ.6) ) THEN
253	adcroft	1.3	calc_fluxes_Y=.TRUE.
254	dimitri	1.24	ELSEIF (ipass.EQ.2 .AND. (nCFace.EQ.3 .OR. nCFace.EQ.4) ) THEN
255	adcroft	1.3	calc_fluxes_X=.TRUE.
256	dimitri	1.24	ELSEIF (ipass.EQ.3 .AND. (nCFace.EQ.2 .OR. nCFace.EQ.3) ) THEN
257	adcroft	1.3	calc_fluxes_Y=.TRUE.
258	dimitri	1.24	ELSEIF (ipass.EQ.3 .AND. (nCFace.EQ.5 .OR. nCFace.EQ.6) ) THEN
259	adcroft	1.3	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	dimitri	1.24	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	cnh	1.26	DO J=1,OLy
310			DO I=1,OLx
311			localTij(1-I, 1-J )= localTij(1-J ,1 )
312			ENDDO
313			ENDDO
314	dimitri	1.24	ENDIF
315			IF ( southEastCorner ) THEN
316	cnh	1.26	DO J=1,OLy
317			DO I=1,OLx
318			localTij(sNx+I, 1-J )=localTij(sNx+J, I )
319			ENDDO
320			ENDDO
321	dimitri	1.24	ENDIF
322			IF ( northWestCorner ) THEN
323	cnh	1.26	DO J=1,OLy
324			DO I=1,OLx
325			localTij( 1-I ,sNy+J)=localTij( 1-J , sNy+1-I )
326			ENDDO
327			ENDDO
328	dimitri	1.24	ENDIF
329			IF ( northEastCorner ) THEN
330	cnh	1.26	DO J=1,OLy
331			DO I=1,OLx
332			localTij(sNx+I,sNy+J)=localTij(sNx+J, sNy+1-I )
333			ENDDO
334			ENDDO
335	dimitri	1.24	ENDIF
336	adcroft	1.3	ENDIF
337
338	adcroft	1.1	C- Advective flux in X
339			DO j=1-Oly,sNy+Oly
340			DO i=1-Olx,sNx+Olx
341			af(i,j) = 0.
342			ENDDO
343			ENDDO
344	heimbach	1.6
345			#ifdef ALLOW_AUTODIFF_TAMC
346	adcroft	1.7	#ifndef DISABLE_MULTIDIM_ADVECTION
347	heimbach	1.14	CADJ STORE localTij(:,:) =
348			CADJ & comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
349	heimbach	1.6	#endif
350			#endif /* ALLOW_AUTODIFF_TAMC */
351
352	adcroft	1.1	IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
353			CALL GAD_FLUXLIMIT_ADV_X(
354			& bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)
355			ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
356			CALL GAD_DST3_ADV_X(
357			& bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)
358			ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
359			CALL GAD_DST3FL_ADV_X(
360			& bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)
361			ELSE
362	adcroft	1.9	STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'
363	adcroft	1.1	ENDIF
364	heimbach	1.6
365	adcroft	1.1	DO j=1-Oly,sNy+Oly
366			DO i=1-Olx,sNx+Olx-1
367			localTij(i,j)=localTij(i,j)-deltaTtracer*
368			& _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
369			& *recip_rA(i,j,bi,bj)
370			& *( af(i+1,j)-af(i,j)
371			& -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
372			& )
373			ENDDO
374			ENDDO
375
376			#ifdef ALLOW_OBCS
377			C-- Apply open boundary conditions
378			IF (useOBCS) THEN
379			IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
380			CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )
381			ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
382			CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
383			END IF
384			END IF
385			#endif /* ALLOW_OBCS */
386
387	adcroft	1.3	C-- End of X direction
388			ENDIF
389
390			C-- Y direction
391			IF (calc_fluxes_Y) THEN
392
393	dimitri	1.24	IF (useCubedSphereExchange) THEN
394	adcroft	1.3	C-- Internal exchange for calculations in Y
395	dimitri	1.24	C-- For cube face corners we need to duplicate the
396			C-- j-1 and j+1 values into the null space as follows:
397			C
398			C o SW corner: copy T(0,1) into T(0,0) e.g.
399			C \|
400			C \| x T(1,1)
401			C \|
402			C ----------------\|-----------
403			C \|
404			C x T(0,0)<====== x T(1,0)
405			C \|
406			C
407			C o NW corner: copy T( 0,sNy ) into T( 0,sNy+1) e.g.
408			C \|
409			C x T(0,sNy+1)<=== x T(1,sNy+1)
410			C \|
411			C ----------------\|-----------
412			C \|
413			C \| x T(1,sNy)
414			C \|
415			C
416			C o NE corner: copy T(sNx+1,sNy ) into T(sNx+1,sNy+1) e.g.
417			C \|
418			C x T(sNx,sNy+1)=====>x T(sNx+1,sNy+1)
419			C \|
420			C ----------------\|-----------
421			C \|
422			C x T(sNx,sNy) \|
423			C \|
424			C o SE corner: copy T(sNx+1,1 ) into T(sNx+1,0 ) e.g.
425			C \|
426			C x T(sNx,1) \|
427			C \|
428			C ----------------\|-----------
429			C \|
430			C x T(sNx,0) =====>x T(sNx+1, 0)
431			IF ( southWestCorner ) THEN
432	cnh	1.26	DO J=1,Oly
433			DO I=1,Olx
434			localTij( 1-i , 1-j ) = localTij(j , 1-i )
435			ENDDO
436			ENDDO
437	dimitri	1.24	ENDIF
438			IF ( southEastCorner ) THEN
439	cnh	1.26	DO J=1,Oly
440			DO I=1,Olx
441			localTij(sNx+i, 1-j ) = localTij(sNx+1-j, 1-i )
442			ENDDO
443			ENDDO
444	dimitri	1.24	ENDIF
445			IF ( northWestCorner ) THEN
446	cnh	1.26	DO J=1,Oly
447			DO I=1,Olx
448			localTij( 1-i ,sNy+j) = localTij(j ,sNy+i)
449			ENDDO
450			ENDDO
451	dimitri	1.24	ENDIF
452			IF ( northEastCorner ) THEN
453	cnh	1.26	DO J=1,Oly
454			DO I=1,Olx
455			localTij(sNx+i,sNy+j) = localTij(sNx+1-j,sNy+i)
456			ENDDO
457			ENDDO
458	dimitri	1.24	ENDIF
459	adcroft	1.3	ENDIF
460
461	adcroft	1.1	C- Advective flux in Y
462			DO j=1-Oly,sNy+Oly
463			DO i=1-Olx,sNx+Olx
464			af(i,j) = 0.
465			ENDDO
466			ENDDO
467	heimbach	1.6
468			#ifdef ALLOW_AUTODIFF_TAMC
469	adcroft	1.7	#ifndef DISABLE_MULTIDIM_ADVECTION
470	heimbach	1.14	CADJ STORE localTij(:,:) =
471			CADJ & comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
472	heimbach	1.6	#endif
473			#endif /* ALLOW_AUTODIFF_TAMC */
474
475	adcroft	1.1	IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
476			CALL GAD_FLUXLIMIT_ADV_Y(
477			& bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid)
478			ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
479			CALL GAD_DST3_ADV_Y(
480			& bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid)
481			ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
482			CALL GAD_DST3FL_ADV_Y(
483			& bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid)
484			ELSE
485			STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
486			ENDIF
487	heimbach	1.6
488	adcroft	1.1	DO j=1-Oly,sNy+Oly-1
489			DO i=1-Olx,sNx+Olx
490			localTij(i,j)=localTij(i,j)-deltaTtracer*
491			& _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
492			& *recip_rA(i,j,bi,bj)
493			& *( af(i,j+1)-af(i,j)
494			& -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
495			& )
496			ENDDO
497			ENDDO
498	adcroft	1.3
499	adcroft	1.1	#ifdef ALLOW_OBCS
500			C-- Apply open boundary conditions
501			IF (useOBCS) THEN
502			IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
503			CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )
504			ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
505			CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
506			END IF
507			END IF
508			#endif /* ALLOW_OBCS */
509	adcroft	1.3
510			C-- End of Y direction
511			ENDIF
512
513	jmc	1.18	C-- End of ipass loop
514	adcroft	1.1	ENDDO
515
516	jmc	1.18	IF ( implicitAdvection ) THEN
517			C- explicit advection is done ; store tendency in gTracer:
518			DO j=1-Oly,sNy+Oly
519			DO i=1-Olx,sNx+Olx
520			gTracer(i,j,k,bi,bj)=
521			& (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer
522			ENDDO
523			ENDDO
524			ELSE
525			C- horizontal advection done; store intermediate result in 3D array:
526			DO j=1-Oly,sNy+Oly
527			DO i=1-Olx,sNx+Olx
528			localTijk(i,j,k)=localTij(i,j)
529			ENDDO
530			ENDDO
531			ENDIF
532	adcroft	1.1
533			C-- End of K loop for horizontal fluxes
534			ENDDO
535
536	jmc	1.18	IF ( .NOT.implicitAdvection ) THEN
537	adcroft	1.1	C-- Start of k loop for vertical flux
538	jmc	1.18	DO k=Nr,1,-1
539	heimbach	1.6	#ifdef ALLOW_AUTODIFF_TAMC
540	heimbach	1.16	kkey = (igadkey-1)*Nr + k
541	heimbach	1.6	#endif /* ALLOW_AUTODIFF_TAMC */
542	adcroft	1.1	C-- kup Cycles through 1,2 to point to w-layer above
543			C-- kDown Cycles through 2,1 to point to w-layer below
544	jmc	1.18	kup = 1+MOD(k+1,2)
545			kDown= 1+MOD(k,2)
546			c kp1=min(Nr,k+1)
547			kp1Msk=1.
548			if (k.EQ.Nr) kp1Msk=0.
549	heimbach	1.6
550	jmc	1.11	C-- Compute Vertical transport
551	jmc	1.22	#ifdef ALLOW_AIM
552			C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
553			IF ( k.EQ.1 .OR.
554			& (useAIM .AND. tracerIdentity.EQ.GAD_SALINITY .AND. k.EQ.Nr)
555			& ) THEN
556			#else
557			IF ( k.EQ.1 ) THEN
558			#endif
559	jmc	1.11
560			C- Surface interface :
561	jmc	1.18	DO j=1-Oly,sNy+Oly
562			DO i=1-Olx,sNx+Olx
563	jmc	1.22	rTransKp1(i,j) = kp1Msk*rTrans(i,j)
564	jmc	1.18	rTrans(i,j) = 0.
565			fVerT(i,j,kUp) = 0.
566			af(i,j) = 0.
567			ENDDO
568			ENDDO
569	jmc	1.11
570	jmc	1.18	ELSE
571			C- Interior interface :
572	jmc	1.11
573	jmc	1.18	DO j=1-Oly,sNy+Oly
574			DO i=1-Olx,sNx+Olx
575			rTransKp1(i,j) = kp1Msk*rTrans(i,j)
576			rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)
577			& *maskC(i,j,k-1,bi,bj)
578			af(i,j) = 0.
579			ENDDO
580			ENDDO
581	jmc	1.11
582			#ifdef ALLOW_GMREDI
583			C-- Residual transp = Bolus transp + Eulerian transp
584	jmc	1.18	IF (useGMRedi)
585	jmc	1.11	& CALL GMREDI_CALC_WFLOW(
586			& rTrans, bi, bj, k, myThid)
587			#endif /* ALLOW_GMREDI */
588
589	heimbach	1.16	#ifdef ALLOW_AUTODIFF_TAMC
590			CADJ STORE localTijk(:,:,k)
591			CADJ & = comlev1_bibj_k_gad, key=kkey, byte=isbyte
592			CADJ STORE rTrans(:,:)
593			CADJ & = comlev1_bibj_k_gad, key=kkey, byte=isbyte
594			#endif /* ALLOW_AUTODIFF_TAMC */
595
596	adcroft	1.1	C- Compute vertical advective flux in the interior:
597	jmc	1.23	IF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
598	jmc	1.18	CALL GAD_FLUXLIMIT_ADV_R(
599	jmc	1.17	& bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
600	jmc	1.23	ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
601	jmc	1.18	CALL GAD_DST3_ADV_R(
602	jmc	1.17	& bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
603	jmc	1.23	ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
604	jmc	1.18	CALL GAD_DST3FL_ADV_R(
605	jmc	1.17	& bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
606	jmc	1.18	ELSE
607			STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
608			ENDIF
609	jmc	1.11	C- add the advective flux to fVerT
610	jmc	1.18	DO j=1-Oly,sNy+Oly
611			DO i=1-Olx,sNx+Olx
612			fVerT(i,j,kUp) = af(i,j)
613			ENDDO
614			ENDDO
615	jmc	1.11
616			C- end Surface/Interior if bloc
617	jmc	1.18	ENDIF
618	heimbach	1.16
619			#ifdef ALLOW_AUTODIFF_TAMC
620			CADJ STORE rTrans(:,:)
621			CADJ & = comlev1_bibj_k_gad, key=kkey, byte=isbyte
622			CADJ STORE rTranskp1(:,:)
623			CADJ & = comlev1_bibj_k_gad, key=kkey, byte=isbyte
624			#endif /* ALLOW_AUTODIFF_TAMC */
625	adcroft	1.1
626	jmc	1.18	C-- Divergence of vertical fluxes
627			DO j=1-Oly,sNy+Oly
628			DO i=1-Olx,sNx+Olx
629			localTij(i,j)=localTijk(i,j,k)-deltaTtracer*
630			& _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
631			& *recip_rA(i,j,bi,bj)
632			& *( fVerT(i,j,kUp)-fVerT(i,j,kDown)
633			& -tracer(i,j,k,bi,bj)*(rTrans(i,j)-rTransKp1(i,j))
634			& )*rkFac
635			gTracer(i,j,k,bi,bj)=
636			& (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer
637			ENDDO
638			ENDDO
639	adcroft	1.1
640			C-- End of K loop for vertical flux
641	jmc	1.18	ENDDO
642			C-- end of if not.implicitAdvection block
643			ENDIF
644	adcroft	1.1
645			RETURN
646			END