C $Header: /u/gcmpack/MITgcm/pkg/cheapaml/gad_2d_calc_rhs.F,v 1.9 2013/02/18 21:16:15 jmc Exp $
C $Name: $
#include "GAD_OPTIONS.h"
CBOP
C !ROUTINE: GAD_2d_CALC_RHS
C !INTERFACE: ==========================================================
SUBROUTINE GAD_2D_CALC_RHS(
I bi,bj,iMin,iMax,jMin,jMax,
I uTrans,vTrans,
I uVel, vVel,
I diffKh, Tracer,
I deltaTtracer, zu,
I useFluxLimit, cheapamlXperiodic, cheapamlYperiodic,
O wVel,
U gTracer,
I myTime, myIter, myThid )
C !DESCRIPTION:
C Calculates the tendency of a tracer due to advection and diffusion.
C Because horizontal velocity field is potential divergent, it is
C necessary to compute the vertical velocity and to compute the
C vertical flux divergence.
C The fluxes in each direction are computed independently and then
C the tendency is set to the divergence of these fluxes.
C In Cheapaml, it is always assumed the boundaries are open, and
C a simple open boundary implementation is used, whereby if the
C transport is outward directed, upwind weighting is used
C for the advective flux and the diffusive flux is shut off.
C If the transport is inward directed, the advective flux is computed
C using the Tr file, as is the diffusive flux.
C
C The tendency is the divergence of the fluxes:
C \begin{equation*}
C G_\theta = G_\theta + \nabla \cdot {\bf F}
C \end{equation*}
C
C The tendency 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"
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 uTrans,vTrans :: 2-D arrays of volume transports at U,V points
C uVel,vVel, :: 2 components of the velcity field (2-D array)
C diffKh :: horizontal diffusion coefficient
C Tracer :: tracer field
C deltaTtracer :: atmospheric tracer time step
C zu ::
C useFluxLimit ::
C cheapamlXperiodic :: cheapaml domain is periodic in X dir
C cheapamlYperiodic :: cheapaml domain is periodic in Y dir
C myTime :: current time
C myIter :: iteration number
C myThid :: thread number
INTEGER bi,bj,iMin,iMax,jMin,jMax
_RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL uVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
_RL vVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
_RL diffKh
_RL Tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
_RL deltaTtracer, zu
LOGICAL useFluxLimit
LOGICAL cheapamlXperiodic, cheapamlYperiodic
_RL myTime
INTEGER myIter, myThid
C !OUTPUT PARAMETERS: ==================================================
C wVel ::
C gTracer :: tendency array
_RL wVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
_RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
C !LOCAL VARIABLES: ====================================================
C i,j :: loop indices
C fZon :: zonal fluxes
C fMer :: meridional fluxes
C fVer :: vertical fluxes
C af :: advective fluxes
C df :: diffusive flux
C localT :: local copy of tracer field
C localW :: local vertical velocity
C maskLocW :: local copy of West Face Land Mask
C maskLocS :: local copy of South Face Land Mask
INTEGER i,j,iG,jG, ii, jj
_RL fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL fVer (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 localW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RS maskLocW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RS maskLocS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
CEOP
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
df(i,j) = 0. _d 0
ENDDO
ENDDO
C-- Make local copy of tracer array
C-- and compute w
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
localT(i,j)=tracer(i,j,bi,bj)
ENDDO
ENDDO
DO j=1-OLy,sNy+OLy-1
DO i=1-OLx,sNx+OLx-1
localW(i,j)=-(recip_dxC(i,j,bi,bj)
& *(uVel(i+1,j,bi,bj)-uVel(i,j,bi,bj))
& + (vVel(i,j+1,bi,bj)-vVel(i,j,bi,bj))*recip_dyC(i,j,bi,bj)
& )*2. _d 0 *zu
wvel(i,j,bi,bj) = localW(i,j)/2. _d 0 /zu
ENDDO
ENDDO
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
af(i,j) = 0. _d 0
ENDDO
ENDDO
C prepare boundary tracer values
IF(.NOT.cheapamlXperiodic)THEN
DO j=1,sNy
DO i=1,sNx
iG=myXGlobalLo-1+(bi-1)*sNx+i
IF (iG.EQ.2) THEN
IF (uVel(i,j,bi,bj).LT.0. _d 0) THEN
DO ii=1-OLx,1
localT(ii,j)=localT(i,j)
ENDDO
ENDIF
ELSEIF (iG.EQ.Nx-1) THEN
IF (uVel(i+1,j,bi,bj).GT.0. _d 0) THEN
DO ii=sNx,sNx+OLx
localT(ii,j)=localT(i,j)
ENDDO
ENDIF
ENDIF
ENDDO
ENDDO
ENDIF
C - Advective flux in X
IF (useFluxLimit) THEN
C make local copy of west land mask
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
c maskLocW(i,j)=maskW(i,j,1,bi,bj)
C change to allow advection on land
maskLocW(i,j)=1. _d 0
ENDDO
ENDDO
CALL GAD_DST3FL_ADV_X(
I bi,bj,1,.TRUE., deltaTtracer,
I uTrans, uTrans,
I maskLocW, localT,
O af,
I myThid )
ELSE
CALL GAD_C2_2D_ADV_X(bi,bj,uVel,localT,af,myThid)
ENDIF
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
fZon(i,j)=fZon(i,j)+af(i,j)
ENDDO
ENDDO
C - Diffusive flux in X
IF (diffKh.NE.0.) THEN
CALL GAD_DIFF_2D_X(bi,bj,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
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
fZon(i,j) = fZon(i,j) + df(i,j)
ENDDO
ENDDO
C repair boundary tracer values
IF(.NOT.cheapamlXperiodic)THEN
DO j=1,sNy
DO i=1,sNx
iG=myXGlobalLo-1+(bi-1)*sNx+i
IF (iG.EQ.2) THEN
IF (uVel(i,j,bi,bj).LT.0. _d 0) THEN
DO ii=1-OLx,1
localT(ii,j)=tracer(ii,j,bi,bj)
ENDDO
ENDIF
ELSEIF (iG.EQ.Nx-1) THEN
IF (uVel(i+1,j,bi,bj).GT.0. _d 0) THEN
DO ii=sNx,sNx+OLx
localT(ii,j)=tracer(ii,j,bi,bj)
ENDDO
ENDIF
ENDIF
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
af(i,j) = 0. _d 0
df(i,j) = 0. _d 0
ENDDO
ENDDO
C - Advective flux in Y
C prepare boundary tracer values
IF(.NOT.cheapamlYperiodic)THEN
DO j=1,sNy
jG = myYGlobalLo-1+(bj-1)*sNy+j
DO i=1,sNx
IF (jG.EQ.2) THEN
IF (vVel(i,j,bi,bj).LT.0. _d 0) THEN
DO jj=1-OLy,1
localT(i,jj)=localT(i,j)
ENDDO
ENDIF
ELSEIF (jG.EQ.Ny-1) THEN
IF (vVel(i,j+1,bi,bj).GT.0. _d 0) THEN
DO jj=sNy,sNy+OLy
localT(i,jj)=localT(i,j)
ENDDO
ENDIF
ENDIF
ENDDO
ENDDO
ENDIF
IF (useFluxLimit) THEN
C make local copy of south land mask
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
c maskLocS(i,j)=maskS(i,j,1,bi,bj)
C Change to allow advection on land
maskLocS(i,j)= 1. _d 0
ENDDO
ENDDO
CALL GAD_DST3FL_ADV_Y(
I bi,bj,1,.TRUE., deltaTtracer,
I vTrans, vTrans,
I maskLocS, localT,
O af,
I myThid )
ELSE
CALL GAD_C2_2D_ADV_Y(bi,bj,vVel,localT,af,myThid)
ENDIF
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
fMer(i,j) = fMer(i,j) + af(i,j)
ENDDO
ENDDO
C - Diffusive flux in Y
IF (diffKh.NE.0.) THEN
CALL GAD_DIFF_2D_Y(bi,bj,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
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
fMer(i,j) = fMer(i,j) + df(i,j)
ENDDO
ENDDO
C-- Initialize net flux in R direction
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
fVer(i,j) = 0. _d 0
af(i,j) = 0. _d 0
ENDDO
ENDDO
C - Advective flux in R
C at top of the cell
IF (useFluxLimit) THEN
CALL CHEAPAML_GAD_DST3FL_ADV_R(
I bi,bj,deltaTtracer,
I localW,
I localT,
O af,
I myThid )
ELSE
CALL CHEAPAML_GAD_C2_ADV_R(
I bi,bj,
I localW,
I localT,
O af,
I myThid )
ENDIF
DO j=1-OLy,sNy+OLy-1
DO i=1-OLx,sNx+OLx-1
fVer(i,j) = fVer(i,j) + af(i,j)
ENDDO
ENDDO
C no need to repair localT
C-- Divergence of fluxes
DO j=1-OLy,sNy+OLy-1
DO i=1-OLx,sNx+OLx-1
gTracer(i,j,bi,bj)=
& -(recip_dxC(i,j,bi,bj)
& *(fZon(i+1,j)-fZon(i,j))
& +recip_dyC(i,j,bi,bj)
& *(fMer(i,j+1)-fMer(i,j))
& +1. _d 0/2. _d 0/zu*fVer(i,j)
& )
ENDDO
ENDDO
RETURN
END