pkg/generic_advdiff/gad_calc_rhs.F

C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_calc_rhs.F,v 1.26 2004/06/26 02:38:54 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,rTransKp1,maskUp,
     I           uVel, vVel, wVel,
     I           diffKh, diffK4, KappaRT, Tracer,
     I           tracerIdentity, advectionScheme, vertAdvecScheme,
     I           calcAdvection, implicitAdvection,
     U           fVerT, gTracer,
     I           myTime, myIter, 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        :: loop range for called routines
C jMin,jMax        :: loop range for called routines
C kup              :: index into 2 1/2D array, toggles between 1|2
C kdown            :: index into 2 1/2D array, toggles between 2|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    :: 2-D arrays of volume transports at U,V points
C rTrans           :: 2-D arrays of volume transports at W points
C rTransKp1        :: 2-D array of volume trans at W pts, interf 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   :: tracer identifier (required for KPP,GM)
C advectionScheme  :: advection scheme to use (Horizontal plane)
C vertAdvecScheme  :: advection scheme to use (Vertical direction)
C calcAdvection    :: =False if Advec computed with multiDim scheme
C implicitAdvection:: =True if vertical Advec computed implicitly
C myTime           :: current time
C myIter           :: iteration number
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, vertAdvecScheme
      LOGICAL calcAdvection
      LOGICAL implicitAdvection
      _RL     myTime
      INTEGER myIter, 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,
     I          uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), localT,
     O          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,
     I          uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), localT,
     O          af, myThid )
      ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
       CALL GAD_DST3FL_ADV_X( bi,bj,k, deltaTtracer,
     I          uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), localT,
     O          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,
     I          vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), localT,
     O          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,
     I          vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), localT,
     O          af, myThid )
      ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
       CALL GAD_DST3FL_ADV_Y( bi,bj,k, deltaTtracer,
     I          vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), localT,
     O          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
#ifdef ALLOW_AIM
C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
      IF (calcAdvection .AND. .NOT.implicitAdvection .AND. K.GE.2 .AND.
     &     (.NOT.useAIM .OR.tracerIdentity.NE.GAD_SALINITY .OR.K.LT.Nr)
     &   ) THEN
#else
      IF (calcAdvection .AND. .NOT.implicitAdvection .AND. K.GE.2) THEN
#endif
C-    Compute vertical advective flux in the interior:
       IF (vertAdvecScheme.EQ.ENUM_CENTERED_2ND) THEN
        CALL GAD_C2_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
       ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
        CALL GAD_FLUXLIMIT_ADV_R(
     &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)
       ELSEIF (vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN
        CALL GAD_U3_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
       ELSEIF (vertAdvecScheme.EQ.ENUM_CENTERED_4TH) THEN
        CALL GAD_C4_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
       ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
        CALL GAD_DST3_ADV_R(
     &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)
       ELSEIF (vertAdvecScheme.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 vertAdvecScheme (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

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