pkg/generic_advdiff/gad_calc_rhs.F

C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_calc_rhs.F,v 1.36 2005/06/22 00:27:47 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, KappaR, 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 KappaR           :: 2-D array for vertical diffusion coefficient, interf k
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 KappaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _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
#ifdef ALLOW_DIAGNOSTICS
      CHARACTER*8 diagName
      CHARACTER*4 GAD_DIAG_SUFX, diagSufx 
      EXTERNAL    GAD_DIAG_SUFX
#endif
      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

#ifdef ALLOW_DIAGNOSTICS
C--   Set diagnostic suffix for the current tracer
      IF ( useDiagnostics ) THEN
        diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )
      ENDIF
#endif

      advFac  = 0. _d 0
      IF (calcAdvection) advFac = 1. _d 0
      rAdvFac = rkSign*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_UPWIND_1RST 
     &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
          CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme, 
     I            dTtracerLev(k), uTrans, uVel, localT,
     O            af, myThid )
        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
          CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, dTtracerLev(k),
     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, dTtracerLev(k),
     I            uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), localT,
     O            af, myThid )
        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
         IF ( inAdMode ) THEN
cph This block is to trick the adjoint:
cph IF inAdExact=.FALSE., we want to use DST3 
cph with limiters in forward, but without limiters in reverse.
          CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k),
     I           uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), localT,
     O           af, myThid )
         ELSE
          CALL GAD_DST3FL_ADV_X( bi,bj,k, dTtracerLev(k),
     I           uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), localT,
     O           af, myThid )
         ENDIF
        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
#ifdef ALLOW_DIAGNOSTICS
        IF ( useDiagnostics ) THEN
          diagName = 'ADVx'//diagSufx
          CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
        ENDIF
#endif
      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

C-    Add bi-harmonic diffusive flux in X
      IF (diffK4 .NE. 0.) THEN
       CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)
      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

#ifdef ALLOW_DIAGNOSTICS
C-    Diagnostics of Tracer flux in X dir (mainly Diffusive term),
C       excluding advective terms:
      IF ( useDiagnostics .AND.
     &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. useGMRedi) ) THEN
          diagName = 'DIFx'//diagSufx
          CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
      ENDIF
#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_UPWIND_1RST 
     &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
          CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, 
     I            dTtracerLev(k), vTrans, vVel, localT,
     O            af, myThid )
        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
          CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, dTtracerLev(k),
     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, dTtracerLev(k),
     I            vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), localT,
     O            af, myThid )
        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
         IF ( inAdMode ) THEN
cph This block is to trick the adjoint:
cph IF inAdExact=.FALSE., we want to use DST3 
cph with limiters in forward, but without limiters in reverse.
          CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),
     I           vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), localT,
     O           af, myThid )
         ELSE
          CALL GAD_DST3FL_ADV_Y( bi,bj,k, dTtracerLev(k),
     I           vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), localT,
     O           af, myThid )
         ENDIF
        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
#ifdef ALLOW_DIAGNOSTICS
        IF ( useDiagnostics ) THEN
          diagName = 'ADVy'//diagSufx
          CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
        ENDIF
#endif
      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

C-    Add bi-harmonic flux in Y
      IF (diffK4 .NE. 0.) THEN
       CALL GAD_BIHARM_Y(bi,bj,k,yA,df4,diffK4,df,myThid)
      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

#ifdef ALLOW_DIAGNOSTICS
C-    Diagnostics of Tracer flux in Y dir (mainly Diffusive terms),
C       excluding advective terms:
      IF ( useDiagnostics .AND.
     &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. useGMRedi) ) THEN
          diagName = 'DIFy'//diagSufx
          CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
      ENDIF
#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_UPWIND_1RST 
     &          .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
          CALL GAD_DST2U1_ADV_R( bi,bj,k, vertAdvecScheme, 
     I         dTtracerLev(k),rTrans,wVel,Tracer(1-Olx,1-Oly,1,bi,bj),
     O         af, myThid )
        ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
          CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k,
     I         dTtracerLev(k),rTrans,wVel,Tracer(1-Olx,1-Oly,1,bi,bj),
     O         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,
     I         dTtracerLev(k),rTrans,wVel,Tracer(1-Olx,1-Oly,1,bi,bj),
     O         af, myThid )
        ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
cph This block is to trick the adjoint:
cph IF inAdExact=.FALSE., we want to use DST3 
cph with limiters in forward, but without limiters in reverse.
          IF ( inAdMode ) THEN
           CALL GAD_DST3_ADV_R( bi,bj,k,
     I         dTtracerLev(k),rTrans,wVel,Tracer(1-Olx,1-Oly,1,bi,bj),
     O         af, myThid )
          ELSE
           CALL GAD_DST3FL_ADV_R( bi,bj,k,
     I         dTtracerLev(k),rTrans,wVel,Tracer(1-Olx,1-Oly,1,bi,bj),
     O         af, myThid )
          ENDIF
        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
#ifdef ALLOW_DIAGNOSTICS
        IF ( useDiagnostics ) THEN
          diagName = 'ADVr'//diagSufx
          CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
C        does it only if k=1 (never the case here)
          IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)
        ENDIF
#endif
      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,KappaR,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_DIAGNOSTICS
C-    Diagnostics of Tracer flux in R dir (mainly Diffusive terms), 
C       Explicit terms only & excluding advective terms:
      IF ( useDiagnostics .AND.
     &    (.NOT.implicitDiffusion .OR. useGMRedi) ) THEN
          diagName = 'DFrE'//diagSufx
          CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
      ENDIF
#endif

#ifdef ALLOW_KPP
C-    Set non local KPP transport term (ghat):
      IF ( useKPP .AND. k.GE.2 ) 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
        CALL KPP_TRANSPORT_T(
     I     iMin,iMax,jMin,jMax,bi,bj,k,km1,
     O     df )
       ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
        CALL KPP_TRANSPORT_S(
     I     iMin,iMax,jMin,jMax,bi,bj,k,km1,
     O     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,
     O     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,kDown)-fVerT(i,j,kUp))*rkSign
     &     -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))
     &                   +(vTrans(i,j+1)-vTrans(i,j))
     &                   +(rTransKp1(i,j)-rTrans(i,j))*rAdvFac
     &                  )*advFac
     &    )
       ENDDO
      ENDDO

#ifdef ALLOW_DEBUG
      IF ( debugLevel .GE. debLevB
     &   .AND. tracerIdentity.EQ.GAD_TEMPERATURE
     &   .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	C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_calc_rhs.F,v 1.36 2005/06/22 00:27:47 jmc Exp $
2	C $Name: $
3
4	#include "GAD_OPTIONS.h"
5
6	CBOP
7	C !ROUTINE: GAD_CALC_RHS
8
9	C !INTERFACE: ==========================================================
10	SUBROUTINE GAD_CALC_RHS(
11	I bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
12	I xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
13	I uVel, vVel, wVel,
14	I diffKh, diffK4, KappaR, Tracer,
15	I tracerIdentity, advectionScheme, vertAdvecScheme,
16	I calcAdvection, implicitAdvection,
17	U fVerT, gTracer,
18	I myTime, myIter, myThid )
19
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	IMPLICIT NONE
41	#include "SIZE.h"
42	#include "EEPARAMS.h"
43	#include "PARAMS.h"
44	#include "GRID.h"
45	#include "SURFACE.h"
46	#include "GAD.h"
47
48	#ifdef ALLOW_AUTODIFF_TAMC
49	#include "tamc.h"
50	#include "tamc_keys.h"
51	#endif /* ALLOW_AUTODIFF_TAMC */
52
53	C !INPUT PARAMETERS: ===================================================
54	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 KappaR :: 2-D array for vertical diffusion coefficient, interf k
69	C Tracer :: tracer field
70	C tracerIdentity :: tracer identifier (required for KPP,GM)
71	C advectionScheme :: advection scheme to use (Horizontal plane)
72	C vertAdvecScheme :: advection scheme to use (Vertical direction)
73	C calcAdvection :: =False if Advec computed with multiDim scheme
74	C implicitAdvection:: =True if vertical Advec computed implicitly
75	C myTime :: current time
76	C myIter :: iteration number
77	C myThid :: thread number
78	INTEGER bi,bj,iMin,iMax,jMin,jMax
79	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	_RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
86	_RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
87	_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	_RL diffKh, diffK4
91	_RL KappaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92	_RL Tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
93	INTEGER tracerIdentity
94	INTEGER advectionScheme, vertAdvecScheme
95	LOGICAL calcAdvection
96	LOGICAL implicitAdvection
97	_RL myTime
98	INTEGER myIter, myThid
99
100	C !OUTPUT PARAMETERS: ==================================================
101	C gTracer :: tendancy array
102	C fVerT :: 2 1/2D arrays for vertical advective flux
103	_RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
104	_RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
105
106	C !LOCAL VARIABLES: ====================================================
107	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	#ifdef ALLOW_DIAGNOSTICS
115	CHARACTER*8 diagName
116	CHARACTER*4 GAD_DIAG_SUFX, diagSufx
117	EXTERNAL GAD_DIAG_SUFX
118	#endif
119	INTEGER i,j
120	_RL df4 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
121	_RL fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
122	_RL fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
123	_RL af (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
124	_RL df (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
125	_RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
126	_RL advFac, rAdvFac
127	CEOP
128
129	#ifdef ALLOW_AUTODIFF_TAMC
130	C-- only the kUp part of fverT is set in this subroutine
131	C-- the kDown is still required
132	fVerT(1,1,kDown) = fVerT(1,1,kDown)
133	#endif
134
135	#ifdef ALLOW_DIAGNOSTICS
136	C-- Set diagnostic suffix for the current tracer
137	IF ( useDiagnostics ) THEN
138	diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )
139	ENDIF
140	#endif
141
142	advFac = 0. _d 0
143	IF (calcAdvection) advFac = 1. _d 0
144	rAdvFac = rkSign*advFac
145	IF (implicitAdvection) rAdvFac = 0. _d 0
146
147	DO j=1-OLy,sNy+OLy
148	DO i=1-OLx,sNx+OLx
149	fZon(i,j) = 0. _d 0
150	fMer(i,j) = 0. _d 0
151	fVerT(i,j,kUp) = 0. _d 0
152	df(i,j) = 0. _d 0
153	df4(i,j) = 0. _d 0
154	ENDDO
155	ENDDO
156
157	C-- Make local copy of tracer array
158	DO j=1-OLy,sNy+OLy
159	DO i=1-OLx,sNx+OLx
160	localT(i,j)=Tracer(i,j,k,bi,bj)
161	ENDDO
162	ENDDO
163
164	C-- Unless we have already calculated the advection terms we initialize
165	C the tendency to zero.
166	C <== now done earlier at the beginning of thermodynamics.
167	c IF (calcAdvection) THEN
168	c DO j=1-Oly,sNy+Oly
169	c DO i=1-Olx,sNx+Olx
170	c gTracer(i,j,k,bi,bj)=0. _d 0
171	c ENDDO
172	c ENDDO
173	c ENDIF
174
175	C-- Pre-calculate del^2 T if bi-harmonic coefficient is non-zero
176	IF (diffK4 .NE. 0.) THEN
177	CALL GAD_GRAD_X(bi,bj,k,xA,localT,fZon,myThid)
178	CALL GAD_GRAD_Y(bi,bj,k,yA,localT,fMer,myThid)
179	CALL GAD_DEL2(bi,bj,k,fZon,fMer,df4,myThid)
180	ENDIF
181
182	C-- Initialize net flux in X direction
183	DO j=1-Oly,sNy+Oly
184	DO i=1-Olx,sNx+Olx
185	fZon(i,j) = 0. _d 0
186	ENDDO
187	ENDDO
188
189	C- Advective flux in X
190	IF (calcAdvection) THEN
191	IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
192	CALL GAD_C2_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
193	ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
194	& .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
195	CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme,
196	I dTtracerLev(k), uTrans, uVel, localT,
197	O af, myThid )
198	ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
199	CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, dTtracerLev(k),
200	I uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), localT,
201	O af, myThid )
202	ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
203	CALL GAD_U3_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
204	ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
205	CALL GAD_C4_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
206	ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
207	CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k),
208	I uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), localT,
209	O af, myThid )
210	ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
211	IF ( inAdMode ) THEN
212	cph This block is to trick the adjoint:
213	cph IF inAdExact=.FALSE., we want to use DST3
214	cph with limiters in forward, but without limiters in reverse.
215	CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k),
216	I uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), localT,
217	O af, myThid )
218	ELSE
219	CALL GAD_DST3FL_ADV_X( bi,bj,k, dTtracerLev(k),
220	I uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), localT,
221	O af, myThid )
222	ENDIF
223	ELSE
224	STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'
225	ENDIF
226	DO j=1-Oly,sNy+Oly
227	DO i=1-Olx,sNx+Olx
228	fZon(i,j) = fZon(i,j) + af(i,j)
229	ENDDO
230	ENDDO
231	#ifdef ALLOW_DIAGNOSTICS
232	IF ( useDiagnostics ) THEN
233	diagName = 'ADVx'//diagSufx
234	CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
235	ENDIF
236	#endif
237	ENDIF
238
239	C- Diffusive flux in X
240	IF (diffKh.NE.0.) THEN
241	CALL GAD_DIFF_X(bi,bj,k,xA,diffKh,localT,df,myThid)
242	ELSE
243	DO j=1-Oly,sNy+Oly
244	DO i=1-Olx,sNx+Olx
245	df(i,j) = 0. _d 0
246	ENDDO
247	ENDDO
248	ENDIF
249
250	C- Add bi-harmonic diffusive flux in X
251	IF (diffK4 .NE. 0.) THEN
252	CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)
253	ENDIF
254
255	#ifdef ALLOW_GMREDI
256	C- GM/Redi flux in X
257	IF (useGMRedi) THEN
258	C note should update GMREDI_XTRANSPORT to use localT and set df aja
259	CALL GMREDI_XTRANSPORT(
260	I iMin,iMax,jMin,jMax,bi,bj,K,
261	I xA,Tracer,tracerIdentity,
262	U df,
263	I myThid)
264	ENDIF
265	#endif
266	DO j=1-Oly,sNy+Oly
267	DO i=1-Olx,sNx+Olx
268	fZon(i,j) = fZon(i,j) + df(i,j)
269	ENDDO
270	ENDDO
271
272	#ifdef ALLOW_DIAGNOSTICS
273	C- Diagnostics of Tracer flux in X dir (mainly Diffusive term),
274	C excluding advective terms:
275	IF ( useDiagnostics .AND.
276	& (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. useGMRedi) ) THEN
277	diagName = 'DIFx'//diagSufx
278	CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
279	ENDIF
280	#endif
281
282	C-- Initialize net flux in Y direction
283	DO j=1-Oly,sNy+Oly
284	DO i=1-Olx,sNx+Olx
285	fMer(i,j) = 0. _d 0
286	ENDDO
287	ENDDO
288
289	C- Advective flux in Y
290	IF (calcAdvection) THEN
291	IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
292	CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
293	ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
294	& .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
295	CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme,
296	I dTtracerLev(k), vTrans, vVel, localT,
297	O af, myThid )
298	ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
299	CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, dTtracerLev(k),
300	I vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), localT,
301	O af, myThid )
302	ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
303	CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
304	ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
305	CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
306	ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
307	CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),
308	I vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), localT,
309	O af, myThid )
310	ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
311	IF ( inAdMode ) THEN
312	cph This block is to trick the adjoint:
313	cph IF inAdExact=.FALSE., we want to use DST3
314	cph with limiters in forward, but without limiters in reverse.
315	CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),
316	I vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), localT,
317	O af, myThid )
318	ELSE
319	CALL GAD_DST3FL_ADV_Y( bi,bj,k, dTtracerLev(k),
320	I vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), localT,
321	O af, myThid )
322	ENDIF
323	ELSE
324	STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'
325	ENDIF
326	DO j=1-Oly,sNy+Oly
327	DO i=1-Olx,sNx+Olx
328	fMer(i,j) = fMer(i,j) + af(i,j)
329	ENDDO
330	ENDDO
331	#ifdef ALLOW_DIAGNOSTICS
332	IF ( useDiagnostics ) THEN
333	diagName = 'ADVy'//diagSufx
334	CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
335	ENDIF
336	#endif
337	ENDIF
338
339	C- Diffusive flux in Y
340	IF (diffKh.NE.0.) THEN
341	CALL GAD_DIFF_Y(bi,bj,k,yA,diffKh,localT,df,myThid)
342	ELSE
343	DO j=1-Oly,sNy+Oly
344	DO i=1-Olx,sNx+Olx
345	df(i,j) = 0. _d 0
346	ENDDO
347	ENDDO
348	ENDIF
349
350	C- Add bi-harmonic flux in Y
351	IF (diffK4 .NE. 0.) THEN
352	CALL GAD_BIHARM_Y(bi,bj,k,yA,df4,diffK4,df,myThid)
353	ENDIF
354
355	#ifdef ALLOW_GMREDI
356	C- GM/Redi flux in Y
357	IF (useGMRedi) THEN
358	C note should update GMREDI_YTRANSPORT to use localT and set df aja
359	CALL GMREDI_YTRANSPORT(
360	I iMin,iMax,jMin,jMax,bi,bj,K,
361	I yA,Tracer,tracerIdentity,
362	U df,
363	I myThid)
364	ENDIF
365	#endif
366	DO j=1-Oly,sNy+Oly
367	DO i=1-Olx,sNx+Olx
368	fMer(i,j) = fMer(i,j) + df(i,j)
369	ENDDO
370	ENDDO
371
372	#ifdef ALLOW_DIAGNOSTICS
373	C- Diagnostics of Tracer flux in Y dir (mainly Diffusive terms),
374	C excluding advective terms:
375	IF ( useDiagnostics .AND.
376	& (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. useGMRedi) ) THEN
377	diagName = 'DIFy'//diagSufx
378	CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
379	ENDIF
380	#endif
381
382	C-- Compute vertical flux fVerT(kUp) at interface k (between k-1 & k):
383	C- Advective flux in R
384	#ifdef ALLOW_AIM
385	C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
386	IF (calcAdvection .AND. .NOT.implicitAdvection .AND. K.GE.2 .AND.
387	& (.NOT.useAIM .OR.tracerIdentity.NE.GAD_SALINITY .OR.K.LT.Nr)
388	& ) THEN
389	#else
390	IF (calcAdvection .AND. .NOT.implicitAdvection .AND. K.GE.2) THEN
391	#endif
392	C- Compute vertical advective flux in the interior:
393	IF (vertAdvecScheme.EQ.ENUM_CENTERED_2ND) THEN
394	CALL GAD_C2_ADV_R(bi,bj,k,rTrans,Tracer,af,myThid)
395	ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
396	& .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
397	CALL GAD_DST2U1_ADV_R( bi,bj,k, vertAdvecScheme,
398	I dTtracerLev(k),rTrans,wVel,Tracer(1-Olx,1-Oly,1,bi,bj),
399	O af, myThid )
400	ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
401	CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k,
402	I dTtracerLev(k),rTrans,wVel,Tracer(1-Olx,1-Oly,1,bi,bj),
403	O af, myThid )
404	ELSEIF (vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN
405	CALL GAD_U3_ADV_R(bi,bj,k,rTrans,Tracer,af,myThid)
406	ELSEIF (vertAdvecScheme.EQ.ENUM_CENTERED_4TH) THEN
407	CALL GAD_C4_ADV_R(bi,bj,k,rTrans,Tracer,af,myThid)
408	ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
409	CALL GAD_DST3_ADV_R( bi,bj,k,
410	I dTtracerLev(k),rTrans,wVel,Tracer(1-Olx,1-Oly,1,bi,bj),
411	O af, myThid )
412	ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
413	cph This block is to trick the adjoint:
414	cph IF inAdExact=.FALSE., we want to use DST3
415	cph with limiters in forward, but without limiters in reverse.
416	IF ( inAdMode ) THEN
417	CALL GAD_DST3_ADV_R( bi,bj,k,
418	I dTtracerLev(k),rTrans,wVel,Tracer(1-Olx,1-Oly,1,bi,bj),
419	O af, myThid )
420	ELSE
421	CALL GAD_DST3FL_ADV_R( bi,bj,k,
422	I dTtracerLev(k),rTrans,wVel,Tracer(1-Olx,1-Oly,1,bi,bj),
423	O af, myThid )
424	ENDIF
425	ELSE
426	STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'
427	ENDIF
428	C- add the advective flux to fVerT
429	DO j=1-Oly,sNy+Oly
430	DO i=1-Olx,sNx+Olx
431	fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)
432	ENDDO
433	ENDDO
434	#ifdef ALLOW_DIAGNOSTICS
435	IF ( useDiagnostics ) THEN
436	diagName = 'ADVr'//diagSufx
437	CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
438	C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
439	C does it only if k=1 (never the case here)
440	IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)
441	ENDIF
442	#endif
443	ENDIF
444
445	C- Diffusive flux in R
446	C Note: For K=1 then KM1=1 and this gives a dT/dr = 0 upper
447	C boundary condition.
448	IF (implicitDiffusion) THEN
449	DO j=1-Oly,sNy+Oly
450	DO i=1-Olx,sNx+Olx
451	df(i,j) = 0. _d 0
452	ENDDO
453	ENDDO
454	ELSE
455	CALL GAD_DIFF_R(bi,bj,k,KappaR,Tracer,df,myThid)
456	ENDIF
457
458	#ifdef ALLOW_GMREDI
459	C- GM/Redi flux in R
460	IF (useGMRedi) THEN
461	C note should update GMREDI_RTRANSPORT to set df aja
462	CALL GMREDI_RTRANSPORT(
463	I iMin,iMax,jMin,jMax,bi,bj,K,
464	I Tracer,tracerIdentity,
465	U df,
466	I myThid)
467	ENDIF
468	#endif
469
470	DO j=1-Oly,sNy+Oly
471	DO i=1-Olx,sNx+Olx
472	fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)
473	ENDDO
474	ENDDO
475
476	#ifdef ALLOW_DIAGNOSTICS
477	C- Diagnostics of Tracer flux in R dir (mainly Diffusive terms),
478	C Explicit terms only & excluding advective terms:
479	IF ( useDiagnostics .AND.
480	& (.NOT.implicitDiffusion .OR. useGMRedi) ) THEN
481	diagName = 'DFrE'//diagSufx
482	CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
483	ENDIF
484	#endif
485
486	#ifdef ALLOW_KPP
487	C- Set non local KPP transport term (ghat):
488	IF ( useKPP .AND. k.GE.2 ) THEN
489	DO j=1-Oly,sNy+Oly
490	DO i=1-Olx,sNx+Olx
491	df(i,j) = 0. _d 0
492	ENDDO
493	ENDDO
494	IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
495	CALL KPP_TRANSPORT_T(
496	I iMin,iMax,jMin,jMax,bi,bj,k,km1,
497	O df )
498	ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
499	CALL KPP_TRANSPORT_S(
500	I iMin,iMax,jMin,jMax,bi,bj,k,km1,
501	O df )
502	#ifdef ALLOW_PTRACERS
503	ELSEIF (tracerIdentity .GE. GAD_TR1) THEN
504	CALL KPP_TRANSPORT_PTR(
505	I iMin,iMax,jMin,jMax,bi,bj,k,km1,
506	I tracerIdentity-GAD_TR1+1,
507	O df )
508	#endif
509	ELSE
510	PRINT*,'invalid tracer indentity: ', tracerIdentity
511	STOP 'GAD_CALC_RHS: Ooops'
512	ENDIF
513	DO j=1-Oly,sNy+Oly
514	DO i=1-Olx,sNx+Olx
515	fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)
516	ENDDO
517	ENDDO
518	ENDIF
519	#endif
520
521	C-- Divergence of fluxes
522	DO j=1-Oly,sNy+Oly-1
523	DO i=1-Olx,sNx+Olx-1
524	gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)
525	& -_recip_hFacC(i,j,k,bi,bj)recip_drF(k)recip_rA(i,j,bi,bj)
526	& *( (fZon(i+1,j)-fZon(i,j))
527	& +(fMer(i,j+1)-fMer(i,j))
528	& +(fVerT(i,j,kDown)-fVerT(i,j,kUp))*rkSign
529	& -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))
530	& +(vTrans(i,j+1)-vTrans(i,j))
531	& +(rTransKp1(i,j)-rTrans(i,j))*rAdvFac
532	& )*advFac
533	& )
534	ENDDO
535	ENDDO
536
537	#ifdef ALLOW_DEBUG
538	IF ( debugLevel .GE. debLevB
539	& .AND. tracerIdentity.EQ.GAD_TEMPERATURE
540	& .AND. k.EQ.2 .AND. myIter.EQ.1+nIter0
541	& .AND. nPx.EQ.1 .AND. nPy.EQ.1
542	& .AND. useCubedSphereExchange ) THEN
543	CALL DEBUG_CS_CORNER_UV( ' fZon,fMer from GAD_CALC_RHS',
544	& fZon,fMer, k, standardMessageUnit,bi,bj,myThid )
545	ENDIF
546	#endif /* ALLOW_DEBUG */
547
548	RETURN
549	END