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
b15f6f1e9f Jean*0176      &                           * implicDiv2DFlow / deltaTMom
                0177      &                           * rAc_3dMean * SQRT( n3dWetPts )
e6e223b277 Jean*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: ',
b15f6f1e9f Jean*0187      &      'cg3dTolerance =', cg3dTolerance,
                0188      &      ' (Area=', rAc_3dMean*SQRT(n3dWetPts), ')'
e6e223b277 Jean*0189         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0190      &                      SQUEEZE_RIGHT, myThid )
                0191        ENDIF
88830be691 Alis*0192        WRITE(msgBuf,*) '                               '
d540b62759 Jean*0193        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0194      &                     SQUEEZE_RIGHT, myThid )
88830be691 Alis*0195       _END_MASTER( myThid )
4606c28752 Jean*0196
88830be691 Alis*0197       DO bj=myByLo(myThid),myByHi(myThid)
                0198        DO bi=myBxLo(myThid),myBxHi(myThid)
e4f5a8bbfa Jean*0199 C-    Set solver main diagonal term
d540b62759 Jean*0200         DO k=1,Nr
                0201          DO j=1,sNy
                0202           DO i=1,sNx
                0203            aW = aW3d( i, j, k, bi,bj)
                0204            aE = aW3d(i+1,j, k, bi,bj)
                0205            aN = aS3d( i,j+1,k, bi,bj)
                0206            aS = aS3d( i, j, k, bi,bj)
                0207            aU = aV3d( i, j, k, bi,bj)
                0208            IF ( k .NE. Nr  ) THEN
                0209             aL = aV3d(i, j,k+1,bi,bj)
e4f5a8bbfa Jean*0210            ELSE
                0211             aL = 0.
                0212            ENDIF
d540b62759 Jean*0213            aC3d(i,j,k,bi,bj) = -aW-aE-aN-aS-aU-aL
e4f5a8bbfa Jean*0214           ENDDO
                0215          ENDDO
                0216         ENDDO
                0217 C-    Add free-surface source term
798efca3d8 Jean*0218         IF ( selectNHfreeSurf.GE.1 ) THEN
                0219          DO j=1,sNy
                0220           DO i=1,sNx
                0221            locGamma = drC(1)*recip_Bo(i,j,bi,bj)
                0222      &              /( deltaTMom*deltaTFreeSurf
                0223      &                *implicitNHPress*implicDiv2DFlow )
                0224            ks = 1
                0225 c          ks = kSurfC(i,j,bi,bj)
                0226 c          IF ( ks.LE.Nr ) THEN
                0227              aC3d(i,j,ks,bi,bj) = aC3d(i,j,ks,bi,bj)
                0228      &         - freeSurfFac*recip_Bo(i,j,bi,bj)
                0229      &          *rA(i,j,bi,bj)*deepFac2F(ks)/deltaTMom/deltaTFreeSurf
                0230      &          / (1. _d 0 + locGamma )
                0231 c          ENDIF
                0232           ENDDO
                0233          ENDDO
                0234         ELSE
d540b62759 Jean*0235          DO j=1,sNy
                0236           DO i=1,sNx
23d1f65433 Jean*0237            ks = kSurfC(i,j,bi,bj)
e4f5a8bbfa Jean*0238            IF ( ks.LE.Nr ) THEN
d540b62759 Jean*0239              aC3d(i,j,ks,bi,bj) = aC3d(i,j,ks,bi,bj)
                0240      &         - freeSurfFac*recip_Bo(i,j,bi,bj)
798efca3d8 Jean*0241      &          *rA(i,j,bi,bj)*deepFac2F(ks)/deltaTMom/deltaTFreeSurf
e4f5a8bbfa Jean*0242            ENDIF
d540b62759 Jean*0243           ENDDO
e4f5a8bbfa Jean*0244          ENDDO
d540b62759 Jean*0245         ENDIF
e4f5a8bbfa Jean*0246 C-    Matrix solver normalisation
d540b62759 Jean*0247         DO k=1,Nr
                0248          DO j=1,sNy
                0249           DO i=1,sNx
                0250            aW3d(i,j,k,bi,bj) = aW3d(i,j,k,bi,bj)*myNorm
                0251            aS3d(i,j,k,bi,bj) = aS3d(i,j,k,bi,bj)*myNorm
                0252            aV3d(i,j,k,bi,bj) = aV3d(i,j,k,bi,bj)*myNorm
                0253            aC3d(i,j,k,bi,bj) = aC3d(i,j,k,bi,bj)*myNorm
88830be691 Alis*0254           ENDDO
                0255          ENDDO
                0256         ENDDO
                0257        ENDDO
                0258       ENDDO
e4f5a8bbfa Jean*0259
88830be691 Alis*0260 C--   Update overlap regions
9e4e1e8791 Jean*0261       CALL EXCH_UV_XYZ_RS(aW3d,aS3d,.FALSE.,myThid)
12c8b75709 Jean*0262       _EXCH_XYZ_RS(aV3d, myThid)
                0263       _EXCH_XYZ_RS(aC3d, myThid)
88830be691 Alis*0264 CcnhDebugStarts
                0265 C     CALL PLOT_FIELD_XYZRS( aW3d, 'AW3D INI_CG3D.1' , Nr, 1, myThid )
                0266 C     CALL PLOT_FIELD_XYZRS( aS3d, 'AS3D INI_CG3D.1' , Nr, 1, myThid )
                0267 CcnhDebugEnds
                0268
                0269 C--   Initialise preconditioner
d540b62759 Jean*0270 C     For now PC is just the identity. Change to
88830be691 Alis*0271 C     be LU factorization of d2/dz2 later. Note
                0272 C     check for consistency with S/R CG3D before
                0273 C     assuming zML is lower and zMU is upper!
                0274       DO bj=myByLo(myThid),myByHi(myThid)
                0275        DO bi=myBxLo(myThid),myBxHi(myThid)
d540b62759 Jean*0276         DO k=1,Nr
                0277          DO j=1,sNy
                0278           DO i=1,sNx
                0279            IF ( aC3d(i,j,k,bi,bj) .NE. 0. ) THEN
                0280             zMC(i,j,k,bi,bj) = aC3d(i,j,k,bi,bj)
                0281             zML(i,j,k,bi,bj) = aV3d(i,j,k,bi,bj)
                0282             IF ( k.NE.Nr ) THEN
                0283              zMU(i,j,k,bi,bj)= aV3d(i,j,k+1,bi,bj)
                0284             ELSE
                0285              zMU(i,j,k,bi,bj)= 0.
                0286             ENDIF
88830be691 Alis*0287 CcnhDebugStarts
                0288 C           zMC(i,j,k,bi,bj) = 1.
                0289 C           zMU(i,j,k,bi,bj) = 0.
                0290 C           zML(i,j,k,bi,bj) = 0.
                0291 CcnhDebugEnds
                0292            ELSE
                0293             zMC(i,j,k,bi,bj) = 1. _d 0
                0294             zMU(i,j,k,bi,bj) = 0.
                0295             zML(i,j,k,bi,bj) = 0.
                0296            ENDIF
                0297           ENDDO
                0298          ENDDO
                0299         ENDDO
d540b62759 Jean*0300         k = 1
                0301          DO j=1,sNy
                0302           DO i=1,sNx
                0303            zMC(i,j,k,bi,bj) = 1. _d 0 / zMC(i,j,k,bi,bj)
                0304            zMU(i,j,k,bi,bj) = zMU(i,j,k,bi,bj)*zMC(i,j,k,bi,bj)
                0305           ENDDO
88830be691 Alis*0306          ENDDO
d540b62759 Jean*0307         DO k=2,Nr
                0308          DO j=1,sNy
                0309           DO i=1,sNx
88830be691 Alis*0310            zMC(i,j,k,bi,bj) = 1. _d 0 /
                0311      &     (zMC(i,j,k,bi,bj)-zML(i,j,k,bi,bj)*zMU(i,j,k-1,bi,bj))
d540b62759 Jean*0312            zMU(i,j,k,bi,bj) = zMU(i,j,k,bi,bj)*zMC(i,j,k,bi,bj)
88830be691 Alis*0313           ENDDO
                0314          ENDDO
                0315         ENDDO
d540b62759 Jean*0316         DO k=1,Nr
                0317          DO j=1,sNy
                0318           DO i=1,sNx
                0319            IF ( aC3d(i,j,k,bi,bj) .EQ. 0. ) THEN
9bec0e0933 Alis*0320             zMC(i,j,k,bi,bj) = 1.
88830be691 Alis*0321             zML(i,j,k,bi,bj) = 0.
                0322             zMU(i,j,k,bi,bj) = 0.
                0323 CcnhDebugStarts
                0324 C          ELSE
                0325 C           zMC(i,j,k,bi,bj) = 1.
                0326 C           zML(i,j,k,bi,bj) = 0.
                0327 C           zMU(i,j,k,bi,bj) = 0.
                0328 CcnhDEbugEnds
                0329            ENDIF
                0330           ENDDO
                0331          ENDDO
                0332         ENDDO
                0333        ENDDO
                0334       ENDDO
                0335 C--   Update overlap regions
12c8b75709 Jean*0336       _EXCH_XYZ_RS(zMC, myThid)
                0337       _EXCH_XYZ_RS(zML, myThid)
                0338       _EXCH_XYZ_RS(zMU, myThid)
88830be691 Alis*0339
23d1f65433 Jean*0340       IF ( debugLevel .GE. debLevC ) THEN
d540b62759 Jean*0341         CALL WRITE_FLD_XYZ_RS( 'zMC',' ',zMC, 0, myThid )
                0342         CALL WRITE_FLD_XYZ_RS( 'zML',' ',zML, 0, myThid )
                0343         CALL WRITE_FLD_XYZ_RS( 'zMU',' ',zMU, 0, myThid )
                0344       ENDIF
88830be691 Alis*0345 CcnhDebugStarts
9e4e1e8791 Jean*0346 c     DO k=1,Nr
                0347 c     DO j=1-OLy,sNy+OLy
                0348 c     DO i=1-OLx,sNx+OLx
                0349 c      phi(i,j,1,1) = zMc(i,j,k,1,1)
                0350 c     ENDDO
                0351 c     ENDDO
88830be691 Alis*0352 C     CALL PLOT_FIELD_XYRS( phi, 'zMC INI_CG3D.1' , 1, myThid )
9e4e1e8791 Jean*0353 c     ENDDO
88830be691 Alis*0354 C     CALL PLOT_FIELD_XYRS( zMU, 'zMU INI_CG3D.1' , Nr, 1, myThid )
                0355 C     CALL PLOT_FIELD_XYRS( zML, 'zML INI_CG3D.1' , Nr, 1, myThid )
                0356 CcnhDebugEnds
                0357
cb7fa97db9 Jean*0358 C--   end if (use3Dsolver)
                0359       ENDIF
                0360
88830be691 Alis*0361 #endif /* ALLOW_NONHYDROSTATIC */
                0362
                0363       RETURN
                0364       END