pkg/generic_advdiff/gad_advection.F

C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_advection.F,v 1.24 2004/06/28 21:10:55 dimitri 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
        localTij(0    ,0    )= localTij(0    ,1  )
       ENDIF
       IF ( southEastCorner ) THEN
        localTij(sNx+1,0    )= localTij(sNx+1,1  )
       ENDIF
       IF ( northWestCorner ) THEN
        localTij(0    ,sNy+1)= localTij(0    ,sNy)
       ENDIF
       IF ( northEastCorner ) THEN
        localTij(sNx+1,sNy+1)= localTij(sNx+1,sNy)
       ENDIF
      ENDIF

C-    Advective flux in X
      DO j=1-Oly,sNy+Oly
       DO i=1-Olx,sNx+Olx
        af(i,j) = 0.
       ENDDO
      ENDDO

#ifdef ALLOW_AUTODIFF_TAMC
#ifndef DISABLE_MULTIDIM_ADVECTION
CADJ STORE localTij(:,:)  = 
CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
#endif
#endif /* ALLOW_AUTODIFF_TAMC */

      IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
       CALL GAD_FLUXLIMIT_ADV_X(
     &      bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)
      ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
       CALL GAD_DST3_ADV_X(
     &       bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)
      ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
       CALL GAD_DST3FL_ADV_X(
     &       bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)
      ELSE
       STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'
      ENDIF

      DO j=1-Oly,sNy+Oly
       DO i=1-Olx,sNx+Olx-1
        localTij(i,j)=localTij(i,j)-deltaTtracer*
     &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
     &    *recip_rA(i,j,bi,bj)
     &    *( af(i+1,j)-af(i,j)
     &      -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
     &     )
       ENDDO
      ENDDO

#ifdef ALLOW_OBCS
C--   Apply open boundary conditions
      IF (useOBCS) THEN
       IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
        CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )
       ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
        CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
       END IF
      END IF
#endif /* ALLOW_OBCS */

C--   End of X direction
      ENDIF

C--   Y direction
      IF (calc_fluxes_Y) THEN

      IF (useCubedSphereExchange) THEN
C--   Internal exchange for calculations in Y
C--    For cube face corners we need to duplicate the
C--    j-1 and j+1 values into the null space as follows:
C
C      o SW corner: copy T(0,1) into T(0,0) e.g.
C                      |
C                      |  x T(1,1)
C                      |
C      ----------------|-----------
C                      |
C         x T(0,0)<====== x T(1,0)
C                      |
C
C      o NW corner: copy T(    0,sNy  ) into T(    0,sNy+1) e.g.
C                      |
C         x T(0,sNy+1)<=== x T(1,sNy+1)
C                      |
C      ----------------|-----------
C                      |
C                      |   x T(1,sNy)
C                      |
C
C      o NE corner: copy T(sNx+1,sNy  ) into T(sNx+1,sNy+1) e.g.
C                      |
C      x T(sNx,sNy+1)=====>x T(sNx+1,sNy+1)
C                      |     
C      ----------------|-----------
C                      |     
C      x T(sNx,sNy)    |                       
C                      |
C      o SE corner: copy T(sNx+1,1    ) into T(sNx+1,0    ) e.g.
C                      |
C         x T(sNx,1)   |                    
C                      |    
C      ----------------|-----------
C                      |     
C         x T(sNx,0) =====>x T(sNx+1,    0)
       IF ( southWestCorner ) THEN
         localTij(    0,0    ) = localTij(  1,0    )
       ENDIF
       IF ( southEastCorner ) THEN
         localTij(sNx+1,0    ) = localTij(sNx,0    )
       ENDIF
       IF ( northWestCorner ) THEN
         localTij(0    ,sNy+1) = localTij(  1,sNy+1)
       ENDIF
       IF ( northEastCorner ) THEN
         localTij(sNx+1,sNy+1) = localTij(sNx,sNy+1)
       ENDIF
      ENDIF

C-    Advective flux in Y
      DO j=1-Oly,sNy+Oly
       DO i=1-Olx,sNx+Olx
        af(i,j) = 0.
       ENDDO
      ENDDO

#ifdef ALLOW_AUTODIFF_TAMC 
#ifndef DISABLE_MULTIDIM_ADVECTION
CADJ STORE localTij(:,:)  = 
CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
#endif
#endif /* ALLOW_AUTODIFF_TAMC */

      IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
       CALL GAD_FLUXLIMIT_ADV_Y(
     &       bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid)
      ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
       CALL GAD_DST3_ADV_Y(
     &       bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid)
      ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
       CALL GAD_DST3FL_ADV_Y(
     &       bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid)
      ELSE
       STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
      ENDIF

      DO j=1-Oly,sNy+Oly-1
       DO i=1-Olx,sNx+Olx
        localTij(i,j)=localTij(i,j)-deltaTtracer*
     &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
     &    *recip_rA(i,j,bi,bj)
     &    *( af(i,j+1)-af(i,j)
     &      -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
     &     )
       ENDDO
      ENDDO

#ifdef ALLOW_OBCS
C--   Apply open boundary conditions
      IF (useOBCS) THEN
       IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
        CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )
       ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
        CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
       END IF
      END IF
#endif /* ALLOW_OBCS */

C--   End of Y direction
      ENDIF

C--   End of ipass loop
      ENDDO

      IF ( implicitAdvection ) THEN
C-    explicit advection is done ; store tendency in gTracer:
        DO j=1-Oly,sNy+Oly
         DO i=1-Olx,sNx+Olx
          gTracer(i,j,k,bi,bj)=
     &     (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer
         ENDDO
        ENDDO
      ELSE
C-    horizontal advection done; store intermediate result in 3D array:
       DO j=1-Oly,sNy+Oly
        DO i=1-Olx,sNx+Olx
         localTijk(i,j,k)=localTij(i,j)
        ENDDO
       ENDDO
      ENDIF

C--   End of K loop for horizontal fluxes
      ENDDO

      IF ( .NOT.implicitAdvection ) THEN
C--   Start of k loop for vertical flux
       DO k=Nr,1,-1
#ifdef ALLOW_AUTODIFF_TAMC 
         kkey = (igadkey-1)*Nr + k
#endif /* ALLOW_AUTODIFF_TAMC */
C--   kup    Cycles through 1,2 to point to w-layer above
C--   kDown  Cycles through 2,1 to point to w-layer below
        kup  = 1+MOD(k+1,2)
        kDown= 1+MOD(k,2)
c       kp1=min(Nr,k+1)
        kp1Msk=1.
        if (k.EQ.Nr) kp1Msk=0.

C-- Compute Vertical transport
#ifdef ALLOW_AIM
C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
        IF ( k.EQ.1 .OR.
     &     (useAIM .AND. tracerIdentity.EQ.GAD_SALINITY .AND. k.EQ.Nr)
     &              ) THEN
#else
        IF ( k.EQ.1 ) THEN
#endif

C- Surface interface :
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
           rTransKp1(i,j) = kp1Msk*rTrans(i,j)
           rTrans(i,j) = 0.
           fVerT(i,j,kUp) = 0.
           af(i,j) = 0.
          ENDDO
         ENDDO

        ELSE
C- Interior interface :

         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
           rTransKp1(i,j) = kp1Msk*rTrans(i,j)
           rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)
     &                 *maskC(i,j,k-1,bi,bj)
           af(i,j) = 0.
          ENDDO
         ENDDO

#ifdef ALLOW_GMREDI
C--   Residual transp = Bolus transp + Eulerian transp
         IF (useGMRedi) 
     &   CALL GMREDI_CALC_WFLOW(
     &                    rTrans, bi, bj, k, myThid)
#endif /* ALLOW_GMREDI */

#ifdef ALLOW_AUTODIFF_TAMC 
CADJ STORE localTijk(:,:,k)  
CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
CADJ STORE rTrans(:,:)  
CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */

C-    Compute vertical advective flux in the interior:
         IF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
          CALL GAD_FLUXLIMIT_ADV_R(
     &        bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
         ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
          CALL GAD_DST3_ADV_R(
     &        bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
         ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
          CALL GAD_DST3FL_ADV_R(
     &        bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
         ELSE
          STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
         ENDIF
C-    add the advective flux to fVerT
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
           fVerT(i,j,kUp) = af(i,j)
          ENDDO
         ENDDO

C- end Surface/Interior if bloc
        ENDIF

#ifdef ALLOW_AUTODIFF_TAMC 
CADJ STORE rTrans(:,:)  
CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
CADJ STORE rTranskp1(:,:)  
CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */

C--   Divergence of vertical fluxes
        DO j=1-Oly,sNy+Oly
         DO i=1-Olx,sNx+Olx
          localTij(i,j)=localTijk(i,j,k)-deltaTtracer*
     &     _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
     &     *recip_rA(i,j,bi,bj)
     &     *( fVerT(i,j,kUp)-fVerT(i,j,kDown)
     &       -tracer(i,j,k,bi,bj)*(rTrans(i,j)-rTransKp1(i,j))
     &      )*rkFac
          gTracer(i,j,k,bi,bj)=
     &     (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer
         ENDDO
        ENDDO
 
C--   End of K loop for vertical flux
       ENDDO
C--   end of if not.implicitAdvection block
      ENDIF 

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