pkg/generic_advdiff/gad_advection.F

C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_advection.F,v 1.16 2003/11/25 23:31:44 heimbach Exp $
C $Name:  $

CBOI
C !TITLE: pkg/generic\_advdiff
C !AUTHORS: adcroft@mit.edu
C !INTRODUCTION: Generic Advection Diffusion Package
C
C Package "generic\_advdiff" provides a common set of routines for calculating
C advective/diffusive fluxes for tracers (cell centered quantities on a C-grid).
C
C Many different advection schemes are available: the standard centered
C second order, centered fourth order and upwind biased third order schemes
C are known as linear methods and require some stable time-stepping method
C such as Adams-Bashforth. Alternatives such as flux-limited schemes are
C stable in the forward sense and are best combined with the multi-dimensional
C method provided in gad\_advection.
C
C There are two high-level routines:
C  \begin{itemize}
C  \item{GAD\_CALC\_RHS} calculates all fluxes at time level "n" and is used
C  for the standard linear schemes. This must be used in conjuction with
C  Adams-Bashforth time-stepping. Diffusive and parameterized fluxes are
C  always calculated here.
C  \item{GAD\_ADVECTION} calculates just the advective fluxes using the
C  non-linear schemes and can not be used in conjuction with Adams-Bashforth
C  time-stepping.
C  \end{itemize}
CEOI

#include "GAD_OPTIONS.h"

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---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

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    :: vertical advection treated implicitly (later on)
C  advectionScheme      :: advection scheme to use
C  tracerIdentity       :: identifier for the tracer (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,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,rTrans :: 2-D arrays of volume transports at U,V and 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 used as temporary local copy of tracer fld
C  localTijk            :: 3-D array used as temporary local copy of tracer fld
C  kp1Msk               :: flag (0,1) to act as 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 to make 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

      DO j=1-Oly,sNy+Oly
       DO i=1-Olx,sNx+Olx
        localTijk(i,j,k)=localTij(i,j)
       ENDDO
      ENDDO

C--   End of ipass loop
      ENDDO

C--   End of K loop for horizontal fluxes
      ENDDO

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
C     Note: wVel needs to be masked 

      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 */

       IF ( .NOT.implicitAdvection ) THEN
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
       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 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

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