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