pkg/generic_advdiff/gad_advection.F

C $Header: /u/gcmpack/models/MITgcmUV/pkg/generic_advdiff/gad_advection.F,v 1.5 2001/09/21 13:11:43 adcroft Exp $
C $Name:  $

CBOI
C !TITLE: pkg/generic\_advdiff
C !AUTHORS: adcroft@mit.edu
C !INTRODUCTION:
C \section{Generica 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(bi,bj,advectionScheme,tracerIdentity,
     U                         Tracer,Gtracer,
     I                         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 "DYNVARS.h"
#include "GRID.h"
#include "GAD.h"
#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
# include "tamc_keys.h"
#endif

C !INPUT PARAMETERS: ===================================================
C  bi,bj                :: tile indices
C  advectionScheme      :: advection scheme to use
C  tracerIdentity       :: identifier for the tracer (required only for OBCS)
C  Tracer               :: tracer field
C  myTime               :: current time
C  myIter               :: iteration number
C  myThid               :: thread number
      INTEGER bi,bj
      INTEGER advectionScheme
      INTEGER tracerIdentity
      _RL Gtracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
      _RL myTime
      INTEGER myIter
      INTEGER myThid

C !OUTPUT PARAMETERS: ==================================================
C  gTracer              :: tendancy array
      _RL Tracer(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  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,kp1
      _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 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
          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
          ikey = (act1 + 1) + act2*max1
     &                      + act3*max1*max2
     &                      + act4*max1*max2*max3
#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
       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 = (ikey-1)*Nr + k
CADJ STORE tracer(:,:,k,bi,bj) = comlev1_bibj_k, 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)

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
      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)   *maxpass
     &                   + (ikey-1)*maxpass*Nr
         IF (nipass .GT. maxpass) THEN
          STOP 'GAD_ADVECTION: nipass > maxpass. 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
#ifdef ALLOW_MULTIDIM_ADVECTION
CADJ STORE localTij(:,:)  = comlev1_bibj_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 mutli-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 
#ifdef ALLOW_MULTIDIM_ADVECTION
CADJ STORE localTij(:,:)  = comlev1_bibj_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 = (ikey-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--   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)

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

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

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