pkg/generic_advdiff/gad_advection.F

C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_advection.F,v 1.28 2004/09/21 12:13:44 jmc 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  maskLocW      :: 2-D array for mask at West points
C  maskLocS      :: 2-D array for mask at South 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  afx           :: 2-D array for horizontal advective flux, x direction
C  afy           :: 2-D array for horizontal advective flux, y direction
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)
      _RS maskLocW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RS maskLocS(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 afx     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL afy     (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, withSigns
      INTEGER nipass,ipass
      INTEGER myTile, nCFace
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 and mask West & South
      DO j=1-OLy,sNy+OLy
       DO i=1-OLx,sNx+OLx
        localTij(i,j)=tracer(i,j,k,bi,bj)
        maskLocW(i,j)=maskW(i,j,k,bi,bj)
        maskLocS(i,j)=maskS(i,j,k,bi,bj)
       ENDDO
      ENDDO

      IF (useCubedSphereExchange) THEN
        withSigns = .FALSE.
        CALL FILL_CS_CORNER_UV_RS( 
     &            withSigns, maskLocW,maskLocS, bi,bj, myThid )
      ENDIF
#ifdef ALLOW_EXCH2
       myTile = W2_myTileList(bi)
       nCFace = exch2_myFace(myTile)
#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

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---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C--   X direction
      IF (calc_fluxes_X) THEN

C--   Internal exchange for calculations in X
      IF (useCubedSphereExchange) THEN
         CALL FILL_CS_CORNER_TR_RL( .TRUE., localTij, bi,bj, myThid )
      ENDIF

C-    Advective flux in X
      DO j=1-Oly,sNy+Oly
       DO i=1-Olx,sNx+Olx
        afx(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,
     I                            uTrans, uVel, maskLocW, localTij,
     O                            afx, myThid )
      ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
        CALL GAD_DST3_ADV_X(      bi,bj,k, deltaTtracer,
     I                            uTrans, uVel, maskLocW, localTij,
     O                            afx, myThid )
      ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
        CALL GAD_DST3FL_ADV_X(    bi,bj,k, deltaTtracer,
     I                            uTrans, uVel, maskLocW, localTij,
     O                            afx, 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)
     &    *( afx(i+1,j)-afx(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---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C--   Y direction
      IF (calc_fluxes_Y) THEN

C--   Internal exchange for calculations in Y
      IF (useCubedSphereExchange) THEN
         CALL FILL_CS_CORNER_TR_RL(.FALSE., localTij, bi,bj, myThid )
      ENDIF

C-    Advective flux in Y
      DO j=1-Oly,sNy+Oly
       DO i=1-Olx,sNx+Olx
        afy(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,
     I                            vTrans, vVel, maskLocS, localTij,
     O                            afy, myThid )
      ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
        CALL GAD_DST3_ADV_Y(      bi,bj,k, deltaTtracer,
     I                            vTrans, vVel, maskLocS, localTij,
     O                            afy, myThid )
      ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
        CALL GAD_DST3FL_ADV_Y(    bi,bj,k, deltaTtracer,
     I                            vTrans, vVel, maskLocS, localTij,
     O                            afy, 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)
     &    *( afy(i,j+1)-afy(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

#ifdef ALLOW_DEBUG
      IF ( debugLevel .GE. debLevB
     &   .AND. k.EQ.3 .AND. myIter.EQ.1+nIter0
     &   .AND. nPx.EQ.1 .AND. nPy.EQ.1
     &   .AND. useCubedSphereExchange ) THEN
        CALL DEBUG_CS_CORNER_UV( ' afx,afy from GAD_ADVECTION',
     &             afx,afy, k, standardMessageUnit,bi,bj,myThid )
      ENDIF
#endif /* ALLOW_DEBUG */

C--   End of K loop for horizontal fluxes
      ENDDO

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

      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.
          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)
           fVerT(i,j,kUp) = 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
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
           CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, deltaTtracer,
     I                               rTrans, wVel, localTijk,
     O                               fVerT(1-Olx,1-Oly,kUp), myThid )
         ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
           CALL GAD_DST3_ADV_R(      bi,bj,k, deltaTtracer,
     I                               rTrans, wVel, localTijk,
     O                               fVerT(1-Olx,1-Oly,kUp), myThid )
         ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
           CALL GAD_DST3FL_ADV_R(    bi,bj,k, deltaTtracer,
     I                               rTrans, wVel, localTijk,
     O                               fVerT(1-Olx,1-Oly,kUp), myThid )
         ELSE
          STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
         ENDIF

C- end Surface/Interior if bloc
        ENDIF

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

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

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