pkg/generic_advdiff/gad_calc_rhs.F

C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_calc_rhs.F,v 1.22 2003/09/25 03:01:59 dimitri 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,rTransKp1,maskUp,
     I           uVel, vVel, wVel,
     I           diffKh, diffK4, KappaRT, Tracer,
     I           tracerIdentity, advectionScheme, 
     I           calcAdvection, implicitAdvection,
     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 "SURFACE.h"
#include "GAD.h"

#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  rTransKp1            :: 2-D array of volume transport at W pt, interface k+1
C  maskUp               :: 2-D array for mask at W points
C  uVel, vVel, wVel     :: 3 components of the velcity field (3-D array)
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  implicitAdvection    :: =True if vertical Advec term is computed implicitly
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)
      _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _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 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
      LOGICAL implicitAdvection
      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)
      _RL advFac, rAdvFac
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

      advFac  = 0. _d 0
      IF (calcAdvection) advFac = 1. _d 0
      rAdvFac = rkFac*advFac
      IF (implicitAdvection) rAdvFac = 0. _d 0

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

C--   Compute vertical flux fVerT(kUp) at interface k (between k-1 & k):
C-    Advective flux in R
      IF (calcAdvection .AND. .NOT.implicitAdvection .AND. 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-     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) 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
     &     -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))
     &                   +(vTrans(i,j+1)-vTrans(i,j))
     &                   +(rTrans(i,j)-rTransKp1(i,j))*rAdvFac
     &                  )*advFac
     &    )
       ENDDO
      ENDDO

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