Back to home page

MITgcm
 
 

    

File indexing completed on 2018-03-02 18:36:46 UTC

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

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!