pkg/generic_advdiff/gad_calc_rhs.F

C $Header: /u/gcmpack/models/MITgcmUV/pkg/generic_advdiff/gad_calc_rhs.F,v 1.11 2001/09/19 20:45:09 adcroft Exp $
C $Name:  $

#include "GAD_OPTIONS.h"

CBOP
C !ROUTINE: GAD_CALC_RHS

C !INTERFACE: ==========================================================
      SUBROUTINE GAD_CALC_RHS( 
     I           bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
     I           xA,yA,uTrans,vTrans,rTrans,maskUp,
     I           diffKh, diffK4, KappaRT, Tracer,
     I           tracerIdentity, advectionScheme,
     U           fVerT, gTracer,
     I           myThid )

C !DESCRIPTION:
C Calculates the tendancy of a tracer due to advection and diffusion.
C It calculates the fluxes in each direction indepentently and then
C sets the tendancy to the divergence of these fluxes. The advective
C fluxes are only calculated here when using the linear advection schemes
C otherwise only the diffusive and parameterized fluxes are calculated.
C
C Contributions to the flux are calculated and added:
C \begin{equation*}
C {\bf F} = {\bf F}_{adv} + {\bf F}_{diff} +{\bf F}_{GM} + {\bf F}_{KPP}
C \end{equation*}
C
C The tendancy is the divergence of the fluxes:
C \begin{equation*}
C G_\theta = G_\theta + \nabla \cdot {\bf F}
C \end{equation*}
C
C The tendancy is assumed to contain data on entry.

C !USES: ===============================================================
      IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DYNVARS.h"
#include "GAD.h"

C !INPUT PARAMETERS: ===================================================
C  bi,bj                :: tile indices
C  iMin,iMax,jMin,jMax  :: loop range for called routines
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  maskUp               :: 2-D array for mask at W points
C  diffKh               :: horizontal diffusion coefficient
C  diffK4               :: bi-harmonic diffusion coefficient
C  KappaRT              :: 3-D array for vertical diffusion coefficient
C  Tracer               :: tracer field
C  tracerIdentity       :: identifier for the tracer (required only for KPP)
C  advectionScheme      :: advection scheme to use
C  myThid               :: thread number
      INTEGER bi,bj,iMin,iMax,jMin,jMax
      INTEGER k,kUp,kDown,kM1
      _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)
      _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL diffKh, diffK4
      _RL KappaRT(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL Tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
      INTEGER tracerIdentity
      INTEGER advectionScheme
      INTEGER myThid

C !OUTPUT PARAMETERS: ==================================================
C  gTracer              :: tendancy array
C  fVerT                :: 2 1/2D arrays for vertical advective flux
      _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
      _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)

C !LOCAL VARIABLES: ====================================================
C  i,j                  :: loop indices
C  df4                  :: used for storing del^2 T for bi-harmonic term
C  fZon                 :: zonal flux
C  fmer                 :: meridional flux
C  af                   :: advective flux
C  df                   :: diffusive flux
C  localT               :: local copy of tracer field
      INTEGER i,j
      _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
CEOP

#ifdef ALLOW_AUTODIFF_TAMC
C--   only the kUp part of fverT is set in this subroutine
C--   the kDown is still required
      fVerT(1,1,kDown) = fVerT(1,1,kDown)
#endif
      DO j=1-OLy,sNy+OLy
       DO i=1-OLx,sNx+OLx
        fZon(i,j)      = 0. _d 0
        fMer(i,j)      = 0. _d 0
        fVerT(i,j,kUp) = 0. _d 0
       ENDDO
      ENDDO

C--   Make local copy of tracer array
      DO j=1-OLy,sNy+OLy
       DO i=1-OLx,sNx+OLx
        localT(i,j)=tracer(i,j,k,bi,bj)
       ENDDO
      ENDDO

C--   Unless we have already calculated the advection terms we initialize
C     the tendency to zero.
      IF (.NOT. multiDimAdvection .OR.
     &    advectionScheme.EQ.ENUM_CENTERED_2ND .OR.
     &    advectionScheme.EQ.ENUM_UPWIND_3RD .OR.
     &    advectionScheme.EQ.ENUM_CENTERED_4TH ) THEN
       DO j=1-Oly,sNy+Oly
        DO i=1-Olx,sNx+Olx
         gTracer(i,j,k,bi,bj)=0. _d 0
        ENDDO
       ENDDO
      ENDIF

C--   Pre-calculate del^2 T if bi-harmonic coefficient is non-zero
      IF (diffK4 .NE. 0.) THEN
       CALL GAD_GRAD_X(bi,bj,k,xA,localT,fZon,myThid)
       CALL GAD_GRAD_Y(bi,bj,k,yA,localT,fMer,myThid)
       CALL GAD_DEL2(bi,bj,k,fZon,fMer,df4,myThid)
      ENDIF

C--   Initialize net flux in X direction
      DO j=1-Oly,sNy+Oly
       DO i=1-Olx,sNx+Olx
        fZon(i,j) = 0. _d 0
       ENDDO
      ENDDO

C-    Advective flux in X
      IF (.NOT. multiDimAdvection .OR.
     &    advectionScheme.EQ.ENUM_CENTERED_2ND .OR.
     &    advectionScheme.EQ.ENUM_UPWIND_3RD .OR.
     &    advectionScheme.EQ.ENUM_CENTERED_4TH ) THEN
      IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
       CALL GAD_C2_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
      ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
       CALL GAD_FLUXLIMIT_ADV_X(
     &      bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)
      ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
       CALL GAD_U3_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
      ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
       CALL GAD_C4_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
      ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
       CALL GAD_DST3_ADV_X(
     &       bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)
      ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
       CALL GAD_DST3FL_ADV_X(
     &       bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)
      ELSE
       STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'
      ENDIF
      DO j=1-Oly,sNy+Oly
       DO i=1-Olx,sNx+Olx
        fZon(i,j) = fZon(i,j) + af(i,j)
       ENDDO
      ENDDO
      ENDIF

C-    Diffusive flux in X
      IF (diffKh.NE.0.) THEN
       CALL GAD_DIFF_X(bi,bj,k,xA,diffKh,localT,df,myThid)
      ELSE
       DO j=1-Oly,sNy+Oly
        DO i=1-Olx,sNx+Olx
         df(i,j) = 0. _d 0
        ENDDO
       ENDDO
      ENDIF

#ifdef ALLOW_GMREDI
C-    GM/Redi flux in X
      IF (useGMRedi) THEN
C *note* should update GMREDI_XTRANSPORT to use localT and set df  *aja*
        CALL GMREDI_XTRANSPORT(
     I     iMin,iMax,jMin,jMax,bi,bj,K,
     I     xA,Tracer,
     U     df,
     I     myThid)
      ENDIF
#endif
      DO j=1-Oly,sNy+Oly
       DO i=1-Olx,sNx+Olx
        fZon(i,j) = fZon(i,j) + df(i,j)
       ENDDO
      ENDDO

C-    Bi-harmonic duffusive flux in X
      IF (diffK4 .NE. 0.) THEN
       CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)
       DO j=1-Oly,sNy+Oly
        DO i=1-Olx,sNx+Olx
         fZon(i,j) = fZon(i,j) + df(i,j)
        ENDDO
       ENDDO
      ENDIF

C--   Initialize net flux in Y direction
      DO j=1-Oly,sNy+Oly
       DO i=1-Olx,sNx+Olx
        fMer(i,j) = 0. _d 0
       ENDDO
      ENDDO

C-    Advective flux in Y
      IF (.NOT. multiDimAdvection .OR.
     &    advectionScheme.EQ.ENUM_CENTERED_2ND .OR.
     &    advectionScheme.EQ.ENUM_UPWIND_3RD .OR.
     &    advectionScheme.EQ.ENUM_CENTERED_4TH ) THEN
      IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
       CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
      ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
       CALL GAD_FLUXLIMIT_ADV_Y(
     &       bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)
      ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
       CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
      ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
       CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
      ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
       CALL GAD_DST3_ADV_Y(
     &       bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)
      ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
       CALL GAD_DST3FL_ADV_Y(
     &       bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)
      ELSE
       STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'
      ENDIF
      DO j=1-Oly,sNy+Oly
       DO i=1-Olx,sNx+Olx
        fMer(i,j) = fMer(i,j) + af(i,j)
       ENDDO
      ENDDO
      ENDIF

C-    Diffusive flux in Y
      IF (diffKh.NE.0.) THEN
       CALL GAD_DIFF_Y(bi,bj,k,yA,diffKh,localT,df,myThid)
      ELSE
       DO j=1-Oly,sNy+Oly
        DO i=1-Olx,sNx+Olx
         df(i,j) = 0. _d 0
        ENDDO
       ENDDO
      ENDIF

#ifdef ALLOW_GMREDI
C-    GM/Redi flux in Y
      IF (useGMRedi) THEN
C *note* should update GMREDI_YTRANSPORT to use localT and set df  *aja*
       CALL GMREDI_YTRANSPORT(
     I     iMin,iMax,jMin,jMax,bi,bj,K,
     I     yA,Tracer,
     U     df,
     I     myThid)
      ENDIF
#endif
      DO j=1-Oly,sNy+Oly
       DO i=1-Olx,sNx+Olx
        fMer(i,j) = fMer(i,j) + df(i,j)
       ENDDO
      ENDDO

C-    Bi-harmonic flux in Y
      IF (diffK4 .NE. 0.) THEN
       CALL GAD_BIHARM_Y(bi,bj,k,yA,df4,diffK4,df,myThid)
       DO j=1-Oly,sNy+Oly
        DO i=1-Olx,sNx+Olx
         fMer(i,j) = fMer(i,j) + df(i,j)
        ENDDO
       ENDDO
      ENDIF

C-    Advective flux in R
      IF (.NOT. multiDimAdvection .OR.
     &    advectionScheme.EQ.ENUM_CENTERED_2ND .OR.
     &    advectionScheme.EQ.ENUM_UPWIND_3RD .OR.
     &    advectionScheme.EQ.ENUM_CENTERED_4TH ) THEN
C     Note: wVel needs to be masked 
      IF (K.GE.2) THEN
C-    Compute vertical advective flux in the interior:
       IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
        CALL GAD_C2_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
       ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
        CALL GAD_FLUXLIMIT_ADV_R(
     &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)
       ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
        CALL GAD_U3_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
       ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
        CALL GAD_C4_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
       ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
        CALL GAD_DST3_ADV_R(
     &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)
       ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
        CALL GAD_DST3FL_ADV_R(
     &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)
       ELSE
        STOP 'GAD_CALC_RHS: Bad advectionScheme (R)'
       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)*Tracer(i,j,k,bi,bj)
        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)*Tracer(i,j,k,bi,bj)
        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) = fVerT(i,j,kUp) + af(i,j)
       ENDDO
      ENDDO
      ENDIF

C-    Diffusive flux in R
C     Note: For K=1 then KM1=1 and this gives a dT/dr = 0 upper
C           boundary condition.
      IF (implicitDiffusion) THEN
       DO j=1-Oly,sNy+Oly
        DO i=1-Olx,sNx+Olx
         df(i,j) = 0. _d 0
        ENDDO
       ENDDO
      ELSE
       CALL GAD_DIFF_R(bi,bj,k,KappaRT,tracer,df,myThid)
      ENDIF
c     DO j=1-Oly,sNy+Oly
c      DO i=1-Olx,sNx+Olx
c       fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)
c      ENDDO
c     ENDDO

#ifdef ALLOW_GMREDI
C-    GM/Redi flux in R
      IF (useGMRedi) THEN
C *note* should update GMREDI_RTRANSPORT to set df  *aja*
       CALL GMREDI_RTRANSPORT(
     I     iMin,iMax,jMin,jMax,bi,bj,K,
     I     Tracer,
     U     df,
     I     myThid)
c      DO j=1-Oly,sNy+Oly
c       DO i=1-Olx,sNx+Olx
c        fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)
c       ENDDO
c      ENDDO
      ENDIF
#endif

      DO j=1-Oly,sNy+Oly
       DO i=1-Olx,sNx+Olx
        fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)
       ENDDO
      ENDDO

#ifdef ALLOW_KPP
C-    Add non local KPP transport term (ghat) to diffusive T flux.
      IF (useKPP) THEN
       DO j=1-Oly,sNy+Oly
        DO i=1-Olx,sNx+Olx
         df(i,j) = 0. _d 0
        ENDDO
       ENDDO
       IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
C *note* should update KPP_TRANSPORT_T to set df  *aja*
        CALL KPP_TRANSPORT_T(
     I     iMin,iMax,jMin,jMax,bi,bj,k,km1,
     I     KappaRT,
     U     df )
       ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
        CALL KPP_TRANSPORT_S(
     I     iMin,iMax,jMin,jMax,bi,bj,k,km1,
     I     KappaRT,
     U     df )
       ELSE
        STOP 'GAD_CALC_RHS: Ooops'
       ENDIF
       DO j=1-Oly,sNy+Oly
        DO i=1-Olx,sNx+Olx
         fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)
        ENDDO
       ENDDO
      ENDIF
#endif

C--   Divergence of fluxes
      DO j=1-Oly,sNy+Oly-1
       DO i=1-Olx,sNx+Olx-1
        gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)
     &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
     &    *recip_rA(i,j,bi,bj)
     &    *(
     &    +( fZon(i+1,j)-fZon(i,j) )
     &    +( fMer(i,j+1)-fMer(i,j) )
     &    +( fVerT(i,j,kUp)-fVerT(i,j,kDown) )*rkFac
     &    )
       ENDDO
      ENDDO

      RETURN
      END
1	heimbach	1.12	C $Header: /u/gcmpack/models/MITgcmUV/pkg/generic_advdiff/gad_calc_rhs.F,v 1.11 2001/09/19 20:45:09 adcroft Exp $
2	jmc	1.2	C $Name: $
3	adcroft	1.1
4			#include "GAD_OPTIONS.h"
5
6	adcroft	1.11	CBOP
7			C !ROUTINE: GAD_CALC_RHS
8
9			C !INTERFACE: ==========================================================
10	adcroft	1.1	SUBROUTINE GAD_CALC_RHS(
11			I bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
12			I xA,yA,uTrans,vTrans,rTrans,maskUp,
13			I diffKh, diffK4, KappaRT, Tracer,
14	adcroft	1.3	I tracerIdentity, advectionScheme,
15	adcroft	1.1	U fVerT, gTracer,
16			I myThid )
17	adcroft	1.11
18			C !DESCRIPTION:
19			C Calculates the tendancy of a tracer due to advection and diffusion.
20			C It calculates the fluxes in each direction indepentently and then
21			C sets the tendancy to the divergence of these fluxes. The advective
22			C fluxes are only calculated here when using the linear advection schemes
23			C otherwise only the diffusive and parameterized fluxes are calculated.
24			C
25			C Contributions to the flux are calculated and added:
26			C \begin{equation*}
27			C {\bf F} = {\bf F}_{adv} + {\bf F}_{diff} +{\bf F}_{GM} + {\bf F}_{KPP}
28			C \end{equation*}
29			C
30			C The tendancy is the divergence of the fluxes:
31			C \begin{equation*}
32			C G_\theta = G_\theta + \nabla \cdot {\bf F}
33			C \end{equation*}
34			C
35			C The tendancy is assumed to contain data on entry.
36
37			C !USES: ===============================================================
38	adcroft	1.1	IMPLICIT NONE
39			#include "SIZE.h"
40			#include "EEPARAMS.h"
41			#include "PARAMS.h"
42			#include "GRID.h"
43			#include "DYNVARS.h"
44			#include "GAD.h"
45
46	adcroft	1.11	C !INPUT PARAMETERS: ===================================================
47			C bi,bj :: tile indices
48			C iMin,iMax,jMin,jMax :: loop range for called routines
49			C kup :: index into 2 1/2D array, toggles between 1 and 2
50			C kdown :: index into 2 1/2D array, toggles between 2 and 1
51			C kp1 :: =k+1 for k<Nr, =Nr for k=Nr
52			C xA,yA :: areas of X and Y face of tracer cells
53			C uTrans,vTrans,rTrans :: 2-D arrays of volume transports at U,V and W points
54			C maskUp :: 2-D array for mask at W points
55			C diffKh :: horizontal diffusion coefficient
56			C diffK4 :: bi-harmonic diffusion coefficient
57			C KappaRT :: 3-D array for vertical diffusion coefficient
58			C Tracer :: tracer field
59			C tracerIdentity :: identifier for the tracer (required only for KPP)
60			C advectionScheme :: advection scheme to use
61			C myThid :: thread number
62			INTEGER bi,bj,iMin,iMax,jMin,jMax
63	adcroft	1.1	INTEGER k,kUp,kDown,kM1
64			_RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
65			_RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
66			_RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
67			_RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68			_RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
69			_RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70			_RL diffKh, diffK4
71			_RL KappaRT(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
72			_RL Tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
73			INTEGER tracerIdentity
74	adcroft	1.3	INTEGER advectionScheme
75	adcroft	1.11	INTEGER myThid
76
77			C !OUTPUT PARAMETERS: ==================================================
78			C gTracer :: tendancy array
79			C fVerT :: 2 1/2D arrays for vertical advective flux
80			_RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
81	adcroft	1.1	_RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
82
83	adcroft	1.11	C !LOCAL VARIABLES: ====================================================
84			C i,j :: loop indices
85			C df4 :: used for storing del^2 T for bi-harmonic term
86			C fZon :: zonal flux
87			C fmer :: meridional flux
88			C af :: advective flux
89			C df :: diffusive flux
90			C localT :: local copy of tracer field
91	adcroft	1.1	INTEGER i,j
92			_RL df4 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93			_RL fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94			_RL fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
95			_RL af (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
96			_RL df (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97			_RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98	adcroft	1.11	CEOP
99	adcroft	1.1
100			#ifdef ALLOW_AUTODIFF_TAMC
101			C-- only the kUp part of fverT is set in this subroutine
102			C-- the kDown is still required
103			fVerT(1,1,kDown) = fVerT(1,1,kDown)
104			#endif
105			DO j=1-OLy,sNy+OLy
106			DO i=1-OLx,sNx+OLx
107	heimbach	1.12	fZon(i,j) = 0. _d 0
108			fMer(i,j) = 0. _d 0
109			fVerT(i,j,kUp) = 0. _d 0
110	adcroft	1.1	ENDDO
111			ENDDO
112
113			C-- Make local copy of tracer array
114			DO j=1-OLy,sNy+OLy
115			DO i=1-OLx,sNx+OLx
116			localT(i,j)=tracer(i,j,k,bi,bj)
117			ENDDO
118			ENDDO
119
120	adcroft	1.8	C-- Unless we have already calculated the advection terms we initialize
121			C the tendency to zero.
122			IF (.NOT. multiDimAdvection .OR.
123			& advectionScheme.EQ.ENUM_CENTERED_2ND .OR.
124			& advectionScheme.EQ.ENUM_UPWIND_3RD .OR.
125			& advectionScheme.EQ.ENUM_CENTERED_4TH ) THEN
126			DO j=1-Oly,sNy+Oly
127			DO i=1-Olx,sNx+Olx
128	heimbach	1.12	gTracer(i,j,k,bi,bj)=0. _d 0
129	adcroft	1.8	ENDDO
130			ENDDO
131			ENDIF
132	adcroft	1.1
133			C-- Pre-calculate del^2 T if bi-harmonic coefficient is non-zero
134			IF (diffK4 .NE. 0.) THEN
135			CALL GAD_GRAD_X(bi,bj,k,xA,localT,fZon,myThid)
136			CALL GAD_GRAD_Y(bi,bj,k,yA,localT,fMer,myThid)
137			CALL GAD_DEL2(bi,bj,k,fZon,fMer,df4,myThid)
138			ENDIF
139
140			C-- Initialize net flux in X direction
141			DO j=1-Oly,sNy+Oly
142			DO i=1-Olx,sNx+Olx
143	heimbach	1.12	fZon(i,j) = 0. _d 0
144	adcroft	1.1	ENDDO
145			ENDDO
146
147			C- Advective flux in X
148	adcroft	1.8	IF (.NOT. multiDimAdvection .OR.
149			& advectionScheme.EQ.ENUM_CENTERED_2ND .OR.
150			& advectionScheme.EQ.ENUM_UPWIND_3RD .OR.
151			& advectionScheme.EQ.ENUM_CENTERED_4TH ) THEN
152	adcroft	1.3	IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
153	adcroft	1.1	CALL GAD_C2_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
154	adcroft	1.3	ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
155	adcroft	1.1	CALL GAD_FLUXLIMIT_ADV_X(
156			& bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)
157	adcroft	1.3	ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
158	jmc	1.2	CALL GAD_U3_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
159	adcroft	1.3	ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
160	adcroft	1.1	CALL GAD_C4_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
161	adcroft	1.4	ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
162			CALL GAD_DST3_ADV_X(
163			& bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)
164			ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
165			CALL GAD_DST3FL_ADV_X(
166			& bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)
167	adcroft	1.1	ELSE
168	adcroft	1.3	STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'
169	adcroft	1.1	ENDIF
170	adcroft	1.5	DO j=1-Oly,sNy+Oly
171			DO i=1-Olx,sNx+Olx
172	adcroft	1.1	fZon(i,j) = fZon(i,j) + af(i,j)
173			ENDDO
174			ENDDO
175	adcroft	1.8	ENDIF
176	adcroft	1.1
177			C- Diffusive flux in X
178			IF (diffKh.NE.0.) THEN
179			CALL GAD_DIFF_X(bi,bj,k,xA,diffKh,localT,df,myThid)
180			ELSE
181	adcroft	1.5	DO j=1-Oly,sNy+Oly
182			DO i=1-Olx,sNx+Olx
183	heimbach	1.12	df(i,j) = 0. _d 0
184	adcroft	1.1	ENDDO
185			ENDDO
186			ENDIF
187
188			#ifdef ALLOW_GMREDI
189			C- GM/Redi flux in X
190			IF (useGMRedi) THEN
191			C note should update GMREDI_XTRANSPORT to use localT and set df aja
192			CALL GMREDI_XTRANSPORT(
193			I iMin,iMax,jMin,jMax,bi,bj,K,
194			I xA,Tracer,
195			U df,
196			I myThid)
197			ENDIF
198			#endif
199	adcroft	1.5	DO j=1-Oly,sNy+Oly
200			DO i=1-Olx,sNx+Olx
201	adcroft	1.1	fZon(i,j) = fZon(i,j) + df(i,j)
202			ENDDO
203			ENDDO
204
205			C- Bi-harmonic duffusive flux in X
206			IF (diffK4 .NE. 0.) THEN
207			CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)
208	adcroft	1.5	DO j=1-Oly,sNy+Oly
209			DO i=1-Olx,sNx+Olx
210	adcroft	1.1	fZon(i,j) = fZon(i,j) + df(i,j)
211			ENDDO
212			ENDDO
213			ENDIF
214
215			C-- Initialize net flux in Y direction
216			DO j=1-Oly,sNy+Oly
217			DO i=1-Olx,sNx+Olx
218	heimbach	1.12	fMer(i,j) = 0. _d 0
219	adcroft	1.1	ENDDO
220			ENDDO
221
222			C- Advective flux in Y
223	adcroft	1.8	IF (.NOT. multiDimAdvection .OR.
224			& advectionScheme.EQ.ENUM_CENTERED_2ND .OR.
225			& advectionScheme.EQ.ENUM_UPWIND_3RD .OR.
226			& advectionScheme.EQ.ENUM_CENTERED_4TH ) THEN
227	adcroft	1.3	IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
228	adcroft	1.1	CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
229	adcroft	1.3	ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
230	adcroft	1.1	CALL GAD_FLUXLIMIT_ADV_Y(
231			& bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)
232	adcroft	1.3	ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
233	jmc	1.2	CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
234	adcroft	1.3	ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
235	adcroft	1.1	CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
236	adcroft	1.4	ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
237			CALL GAD_DST3_ADV_Y(
238			& bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)
239			ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
240			CALL GAD_DST3FL_ADV_Y(
241			& bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)
242	adcroft	1.1	ELSE
243	adcroft	1.3	STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'
244	adcroft	1.1	ENDIF
245			DO j=1-Oly,sNy+Oly
246			DO i=1-Olx,sNx+Olx
247			fMer(i,j) = fMer(i,j) + af(i,j)
248			ENDDO
249			ENDDO
250	adcroft	1.8	ENDIF
251	adcroft	1.1
252			C- Diffusive flux in Y
253			IF (diffKh.NE.0.) THEN
254			CALL GAD_DIFF_Y(bi,bj,k,yA,diffKh,localT,df,myThid)
255			ELSE
256			DO j=1-Oly,sNy+Oly
257			DO i=1-Olx,sNx+Olx
258	heimbach	1.12	df(i,j) = 0. _d 0
259	adcroft	1.1	ENDDO
260			ENDDO
261			ENDIF
262
263			#ifdef ALLOW_GMREDI
264			C- GM/Redi flux in Y
265			IF (useGMRedi) THEN
266	heimbach	1.7	C note should update GMREDI_YTRANSPORT to use localT and set df aja
267	adcroft	1.1	CALL GMREDI_YTRANSPORT(
268			I iMin,iMax,jMin,jMax,bi,bj,K,
269			I yA,Tracer,
270			U df,
271			I myThid)
272			ENDIF
273			#endif
274			DO j=1-Oly,sNy+Oly
275			DO i=1-Olx,sNx+Olx
276			fMer(i,j) = fMer(i,j) + df(i,j)
277			ENDDO
278			ENDDO
279
280			C- Bi-harmonic flux in Y
281			IF (diffK4 .NE. 0.) THEN
282			CALL GAD_BIHARM_Y(bi,bj,k,yA,df4,diffK4,df,myThid)
283			DO j=1-Oly,sNy+Oly
284			DO i=1-Olx,sNx+Olx
285			fMer(i,j) = fMer(i,j) + df(i,j)
286			ENDDO
287			ENDDO
288			ENDIF
289
290			C- Advective flux in R
291	adcroft	1.8	IF (.NOT. multiDimAdvection .OR.
292			& advectionScheme.EQ.ENUM_CENTERED_2ND .OR.
293			& advectionScheme.EQ.ENUM_UPWIND_3RD .OR.
294			& advectionScheme.EQ.ENUM_CENTERED_4TH ) THEN
295	jmc	1.2	C Note: wVel needs to be masked
296			IF (K.GE.2) THEN
297			C- Compute vertical advective flux in the interior:
298	adcroft	1.3	IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
299	jmc	1.2	CALL GAD_C2_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
300	adcroft	1.3	ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
301	jmc	1.2	CALL GAD_FLUXLIMIT_ADV_R(
302	adcroft	1.1	& bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)
303	adcroft	1.3	ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
304	jmc	1.2	CALL GAD_U3_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
305	adcroft	1.3	ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
306	jmc	1.2	CALL GAD_C4_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
307	adcroft	1.4	ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
308	adcroft	1.9	CALL GAD_DST3_ADV_R(
309			& bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)
310	adcroft	1.4	ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
311	adcroft	1.9	CALL GAD_DST3FL_ADV_R(
312			& bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)
313	jmc	1.2	ELSE
314	adcroft	1.3	STOP 'GAD_CALC_RHS: Bad advectionScheme (R)'
315	jmc	1.2	ENDIF
316			C- Surface "correction" term at k>1 :
317			DO j=1-Oly,sNy+Oly
318			DO i=1-Olx,sNx+Olx
319			af(i,j) = af(i,j)
320			& + (maskC(i,j,k,bi,bj)-maskC(i,j,k-1,bi,bj))*
321			& rTrans(i,j)*Tracer(i,j,k,bi,bj)
322			ENDDO
323			ENDDO
324	adcroft	1.1	ELSE
325	jmc	1.2	C- Surface "correction" term at k=1 :
326			DO j=1-Oly,sNy+Oly
327			DO i=1-Olx,sNx+Olx
328			af(i,j) = rTrans(i,j)*Tracer(i,j,k,bi,bj)
329			ENDDO
330			ENDDO
331	adcroft	1.1	ENDIF
332	jmc	1.2	C- add the advective flux to fVerT
333	adcroft	1.5	DO j=1-Oly,sNy+Oly
334			DO i=1-Olx,sNx+Olx
335	adcroft	1.11	fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)
336	adcroft	1.1	ENDDO
337			ENDDO
338	adcroft	1.8	ENDIF
339	adcroft	1.1
340			C- Diffusive flux in R
341			C Note: For K=1 then KM1=1 and this gives a dT/dr = 0 upper
342			C boundary condition.
343			IF (implicitDiffusion) THEN
344	adcroft	1.5	DO j=1-Oly,sNy+Oly
345			DO i=1-Olx,sNx+Olx
346	heimbach	1.12	df(i,j) = 0. _d 0
347	adcroft	1.1	ENDDO
348			ENDDO
349			ELSE
350			CALL GAD_DIFF_R(bi,bj,k,KappaRT,tracer,df,myThid)
351			ENDIF
352	adcroft	1.5	c DO j=1-Oly,sNy+Oly
353			c DO i=1-Olx,sNx+Olx
354	adcroft	1.11	c fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)
355	adcroft	1.1	c ENDDO
356			c ENDDO
357
358			#ifdef ALLOW_GMREDI
359			C- GM/Redi flux in R
360			IF (useGMRedi) THEN
361			C note should update GMREDI_RTRANSPORT to set df aja
362			CALL GMREDI_RTRANSPORT(
363			I iMin,iMax,jMin,jMax,bi,bj,K,
364	adcroft	1.6	I Tracer,
365	adcroft	1.1	U df,
366			I myThid)
367	adcroft	1.5	c DO j=1-Oly,sNy+Oly
368			c DO i=1-Olx,sNx+Olx
369	adcroft	1.11	c fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)
370	adcroft	1.1	c ENDDO
371			c ENDDO
372			ENDIF
373			#endif
374
375	adcroft	1.5	DO j=1-Oly,sNy+Oly
376			DO i=1-Olx,sNx+Olx
377	adcroft	1.11	fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)
378	adcroft	1.1	ENDDO
379			ENDDO
380
381			#ifdef ALLOW_KPP
382			C- Add non local KPP transport term (ghat) to diffusive T flux.
383			IF (useKPP) THEN
384	adcroft	1.5	DO j=1-Oly,sNy+Oly
385			DO i=1-Olx,sNx+Olx
386	heimbach	1.12	df(i,j) = 0. _d 0
387	adcroft	1.1	ENDDO
388			ENDDO
389			IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
390			C note should update KPP_TRANSPORT_T to set df aja
391			CALL KPP_TRANSPORT_T(
392			I iMin,iMax,jMin,jMax,bi,bj,k,km1,
393			I KappaRT,
394			U df )
395			ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
396			CALL KPP_TRANSPORT_S(
397			I iMin,iMax,jMin,jMax,bi,bj,k,km1,
398			I KappaRT,
399			U df )
400			ELSE
401			STOP 'GAD_CALC_RHS: Ooops'
402			ENDIF
403	adcroft	1.5	DO j=1-Oly,sNy+Oly
404			DO i=1-Olx,sNx+Olx
405	adcroft	1.11	fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)
406	adcroft	1.1	ENDDO
407			ENDDO
408			ENDIF
409			#endif
410
411			C-- Divergence of fluxes
412	adcroft	1.10	DO j=1-Oly,sNy+Oly-1
413			DO i=1-Olx,sNx+Olx-1
414	adcroft	1.8	gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)
415	adcroft	1.1	& -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
416			& *recip_rA(i,j,bi,bj)
417			& *(
418			& +( fZon(i+1,j)-fZon(i,j) )
419			& +( fMer(i,j+1)-fMer(i,j) )
420			& +( fVerT(i,j,kUp)-fVerT(i,j,kDown) )*rkFac
421			& )
422			ENDDO
423			ENDDO
424
425			RETURN
426			END