model/src/ini_cg3d.F

f293c3e8bc Ed H*0001 #include "PACKAGES_CONFIG.h"
88830be691 Alis*0002 #include "CPP_OPTIONS.h"
                0003
9366854e02 Chri*0004 CBOP
9e4e1e8791 Jean*0005 C     !ROUTINE: INI_CG3D
9366854e02 Chri*0006 C     !INTERFACE:
88830be691 Alis*0007       SUBROUTINE INI_CG3D( myThid )
9366854e02 Chri*0008 C     !DESCRIPTION: \bv
                0009 C     *==========================================================*
d540b62759 Jean*0010 C     | SUBROUTINE INI_CG3D
                0011 C     | o Initialise 3d conjugate gradient solver operators.
9366854e02 Chri*0012 C     *==========================================================*
d540b62759 Jean*0013 C     | These arrays are purely a function of the basin geom.
                0014 C     | We set then here once and them use then repeatedly.
9366854e02 Chri*0015 C     *==========================================================*
                0016 C     \ev
88830be691 Alis*0017
9366854e02 Chri*0018 C     !USES:
                0019       IMPLICIT NONE
88830be691 Alis*0020 C     === Global variables ===
                0021 #include "SIZE.h"
                0022 #include "EEPARAMS.h"
e6e223b277 Jean*0023 #ifdef ALLOW_NONHYDROSTATIC
                0024 # include "PARAMS.h"
                0025 # include "GRID.h"
                0026 # include "SURFACE.h"
                0027 # include "CG3D.h"
                0028 #endif /* ALLOW_NONHYDROSTATIC */
88830be691 Alis*0029
9366854e02 Chri*0030 C     !INPUT/OUTPUT PARAMETERS:
d540b62759 Jean*0031 C     myThid   :: My Thread Id number
88830be691 Alis*0032       INTEGER myThid
                0033
87f515096d Alis*0034 #ifdef ALLOW_NONHYDROSTATIC
9366854e02 Chri*0035 C     !LOCAL VARIABLES:
d540b62759 Jean*0036 C     bi,bj    :: tile indices
                0037 C     i,j,k    :: Loop counters
                0038 C     faceArea :: Temporary used to hold cell face areas.
                0039 C     myNorm   :: Work variable used in clculating normalisation factor
88830be691 Alis*0040       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0041       INTEGER bi, bj
d540b62759 Jean*0042       INTEGER i, j, k, ks
88830be691 Alis*0043       _RL     faceArea
                0044       _RS     myNorm
                0045       _RL     theRecip_Dr
e6e223b277 Jean*0046       _RL     cg3dTolerance
d540b62759 Jean*0047       _RL     aU, aL, aW, aE, aN, aS
cb7fa97db9 Jean*0048       _RL     tmpFac, nh_Fac, igwFac
798efca3d8 Jean*0049       _RL     locGamma
9366854e02 Chri*0050 CEOP
88830be691 Alis*0051
                0052 CcnhDebugStarts
9e4e1e8791 Jean*0053 c     _RL    phi(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
88830be691 Alis*0054 CcnhDebugEnds
                0055
9e4e1e8791 Jean*0056 C--   Initialise to zero over the full range of indices
                0057       DO bj=myByLo(myThid),myByHi(myThid)
                0058        DO bi=myBxLo(myThid),myBxHi(myThid)
d540b62759 Jean*0059         DO k=1,Nr
9e4e1e8791 Jean*0060 C-       From common bloc CG3D_R:
d540b62759 Jean*0061          DO j=1-OLy,sNy+OLy
                0062           DO i=1-OLx,sNx+OLx
                0063            aW3d(i,j,k,bi,bj) = 0.
                0064            aS3d(i,j,k,bi,bj) = 0.
                0065            aV3d(i,j,k,bi,bj) = 0.
                0066            aC3d(i,j,k,bi,bj) = 0.
                0067            zMC (i,j,k,bi,bj) = 0.
                0068            zML (i,j,k,bi,bj) = 0.
                0069            zMU (i,j,k,bi,bj) = 0.
9e4e1e8791 Jean*0070           ENDDO
                0071          ENDDO
                0072 C-       From common bloc CG3D_WK_R:
d540b62759 Jean*0073          DO j=0,sNy+1
                0074           DO i=0,sNx+1
                0075            cg3d_q(i,j,k,bi,bj) = 0.
                0076            cg3d_r(i,j,k,bi,bj) = 0.
                0077            cg3d_s(i,j,k,bi,bj) = 0.
9e4e1e8791 Jean*0078           ENDDO
                0079          ENDDO
                0080         ENDDO
                0081        ENDDO
                0082       ENDDO
                0083
cb7fa97db9 Jean*0084       nh_Fac = 0.
                0085       igwFac = 0.
                0086       IF ( nonHydrostatic
                0087      &      .AND. nh_Am2.NE.0. ) nh_Fac = 1. _d 0 / nh_Am2
                0088       IF ( implicitIntGravWave ) igwFac = 1. _d 0
                0089
                0090       IF ( use3Dsolver ) THEN
88830be691 Alis*0091 C--   Initialise laplace operator
                0092 C     aW3d: Ax/dX
                0093 C     aS3d: Ay/dY
                0094 C     aV3d: Ar/dR
                0095       myNorm = 0. _d 0
                0096       DO bj=myByLo(myThid),myByHi(myThid)
                0097        DO bi=myBxLo(myThid),myBxHi(myThid)
d540b62759 Jean*0098         DO k=1,Nr
                0099          DO j=1,sNy
                0100           DO i=1,sNx+1
                0101            faceArea = _dyG(i,j,bi,bj)*drF(k)
6b917c5a38 Jean*0102      &               *_hFacW(i,j,k,bi,bj)
                0103 #ifdef ALLOW_OBCS
                0104      &               *maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj)
                0105 #endif
d540b62759 Jean*0106            aW3d(i,j,k,bi,bj) = faceArea*recip_dxC(i,j,bi,bj)
                0107      &                        *implicitNHPress*implicDiv2DFlow
                0108            myNorm = MAX(ABS(aW3d(i,j,k,bi,bj)),myNorm)
e4f5a8bbfa Jean*0109           ENDDO
                0110          ENDDO
4606c28752 Jean*0111 C  deep-model: *deepFacC (faceArea), /deepFacC (recip_dx,y): => no net effect
d540b62759 Jean*0112          DO j=1,sNy+1
                0113           DO i=1,sNx
e6e223b277 Jean*0114            faceArea = _dxG(i,j,bi,bj)*drF(k)
6b917c5a38 Jean*0115      &               *_hFacS(i,j,k,bi,bj)
                0116 #ifdef ALLOW_OBCS
                0117      &               *maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj)
                0118 #endif
d540b62759 Jean*0119            aS3d(i,j,k,bi,bj) = faceArea*recip_dyC(i,j,bi,bj)
                0120      &                        *implicitNHPress*implicDiv2DFlow
                0121            myNorm = MAX(ABS(aS3d(i,j,k,bi,bj)),myNorm)
88830be691 Alis*0122           ENDDO
                0123          ENDDO
                0124         ENDDO
d540b62759 Jean*0125         DO k=1,1
                0126          DO j=1,sNy
                0127           DO i=1,sNx
                0128            aV3d(i,j,k,bi,bj) =  0.
88830be691 Alis*0129           ENDDO
                0130          ENDDO
                0131         ENDDO
d540b62759 Jean*0132         DO k=2,Nr
f056be794c Jean*0133          tmpFac = nh_Fac*rVel2wUnit(k)*rVel2wUnit(k)
cb7fa97db9 Jean*0134      &          + igwFac*dBdrRef(k)*deltaTMom*dTtracerLev(k)
                0135          IF (tmpFac.GT.0. ) tmpFac = 1. _d 0 / tmpFac
d540b62759 Jean*0136          DO j=1,sNy
                0137           DO i=1,sNx
                0138            faceArea = _rA(i,j,bi,bj)*maskC(i,j, k ,bi,bj)
                0139      &                              *maskC(i,j,k-1,bi,bj)
4606c28752 Jean*0140      &                              *deepFac2F(k)
6b917c5a38 Jean*0141 #ifdef ALLOW_OBCS
                0142      &                              *maskInC(i,j,bi,bj)
                0143 #endif
d540b62759 Jean*0144            theRecip_Dr = recip_drC(k)
4606c28752 Jean*0145 c          theRecip_Dr =
d540b62759 Jean*0146 caja &      drF(k  )*_hFacC(i,j,k  ,bi,bj)*0.5
                0147 caja &     +drF(k-1)*_hFacC(i,j,k-1,bi,bj)*0.5
4606c28752 Jean*0148 c          IF ( theRecip_Dr .NE. 0. )
9e4e1e8791 Jean*0149 c    &      theRecip_Dr = 1. _d 0/theRecip_Dr
d540b62759 Jean*0150            aV3d(i,j,k,bi,bj) = faceArea*theRecip_Dr*tmpFac
                0151      &                        *implicitNHPress*implicDiv2DFlow
                0152            myNorm = MAX(ABS(aV3d(i,j,k,bi,bj)),myNorm)
88830be691 Alis*0153           ENDDO
                0154          ENDDO
9bec0e0933 Alis*0155         ENDDO
88830be691 Alis*0156        ENDDO
                0157       ENDDO
12c8b75709 Jean*0158       _GLOBAL_MAX_RS( myNorm, myThid )
78c96bee2f Alis*0159       IF ( myNorm .NE. 0. _d 0 ) THEN
                0160        myNorm = 1. _d 0/myNorm
88830be691 Alis*0161       ELSE
                0162        myNorm = 1. _d 0
                0163       ENDIF
4606c28752 Jean*0164
88830be691 Alis*0165       _BEGIN_MASTER( myThid )
4606c28752 Jean*0166 C-- set global parameter in common block:
06bb0cec77 Jean*0167        cg3dNorm = myNorm
e6e223b277 Jean*0168 C-- Define the solver tolerance in the appropriate Unit :
                0169       cg3dNormaliseRHS = cg3dTargetResWunit.LE.zeroRL
                0170       IF (cg3dNormaliseRHS) THEN
                0171 C-  when using a normalisation of RHS, tolerance has no unit => no conversion
                0172         cg3dTolerance = cg3dTargetResidual
                0173       ELSE
                0174 C-  convert Target-Residual (in W unit) to cg3d-solver residual unit [m^2/s^2]
                0175         cg3dTolerance = cg3dNorm * cg3dTargetResWunit
79b5d5775c Jean*0176      &                           * implicDiv2DFlow
e6e223b277 Jean*0177      &                           * globalArea / deltaTMom
                0178       ENDIF
                0179       cg3dTolerance_sq = cg3dTolerance*cg3dTolerance
                0180 C-- print param values:
4606c28752 Jean*0181        WRITE(msgBuf,'(2A,1PE23.16)') 'INI_CG3D: ',
                0182      &      'CG3D normalisation factor = ', cg3dNorm
d540b62759 Jean*0183        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0184      &                     SQUEEZE_RIGHT, myThid )
e6e223b277 Jean*0185        IF (.NOT.cg3dNormaliseRHS) THEN
                0186         WRITE(msgBuf,'(2A,1PE22.15,A,1PE16.10,A)') 'INI_CG3D: ',
                0187      &      'cg3dTolerance =', cg3dTolerance, ' (Area=',globalArea,')'
                0188         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0189      &                      SQUEEZE_RIGHT, myThid )
                0190        ENDIF
88830be691 Alis*0191        WRITE(msgBuf,*) '                               '
d540b62759 Jean*0192        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0193      &                     SQUEEZE_RIGHT, myThid )
88830be691 Alis*0194       _END_MASTER( myThid )
4606c28752 Jean*0195
88830be691 Alis*0196       DO bj=myByLo(myThid),myByHi(myThid)
                0197        DO bi=myBxLo(myThid),myBxHi(myThid)
e4f5a8bbfa Jean*0198 C-    Set solver main diagonal term
d540b62759 Jean*0199         DO k=1,Nr
                0200          DO j=1,sNy
                0201           DO i=1,sNx
                0202            aW = aW3d( i, j, k, bi,bj)
                0203            aE = aW3d(i+1,j, k, bi,bj)
                0204            aN = aS3d( i,j+1,k, bi,bj)
                0205            aS = aS3d( i, j, k, bi,bj)
                0206            aU = aV3d( i, j, k, bi,bj)
                0207            IF ( k .NE. Nr  ) THEN
                0208             aL = aV3d(i, j,k+1,bi,bj)
e4f5a8bbfa Jean*0209            ELSE
                0210             aL = 0.
                0211            ENDIF
d540b62759 Jean*0212            aC3d(i,j,k,bi,bj) = -aW-aE-aN-aS-aU-aL
e4f5a8bbfa Jean*0213           ENDDO
                0214          ENDDO
                0215         ENDDO
                0216 C-    Add free-surface source term
798efca3d8 Jean*0217         IF ( selectNHfreeSurf.GE.1 ) THEN
                0218          DO j=1,sNy
                0219           DO i=1,sNx
                0220            locGamma = drC(1)*recip_Bo(i,j,bi,bj)
                0221      &              /( deltaTMom*deltaTFreeSurf
                0222      &                *implicitNHPress*implicDiv2DFlow )
                0223            ks = 1
                0224 c          ks = kSurfC(i,j,bi,bj)
                0225 c          IF ( ks.LE.Nr ) THEN
                0226              aC3d(i,j,ks,bi,bj) = aC3d(i,j,ks,bi,bj)
                0227      &         - freeSurfFac*recip_Bo(i,j,bi,bj)
                0228      &          *rA(i,j,bi,bj)*deepFac2F(ks)/deltaTMom/deltaTFreeSurf
                0229      &          / (1. _d 0 + locGamma )
                0230 c          ENDIF
                0231           ENDDO
                0232          ENDDO
                0233         ELSE
d540b62759 Jean*0234          DO j=1,sNy
                0235           DO i=1,sNx
23d1f65433 Jean*0236            ks = kSurfC(i,j,bi,bj)
e4f5a8bbfa Jean*0237            IF ( ks.LE.Nr ) THEN
d540b62759 Jean*0238              aC3d(i,j,ks,bi,bj) = aC3d(i,j,ks,bi,bj)
                0239      &         - freeSurfFac*recip_Bo(i,j,bi,bj)
798efca3d8 Jean*0240      &          *rA(i,j,bi,bj)*deepFac2F(ks)/deltaTMom/deltaTFreeSurf
e4f5a8bbfa Jean*0241            ENDIF
d540b62759 Jean*0242           ENDDO
e4f5a8bbfa Jean*0243          ENDDO
d540b62759 Jean*0244         ENDIF
e4f5a8bbfa Jean*0245 C-    Matrix solver normalisation
d540b62759 Jean*0246         DO k=1,Nr
                0247          DO j=1,sNy
                0248           DO i=1,sNx
                0249            aW3d(i,j,k,bi,bj) = aW3d(i,j,k,bi,bj)*myNorm
                0250            aS3d(i,j,k,bi,bj) = aS3d(i,j,k,bi,bj)*myNorm
                0251            aV3d(i,j,k,bi,bj) = aV3d(i,j,k,bi,bj)*myNorm
                0252            aC3d(i,j,k,bi,bj) = aC3d(i,j,k,bi,bj)*myNorm
88830be691 Alis*0253           ENDDO
                0254          ENDDO
                0255         ENDDO
                0256        ENDDO
                0257       ENDDO
e4f5a8bbfa Jean*0258
88830be691 Alis*0259 C--   Update overlap regions
9e4e1e8791 Jean*0260       CALL EXCH_UV_XYZ_RS(aW3d,aS3d,.FALSE.,myThid)
12c8b75709 Jean*0261       _EXCH_XYZ_RS(aV3d, myThid)
                0262       _EXCH_XYZ_RS(aC3d, myThid)
88830be691 Alis*0263 CcnhDebugStarts
                0264 C     CALL PLOT_FIELD_XYZRS( aW3d, 'AW3D INI_CG3D.1' , Nr, 1, myThid )
                0265 C     CALL PLOT_FIELD_XYZRS( aS3d, 'AS3D INI_CG3D.1' , Nr, 1, myThid )
                0266 CcnhDebugEnds
                0267
                0268 C--   Initialise preconditioner
d540b62759 Jean*0269 C     For now PC is just the identity. Change to
88830be691 Alis*0270 C     be LU factorization of d2/dz2 later. Note
                0271 C     check for consistency with S/R CG3D before
                0272 C     assuming zML is lower and zMU is upper!
                0273       DO bj=myByLo(myThid),myByHi(myThid)
                0274        DO bi=myBxLo(myThid),myBxHi(myThid)
d540b62759 Jean*0275         DO k=1,Nr
                0276          DO j=1,sNy
                0277           DO i=1,sNx
                0278            IF ( aC3d(i,j,k,bi,bj) .NE. 0. ) THEN
                0279             zMC(i,j,k,bi,bj) = aC3d(i,j,k,bi,bj)
                0280             zML(i,j,k,bi,bj) = aV3d(i,j,k,bi,bj)
                0281             IF ( k.NE.Nr ) THEN
                0282              zMU(i,j,k,bi,bj)= aV3d(i,j,k+1,bi,bj)
                0283             ELSE
                0284              zMU(i,j,k,bi,bj)= 0.
                0285             ENDIF
88830be691 Alis*0286 CcnhDebugStarts
                0287 C           zMC(i,j,k,bi,bj) = 1.
                0288 C           zMU(i,j,k,bi,bj) = 0.
                0289 C           zML(i,j,k,bi,bj) = 0.
                0290 CcnhDebugEnds
                0291            ELSE
                0292             zMC(i,j,k,bi,bj) = 1. _d 0
                0293             zMU(i,j,k,bi,bj) = 0.
                0294             zML(i,j,k,bi,bj) = 0.
                0295            ENDIF
                0296           ENDDO
                0297          ENDDO
                0298         ENDDO
d540b62759 Jean*0299         k = 1
                0300          DO j=1,sNy
                0301           DO i=1,sNx
                0302            zMC(i,j,k,bi,bj) = 1. _d 0 / zMC(i,j,k,bi,bj)
                0303            zMU(i,j,k,bi,bj) = zMU(i,j,k,bi,bj)*zMC(i,j,k,bi,bj)
                0304           ENDDO
88830be691 Alis*0305          ENDDO
d540b62759 Jean*0306         DO k=2,Nr
                0307          DO j=1,sNy
                0308           DO i=1,sNx
88830be691 Alis*0309            zMC(i,j,k,bi,bj) = 1. _d 0 /
                0310      &     (zMC(i,j,k,bi,bj)-zML(i,j,k,bi,bj)*zMU(i,j,k-1,bi,bj))
d540b62759 Jean*0311            zMU(i,j,k,bi,bj) = zMU(i,j,k,bi,bj)*zMC(i,j,k,bi,bj)
88830be691 Alis*0312           ENDDO
                0313          ENDDO
                0314         ENDDO
d540b62759 Jean*0315         DO k=1,Nr
                0316          DO j=1,sNy
                0317           DO i=1,sNx
                0318            IF ( aC3d(i,j,k,bi,bj) .EQ. 0. ) THEN
9bec0e0933 Alis*0319             zMC(i,j,k,bi,bj) = 1.
88830be691 Alis*0320             zML(i,j,k,bi,bj) = 0.
                0321             zMU(i,j,k,bi,bj) = 0.
                0322 CcnhDebugStarts
                0323 C          ELSE
                0324 C           zMC(i,j,k,bi,bj) = 1.
                0325 C           zML(i,j,k,bi,bj) = 0.
                0326 C           zMU(i,j,k,bi,bj) = 0.
                0327 CcnhDEbugEnds
                0328            ENDIF
                0329           ENDDO
                0330          ENDDO
                0331         ENDDO
                0332        ENDDO
                0333       ENDDO
                0334 C--   Update overlap regions
12c8b75709 Jean*0335       _EXCH_XYZ_RS(zMC, myThid)
                0336       _EXCH_XYZ_RS(zML, myThid)
                0337       _EXCH_XYZ_RS(zMU, myThid)
88830be691 Alis*0338
23d1f65433 Jean*0339       IF ( debugLevel .GE. debLevC ) THEN
d540b62759 Jean*0340         CALL WRITE_FLD_XYZ_RS( 'zMC',' ',zMC, 0, myThid )
                0341         CALL WRITE_FLD_XYZ_RS( 'zML',' ',zML, 0, myThid )
                0342         CALL WRITE_FLD_XYZ_RS( 'zMU',' ',zMU, 0, myThid )
                0343       ENDIF
88830be691 Alis*0344 CcnhDebugStarts
9e4e1e8791 Jean*0345 c     DO k=1,Nr
                0346 c     DO j=1-OLy,sNy+OLy
                0347 c     DO i=1-OLx,sNx+OLx
                0348 c      phi(i,j,1,1) = zMc(i,j,k,1,1)
                0349 c     ENDDO
                0350 c     ENDDO
88830be691 Alis*0351 C     CALL PLOT_FIELD_XYRS( phi, 'zMC INI_CG3D.1' , 1, myThid )
9e4e1e8791 Jean*0352 c     ENDDO
88830be691 Alis*0353 C     CALL PLOT_FIELD_XYRS( zMU, 'zMU INI_CG3D.1' , Nr, 1, myThid )
                0354 C     CALL PLOT_FIELD_XYRS( zML, 'zML INI_CG3D.1' , Nr, 1, myThid )
                0355 CcnhDebugEnds
                0356
cb7fa97db9 Jean*0357 C--   end if (use3Dsolver)
                0358       ENDIF
                0359
88830be691 Alis*0360 #endif /* ALLOW_NONHYDROSTATIC */
                0361
                0362       RETURN
                0363       END