pkg/generic_advdiff/gad_advection.F

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