pkg/generic_advdiff/gad_advection.F

C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_advection.F,v 1.20 2004/03/28 19:28:34 edhill 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, 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

C !INPUT PARAMETERS: ===================================================
C  implicitAdvection :: implicit vertical advection (later on)
C  advectionScheme   :: advection scheme to use
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
      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
      _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
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

      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
      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. (bi.EQ.1 .OR. bi.EQ.2) ) THEN
        calc_fluxes_X=.TRUE.
       ELSEIF (ipass.EQ.1 .AND. (bi.EQ.4 .OR. bi.EQ.5) ) THEN
        calc_fluxes_Y=.TRUE.
       ELSEIF (ipass.EQ.2 .AND. (bi.EQ.1 .OR. bi.EQ.6) ) THEN
        calc_fluxes_Y=.TRUE.
       ELSEIF (ipass.EQ.2 .AND. (bi.EQ.3 .OR. bi.EQ.4) ) THEN
        calc_fluxes_X=.TRUE.
       ELSEIF (ipass.EQ.3 .AND. (bi.EQ.2 .OR. bi.EQ.3) ) THEN
        calc_fluxes_Y=.TRUE.
       ELSEIF (ipass.EQ.3 .AND. (bi.EQ.5 .OR. bi.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
       DO j=1,Oly
        DO i=1,Olx
         localTij( 1-i , 1-j )=localTij( 1-j ,    i    )
         localTij( 1-i ,sNy+j)=localTij( 1-j , sNy+1-i )
         localTij(sNx+i, 1-j )=localTij(sNx+j,    i    )
         localTij(sNx+i,sNy+j)=localTij(sNx+j, sNy+1-i )
        ENDDO
       ENDDO
      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

C--   Internal exchange for calculations in Y
      IF (useCubedSphereExchange) THEN
       DO j=1,Oly
        DO i=1,Olx
         localTij( 1-i , 1-j )=localTij(   j   , 1-i )
         localTij( 1-i ,sNy+j)=localTij(   j   ,sNy+i)
         localTij(sNx+i, 1-j )=localTij(sNx+1-j, 1-i )
         localTij(sNx+i,sNy+j)=localTij(sNx+1-j,sNy+i)
        ENDDO
       ENDDO
      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
        IF (k.EQ.1) THEN

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