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	C $Header: /u/u0/gcmpack/MITgcm/pkg/generic_advdiff/gad_advection.F,v 1.25 2004/06/30 23:45:35 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	#endif
48	#ifdef ALLOW_EXCH2
49	#include "W2_EXCH2_TOPOLOGY.h"
50	#include "W2_EXCH2_PARAMS.h"
51	#endif /* ALLOW_EXCH2 */
52
53	C !INPUT PARAMETERS: ===================================================
54	C implicitAdvection :: implicit vertical advection (later on)
55	C advectionScheme :: advection scheme to use (Horizontal plane)
56	C vertAdvecScheme :: advection scheme to use (vertical direction)
57	C tracerIdentity :: tracer identifier (required only for OBCS)
58	C uVel :: velocity, zonal component
59	C vVel :: velocity, meridional component
60	C wVel :: velocity, vertical component
61	C tracer :: tracer field
62	C bi,bj :: tile indices
63	C myTime :: current time
64	C myIter :: iteration number
65	C myThid :: thread number
66	LOGICAL implicitAdvection
67	INTEGER advectionScheme, vertAdvecScheme
68	INTEGER tracerIdentity
69	_RL uVel (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
70	_RL vVel (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
71	_RL wVel (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
72	_RL tracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
73	INTEGER bi,bj
74	_RL myTime
75	INTEGER myIter
76	INTEGER myThid
77
78	C !OUTPUT PARAMETERS: ==================================================
79	C gTracer :: tendancy array
80	_RL gTracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
81
82	C !LOCAL VARIABLES: ====================================================
83	C maskUp :: 2-D array for mask at W points
84	C iMin,iMax, :: loop range for called routines
85	C jMin,jMax :: loop range for called routines
86	C i,j,k :: loop indices
87	C kup :: index into 2 1/2D array, toggles between 1 and 2
88	C kdown :: index into 2 1/2D array, toggles between 2 and 1
89	C kp1 :: =k+1 for k<Nr, =Nr for k=Nr
90	C xA,yA :: areas of X and Y face of tracer cells
91	C uTrans,vTrans :: 2-D arrays of volume transports at U,V points
92	C rTrans :: 2-D arrays of volume transports at W points
93	C rTransKp1 :: vertical volume transport at interface k+1
94	C af :: 2-D array for horizontal advective flux
95	C fVerT :: 2 1/2D arrays for vertical advective flux
96	C localTij :: 2-D array, temporary local copy of tracer fld
97	C localTijk :: 3-D array, temporary local copy of tracer fld
98	C kp1Msk :: flag (0,1) for over-riding mask for W levels
99	C calc_fluxes_X :: logical to indicate to calculate fluxes in X dir
100	C calc_fluxes_Y :: logical to indicate to calculate fluxes in Y dir
101	C nipass :: number of passes in multi-dimensional method
102	C ipass :: number of the current pass being made
103	C myTile :: variables used to determine which cube face
104	C nCFace :: owns a tile for cube grid runs using
105	C :: multi-dim advection.
106	_RS maskUp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
107	INTEGER iMin,iMax,jMin,jMax
108	INTEGER i,j,k,kup,kDown
109	_RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
110	_RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
111	_RL uTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
112	_RL vTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
113	_RL rTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
114	_RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
115	_RL af (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116	_RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
117	_RL localTij(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
118	_RL localTijk(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
119	_RL kp1Msk
120	LOGICAL calc_fluxes_X,calc_fluxes_Y
121	INTEGER nipass,ipass
122	INTEGER myTile, nCFace
123	LOGICAL southWestCorner
124	LOGICAL southEastCorner
125	LOGICAL northWestCorner
126	LOGICAL northEastCorner
127	CEOP
128
129	#ifdef ALLOW_AUTODIFF_TAMC
130	act0 = tracerIdentity - 1
131	max0 = maxpass
132	act1 = bi - myBxLo(myThid)
133	max1 = myBxHi(myThid) - myBxLo(myThid) + 1
134	act2 = bj - myByLo(myThid)
135	max2 = myByHi(myThid) - myByLo(myThid) + 1
136	act3 = myThid - 1
137	max3 = nTx*nTy
138	act4 = ikey_dynamics - 1
139	igadkey = (act0 + 1)
140	& + act1*max0
141	& + act2max0max1
142	& + act3max0max1*max2
143	& + act4max0max1max2max3
144	if (tracerIdentity.GT.maxpass) then
145	print *, 'ph-pass gad_advection ', maxpass, tracerIdentity
146	STOP 'maxpass seems smaller than tracerIdentity'
147	endif
148	#endif /* ALLOW_AUTODIFF_TAMC */
149
150	C-- Set up work arrays with valid (i.e. not NaN) values
151	C These inital values do not alter the numerical results. They
152	C just ensure that all memory references are to valid floating
153	C point numbers. This prevents spurious hardware signals due to
154	C uninitialised but inert locations.
155	DO j=1-OLy,sNy+OLy
156	DO i=1-OLx,sNx+OLx
157	xA(i,j) = 0. _d 0
158	yA(i,j) = 0. _d 0
159	uTrans(i,j) = 0. _d 0
160	vTrans(i,j) = 0. _d 0
161	rTrans(i,j) = 0. _d 0
162	fVerT(i,j,1) = 0. _d 0
163	fVerT(i,j,2) = 0. _d 0
164	rTransKp1(i,j)= 0. _d 0
165	ENDDO
166	ENDDO
167
168	iMin = 1-OLx
169	iMax = sNx+OLx
170	jMin = 1-OLy
171	jMax = sNy+OLy
172
173	C-- Start of k loop for horizontal fluxes
174	DO k=1,Nr
175	#ifdef ALLOW_AUTODIFF_TAMC
176	kkey = (igadkey-1)*Nr + k
177	CADJ STORE tracer(:,:,k,bi,bj) =
178	CADJ & comlev1_bibj_k_gad, key=kkey, byte=isbyte
179	#endif /* ALLOW_AUTODIFF_TAMC */
180
181	C-- Get temporary terms used by tendency routines
182	CALL CALC_COMMON_FACTORS (
183	I bi,bj,iMin,iMax,jMin,jMax,k,
184	O xA,yA,uTrans,vTrans,rTrans,maskUp,
185	I myThid)
186
187	#ifdef ALLOW_GMREDI
188	C-- Residual transp = Bolus transp + Eulerian transp
189	IF (useGMRedi)
190	& CALL GMREDI_CALC_UVFLOW(
191	& uTrans, vTrans, bi, bj, k, myThid)
192	#endif /* ALLOW_GMREDI */
193
194	C-- Make local copy of tracer array
195	DO j=1-OLy,sNy+OLy
196	DO i=1-OLx,sNx+OLx
197	localTij(i,j)=tracer(i,j,k,bi,bj)
198	ENDDO
199	ENDDO
200
201	cph The following block is needed for useCubedSphereExchange only,
202	cph but needs to be set for all cases to avoid spurious
203	cph TAF dependencies
204	southWestCorner = .TRUE.
205	southEastCorner = .TRUE.
206	northWestCorner = .TRUE.
207	northEastCorner = .TRUE.
208	#ifdef ALLOW_EXCH2
209	myTile = W2_myTileList(bi)
210	nCFace = exch2_myFace(myTile)
211	southWestCorner = exch2_isWedge(myTile).EQ.1
212	& .AND. exch2_isSedge(myTile).EQ.1
213	southEastCorner = exch2_isEedge(myTile).EQ.1
214	& .AND. exch2_isSedge(myTile).EQ.1
215	northEastCorner = exch2_isEedge(myTile).EQ.1
216	& .AND. exch2_isNedge(myTile).EQ.1
217	northWestCorner = exch2_isWedge(myTile).EQ.1
218	& .AND. exch2_isNedge(myTile).EQ.1
219	#else
220	nCFace = bi
221	#endif
222	IF (useCubedSphereExchange) THEN
223
224	nipass=3
225	#ifdef ALLOW_AUTODIFF_TAMC
226	if ( nipass.GT.maxcube )
227	& STOP 'maxcube needs to be = 3'
228	#endif
229	ELSE
230	nipass=1
231	ENDIF
232	cph nipass=1
233
234	C-- Multiple passes for different directions on different tiles
235	C-- For cube need one pass for each of red, green and blue axes.
236	DO ipass=1,nipass
237	#ifdef ALLOW_AUTODIFF_TAMC
238	passkey = ipass + (k-1) *maxcube
239	& + (igadkey-1)maxcubeNr
240	IF (nipass .GT. maxpass) THEN
241	STOP 'GAD_ADVECTION: nipass > maxcube. check tamc.h'
242	ENDIF
243	#endif /* ALLOW_AUTODIFF_TAMC */
244
245	IF (nipass.EQ.3) THEN
246	calc_fluxes_X=.FALSE.
247	calc_fluxes_Y=.FALSE.
248	IF (ipass.EQ.1 .AND. (nCFace.EQ.1 .OR. nCFace.EQ.2) ) THEN
249	calc_fluxes_X=.TRUE.
250	ELSEIF (ipass.EQ.1 .AND. (nCFace.EQ.4 .OR. nCFace.EQ.5) ) THEN
251	calc_fluxes_Y=.TRUE.
252	ELSEIF (ipass.EQ.2 .AND. (nCFace.EQ.1 .OR. nCFace.EQ.6) ) THEN
253	calc_fluxes_Y=.TRUE.
254	ELSEIF (ipass.EQ.2 .AND. (nCFace.EQ.3 .OR. nCFace.EQ.4) ) THEN
255	calc_fluxes_X=.TRUE.
256	ELSEIF (ipass.EQ.3 .AND. (nCFace.EQ.2 .OR. nCFace.EQ.3) ) THEN
257	calc_fluxes_Y=.TRUE.
258	ELSEIF (ipass.EQ.3 .AND. (nCFace.EQ.5 .OR. nCFace.EQ.6) ) THEN
259	calc_fluxes_X=.TRUE.
260	ENDIF
261	ELSE
262	calc_fluxes_X=.TRUE.
263	calc_fluxes_Y=.TRUE.
264	ENDIF
265
266	C-- X direction
267	IF (calc_fluxes_X) THEN
268
269	C-- Internal exchange for calculations in X
270	IF (useCubedSphereExchange) THEN
271	C-- For cube face corners we need to duplicate the
272	C-- i-1 and i+1 values into the null space as follows:
273	C
274	C
275	C o NW corner: copy T( 0,sNy ) into T( 0,sNy+1) e.g.
276	C \|
277	C x T(0,sNy+1) \|
278	C /\ \|
279	C --\|\|------------\|-----------
280	C \|\| \|
281	C x T(0,sNy) \| x T(1,sNy)
282	C \|
283	C
284	C o SW corner: copy T(0,1) into T(0,0) e.g.
285	C \|
286	C x T(0,1) \| x T(1,1)
287	C \|\| \|
288	C --\|\|------------\|-----------
289	C \/ \|
290	C x T(0,0) \|
291	C \|
292	C
293	C o NE corner: copy T(sNx+1,sNy ) into T(sNx+1,sNy+1) e.g.
294	C \|
295	C \| x T(sNx+1,sNy+1)
296	C \| /\
297	C ----------------\|--\|\|-------
298	C \| \|\|
299	C x T(sNx,sNy) \| x T(sNx+1,sNy )
300	C \|
301	C o SE corner: copy T(sNx+1,1 ) into T(sNx+1,0 ) e.g.
302	C \|
303	C x T(sNx,1) \| x T(sNx+1, 1)
304	C \| \|\|
305	C ----------------\|--\|\|-------
306	C \| \/
307	C \| x T(sNx+1, 0)
308	IF ( southWestCorner ) THEN
309	DO J=1,OLy
310	DO I=1,OLx
311	localTij(1-I, 1-J )= localTij(1-J ,1 )
312	ENDDO
313	ENDDO
314	ENDIF
315	IF ( southEastCorner ) THEN
316	DO J=1,OLy
317	DO I=1,OLx
318	localTij(sNx+I, 1-J )=localTij(sNx+J, I )
319	ENDDO
320	ENDDO
321	ENDIF
322	IF ( northWestCorner ) THEN
323	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	ENDIF
329	IF ( northEastCorner ) THEN
330	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	ENDIF
336	ENDIF
337
338	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
345	#ifdef ALLOW_AUTODIFF_TAMC
346	#ifndef DISABLE_MULTIDIM_ADVECTION
347	CADJ STORE localTij(:,:) =
348	CADJ & comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
349	#endif
350	#endif /* ALLOW_AUTODIFF_TAMC */
351
352	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	STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'
363	ENDIF
364
365	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	C-- End of X direction
388	ENDIF
389
390	C-- Y direction
391	IF (calc_fluxes_Y) THEN
392
393	IF (useCubedSphereExchange) THEN
394	C-- Internal exchange for calculations in Y
395	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	DO J=1,Oly
433	DO I=1,Olx
434	localTij( 1-i , 1-j ) = localTij(j , 1-i )
435	ENDDO
436	ENDDO
437	ENDIF
438	IF ( southEastCorner ) THEN
439	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	ENDIF
445	IF ( northWestCorner ) THEN
446	DO J=1,Oly
447	DO I=1,Olx
448	localTij( 1-i ,sNy+j) = localTij(j ,sNy+i)
449	ENDDO
450	ENDDO
451	ENDIF
452	IF ( northEastCorner ) THEN
453	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	ENDIF
459	ENDIF
460
461	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
468	#ifdef ALLOW_AUTODIFF_TAMC
469	#ifndef DISABLE_MULTIDIM_ADVECTION
470	CADJ STORE localTij(:,:) =
471	CADJ & comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
472	#endif
473	#endif /* ALLOW_AUTODIFF_TAMC */
474
475	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
488	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
499	#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
510	C-- End of Y direction
511	ENDIF
512
513	C-- End of ipass loop
514	ENDDO
515
516	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
533	C-- End of K loop for horizontal fluxes
534	ENDDO
535
536	IF ( .NOT.implicitAdvection ) THEN
537	C-- Start of k loop for vertical flux
538	DO k=Nr,1,-1
539	#ifdef ALLOW_AUTODIFF_TAMC
540	kkey = (igadkey-1)*Nr + k
541	#endif /* ALLOW_AUTODIFF_TAMC */
542	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	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
550	C-- Compute Vertical transport
551	#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
560	C- Surface interface :
561	DO j=1-Oly,sNy+Oly
562	DO i=1-Olx,sNx+Olx
563	rTransKp1(i,j) = kp1Msk*rTrans(i,j)
564	rTrans(i,j) = 0.
565	fVerT(i,j,kUp) = 0.
566	af(i,j) = 0.
567	ENDDO
568	ENDDO
569
570	ELSE
571	C- Interior interface :
572
573	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
582	#ifdef ALLOW_GMREDI
583	C-- Residual transp = Bolus transp + Eulerian transp
584	IF (useGMRedi)
585	& CALL GMREDI_CALC_WFLOW(
586	& rTrans, bi, bj, k, myThid)
587	#endif /* ALLOW_GMREDI */
588
589	#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	C- Compute vertical advective flux in the interior:
597	IF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
598	CALL GAD_FLUXLIMIT_ADV_R(
599	& bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
600	ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
601	CALL GAD_DST3_ADV_R(
602	& bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
603	ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
604	CALL GAD_DST3FL_ADV_R(
605	& bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
606	ELSE
607	STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
608	ENDIF
609	C- add the advective flux to fVerT
610	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
616	C- end Surface/Interior if bloc
617	ENDIF
618
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
626	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
640	C-- End of K loop for vertical flux
641	ENDDO
642	C-- end of if not.implicitAdvection block
643	ENDIF
644
645	RETURN
646	END