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)
9366854e02 Chri*0059 CEOP
8ae6a40c09 Alis*0060
97c1a61dc2 Jean*0061 C     Calculate vertical face areas
                0062       DO j=1,sNy+1
                0063         DO i=1,sNx+1
                0064           xA(i,j) = _dyG(i,j,bi,bj)*deepFacC(k)
                0065      &              *drF(k)*_hFacW(i,j,k,bi,bj)*rhoFacC(k)
                0066           yA(i,j) = _dxG(i,j,bi,bj)*deepFacC(k)
                0067      &              *drF(k)*_hFacS(i,j,k,bi,bj)*rhoFacC(k)
                0068         ENDDO
                0069       ENDDO
                0070
8ae6a40c09 Alis*0071 C--   Pressure equation source term
cb7fa97db9 Jean*0072 C     Term is the vertical integral of the divergence of the
8ae6a40c09 Alis*0073 C     time tendency terms along with a relaxation term that
                0074 C     pulls div(U) + dh/dt back toward zero.
                0075
9fdcce8056 Jean*0076       IF (implicDiv2Dflow.EQ.1.) THEN
678051767f Jean*0077 C     Fully Implicit treatment of the Barotropic Flow Divergence
                0078         DO j=1,sNy
                0079          DO i=1,sNx+1
229ac9feb6 Jean*0080           pf(i,j) = xA(i,j)*gU(i,j,k,bi,bj) / deltaTMom
678051767f Jean*0081          ENDDO
                0082         ENDDO
9fdcce8056 Jean*0083       ELSEIF (exactConserv) THEN
                0084 c     ELSEIF (nonlinFreeSurf.GT.0) THEN
                0085 C     Implicit treatment of the Barotropic Flow Divergence
                0086         DO j=1,sNy
                0087          DO i=1,sNx+1
cb7fa97db9 Jean*0088           pf(i,j) = implicDiv2Dflow
229ac9feb6 Jean*0089      &             *xA(i,j)*gU(i,j,k,bi,bj) / deltaTMom
9fdcce8056 Jean*0090          ENDDO
                0091         ENDDO
678051767f Jean*0092       ELSE
                0093 C     Explicit+Implicit part of the Barotropic Flow Divergence
                0094 C      => Filtering of uVel,vVel is necessary
9fdcce8056 Jean*0095 C-- Now the filter are applied in the_correction_step().
                0096 C   We have left this code here to indicate where the filters used to be
cb7fa97db9 Jean*0097 C   in the algorithm before JMC moved them to after the pressure solver.
e4b65e705c Jean*0098 c#ifdef ALLOW_ZONAL_FILT
                0099 c        IF (zonal_filt_lat.LT.90.) THEN
                0100 c          CALL ZONAL_FILTER(
                0101 c    U                     uVel( 1-OLx,1-OLy,k,bi,bj),
                0102 c    I                     hFacW(1-OLx,1-OLy,k,bi,bj),
                0103 c    I                     0, sNy+1, 1, bi, bj, 1, myThid )
                0104 c          CALL ZONAL_FILTER(
                0105 c    U                     vVel( 1-OLx,1-OLy,k,bi,bj),
                0106 c    I                     hFacS(1-OLx,1-OLy,k,bi,bj),
                0107 c    I                     0, sNy+1, 1, bi, bj, 2, myThid )
                0108 c        ENDIF
                0109 c#endif
678051767f Jean*0110         DO j=1,sNy
                0111          DO i=1,sNx+1
9f9dfaf749 Jean*0112           pf(i,j) = ( implicDiv2Dflow * gU(i,j,k,bi,bj)
97c1a61dc2 Jean*0113      &     + (1. _d 0-implicDiv2Dflow)* uVel(i,j,k,bi,bj)
229ac9feb6 Jean*0114      &               ) * xA(i,j) / deltaTMom
678051767f Jean*0115          ENDDO
                0116         ENDDO
                0117       ENDIF
8ae6a40c09 Alis*0118       DO j=1,sNy
                0119        DO i=1,sNx
88830be691 Alis*0120         cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj) +
8ae6a40c09 Alis*0121      &   pf(i+1,j)-pf(i,j)
                0122        ENDDO
                0123       ENDDO
                0124
88830be691 Alis*0125 #ifdef ALLOW_NONHYDROSTATIC
cb7fa97db9 Jean*0126       IF (use3Dsolver) THEN
88830be691 Alis*0127        DO j=1,sNy
                0128         DO i=1,sNx
a7cb4644df Jean*0129          cg3d_b(i,j,k,bi,bj) = ( pf(i+1,j)-pf(i,j) )
88830be691 Alis*0130         ENDDO
                0131        ENDDO
                0132       ENDIF
                0133 #endif
                0134
9fdcce8056 Jean*0135       IF (implicDiv2Dflow.EQ.1.) THEN
678051767f Jean*0136 C     Fully Implicit treatment of the Barotropic Flow Divergence
                0137         DO j=1,sNy+1
                0138          DO i=1,sNx
9f9dfaf749 Jean*0139           pf(i,j) = yA(i,j)*gV(i,j,k,bi,bj) / deltatmom
678051767f Jean*0140          ENDDO
                0141         ENDDO
9fdcce8056 Jean*0142       ELSEIF (exactConserv) THEN
                0143 c     ELSEIF (nonlinFreeSurf.GT.0) THEN
                0144 C     Implicit treatment of the Barotropic Flow Divergence
                0145         DO j=1,sNy+1
                0146          DO i=1,sNx
                0147           pf(i,j) = implicDiv2Dflow
9f9dfaf749 Jean*0148      &             *yA(i,j)*gV(i,j,k,bi,bj) / deltatmom
9fdcce8056 Jean*0149          ENDDO
                0150         ENDDO
678051767f Jean*0151       ELSE
                0152 C     Explicit+Implicit part of the Barotropic Flow Divergence
                0153         DO j=1,sNy+1
                0154          DO i=1,sNx
9f9dfaf749 Jean*0155           pf(i,j) = ( implicDiv2Dflow * gV(i,j,k,bi,bj)
97c1a61dc2 Jean*0156      &     + (1. _d 0-implicDiv2Dflow)* vVel(i,j,k,bi,bj)
229ac9feb6 Jean*0157      &               ) * yA(i,j) / deltaTMom
678051767f Jean*0158          ENDDO
                0159         ENDDO
                0160       ENDIF
8ae6a40c09 Alis*0161       DO j=1,sNy
                0162        DO i=1,sNx
88830be691 Alis*0163         cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj) +
8ae6a40c09 Alis*0164      &   pf(i,j+1)-pf(i,j)
                0165        ENDDO
                0166       ENDDO
                0167
88830be691 Alis*0168 #ifdef ALLOW_NONHYDROSTATIC
cb7fa97db9 Jean*0169       IF (use3Dsolver) THEN
88830be691 Alis*0170        DO j=1,sNy
                0171         DO i=1,sNx
a7cb4644df Jean*0172          cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
                0173      &                       + ( pf(i,j+1)-pf(i,j) )
88830be691 Alis*0174         ENDDO
                0175        ENDDO
                0176       ENDIF
                0177 #endif
c712b7e80b Chri*0178
d2a11ab670 Jean*0179 #ifdef ALLOW_ADDFLUID
                0180       IF ( selectAddFluid.GE.1 ) THEN
                0181         DO j=1,sNy
                0182          DO i=1,sNx
                0183           cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
229ac9feb6 Jean*0184      &         - addMass(i,j,k,bi,bj)*mass2rUnit/deltaTMom
d2a11ab670 Jean*0185          ENDDO
                0186         ENDDO
                0187 #ifdef ALLOW_NONHYDROSTATIC
                0188        IF (use3Dsolver) THEN
                0189         DO j=1,sNy
                0190          DO i=1,sNx
                0191           cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
229ac9feb6 Jean*0192      &         - addMass(i,j,k,bi,bj)*mass2rUnit/deltaTMom
d2a11ab670 Jean*0193          ENDDO
                0194         ENDDO
                0195        ENDIF
                0196 #endif
                0197       ENDIF
                0198 #endif /* ALLOW_ADDFLUID */
                0199
8ae6a40c09 Alis*0200       RETURN
                0201       END