Back to home page

MITgcm
 
 

    

File indexing completed on 2021-08-25 05:10:56 UTC

view on githubraw file Latest commit aecc8b0f on 2021-08-25 02:39:15 UTC
6d54cf9ca1 Ed H*0001 #include "PACKAGES_CONFIG.h"
73470161ee Alis*0002 #include "CPP_OPTIONS.h"
924557e60a Chri*0003 
9366854e02 Chri*0004 CBOP
                0005 C     !ROUTINE: INI_CG2D
                0006 C     !INTERFACE:
924557e60a Chri*0007       SUBROUTINE INI_CG2D( myThid )
                0008 
9366854e02 Chri*0009 C     !DESCRIPTION: \bv
                0010 C     *==========================================================*
9155c64ca4 Jean*0011 C     | SUBROUTINE INI_CG2D
                0012 C     | o Initialise 2d conjugate gradient solver operators.
9366854e02 Chri*0013 C     *==========================================================*
9155c64ca4 Jean*0014 C     | These arrays are purely a function of the basin geom.
                0015 C     | We set then here once and them use then repeatedly.
9366854e02 Chri*0016 C     *==========================================================*
                0017 C     \ev
                0018 
                0019 C     !USES:
                0020       IMPLICIT NONE
924557e60a Chri*0021 C     === Global variables ===
                0022 #include "SIZE.h"
                0023 #include "EEPARAMS.h"
                0024 #include "PARAMS.h"
                0025 #include "GRID.h"
71c6a09c16 Jean*0026 #include "SURFACE.h"
aea29c8517 Alis*0027 #include "CG2D.h"
924557e60a Chri*0028 
9366854e02 Chri*0029 C     !INPUT/OUTPUT PARAMETERS:
924557e60a Chri*0030 C     === Routine arguments ===
                0031 C     myThid - Thread no. that called this routine.
                0032       INTEGER myThid
                0033 
9366854e02 Chri*0034 C     !LOCAL VARIABLES:
924557e60a Chri*0035 C     === Local variables ===
4606c28752 Jean*0036 C     bi,bj  :: tile indices
63f3ae4253 Jean*0037 C     i,j,k  :: Loop counters
7cf1bc4c9a Jean*0038 C     faceArea :: Temporary used to hold cell face areas.
                0039 C     myNorm   :: Work variable used in calculating normalisation factor
924557e60a Chri*0040       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0041       INTEGER bi, bj
63f3ae4253 Jean*0042       INTEGER i, j, k, ks
5c76561bf0 Chri*0043       _RL     faceArea
e7217683b5 Chri*0044       _RS     myNorm
4606c28752 Jean*0045       _RS     aC, aCw, aCs
9366854e02 Chri*0046 CEOP
924557e60a Chri*0047 
4606c28752 Jean*0048 C--   Initialize arrays in common blocs (CG2D.h) ; not really necessary
9e9a4cf401 Jean*0049 C     but safer when EXCH do not fill all the overlap regions.
                0050       DO bj=myByLo(myThid),myByHi(myThid)
                0051        DO bi=myBxLo(myThid),myBxHi(myThid)
63f3ae4253 Jean*0052         DO j=1-OLy,sNy+OLy
                0053          DO i=1-OLx,sNx+OLx
                0054           aW2d(i,j,bi,bj) = 0. _d 0
                0055           aS2d(i,j,bi,bj) = 0. _d 0
                0056           aC2d(i,j,bi,bj) = 0. _d 0
                0057           pW(i,j,bi,bj) = 0. _d 0
                0058           pS(i,j,bi,bj) = 0. _d 0
                0059           pC(i,j,bi,bj) = 0. _d 0
9e9a4cf401 Jean*0060          ENDDO
                0061         ENDDO
                0062        ENDDO
                0063       ENDDO
                0064 
81f3a86735 Patr*0065 C--   Init. scalars
                0066       cg2dNorm = 0. _d 0
                0067       cg2dNormaliseRHS = .FALSE.
                0068       cg2dtolerance = 0. _d 0
                0069 
924557e60a Chri*0070 C--   Initialise laplace operator
                0071 C     aW2d: integral in Z Ax/dX
                0072 C     aS2d: integral in Z Ay/dY
                0073       myNorm = 0. _d 0
                0074       DO bj=myByLo(myThid),myByHi(myThid)
                0075        DO bi=myBxLo(myThid),myBxHi(myThid)
63f3ae4253 Jean*0076         DO j=1,sNy
                0077          DO i=1,sNx
                0078           aW2d(i,j,bi,bj) = 0. _d 0
                0079           aS2d(i,j,bi,bj) = 0. _d 0
924557e60a Chri*0080          ENDDO
                0081         ENDDO
63f3ae4253 Jean*0082         DO k=1,Nr
                0083          DO j=1,sNy
                0084           DO i=1,sNx
4606c28752 Jean*0085 C  deep-model: *deepFacC (faceArea), /deepFacC (recip_dx,y): => no net effect
63f3ae4253 Jean*0086            faceArea = _dyG(i,j,bi,bj)*drF(k)
                0087      &               *_hFacW(i,j,k,bi,bj)
                0088            aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)
4606c28752 Jean*0089      &              + implicSurfPress*implicDiv2DFlow
63f3ae4253 Jean*0090      &               *faceArea*recip_dxC(i,j,bi,bj)
                0091            faceArea = _dxG(i,j,bi,bj)*drF(k)
                0092      &               *_hFacS(i,j,k,bi,bj)
                0093            aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)
4606c28752 Jean*0094      &              + implicSurfPress*implicDiv2DFlow
63f3ae4253 Jean*0095      &               *faceArea*recip_dyC(i,j,bi,bj)
924557e60a Chri*0096           ENDDO
                0097          ENDDO
                0098         ENDDO
63f3ae4253 Jean*0099         DO j=1,sNy
                0100          DO i=1,sNx
7cf1bc4c9a Jean*0101 #ifdef ALLOW_OBCS
                0102            aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)
                0103      &                   *maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj)
                0104            aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)
                0105      &                   *maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj)
                0106 #endif /* ALLOW_OBCS */
                0107            myNorm = MAX(ABS(aW2d(i,j,bi,bj)),myNorm)
                0108            myNorm = MAX(ABS(aS2d(i,j,bi,bj)),myNorm)
924557e60a Chri*0109          ENDDO
                0110         ENDDO
                0111        ENDDO
                0112       ENDDO
12c8b75709 Jean*0113       _GLOBAL_MAX_RS( myNorm, myThid )
78c96bee2f Alis*0114       IF ( myNorm .NE. 0. _d 0 ) THEN
                0115        myNorm = 1. _d 0/myNorm
924557e60a Chri*0116       ELSE
                0117        myNorm = 1. _d 0
                0118       ENDIF
                0119       DO bj=myByLo(myThid),myByHi(myThid)
                0120        DO bi=myBxLo(myThid),myBxHi(myThid)
63f3ae4253 Jean*0121         DO j=1,sNy
                0122          DO i=1,sNx
                0123           aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)*myNorm
                0124           aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)*myNorm
924557e60a Chri*0125          ENDDO
                0126         ENDDO
                0127        ENDDO
                0128       ENDDO
4606c28752 Jean*0129 
924557e60a Chri*0130 C--   Update overlap regions
                0131 CcnhDebugStarts
42bd47f06f Chri*0132 C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.1' , 1, myThid )
                0133 C     CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.1' , 1, myThid )
924557e60a Chri*0134 CcnhDebugEnds
9155c64ca4 Jean*0135       CALL EXCH_UV_XY_RS( aW2d, aS2d, .FALSE., myThid )
924557e60a Chri*0136 CcnhDebugStarts
88830be691 Alis*0137 C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )
                0138 C     CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.2' , 1, myThid )
924557e60a Chri*0139 CcnhDebugEnds
                0140 
06bb0cec77 Jean*0141       _BEGIN_MASTER(myThid)
                0142 C-- set global parameter in common block:
                0143       cg2dNorm = myNorm
aea29c8517 Alis*0144 C-- Define the solver tolerance in the appropriate Unit :
0417f72e5a Jean*0145       cg2dNormaliseRHS = cg2dTargetResWunit.LE.0.
aea29c8517 Alis*0146       IF (cg2dNormaliseRHS) THEN
                0147 C-  when using a normalisation of RHS, tolerance has no unit => no conversion
                0148         cg2dTolerance = cg2dTargetResidual
                0149       ELSE
                0150 C-  convert Target-Residual (in W unit) to cg2d-solver residual unit [m^2/s^2]
                0151         cg2dTolerance = cg2dNorm * cg2dTargetResWunit
aa6b2555c8 Jean*0152      &                           * globalArea / deltaTMom
aea29c8517 Alis*0153       ENDIF
06bb0cec77 Jean*0154       _END_MASTER(myThid)
0417f72e5a Jean*0155 
                0156 CcnhDebugStarts
                0157       _BEGIN_MASTER( myThid )
                0158        WRITE(msgBuf,'(2A,1PE23.16)') 'INI_CG2D: ',
                0159      &      'CG2D normalisation factor = ', cg2dNorm
                0160        CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
                0161        IF (.NOT.cg2dNormaliseRHS) THEN
                0162         WRITE(msgBuf,'(2A,1PE22.15,A,1PE16.10,A)') 'INI_CG2D: ',
                0163      &      'cg2dTolerance =', cg2dTolerance, ' (Area=',globalArea,')'
                0164         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
                0165        ENDIF
                0166        WRITE(msgBuf,*) ' '
                0167        CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
                0168       _END_MASTER( myThid )
                0169 CcnhDebugEnds
4606c28752 Jean*0170 
924557e60a Chri*0171 C--   Initialise preconditioner
181bb5586b Chri*0172 C     Note. 20th May 1998
                0173 C           I made a weird discovery! In the model paper we argue
                0174 C           for the form of the preconditioner used here ( see
                0175 C           A Finite-volume, Incompressible Navier-Stokes Model
                0176 C           ...., Marshall et. al ). The algebra gives a simple
                0177 C           0.5 factor for the averaging of ac and aCw to get a
                0178 C           symmettric pre-conditioner. By using a factor of 0.51
                0179 C           i.e. scaling the off-diagonal terms in the
                0180 C           preconditioner down slightly I managed to get the
                0181 C           number of iterations for convergence in a test case to
                0182 C           drop form 192 -> 134! Need to investigate this further!
                0183 C           For now I have introduced a parameter cg2dpcOffDFac which
                0184 C           defaults to 0.51 but can be set at runtime.
924557e60a Chri*0185       DO bj=myByLo(myThid),myByHi(myThid)
                0186        DO bi=myBxLo(myThid),myBxHi(myThid)
63f3ae4253 Jean*0187 C-    calculate and store solver main diagonal :
aa6b2555c8 Jean*0188         DO j=0,sNy
                0189          DO i=0,sNx
                0190            ks = kSurfC(i,j,bi,bj)
63f3ae4253 Jean*0191            aC2d(i,j,bi,bj) = -(
                0192      &       aW2d(i,j,bi,bj) + aW2d(i+1,j, bi,bj)
                0193      &      +aS2d(i,j,bi,bj) + aS2d( i,j+1,bi,bj)
                0194      &      +freeSurfFac*myNorm*recip_Bo(i,j,bi,bj)*deepFac2F(ks)
aa6b2555c8 Jean*0195      &                  *rA(i,j,bi,bj)/deltaTMom/deltaTFreeSurf
63f3ae4253 Jean*0196      &                        )
                0197          ENDDO
                0198         ENDDO
                0199         DO j=1,sNy
                0200          DO i=1,sNx
                0201            aC  = aC2d( i, j, bi,bj)
                0202            aCs = aC2d( i,j-1,bi,bj)
                0203            aCw = aC2d(i-1,j, bi,bj)
                0204            IF ( aC .EQ. 0. ) THEN
                0205             pC(i,j,bi,bj) = 1. _d 0
                0206            ELSE
                0207             pC(i,j,bi,bj) =  1. _d 0 / aC
                0208            ENDIF
                0209            IF ( aC + aCw .EQ. 0. ) THEN
                0210             pW(i,j,bi,bj) = 0.
                0211            ELSE
                0212             pW(i,j,bi,bj) =
                0213      &         -aW2d(i,j,bi,bj)/((cg2dpcOffDFac *(aCw+aC))**2 )
                0214            ENDIF
                0215            IF ( aC + aCs .EQ. 0. ) THEN
                0216             pS(i,j,bi,bj) = 0.
                0217            ELSE
                0218             pS(i,j,bi,bj) =
                0219      &         -aS2d(i,j,bi,bj)/((cg2dpcOffDFac *(aCs+aC))**2 )
                0220            ENDIF
                0221 C          pC(i,j,bi,bj) = 1.
                0222 C          pW(i,j,bi,bj) = 0.
                0223 C          pS(i,j,bi,bj) = 0.
924557e60a Chri*0224          ENDDO
                0225         ENDDO
                0226        ENDDO
                0227       ENDDO
                0228 C--   Update overlap regions
9155c64ca4 Jean*0229       CALL EXCH_XY_RS( pC, myThid )
                0230       CALL EXCH_UV_XY_RS( pW, pS, .FALSE., myThid )
255afb6daa Chri*0231 CcnhDebugStarts
88830be691 Alis*0232 C     CALL PLOT_FIELD_XYRS( pC, 'pC   INI_CG2D.2' , 1, myThid )
                0233 C     CALL PLOT_FIELD_XYRS( pW, 'pW   INI_CG2D.2' , 1, myThid )
                0234 C     CALL PLOT_FIELD_XYRS( pS, 'pS   INI_CG2D.2' , 1, myThid )
255afb6daa Chri*0235 CcnhDebugEnds
924557e60a Chri*0236 
                0237       RETURN
                0238       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!