pkg/generic_advdiff/gad_calc_rhs.F

C $Header: /usr/local/gcmpack/MITgcm/pkg/generic_advdiff/gad_calc_rhs.F,v 1.20 2003/09/22 22:13:11 jmc 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, calcAdvection,
     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 "SURFACE.h"
#include "GAD.h"
#ifdef ALLOW_PTRACERS
#include "PTRACERS_OPTIONS.h"
#include "PTRACERS.h"
#endif

#ifdef ALLOW_AUTODIFF_TAMC
#include "tamc.h"
#include "tamc_keys.h"
#endif /* ALLOW_AUTODIFF_TAMC */

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 for KPP and GM)
C  advectionScheme      :: advection scheme to use
C  calcAdvection        :: =False if Advec terms computed with multiDim scheme
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
      LOGICAL calcAdvection
      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
        df(i,j)        = 0. _d 0
        df4(i,j)       = 0. _d 0
        localT(i,j)    = 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.
C     <== now done earlier at the beginning of thermodynamics.
c     IF (calcAdvection) THEN
c      DO j=1-Oly,sNy+Oly
c       DO i=1-Olx,sNx+Olx
c        gTracer(i,j,k,bi,bj)=0. _d 0
c       ENDDO
c      ENDDO
c     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 (calcAdvection) 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,tracerIdentity,
     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 (calcAdvection) 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,tracerIdentity,
     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

#ifdef NONLIN_FRSURF
C--   Compute vertical flux fVerT(kDown) at interface k+1 (between k & k+1):
      IF ( calcAdvection .AND. K.EQ.Nr .AND.
     &     useRealFreshWaterFlux .AND.
     &     buoyancyRelation .EQ. 'OCEANICP' ) THEN   
       DO j=1-Oly,sNy+Oly
        DO i=1-Olx,sNx+Olx
         fVerT(i,j,kDown) = convertEmP2rUnit*PmEpR(i,j,bi,bj)
     &     *rA(i,j,bi,bj)*maskC(i,j,k,bi,bj)*Tracer(i,j,k,bi,bj)
        ENDDO
       ENDDO 
      ENDIF
#endif /* NONLIN_FRSURF */

C--   Compute vertical flux fVerT(kUp) at interface k (between k-1 & k):
C-    Advective flux in R
      IF (calcAdvection) 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

#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,tracerIdentity,
     U     df,
     I     myThid)
      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 )
#ifdef ALLOW_PTRACERS
       ELSEIF (tracerIdentity .GE. GAD_TR1 .AND.
     &         tracerIdentity .LE. (GAD_TR1+PTRACERS_numInUse-1)) THEN
        CALL KPP_TRANSPORT_PTR(
     I     iMin,iMax,jMin,jMax,bi,bj,k,km1,
     I     tracerIdentity-GAD_TR1+1,KappaRT,
     U     df )
#endif
       ELSE
        PRINT*,'invalid tracer indentity: ', tracerIdentity
        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

#ifdef NONLIN_FRSURF
C-- account for 3.D divergence of the flow in rStar coordinate:
      IF (calcAdvection .AND. select_rStar.GT.0) THEN
       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)
     &     - (rStarExpC(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf
     &       *tracer(i,j,k,bi,bj)*maskC(i,j,k,bi,bj)
        ENDDO
       ENDDO
      ENDIF
      IF (calcAdvection .AND. select_rStar.LT.0) THEN
       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)
     &     - rStarDhCDt(i,j,bi,bj)
     &       *tracer(i,j,k,bi,bj)*maskC(i,j,k,bi,bj)
        ENDDO
       ENDDO
      ENDIF
#endif /* NONLIN_FRSURF */
      

      RETURN
      END
1	dimitri	1.21	C $Header: /usr/local/gcmpack/MITgcm/pkg/generic_advdiff/gad_calc_rhs.F,v 1.20 2003/09/22 22:13:11 jmc 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	jmc	1.14	I tracerIdentity, advectionScheme, calcAdvection,
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	jmc	1.16	#include "SURFACE.h"
45	adcroft	1.1	#include "GAD.h"
46	dimitri	1.19	#ifdef ALLOW_PTRACERS
47			#include "PTRACERS_OPTIONS.h"
48			#include "PTRACERS.h"
49			#endif
50	adcroft	1.1
51	heimbach	1.13	#ifdef ALLOW_AUTODIFF_TAMC
52			#include "tamc.h"
53			#include "tamc_keys.h"
54			#endif /* ALLOW_AUTODIFF_TAMC */
55
56	adcroft	1.11	C !INPUT PARAMETERS: ===================================================
57			C bi,bj :: tile indices
58			C iMin,iMax,jMin,jMax :: loop range for called routines
59			C kup :: index into 2 1/2D array, toggles between 1 and 2
60			C kdown :: index into 2 1/2D array, toggles between 2 and 1
61			C kp1 :: =k+1 for k<Nr, =Nr for k=Nr
62			C xA,yA :: areas of X and Y face of tracer cells
63			C uTrans,vTrans,rTrans :: 2-D arrays of volume transports at U,V and W points
64			C maskUp :: 2-D array for mask at W points
65			C diffKh :: horizontal diffusion coefficient
66			C diffK4 :: bi-harmonic diffusion coefficient
67			C KappaRT :: 3-D array for vertical diffusion coefficient
68			C Tracer :: tracer field
69	dimitri	1.19	C tracerIdentity :: identifier for the tracer (required for KPP and GM)
70	adcroft	1.11	C advectionScheme :: advection scheme to use
71	jmc	1.14	C calcAdvection :: =False if Advec terms computed with multiDim scheme
72	adcroft	1.11	C myThid :: thread number
73			INTEGER bi,bj,iMin,iMax,jMin,jMax
74	adcroft	1.1	INTEGER k,kUp,kDown,kM1
75			_RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76			_RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
77			_RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78			_RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79			_RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
80			_RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
81			_RL diffKh, diffK4
82			_RL KappaRT(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
83			_RL Tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
84			INTEGER tracerIdentity
85	adcroft	1.3	INTEGER advectionScheme
86	jmc	1.14	LOGICAL calcAdvection
87	adcroft	1.11	INTEGER myThid
88
89			C !OUTPUT PARAMETERS: ==================================================
90			C gTracer :: tendancy array
91			C fVerT :: 2 1/2D arrays for vertical advective flux
92			_RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
93	adcroft	1.1	_RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
94
95	adcroft	1.11	C !LOCAL VARIABLES: ====================================================
96			C i,j :: loop indices
97			C df4 :: used for storing del^2 T for bi-harmonic term
98			C fZon :: zonal flux
99			C fmer :: meridional flux
100			C af :: advective flux
101			C df :: diffusive flux
102			C localT :: local copy of tracer field
103	adcroft	1.1	INTEGER i,j
104			_RL df4 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105			_RL fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
106			_RL fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
107			_RL af (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
108			_RL df (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
109			_RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
110	adcroft	1.11	CEOP
111	adcroft	1.1
112			#ifdef ALLOW_AUTODIFF_TAMC
113			C-- only the kUp part of fverT is set in this subroutine
114			C-- the kDown is still required
115			fVerT(1,1,kDown) = fVerT(1,1,kDown)
116			#endif
117	heimbach	1.13
118	adcroft	1.1	DO j=1-OLy,sNy+OLy
119			DO i=1-OLx,sNx+OLx
120	heimbach	1.12	fZon(i,j) = 0. _d 0
121			fMer(i,j) = 0. _d 0
122			fVerT(i,j,kUp) = 0. _d 0
123	heimbach	1.13	df(i,j) = 0. _d 0
124			df4(i,j) = 0. _d 0
125			localT(i,j) = 0. _d 0
126	adcroft	1.1	ENDDO
127			ENDDO
128
129			C-- Make local copy of tracer array
130			DO j=1-OLy,sNy+OLy
131			DO i=1-OLx,sNx+OLx
132			localT(i,j)=tracer(i,j,k,bi,bj)
133			ENDDO
134			ENDDO
135
136	adcroft	1.8	C-- Unless we have already calculated the advection terms we initialize
137			C the tendency to zero.
138	jmc	1.20	C <== now done earlier at the beginning of thermodynamics.
139			c IF (calcAdvection) THEN
140			c DO j=1-Oly,sNy+Oly
141			c DO i=1-Olx,sNx+Olx
142			c gTracer(i,j,k,bi,bj)=0. _d 0
143			c ENDDO
144			c ENDDO
145			c ENDIF
146	adcroft	1.1
147			C-- Pre-calculate del^2 T if bi-harmonic coefficient is non-zero
148			IF (diffK4 .NE. 0.) THEN
149			CALL GAD_GRAD_X(bi,bj,k,xA,localT,fZon,myThid)
150			CALL GAD_GRAD_Y(bi,bj,k,yA,localT,fMer,myThid)
151			CALL GAD_DEL2(bi,bj,k,fZon,fMer,df4,myThid)
152			ENDIF
153
154			C-- Initialize net flux in X direction
155			DO j=1-Oly,sNy+Oly
156			DO i=1-Olx,sNx+Olx
157	heimbach	1.12	fZon(i,j) = 0. _d 0
158	adcroft	1.1	ENDDO
159			ENDDO
160
161			C- Advective flux in X
162	jmc	1.14	IF (calcAdvection) THEN
163	adcroft	1.3	IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
164	adcroft	1.1	CALL GAD_C2_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
165	adcroft	1.3	ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
166	adcroft	1.1	CALL GAD_FLUXLIMIT_ADV_X(
167			& bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)
168	adcroft	1.3	ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
169	jmc	1.2	CALL GAD_U3_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
170	adcroft	1.3	ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
171	adcroft	1.1	CALL GAD_C4_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
172	adcroft	1.4	ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
173			CALL GAD_DST3_ADV_X(
174			& bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)
175			ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
176			CALL GAD_DST3FL_ADV_X(
177			& bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)
178	adcroft	1.1	ELSE
179	adcroft	1.3	STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'
180	adcroft	1.1	ENDIF
181	adcroft	1.5	DO j=1-Oly,sNy+Oly
182			DO i=1-Olx,sNx+Olx
183	adcroft	1.1	fZon(i,j) = fZon(i,j) + af(i,j)
184			ENDDO
185			ENDDO
186	adcroft	1.8	ENDIF
187	adcroft	1.1
188			C- Diffusive flux in X
189			IF (diffKh.NE.0.) THEN
190			CALL GAD_DIFF_X(bi,bj,k,xA,diffKh,localT,df,myThid)
191			ELSE
192	adcroft	1.5	DO j=1-Oly,sNy+Oly
193			DO i=1-Olx,sNx+Olx
194	heimbach	1.12	df(i,j) = 0. _d 0
195	adcroft	1.1	ENDDO
196			ENDDO
197			ENDIF
198
199			#ifdef ALLOW_GMREDI
200			C- GM/Redi flux in X
201			IF (useGMRedi) THEN
202			C note should update GMREDI_XTRANSPORT to use localT and set df aja
203			CALL GMREDI_XTRANSPORT(
204			I iMin,iMax,jMin,jMax,bi,bj,K,
205	heimbach	1.15	I xA,Tracer,tracerIdentity,
206	adcroft	1.1	U df,
207			I myThid)
208			ENDIF
209			#endif
210	adcroft	1.5	DO j=1-Oly,sNy+Oly
211			DO i=1-Olx,sNx+Olx
212	adcroft	1.1	fZon(i,j) = fZon(i,j) + df(i,j)
213			ENDDO
214			ENDDO
215
216			C- Bi-harmonic duffusive flux in X
217			IF (diffK4 .NE. 0.) THEN
218			CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)
219	adcroft	1.5	DO j=1-Oly,sNy+Oly
220			DO i=1-Olx,sNx+Olx
221	adcroft	1.1	fZon(i,j) = fZon(i,j) + df(i,j)
222			ENDDO
223			ENDDO
224			ENDIF
225
226			C-- Initialize net flux in Y direction
227			DO j=1-Oly,sNy+Oly
228			DO i=1-Olx,sNx+Olx
229	heimbach	1.12	fMer(i,j) = 0. _d 0
230	adcroft	1.1	ENDDO
231			ENDDO
232
233			C- Advective flux in Y
234	jmc	1.14	IF (calcAdvection) THEN
235	adcroft	1.3	IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
236	adcroft	1.1	CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
237	adcroft	1.3	ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
238	adcroft	1.1	CALL GAD_FLUXLIMIT_ADV_Y(
239			& bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)
240	adcroft	1.3	ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
241	jmc	1.2	CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
242	adcroft	1.3	ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
243	adcroft	1.1	CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
244	adcroft	1.4	ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
245			CALL GAD_DST3_ADV_Y(
246			& bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)
247			ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
248			CALL GAD_DST3FL_ADV_Y(
249			& bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)
250	adcroft	1.1	ELSE
251	adcroft	1.3	STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'
252	adcroft	1.1	ENDIF
253			DO j=1-Oly,sNy+Oly
254			DO i=1-Olx,sNx+Olx
255			fMer(i,j) = fMer(i,j) + af(i,j)
256			ENDDO
257			ENDDO
258	adcroft	1.8	ENDIF
259	adcroft	1.1
260			C- Diffusive flux in Y
261			IF (diffKh.NE.0.) THEN
262			CALL GAD_DIFF_Y(bi,bj,k,yA,diffKh,localT,df,myThid)
263			ELSE
264			DO j=1-Oly,sNy+Oly
265			DO i=1-Olx,sNx+Olx
266	heimbach	1.12	df(i,j) = 0. _d 0
267	adcroft	1.1	ENDDO
268			ENDDO
269			ENDIF
270
271			#ifdef ALLOW_GMREDI
272			C- GM/Redi flux in Y
273			IF (useGMRedi) THEN
274	heimbach	1.7	C note should update GMREDI_YTRANSPORT to use localT and set df aja
275	adcroft	1.1	CALL GMREDI_YTRANSPORT(
276			I iMin,iMax,jMin,jMax,bi,bj,K,
277	heimbach	1.15	I yA,Tracer,tracerIdentity,
278	adcroft	1.1	U df,
279			I myThid)
280			ENDIF
281			#endif
282			DO j=1-Oly,sNy+Oly
283			DO i=1-Olx,sNx+Olx
284			fMer(i,j) = fMer(i,j) + df(i,j)
285			ENDDO
286			ENDDO
287
288			C- Bi-harmonic flux in Y
289			IF (diffK4 .NE. 0.) THEN
290			CALL GAD_BIHARM_Y(bi,bj,k,yA,df4,diffK4,df,myThid)
291			DO j=1-Oly,sNy+Oly
292			DO i=1-Olx,sNx+Olx
293			fMer(i,j) = fMer(i,j) + df(i,j)
294			ENDDO
295			ENDDO
296			ENDIF
297
298	jmc	1.16	#ifdef NONLIN_FRSURF
299			C-- Compute vertical flux fVerT(kDown) at interface k+1 (between k & k+1):
300			IF ( calcAdvection .AND. K.EQ.Nr .AND.
301			& useRealFreshWaterFlux .AND.
302			& buoyancyRelation .EQ. 'OCEANICP' ) THEN
303			DO j=1-Oly,sNy+Oly
304			DO i=1-Olx,sNx+Olx
305			fVerT(i,j,kDown) = convertEmP2rUnit*PmEpR(i,j,bi,bj)
306			& rA(i,j,bi,bj)maskC(i,j,k,bi,bj)*Tracer(i,j,k,bi,bj)
307			ENDDO
308			ENDDO
309			ENDIF
310			#endif /* NONLIN_FRSURF */
311
312			C-- Compute vertical flux fVerT(kUp) at interface k (between k-1 & k):
313	adcroft	1.1	C- Advective flux in R
314	jmc	1.14	IF (calcAdvection) THEN
315	jmc	1.2	C Note: wVel needs to be masked
316			IF (K.GE.2) THEN
317			C- Compute vertical advective flux in the interior:
318	adcroft	1.3	IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
319	jmc	1.2	CALL GAD_C2_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
320	adcroft	1.3	ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
321	jmc	1.2	CALL GAD_FLUXLIMIT_ADV_R(
322	adcroft	1.1	& bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)
323	adcroft	1.3	ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
324	jmc	1.2	CALL GAD_U3_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
325	adcroft	1.3	ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
326	jmc	1.2	CALL GAD_C4_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
327	adcroft	1.4	ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
328	adcroft	1.9	CALL GAD_DST3_ADV_R(
329			& bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)
330	adcroft	1.4	ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
331	adcroft	1.9	CALL GAD_DST3FL_ADV_R(
332			& bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)
333	jmc	1.2	ELSE
334	adcroft	1.3	STOP 'GAD_CALC_RHS: Bad advectionScheme (R)'
335	jmc	1.2	ENDIF
336			C- Surface "correction" term at k>1 :
337			DO j=1-Oly,sNy+Oly
338			DO i=1-Olx,sNx+Olx
339			af(i,j) = af(i,j)
340			& + (maskC(i,j,k,bi,bj)-maskC(i,j,k-1,bi,bj))*
341			& rTrans(i,j)*Tracer(i,j,k,bi,bj)
342			ENDDO
343			ENDDO
344	adcroft	1.1	ELSE
345	jmc	1.2	C- Surface "correction" term at k=1 :
346			DO j=1-Oly,sNy+Oly
347			DO i=1-Olx,sNx+Olx
348			af(i,j) = rTrans(i,j)*Tracer(i,j,k,bi,bj)
349			ENDDO
350			ENDDO
351	adcroft	1.1	ENDIF
352	jmc	1.2	C- add the advective flux to fVerT
353	adcroft	1.5	DO j=1-Oly,sNy+Oly
354			DO i=1-Olx,sNx+Olx
355	adcroft	1.11	fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)
356	adcroft	1.1	ENDDO
357			ENDDO
358	adcroft	1.8	ENDIF
359	adcroft	1.1
360			C- Diffusive flux in R
361			C Note: For K=1 then KM1=1 and this gives a dT/dr = 0 upper
362			C boundary condition.
363			IF (implicitDiffusion) THEN
364	adcroft	1.5	DO j=1-Oly,sNy+Oly
365			DO i=1-Olx,sNx+Olx
366	heimbach	1.12	df(i,j) = 0. _d 0
367	adcroft	1.1	ENDDO
368			ENDDO
369			ELSE
370			CALL GAD_DIFF_R(bi,bj,k,KappaRT,tracer,df,myThid)
371			ENDIF
372
373			#ifdef ALLOW_GMREDI
374			C- GM/Redi flux in R
375			IF (useGMRedi) THEN
376			C note should update GMREDI_RTRANSPORT to set df aja
377			CALL GMREDI_RTRANSPORT(
378			I iMin,iMax,jMin,jMax,bi,bj,K,
379	heimbach	1.15	I Tracer,tracerIdentity,
380	adcroft	1.1	U df,
381			I myThid)
382			ENDIF
383			#endif
384
385	adcroft	1.5	DO j=1-Oly,sNy+Oly
386			DO i=1-Olx,sNx+Olx
387	adcroft	1.11	fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)
388	adcroft	1.1	ENDDO
389			ENDDO
390
391			#ifdef ALLOW_KPP
392			C- Add non local KPP transport term (ghat) to diffusive T flux.
393			IF (useKPP) THEN
394	adcroft	1.5	DO j=1-Oly,sNy+Oly
395			DO i=1-Olx,sNx+Olx
396	heimbach	1.12	df(i,j) = 0. _d 0
397	adcroft	1.1	ENDDO
398			ENDDO
399			IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
400			C note should update KPP_TRANSPORT_T to set df aja
401			CALL KPP_TRANSPORT_T(
402			I iMin,iMax,jMin,jMax,bi,bj,k,km1,
403			I KappaRT,
404			U df )
405			ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
406			CALL KPP_TRANSPORT_S(
407			I iMin,iMax,jMin,jMax,bi,bj,k,km1,
408			I KappaRT,
409			U df )
410	mlosch	1.18	#ifdef ALLOW_PTRACERS
411	dimitri	1.19	ELSEIF (tracerIdentity .GE. GAD_TR1 .AND.
412			& tracerIdentity .LE. (GAD_TR1+PTRACERS_numInUse-1)) THEN
413	mlosch	1.18	CALL KPP_TRANSPORT_PTR(
414			I iMin,iMax,jMin,jMax,bi,bj,k,km1,
415	dimitri	1.21	I tracerIdentity-GAD_TR1+1,KappaRT,
416	mlosch	1.18	U df )
417			#endif
418	adcroft	1.1	ELSE
419	mlosch	1.18	PRINT*,'invalid tracer indentity: ', tracerIdentity
420	adcroft	1.1	STOP 'GAD_CALC_RHS: Ooops'
421			ENDIF
422	adcroft	1.5	DO j=1-Oly,sNy+Oly
423			DO i=1-Olx,sNx+Olx
424	adcroft	1.11	fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)
425	adcroft	1.1	ENDDO
426			ENDDO
427			ENDIF
428			#endif
429
430			C-- Divergence of fluxes
431	adcroft	1.10	DO j=1-Oly,sNy+Oly-1
432			DO i=1-Olx,sNx+Olx-1
433	adcroft	1.8	gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)
434	adcroft	1.1	& -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
435			& *recip_rA(i,j,bi,bj)
436			& *(
437			& +( fZon(i+1,j)-fZon(i,j) )
438			& +( fMer(i,j+1)-fMer(i,j) )
439			& +( fVerT(i,j,kUp)-fVerT(i,j,kDown) )*rkFac
440			& )
441			ENDDO
442			ENDDO
443	jmc	1.17
444			#ifdef NONLIN_FRSURF
445			C-- account for 3.D divergence of the flow in rStar coordinate:
446			IF (calcAdvection .AND. select_rStar.GT.0) THEN
447			DO j=1-Oly,sNy+Oly-1
448			DO i=1-Olx,sNx+Olx-1
449			gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)
450			& - (rStarExpC(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf
451			& tracer(i,j,k,bi,bj)maskC(i,j,k,bi,bj)
452			ENDDO
453			ENDDO
454			ENDIF
455			IF (calcAdvection .AND. select_rStar.LT.0) THEN
456			DO j=1-Oly,sNy+Oly-1
457			DO i=1-Olx,sNx+Olx-1
458			gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)
459			& - rStarDhCDt(i,j,bi,bj)
460			& tracer(i,j,k,bi,bj)maskC(i,j,k,bi,bj)
461			ENDDO
462			ENDDO
463			ENDIF
464			#endif /* NONLIN_FRSURF */
465
466	adcroft	1.1
467			RETURN
468			END