pkg/cheapaml/cheapaml_calc_rhs.F

cf5b5345a0 Jean*0001 #include "GAD_OPTIONS.h"
                0002
                0003 CBOP
2b5bd8961b Jean*0004 C !ROUTINE: CHEAPAML_CALC_RHS
cf5b5345a0 Jean*0005
                0006 C !INTERFACE: ==========================================================
2b5bd8961b Jean*0007       SUBROUTINE CHEAPAML_CALC_RHS(
58426debb4 Jean*0008      I           bi,bj, iMin,iMax,jMin,jMax,
                0009      I           uTrans, vTrans,
364bacbe7c Jean*0010      I           uVel, vVel,
76b059cef2 Jean*0011      I           diffKh, Tracer,
                0012      I           deltaTtracer, zu,
                0013      I           useFluxLimit, cheapamlXperiodic, cheapamlYperiodic,
                0014      O           wVel,
                0015      U           gTracer,
cf5b5345a0 Jean*0016      I           myTime, myIter, myThid )
                0017
                0018 C !DESCRIPTION:
ced0783fba Jean*0019 C Calculates the tendency of a tracer due to advection and diffusion.
                0020 C Because horizontal velocity field is potential divergent, it is
                0021 C necessary to compute the vertical velocity and to compute the
364bacbe7c Jean*0022 C vertical flux divergence.
ced0783fba Jean*0023 C The fluxes in each direction are computed independently and then
364bacbe7c Jean*0024 C the tendency is set to the divergence of these fluxes.
ced0783fba Jean*0025 C In Cheapaml, it is always assumed the boundaries are open, and
364bacbe7c Jean*0026 C a simple open boundary implementation is used, whereby if the
ced0783fba Jean*0027 C transport is outward directed, upwind weighting is used
                0028 C for the advective flux and the diffusive flux is shut off.
                0029 C If the transport is inward directed, the advective flux is computed
                0030 C using the Tr file, as is the diffusive flux.
cf5b5345a0 Jean*0031 C
                0032 C The tendency is the divergence of the fluxes:
                0033 C \begin{equation*}
                0034 C G_\theta = G_\theta + \nabla \cdot {\bf F}
                0035 C \end{equation*}
                0036 C
                0037 C The tendency is assumed to contain data on entry.
                0038
                0039 C !USES: ===============================================================
                0040       IMPLICIT NONE
                0041 #include "SIZE.h"
                0042 #include "EEPARAMS.h"
                0043 #include "PARAMS.h"
                0044 #include "GRID.h"
                0045 #include "GAD.h"
                0046
                0047 C !INPUT PARAMETERS: ===================================================
                0048 C bi,bj            :: tile indices
                0049 C iMin,iMax        :: loop range for called routines
                0050 C jMin,jMax        :: loop range for called routines
                0051 C uTrans,vTrans    :: 2-D arrays of volume transports at U,V points
                0052 C uVel,vVel,       :: 2 components of the velcity field (2-D array)
                0053 C diffKh           :: horizontal diffusion coefficient
                0054 C Tracer           :: tracer field
364bacbe7c Jean*0055 C deltaTtracer     :: atmospheric tracer time step
76b059cef2 Jean*0056 C zu               ::
                0057 C useFluxLimit     ::
                0058 C cheapamlXperiodic :: cheapaml domain is periodic in X dir
                0059 C cheapamlYperiodic :: cheapaml domain is periodic in Y dir
cf5b5345a0 Jean*0060 C myTime           :: current time
                0061 C myIter           :: iteration number
                0062 C myThid           :: thread number
                0063       INTEGER bi,bj,iMin,iMax,jMin,jMax
                0064       _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0065       _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0066       _RL uVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0067       _RL vVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
76b059cef2 Jean*0068       _RL diffKh
cf5b5345a0 Jean*0069       _RL Tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
76b059cef2 Jean*0070       _RL deltaTtracer, zu
364bacbe7c Jean*0071       LOGICAL useFluxLimit
76b059cef2 Jean*0072       LOGICAL cheapamlXperiodic, cheapamlYperiodic
cf5b5345a0 Jean*0073       _RL     myTime
                0074       INTEGER myIter, myThid
                0075
                0076 C !OUTPUT PARAMETERS: ==================================================
76b059cef2 Jean*0077 C wVel             ::
ced0783fba Jean*0078 C gTracer          :: tendency array
76b059cef2 Jean*0079       _RL wVel   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
cf5b5345a0 Jean*0080       _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0081
                0082 C !LOCAL VARIABLES: ====================================================
                0083 C i,j              :: loop indices
ced0783fba Jean*0084 C fZon             :: zonal fluxes
                0085 C fMer             :: meridional fluxes
364bacbe7c Jean*0086 C fVer             :: vertical fluxes
ced0783fba Jean*0087 C af               :: advective fluxes
cf5b5345a0 Jean*0088 C df               :: diffusive flux
                0089 C localT           :: local copy of tracer field
58426debb4 Jean*0090 C horizDiv         :: divergence of the horizontal wind
364bacbe7c Jean*0091 C maskLocW         :: local copy of West Face Land Mask
                0092 C maskLocS         :: local copy of South Face Land Mask
                0093
58426debb4 Jean*0094       INTEGER i,j, iG,jG, ii, jj
cf5b5345a0 Jean*0095       _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0096       _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0097       _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0098       _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
58426debb4 Jean*0099       _RL horizDiv(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
364bacbe7c Jean*0100       _RS maskLocW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0101       _RS maskLocS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
cf5b5345a0 Jean*0102 CEOP
                0103
58426debb4 Jean*0104 C--   Initialize fluxes
cf5b5345a0 Jean*0105       DO j=1-OLy,sNy+OLy
ced0783fba Jean*0106         DO i=1-OLx,sNx+OLx
58426debb4 Jean*0107           df(i,j)   = 0. _d 0
                0108           fZon(i,j) = 0. _d 0
                0109           fMer(i,j) = 0. _d 0
ced0783fba Jean*0110         ENDDO
cf5b5345a0 Jean*0111       ENDDO
364bacbe7c Jean*0112
cf5b5345a0 Jean*0113 C--   Make local copy of tracer array
364bacbe7c Jean*0114 C--   and compute w
cf5b5345a0 Jean*0115       DO j=1-OLy,sNy+OLy
ced0783fba Jean*0116         DO i=1-OLx,sNx+OLx
                0117           localT(i,j)=tracer(i,j,bi,bj)
4fa4901be6 Nico*0118         ENDDO
                0119       ENDDO
                0120       DO j=1-OLy,sNy+OLy-1
                0121         DO i=1-OLx,sNx+OLx-1
58426debb4 Jean*0122           horizDiv(i,j) = (
                0123      &       ( uTrans(i+1,j) - uTrans(i,j) )
                0124      &     + ( vTrans(i,j+1) - vTrans(i,j) )
                0125      &                    )*recip_rA(i,j,bi,bj)
                0126           wVel(i,j,bi,bj) = -horizDiv(i,j)
ced0783fba Jean*0127         ENDDO
cf5b5345a0 Jean*0128       ENDDO
                0129
364bacbe7c Jean*0130 C     prepare boundary tracer values
76b059cef2 Jean*0131       IF(.NOT.cheapamlXperiodic)THEN
                0132        DO j=1,sNy
                0133         DO i=1,sNx
364bacbe7c Jean*0134          iG=myXGlobalLo-1+(bi-1)*sNx+i
                0135          IF (iG.EQ.2) THEN
                0136            IF (uVel(i,j,bi,bj).LT.0. _d 0) THEN
                0137              DO ii=1-OLx,1
                0138               localT(ii,j)=localT(i,j)
                0139              ENDDO
                0140            ENDIF
                0141          ELSEIF (iG.EQ.Nx-1) THEN
                0142            IF (uVel(i+1,j,bi,bj).GT.0. _d 0) THEN
                0143              DO ii=sNx,sNx+OLx
                0144               localT(ii,j)=localT(i,j)
                0145              ENDDO
                0146            ENDIF
                0147          ENDIF
76b059cef2 Jean*0148         ENDDO
364bacbe7c Jean*0149        ENDDO
76b059cef2 Jean*0150       ENDIF
ced0783fba Jean*0151 C     -    Advective flux in X
364bacbe7c Jean*0152       IF (useFluxLimit) THEN
                0153 C     make local copy of west land mask
ced0783fba Jean*0154         DO j=1-OLy,sNy+OLy
                0155           DO i=1-OLx,sNx+OLx
364bacbe7c Jean*0156 c            maskLocW(i,j)=maskW(i,j,1,bi,bj)
                0157 C     change to allow advection on land
                0158              maskLocW(i,j)=1. _d 0
ced0783fba Jean*0159           ENDDO
                0160         ENDDO
364bacbe7c Jean*0161         CALL GAD_DST3FL_ADV_X(
ced0783fba Jean*0162      I     bi,bj,1,.TRUE., deltaTtracer,
58426debb4 Jean*0163      I     uTrans, uVel(1-OLx,1-OLy,bi,bj),
ced0783fba Jean*0164      I     maskLocW, localT,
58426debb4 Jean*0165      O     fZon,
ced0783fba Jean*0166      I     myThid )
364bacbe7c Jean*0167       ELSE
58426debb4 Jean*0168         CALL GAD_C2_ADV_X( bi,bj,1,uTrans,localT,fZon,myThid )
364bacbe7c Jean*0169       ENDIF
cf5b5345a0 Jean*0170
ced0783fba Jean*0171 C     -    Diffusive flux in X
58426debb4 Jean*0172       IF ( diffKh.NE.zeroRL ) THEN
                0173         CALL GAD_DIFF_X( bi,bj,1, dyG(1-OLx,1-OLy,bi,bj),
                0174      I                   diffKh, localT,
                0175      O                   df, myThid )
364bacbe7c Jean*0176         DO j=1-OLy,sNy+OLy
                0177           DO i=1-OLx,sNx+OLx
58426debb4 Jean*0178             fZon(i,j) = fZon(i,j) + df(i,j)
ced0783fba Jean*0179           ENDDO
cf5b5345a0 Jean*0180         ENDDO
                0181       ENDIF
364bacbe7c Jean*0182
                0183 C     repair boundary tracer values
76b059cef2 Jean*0184       IF(.NOT.cheapamlXperiodic)THEN
                0185        DO j=1,sNy
                0186         DO i=1,sNx
364bacbe7c Jean*0187          iG=myXGlobalLo-1+(bi-1)*sNx+i
                0188          IF (iG.EQ.2) THEN
                0189            IF (uVel(i,j,bi,bj).LT.0. _d 0) THEN
                0190              DO ii=1-OLx,1
                0191               localT(ii,j)=tracer(ii,j,bi,bj)
                0192              ENDDO
                0193            ENDIF
                0194           ELSEIF (iG.EQ.Nx-1) THEN
                0195            IF (uVel(i+1,j,bi,bj).GT.0. _d 0) THEN
                0196              DO ii=sNx,sNx+OLx
                0197               localT(ii,j)=tracer(ii,j,bi,bj)
                0198              ENDDO
                0199            ENDIF
                0200          ENDIF
76b059cef2 Jean*0201         ENDDO
364bacbe7c Jean*0202        ENDDO
76b059cef2 Jean*0203       ENDIF
ced0783fba Jean*0204
                0205 C     -    Advective flux in Y
364bacbe7c Jean*0206
                0207 C     prepare boundary tracer values
76b059cef2 Jean*0208       IF(.NOT.cheapamlYperiodic)THEN
                0209        DO j=1,sNy
                0210         jG = myYGlobalLo-1+(bj-1)*sNy+j
                0211         DO i=1,sNx
364bacbe7c Jean*0212          IF (jG.EQ.2) THEN
                0213            IF (vVel(i,j,bi,bj).LT.0. _d 0) THEN
                0214              DO jj=1-OLy,1
                0215               localT(i,jj)=localT(i,j)
                0216              ENDDO
                0217            ENDIF
                0218          ELSEIF (jG.EQ.Ny-1) THEN
                0219            IF (vVel(i,j+1,bi,bj).GT.0. _d 0) THEN
                0220              DO jj=sNy,sNy+OLy
                0221               localT(i,jj)=localT(i,j)
                0222              ENDDO
                0223            ENDIF
                0224          ENDIF
76b059cef2 Jean*0225         ENDDO
364bacbe7c Jean*0226        ENDDO
76b059cef2 Jean*0227       ENDIF
364bacbe7c Jean*0228
                0229       IF (useFluxLimit) THEN
                0230 C     make local copy of south land mask
ced0783fba Jean*0231         DO j=1-OLy,sNy+OLy
                0232           DO i=1-OLx,sNx+OLx
364bacbe7c Jean*0233 c            maskLocS(i,j)=maskS(i,j,1,bi,bj)
                0234 C Change to allow advection on land
                0235             maskLocS(i,j)= 1. _d 0
ced0783fba Jean*0236           ENDDO
                0237         ENDDO
364bacbe7c Jean*0238         CALL GAD_DST3FL_ADV_Y(
ced0783fba Jean*0239      I     bi,bj,1,.TRUE., deltaTtracer,
58426debb4 Jean*0240      I     vTrans, vVel(1-OLx,1-OLy,bi,bj),
ced0783fba Jean*0241      I     maskLocS, localT,
58426debb4 Jean*0242      O     fMer,
ced0783fba Jean*0243      I     myThid )
364bacbe7c Jean*0244       ELSE
58426debb4 Jean*0245         CALL GAD_C2_ADV_Y( bi,bj,1,vTrans,localT,fMer,myThid )
364bacbe7c Jean*0246       ENDIF
                0247
ced0783fba Jean*0248 C     -    Diffusive flux in Y
58426debb4 Jean*0249       IF ( diffKh.NE.zeroRL ) THEN
                0250         CALL GAD_DIFF_Y( bi,bj,1, dxG(1-OLx,1-OLy,bi,bj),
                0251      I                   diffKh, localT,
                0252      O                   df, myThid )
                0253         IF ( .NOT.cheapamlYperiodic ) THEN
                0254          DO j=1,sNy+1
                0255           jG = myYGlobalLo-1+(bj-1)*sNy+j
                0256           IF ( jG.EQ.1 .OR. jG.EQ.Ny ) THEN
                0257 C- Note: conditions above ire strange. Instead, should it be:
                0258 c         IF ( jG.EQ.1 .OR. jG.EQ.Ny+1 ) THEN
                0259            DO i=1-OLx,sNx+OLx
                0260             df(i,j) = 0. _d 0
                0261            ENDDO
                0262           ENDIF
                0263          ENDDO
                0264         ENDIF
364bacbe7c Jean*0265         DO j=1-OLy,sNy+OLy
                0266           DO i=1-OLx,sNx+OLx
58426debb4 Jean*0267             fMer(i,j) = fMer(i,j) + df(i,j)
ced0783fba Jean*0268           ENDDO
cf5b5345a0 Jean*0269         ENDDO
ced0783fba Jean*0270       ENDIF
364bacbe7c Jean*0271
cf5b5345a0 Jean*0272 C--   Divergence of fluxes
364bacbe7c Jean*0273       DO j=1-OLy,sNy+OLy-1
                0274         DO i=1-OLx,sNx+OLx-1
58426debb4 Jean*0275           gTracer(i,j,bi,bj) =
                0276      &       -( ( fZon(i+1,j) - fZon(i,j) )
                0277      &        + ( fMer(i,j+1) - fMer(i,j) )
                0278      &        )*recip_rA(i,j,bi,bj)
                0279      &       + horizDiv(i,j)*localT(i,j)
ced0783fba Jean*0280         ENDDO
cf5b5345a0 Jean*0281       ENDDO
364bacbe7c Jean*0282
cf5b5345a0 Jean*0283       RETURN
                0284       END