model/src/calc_div_ghat.F

1dbaea09ee Chri*0001 #include "CPP_OPTIONS.h"
6d54cf9ca1 Ed H*0002
9366854e02 Chri*0003 CBOP
                0004 C     !ROUTINE: CALC_DIV_GHAT
                0005 C     !INTERFACE:
                0006       SUBROUTINE CALC_DIV_GHAT(
97c1a61dc2 Jean*0007      I                bi,bj,k,
                0008      U                cg2d_b, cg3d_b,
                0009      I                myThid)
9366854e02 Chri*0010 C     !DESCRIPTION: \bv
                0011 C     *==========================================================*
cb7fa97db9 Jean*0012 C     | S/R CALC_DIV_GHAT
                0013 C     | o Form the right hand-side of the surface pressure eqn.
9366854e02 Chri*0014 C     *==========================================================*
                0015 C     | Right hand side of pressure equation is divergence
                0016 C     | of veclocity tendency (GHAT) term along with a relaxation
                0017 C     | term equal to the barotropic flow field divergence.
                0018 C     *==========================================================*
                0019 C     \ev
8ae6a40c09 Alis*0020
9366854e02 Chri*0021 C     !USES:
8ae6a40c09 Alis*0022       IMPLICIT NONE
                0023 C     == Global variables ==
                0024 #include "SIZE.h"
81bc00c2f0 Chri*0025 #include "EEPARAMS.h"
8ae6a40c09 Alis*0026 #include "PARAMS.h"
                0027 #include "GRID.h"
97c1a61dc2 Jean*0028 #include "DYNVARS.h"
229ac9feb6 Jean*0029 #ifdef ALLOW_ADDFLUID
                0030 # include "FFIELDS.h"
                0031 #endif
8ae6a40c09 Alis*0032
9366854e02 Chri*0033 C     !INPUT/OUTPUT PARAMETERS:
8ae6a40c09 Alis*0034 C     == Routine arguments ==
cb7fa97db9 Jean*0035 C     bi, bj  :: tile indices
                0036 C     k       :: Index of layer.
                0037 C     cg2d_b  :: Conjugate Gradient 2-D solver : Right-hand side vector
97c1a61dc2 Jean*0038 C     cg3d_b  :: Conjugate Gradient 3-D solver : Right-hand side vector
cb7fa97db9 Jean*0039 C     myThid  :: Instance number for this call of CALC_DIV_GHAT
97c1a61dc2 Jean*0040       INTEGER bi,bj
cb7fa97db9 Jean*0041       INTEGER k
71c6a09c16 Jean*0042       _RL cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
08a6f65fd0 Jean*0043 #ifdef ALLOW_NONHYDROSTATIC
97c1a61dc2 Jean*0044       _RL cg3d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
08a6f65fd0 Jean*0045 #else
                0046       _RL cg3d_b(1)
                0047 #endif
8ae6a40c09 Alis*0048       INTEGER myThid
                0049
9366854e02 Chri*0050 C     !LOCAL VARIABLES:
8ae6a40c09 Alis*0051 C     == Local variables ==
97c1a61dc2 Jean*0052 C     i,j    :: Loop counters
                0053 C     xA, yA :: Cell vertical face areas
                0054 C     pf     :: Intermediate array for building RHS source term.
8ae6a40c09 Alis*0055       INTEGER i,j
97c1a61dc2 Jean*0056       _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0057       _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
8ae6a40c09 Alis*0058       _RL pf (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79b5d5775c Jean*0059 #ifdef ALLOW_ADDFLUID
                0060       _RL unitsFac
                0061 #endif
9366854e02 Chri*0062 CEOP
8ae6a40c09 Alis*0063
97c1a61dc2 Jean*0064 C     Calculate vertical face areas
                0065       DO j=1,sNy+1
                0066         DO i=1,sNx+1
                0067           xA(i,j) = _dyG(i,j,bi,bj)*deepFacC(k)
                0068      &              *drF(k)*_hFacW(i,j,k,bi,bj)*rhoFacC(k)
                0069           yA(i,j) = _dxG(i,j,bi,bj)*deepFacC(k)
                0070      &              *drF(k)*_hFacS(i,j,k,bi,bj)*rhoFacC(k)
                0071         ENDDO
                0072       ENDDO
                0073
8ae6a40c09 Alis*0074 C--   Pressure equation source term
cb7fa97db9 Jean*0075 C     Term is the vertical integral of the divergence of the
8ae6a40c09 Alis*0076 C     time tendency terms along with a relaxation term that
                0077 C     pulls div(U) + dh/dt back toward zero.
                0078
9fdcce8056 Jean*0079 C     Implicit treatment of the Barotropic Flow Divergence
79b5d5775c Jean*0080       DO j=1,sNy
                0081        DO i=1,sNx+1
                0082         pf(i,j) = implicDiv2Dflow*xA(i,j)*gU(i,j,k,bi,bj) / deltaTMom
                0083        ENDDO
                0084       ENDDO
8ae6a40c09 Alis*0085       DO j=1,sNy
                0086        DO i=1,sNx
88830be691 Alis*0087         cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj) +
8ae6a40c09 Alis*0088      &   pf(i+1,j)-pf(i,j)
                0089        ENDDO
                0090       ENDDO
                0091
88830be691 Alis*0092 #ifdef ALLOW_NONHYDROSTATIC
cb7fa97db9 Jean*0093       IF (use3Dsolver) THEN
88830be691 Alis*0094        DO j=1,sNy
                0095         DO i=1,sNx
a7cb4644df Jean*0096          cg3d_b(i,j,k,bi,bj) = ( pf(i+1,j)-pf(i,j) )
88830be691 Alis*0097         ENDDO
                0098        ENDDO
                0099       ENDIF
                0100 #endif
                0101
9fdcce8056 Jean*0102 C     Implicit treatment of the Barotropic Flow Divergence
79b5d5775c Jean*0103       DO j=1,sNy+1
                0104        DO i=1,sNx
                0105         pf(i,j) = implicDiv2Dflow*yA(i,j)*gV(i,j,k,bi,bj) / deltaTMom
                0106        ENDDO
                0107       ENDDO
8ae6a40c09 Alis*0108       DO j=1,sNy
                0109        DO i=1,sNx
88830be691 Alis*0110         cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj) +
8ae6a40c09 Alis*0111      &   pf(i,j+1)-pf(i,j)
                0112        ENDDO
                0113       ENDDO
                0114
88830be691 Alis*0115 #ifdef ALLOW_NONHYDROSTATIC
cb7fa97db9 Jean*0116       IF (use3Dsolver) THEN
88830be691 Alis*0117        DO j=1,sNy
                0118         DO i=1,sNx
a7cb4644df Jean*0119          cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
                0120      &                       + ( pf(i,j+1)-pf(i,j) )
88830be691 Alis*0121         ENDDO
                0122        ENDDO
                0123       ENDIF
                0124 #endif
c712b7e80b Chri*0125
d2a11ab670 Jean*0126 #ifdef ALLOW_ADDFLUID
79b5d5775c Jean*0127       unitsFac = mass2rUnit/deltaTMom
                0128       IF ( exactConserv ) unitsFac = unitsFac*implicDiv2Dflow
d2a11ab670 Jean*0129       IF ( selectAddFluid.GE.1 ) THEN
                0130         DO j=1,sNy
                0131          DO i=1,sNx
                0132           cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
79b5d5775c Jean*0133      &         - addMass(i,j,k,bi,bj)*unitsFac
d2a11ab670 Jean*0134          ENDDO
                0135         ENDDO
                0136 #ifdef ALLOW_NONHYDROSTATIC
                0137        IF (use3Dsolver) THEN
                0138         DO j=1,sNy
                0139          DO i=1,sNx
                0140           cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
79b5d5775c Jean*0141      &         - addMass(i,j,k,bi,bj)*unitsFac
d2a11ab670 Jean*0142          ENDDO
                0143         ENDDO
                0144        ENDIF
                0145 #endif
                0146       ENDIF
                0147 #endif /* ALLOW_ADDFLUID */
                0148
8ae6a40c09 Alis*0149       RETURN
                0150       END