model/src/pre_cg3d.F

8991be1c08 Jean*0001 #include "PACKAGES_CONFIG.h"
                0002 #include "CPP_OPTIONS.h"
                0003
                0004 CBOP
                0005 C     !ROUTINE: PRE_CG3D
                0006 C     !INTERFACE:
                0007       SUBROUTINE PRE_CG3D(
                0008      I                     oldFreeSurfTerm,
                0009      I                     cg2d_x,
                0010      U                     cg3d_b,
                0011      I                     myTime, myIter, myThid )
                0012
                0013 C     !DESCRIPTION:
                0014 C     Called from SOLVE_FOR_PRESSURE, before 3-D solver (cg3d):
                0015 C     Finish calculation of 3-D RHS after 2-D inversionis done.
                0016
                0017 C     !USES:
                0018       IMPLICIT NONE
                0019 C     == Global variables
                0020 #include "SIZE.h"
                0021 #include "EEPARAMS.h"
                0022 #include "PARAMS.h"
                0023 #include "GRID.h"
                0024 #include "SURFACE.h"
                0025 #include "FFIELDS.h"
                0026 #include "DYNVARS.h"
                0027 #ifdef ALLOW_NONHYDROSTATIC
                0028 #include "NH_VARS.h"
                0029 #endif
                0030
                0031 C     === Functions ====
                0032 c     LOGICAL  DIFFERENT_MULTIPLE
                0033 c     EXTERNAL DIFFERENT_MULTIPLE
                0034
                0035 C     !INPUT/OUTPUT PARAMETERS:
                0036 C     == Routine arguments ==
                0037 C     oldFreeSurfTerm :: Treat free-surface term in the old way (no exactConserv)
                0038 C     cg2d_x          :: Solution vector of the 2-D solver equation a.x=b
                0039 C     cg3d_b          :: Right Hand side vector of the 3-D solver equation A.X=B
                0040 C     myTime          :: Current time in simulation
                0041 C     myIter          :: Current iteration number in simulation
                0042 C     myThid          :: My Thread Id number
                0043       LOGICAL oldFreeSurfTerm
                0044       _RL     cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0045       _RL     cg3d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0046       _RL     myTime
                0047       INTEGER myIter
                0048       INTEGER myThid
                0049
                0050 #ifdef ALLOW_NONHYDROSTATIC
                0051 C     !LOCAL VARIABLES:
                0052 C     == Local variables ==
8e18cb9146 Jean*0053 C     wSurfP2d :: surface vertical velocity after 2-D solver
8991be1c08 Jean*0054       INTEGER i,j,k,bi,bj
                0055       INTEGER ks, kp1
                0056 c     CHARACTER*10 sufx
                0057 c     CHARACTER*(MAX_LEN_MBUF) msgBuf
8e18cb9146 Jean*0058       _RL     locGamma, surfFac, tmpFac
8991be1c08 Jean*0059       _RL     wFacKm, wFacKp
bf53796fe2 Jean*0060       _RL     uf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0061       _RL     vf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
8e18cb9146 Jean*0062       _RL     wSurfP2d(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
8991be1c08 Jean*0063 CEOP
                0064
                0065 c     IF ( use3Dsolver ) THEN
                0066
                0067 C--   Solve for a three-dimensional pressure term (NH or IGW or both ).
                0068 C     see CG3D.h for the interface to this routine.
                0069        DO bj=myByLo(myThid),myByHi(myThid)
                0070         DO bi=myBxLo(myThid),myBxHi(myThid)
                0071
8e18cb9146 Jean*0072 C--   Calculate updated (after 2-D solver) vertical velocity at the surface
                0073          IF ( oldFreeSurfTerm .OR. implicDiv2DFlow.EQ.zeroRL ) THEN
                0074            DO j=1-OLy,sNy+OLy
                0075             DO i=1-OLx,sNx+OLx
                0076               wSurfP2d(i,j) = 0. _d 0
                0077             ENDDO
                0078            ENDDO
                0079          ELSE
                0080            DO j=1-OLy,sNy+OLy
                0081             DO i=1-OLx,sNx+OLx
                0082                wSurfP2d(i,j) = ( etaN(i,j,bi,bj)-etaH(i,j,bi,bj) )
                0083      &                       / ( implicDiv2DFlow*deltaTFreeSurf )
                0084             ENDDO
                0085            ENDDO
                0086          ENDIF
                0087
                0088 C--   Add EmPmR contribution to top level cg3d_b or to wSurfP2d:
8991be1c08 Jean*0089 C      (has been done for cg2d_b ; and addMass was added by CALC_DIV_GHAT)
8e18cb9146 Jean*0090          IF ( useRealFreshWaterFlux.AND.fluidIsWater ) THEN
                0091           IF ( oldFreeSurfTerm .OR. usingPCoords ) THEN
                0092            tmpFac = freeSurfFac*mass2rUnit*implicDiv2DFlow/deltaTMom
                0093            ks = 1
                0094            IF ( usingPCoords ) ks = Nr
                0095            DO j=1,sNy
                0096             DO i=1,sNx
                0097               cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
                0098      &          + tmpFac*_rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)
                0099      &                                 *maskInC(i,j,bi,bj)
                0100             ENDDO
8991be1c08 Jean*0101            ENDDO
8e18cb9146 Jean*0102           ELSE
                0103            DO j=1-OLy,sNy+OLy
                0104             DO i=1-OLx,sNx+OLx
                0105               wSurfP2d(i,j) = wSurfP2d(i,j)
                0106      &                      + EmPmR(i,j,bi,bj)*mass2rUnit
                0107      &                                 *maskInC(i,j,bi,bj)
                0108             ENDDO
                0109            ENDDO
                0110           ENDIF
                0111          ENDIF
8991be1c08 Jean*0112
                0113 C--   Update or Add free-surface contribution to cg3d_b:
8e18cb9146 Jean*0114          surfFac = 0.
                0115          IF ( selectNHfreeSurf.GE.1 ) THEN
                0116            tmpFac = freeSurfFac*implicDiv2DFlow/deltaTMom
8991be1c08 Jean*0117            DO j=1,sNy
                0118             DO i=1,sNx
8e18cb9146 Jean*0119               locGamma = drC(1)*recip_Bo(i,j,bi,bj)
                0120      &                 /( deltaTMom*deltaTFreeSurf
                0121      &                   *implicitNHPress*implicDiv2DFlow )
                0122               ks = 1
                0123 c             ks = kSurfC(i,j,bi,bj)
                0124 c             IF ( ks.LE.Nr ) THEN
8991be1c08 Jean*0125                cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
8e18cb9146 Jean*0126      &                      + tmpFac*( wSurfP2d(i,j)
                0127      &                               + locGamma*wVel(i,j,ks,bi,bj) )
                0128      &                              /( 1. _d 0 + locGamma )
8991be1c08 Jean*0129      &                              *_rA(i,j,bi,bj)*deepFac2F(ks)
8e18cb9146 Jean*0130 c             ENDIF
                0131 C-    Save wSurfP2d (used in POST_CG3D) into dPhiNH :
                0132               dPhiNH(i,j,bi,bj) = wSurfP2d(i,j)
8991be1c08 Jean*0133             ENDDO
                0134            ENDDO
8e18cb9146 Jean*0135          ELSEIF ( .NOT.oldFreeSurfTerm ) THEN
                0136            tmpFac = freeSurfFac*implicDiv2DFlow/deltaTMom
8991be1c08 Jean*0137            DO j=1,sNy
                0138             DO i=1,sNx
bf53796fe2 Jean*0139               ks = kSurfC(i,j,bi,bj)
8e18cb9146 Jean*0140               IF ( ks.LE.Nr ) THEN
                0141                cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
                0142      &                      + tmpFac*wSurfP2d(i,j)
                0143      &                              *_rA(i,j,bi,bj)*deepFac2F(ks)
                0144               ENDIF
8991be1c08 Jean*0145             ENDDO
                0146            ENDDO
e78345ed77 Jean*0147          ELSEIF ( uniformFreeSurfLev ) THEN
8991be1c08 Jean*0148 C-       Z coordinate: assume surface @ level k=1
8e18cb9146 Jean*0149            surfFac = freeSurfFac*deepFac2F(1)
8991be1c08 Jean*0150          ELSE
                0151 C-       Other than Z coordinate: no assumption on surface level index
                0152            DO j=1,sNy
                0153             DO i=1,sNx
bf53796fe2 Jean*0154               ks = kSurfC(i,j,bi,bj)
8991be1c08 Jean*0155               IF ( ks.LE.Nr ) THEN
                0156                cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
8e18cb9146 Jean*0157      &              +freeSurfFac*etaN(i,j,bi,bj)/deltaTFreeSurf
                0158      &                  *_rA(i,j,bi,bj)*deepFac2F(ks)/deltaTMom
8991be1c08 Jean*0159               ENDIF
                0160             ENDDO
                0161            ENDDO
                0162          ENDIF
                0163
                0164 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0165
                0166 C--   Finish updating cg3d_b: 1) increment in horiz velocity due to new cg2d_x
                0167 C                             2) add vertical velocity contribution.
                0168          DO j=1,sNy+1
                0169           DO i=1,sNx+1
                0170            uf(i,j) = -_recip_dxC(i,j,bi,bj)
                0171      &             * implicSurfPress*implicDiv2DFlow
                0172      &             *(cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
fc5959006e Jean*0173 #ifdef ALLOW_OBCS
                0174      &             *maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj)
                0175 #endif
8991be1c08 Jean*0176            vf(i,j) = -_recip_dyC(i,j,bi,bj)
                0177      &             * implicSurfPress*implicDiv2DFlow
                0178      &             *(cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
                0179 #ifdef ALLOW_OBCS
fc5959006e Jean*0180      &             *maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj)
                0181 #endif
8991be1c08 Jean*0182           ENDDO
fc5959006e Jean*0183          ENDDO
8991be1c08 Jean*0184
                0185 C Note: with implicDiv2DFlow < 1, wVel contribution to cg3d_b is similar to
                0186 C       uVel,vVel contribution to cg2d_b when exactConserv=T, since wVel is
                0187 C       always recomputed from continuity eq (like eta when exactConserv=T)
                0188          k=1
                0189          kp1 = MIN(k+1,Nr)
                0190          wFacKp = implicDiv2DFlow*deepFac2F(kp1)*rhoFacF(kp1)
                0191          IF (k.GE.Nr) wFacKp = 0.
                0192          DO j=1,sNy
                0193           DO i=1,sNx
                0194             cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
                0195      &       +drF(k)*dyG(i+1,j,bi,bj)*_hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
                0196      &       -drF(k)*dyG( i ,j,bi,bj)*_hFacW( i ,j,k,bi,bj)*uf( i ,j)
                0197      &       +drF(k)*dxG(i,j+1,bi,bj)*_hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
                0198      &       -drF(k)*dxG(i, j ,bi,bj)*_hFacS(i, j ,k,bi,bj)*vf(i, j )
8e18cb9146 Jean*0199      &       +( surfFac*etaN(i,j,bi,bj)/deltaTFreeSurf
8991be1c08 Jean*0200      &         -wVel(i,j,kp1,bi,bj)*wFacKp
8e18cb9146 Jean*0201      &        )*_rA(i,j,bi,bj)/deltaTMom
8991be1c08 Jean*0202           ENDDO
                0203          ENDDO
                0204          DO k=2,Nr
                0205           kp1 = MIN(k+1,Nr)
                0206 C-       deepFac & rhoFac cancel with the ones in uf[=del_i(Phi)/dx],vf ;
                0207 C        both appear in wVel term, but at 2 different levels
                0208           wFacKm = implicDiv2DFlow*deepFac2F( k )*rhoFacF( k )
                0209           wFacKp = implicDiv2DFlow*deepFac2F(kp1)*rhoFacF(kp1)
                0210           IF (k.GE.Nr) wFacKp = 0.
                0211           DO j=1,sNy
                0212            DO i=1,sNx
                0213             cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
                0214      &       +drF(k)*dyG(i+1,j,bi,bj)*_hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
                0215      &       -drF(k)*dyG( i ,j,bi,bj)*_hFacW( i ,j,k,bi,bj)*uf( i ,j)
                0216      &       +drF(k)*dxG(i,j+1,bi,bj)*_hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
                0217      &       -drF(k)*dxG(i, j ,bi,bj)*_hFacS(i, j ,k,bi,bj)*vf(i, j )
                0218      &       +( wVel(i,j, k ,bi,bj)*wFacKm*maskC(i,j,k-1,bi,bj)
                0219      &         -wVel(i,j,kp1,bi,bj)*wFacKp
8e18cb9146 Jean*0220      &        )*_rA(i,j,bi,bj)/deltaTMom
8991be1c08 Jean*0221            ENDDO
                0222           ENDDO
                0223          ENDDO
                0224
                0225 #ifdef ALLOW_OBCS
fc5959006e Jean*0226 C- Note: solver matrix is trivial outside OB region (main diagonal only)
                0227 C     => no real need to reset RHS (=cg3d_b) & cg3d_x, except that:
                0228 C    a) normalisation is fct of Max(RHS), which can be large ouside OB region
                0229 C      (would be different if we were solving for increment of phi_nh
                0230 C       instead of directly for phi_nh).
                0231 C       => need to reset RHS to ensure that interior solution does not depend
                0232 C       on ouside OB region.
                0233 C    b) provide directly the trivial solution cg3d_x == 0 for outside OB region
                0234 C      (=> no residual => no effect on solver convergence and interior solution)
8991be1c08 Jean*0235          IF (useOBCS) THEN
                0236           DO k=1,Nr
                0237            DO j=1,sNy
fc5959006e Jean*0238             DO i=1,sNx
                0239               cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
                0240      &                            *maskInC(i,j,bi,bj)
                0241               phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj)
                0242      &                            *maskInC(i,j,bi,bj)
                0243             ENDDO
8991be1c08 Jean*0244            ENDDO
                0245           ENDDO
                0246          ENDIF
                0247 #endif /* ALLOW_OBCS */
                0248
                0249 C-    end bi,bj loops
                0250         ENDDO
                0251        ENDDO
                0252
                0253 c     ENDIF
                0254 #endif /* ALLOW_NONHYDROSTATIC */
                0255
                0256       RETURN
                0257       END