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	C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_advection.F,v 1.27 2004/09/17 23:02:01 heimbach 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	# ifdef ALLOW_PTRACERS
48	# include "PTRACERS_SIZE.h"
49	# endif
50	#endif
51	#ifdef ALLOW_EXCH2
52	#include "W2_EXCH2_TOPOLOGY.h"
53	#include "W2_EXCH2_PARAMS.h"
54	#endif /* ALLOW_EXCH2 */
55
56	C !INPUT PARAMETERS: ===================================================
57	C implicitAdvection :: implicit vertical advection (later on)
58	C advectionScheme :: advection scheme to use (Horizontal plane)
59	C vertAdvecScheme :: advection scheme to use (vertical direction)
60	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	LOGICAL implicitAdvection
70	INTEGER advectionScheme, vertAdvecScheme
71	INTEGER tracerIdentity
72	_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	_RL myTime
78	INTEGER myIter
79	INTEGER myThid
80
81	C !OUTPUT PARAMETERS: ==================================================
82	C gTracer :: tendancy array
83	_RL gTracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
84
85	C !LOCAL VARIABLES: ====================================================
86	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	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	_RS maskUp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
110	INTEGER iMin,iMax,jMin,jMax
111	INTEGER i,j,k,kup,kDown
112	_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	_RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
118	_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	LOGICAL calc_fluxes_X,calc_fluxes_Y
124	INTEGER nipass,ipass
125	INTEGER myTile, nCFace
126	LOGICAL southWestCorner
127	LOGICAL southEastCorner
128	LOGICAL northWestCorner
129	LOGICAL northEastCorner
130	CEOP
131
132	#ifdef ALLOW_AUTODIFF_TAMC
133	act0 = tracerIdentity - 1
134	max0 = maxpass
135	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	igadkey = (act0 + 1)
143	& + act1*max0
144	& + act2max0max1
145	& + act3max0max1*max2
146	& + act4max0max1max2max3
147	if (tracerIdentity.GT.maxpass) then
148	print *, 'ph-pass gad_advection ', maxpass, tracerIdentity
149	STOP 'maxpass seems smaller than tracerIdentity'
150	endif
151	#endif /* ALLOW_AUTODIFF_TAMC */
152
153	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	rTransKp1(i,j)= 0. _d 0
168	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	#ifdef ALLOW_AUTODIFF_TAMC
179	kkey = (igadkey-1)*Nr + k
180	CADJ STORE tracer(:,:,k,bi,bj) =
181	CADJ & comlev1_bibj_k_gad, key=kkey, byte=isbyte
182	#endif /* ALLOW_AUTODIFF_TAMC */
183
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	#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	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	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	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	IF (useCubedSphereExchange) THEN
226
227	nipass=3
228	#ifdef ALLOW_AUTODIFF_TAMC
229	if ( nipass.GT.maxcube )
230	& STOP 'maxcube needs to be = 3'
231	#endif
232	ELSE
233	nipass=1
234	ENDIF
235	cph nipass=1
236
237	C-- Multiple passes for different directions on different tiles
238	C-- For cube need one pass for each of red, green and blue axes.
239	DO ipass=1,nipass
240	#ifdef ALLOW_AUTODIFF_TAMC
241	passkey = ipass + (k-1) *maxcube
242	& + (igadkey-1)maxcubeNr
243	IF (nipass .GT. maxpass) THEN
244	STOP 'GAD_ADVECTION: nipass > maxcube. check tamc.h'
245	ENDIF
246	#endif /* ALLOW_AUTODIFF_TAMC */
247
248	IF (nipass.EQ.3) THEN
249	calc_fluxes_X=.FALSE.
250	calc_fluxes_Y=.FALSE.
251	IF (ipass.EQ.1 .AND. (nCFace.EQ.1 .OR. nCFace.EQ.2) ) THEN
252	calc_fluxes_X=.TRUE.
253	ELSEIF (ipass.EQ.1 .AND. (nCFace.EQ.4 .OR. nCFace.EQ.5) ) THEN
254	calc_fluxes_Y=.TRUE.
255	ELSEIF (ipass.EQ.2 .AND. (nCFace.EQ.1 .OR. nCFace.EQ.6) ) THEN
256	calc_fluxes_Y=.TRUE.
257	ELSEIF (ipass.EQ.2 .AND. (nCFace.EQ.3 .OR. nCFace.EQ.4) ) THEN
258	calc_fluxes_X=.TRUE.
259	ELSEIF (ipass.EQ.3 .AND. (nCFace.EQ.2 .OR. nCFace.EQ.3) ) THEN
260	calc_fluxes_Y=.TRUE.
261	ELSEIF (ipass.EQ.3 .AND. (nCFace.EQ.5 .OR. nCFace.EQ.6) ) THEN
262	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	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	DO j=1,OLy
313	DO i=1,OLx
314	localTij( 1-i , 1-j )=localTij( 1-j , i )
315	ENDDO
316	ENDDO
317	ENDIF
318	IF ( southEastCorner ) THEN
319	DO J=1,OLy
320	DO I=1,OLx
321	localTij(sNx+I, 1-J )=localTij(sNx+J, I )
322	ENDDO
323	ENDDO
324	ENDIF
325	IF ( northWestCorner ) THEN
326	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	ENDIF
332	IF ( northEastCorner ) THEN
333	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	ENDIF
339	ENDIF
340
341	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
348	#ifdef ALLOW_AUTODIFF_TAMC
349	#ifndef DISABLE_MULTIDIM_ADVECTION
350	CADJ STORE localTij(:,:) =
351	CADJ & comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
352	#endif
353	#endif /* ALLOW_AUTODIFF_TAMC */
354
355	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	STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'
366	ENDIF
367
368	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	C-- End of X direction
391	ENDIF
392
393	C-- Y direction
394	IF (calc_fluxes_Y) THEN
395
396	IF (useCubedSphereExchange) THEN
397	C-- Internal exchange for calculations in Y
398	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	DO J=1,Oly
436	DO I=1,Olx
437	localTij( 1-i , 1-j ) = localTij(j , 1-i )
438	ENDDO
439	ENDDO
440	ENDIF
441	IF ( southEastCorner ) THEN
442	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	ENDIF
448	IF ( northWestCorner ) THEN
449	DO J=1,Oly
450	DO I=1,Olx
451	localTij( 1-i ,sNy+j) = localTij(j ,sNy+i)
452	ENDDO
453	ENDDO
454	ENDIF
455	IF ( northEastCorner ) THEN
456	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	ENDIF
462	ENDIF
463
464	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
471	#ifdef ALLOW_AUTODIFF_TAMC
472	#ifndef DISABLE_MULTIDIM_ADVECTION
473	CADJ STORE localTij(:,:) =
474	CADJ & comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
475	#endif
476	#endif /* ALLOW_AUTODIFF_TAMC */
477
478	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
491	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
502	#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
513	C-- End of Y direction
514	ENDIF
515
516	C-- End of ipass loop
517	ENDDO
518
519	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
536	C-- End of K loop for horizontal fluxes
537	ENDDO
538
539	IF ( .NOT.implicitAdvection ) THEN
540	C-- Start of k loop for vertical flux
541	DO k=Nr,1,-1
542	#ifdef ALLOW_AUTODIFF_TAMC
543	kkey = (igadkey-1)*Nr + k
544	#endif /* ALLOW_AUTODIFF_TAMC */
545	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	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
553	C-- Compute Vertical transport
554	#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
563	C- Surface interface :
564	DO j=1-Oly,sNy+Oly
565	DO i=1-Olx,sNx+Olx
566	rTransKp1(i,j) = kp1Msk*rTrans(i,j)
567	rTrans(i,j) = 0.
568	fVerT(i,j,kUp) = 0.
569	af(i,j) = 0.
570	ENDDO
571	ENDDO
572
573	ELSE
574	C- Interior interface :
575
576	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
585	#ifdef ALLOW_GMREDI
586	C-- Residual transp = Bolus transp + Eulerian transp
587	IF (useGMRedi)
588	& CALL GMREDI_CALC_WFLOW(
589	& rTrans, bi, bj, k, myThid)
590	#endif /* ALLOW_GMREDI */
591
592	#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	C- Compute vertical advective flux in the interior:
600	IF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
601	CALL GAD_FLUXLIMIT_ADV_R(
602	& bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
603	ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
604	CALL GAD_DST3_ADV_R(
605	& bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
606	ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
607	CALL GAD_DST3FL_ADV_R(
608	& bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
609	ELSE
610	STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
611	ENDIF
612	C- add the advective flux to fVerT
613	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
619	C- end Surface/Interior if bloc
620	ENDIF
621
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
629	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
643	C-- End of K loop for vertical flux
644	ENDDO
645	C-- end of if not.implicitAdvection block
646	ENDIF
647
648	RETURN
649	END