Back to home page

MITgcm
 
 

    

File indexing completed on 2024-09-20 05:10:19 UTC

view on githubraw file Latest commit e6e223b2 on 2024-09-19 22:01: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
e6e223b277 Jean*0044       _RL     cg2dTolerance
e7217683b5 Chri*0045       _RS     myNorm
4606c28752 Jean*0046       _RS     aC, aCw, aCs
9366854e02 Chri*0047 CEOP
924557e60a Chri*0048 
4606c28752 Jean*0049 C--   Initialize arrays in common blocs (CG2D.h) ; not really necessary
9e9a4cf401 Jean*0050 C     but safer when EXCH do not fill all the overlap regions.
                0051       DO bj=myByLo(myThid),myByHi(myThid)
                0052        DO bi=myBxLo(myThid),myBxHi(myThid)
63f3ae4253 Jean*0053         DO j=1-OLy,sNy+OLy
                0054          DO i=1-OLx,sNx+OLx
                0055           aW2d(i,j,bi,bj) = 0. _d 0
                0056           aS2d(i,j,bi,bj) = 0. _d 0
                0057           aC2d(i,j,bi,bj) = 0. _d 0
                0058           pW(i,j,bi,bj) = 0. _d 0
                0059           pS(i,j,bi,bj) = 0. _d 0
                0060           pC(i,j,bi,bj) = 0. _d 0
9e9a4cf401 Jean*0061          ENDDO
                0062         ENDDO
                0063        ENDDO
                0064       ENDDO
                0065 
81f3a86735 Patr*0066 C--   Init. scalars
e6e223b277 Jean*0067       _BEGIN_MASTER(myThid)
81f3a86735 Patr*0068       cg2dNorm = 0. _d 0
                0069       cg2dNormaliseRHS = .FALSE.
e6e223b277 Jean*0070       cg2dTolerance_sq = 0. _d 0
                0071       _END_MASTER(myThid)
81f3a86735 Patr*0072 
924557e60a Chri*0073 C--   Initialise laplace operator
                0074 C     aW2d: integral in Z Ax/dX
                0075 C     aS2d: integral in Z Ay/dY
                0076       myNorm = 0. _d 0
                0077       DO bj=myByLo(myThid),myByHi(myThid)
                0078        DO bi=myBxLo(myThid),myBxHi(myThid)
63f3ae4253 Jean*0079         DO j=1,sNy
                0080          DO i=1,sNx
                0081           aW2d(i,j,bi,bj) = 0. _d 0
                0082           aS2d(i,j,bi,bj) = 0. _d 0
924557e60a Chri*0083          ENDDO
                0084         ENDDO
63f3ae4253 Jean*0085         DO k=1,Nr
                0086          DO j=1,sNy
                0087           DO i=1,sNx
4606c28752 Jean*0088 C  deep-model: *deepFacC (faceArea), /deepFacC (recip_dx,y): => no net effect
63f3ae4253 Jean*0089            faceArea = _dyG(i,j,bi,bj)*drF(k)
                0090      &               *_hFacW(i,j,k,bi,bj)
                0091            aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)
4606c28752 Jean*0092      &              + implicSurfPress*implicDiv2DFlow
63f3ae4253 Jean*0093      &               *faceArea*recip_dxC(i,j,bi,bj)
                0094            faceArea = _dxG(i,j,bi,bj)*drF(k)
                0095      &               *_hFacS(i,j,k,bi,bj)
                0096            aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)
4606c28752 Jean*0097      &              + implicSurfPress*implicDiv2DFlow
63f3ae4253 Jean*0098      &               *faceArea*recip_dyC(i,j,bi,bj)
924557e60a Chri*0099           ENDDO
                0100          ENDDO
                0101         ENDDO
63f3ae4253 Jean*0102         DO j=1,sNy
                0103          DO i=1,sNx
7cf1bc4c9a Jean*0104 #ifdef ALLOW_OBCS
                0105            aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)
                0106      &                   *maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj)
                0107            aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)
                0108      &                   *maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj)
                0109 #endif /* ALLOW_OBCS */
                0110            myNorm = MAX(ABS(aW2d(i,j,bi,bj)),myNorm)
                0111            myNorm = MAX(ABS(aS2d(i,j,bi,bj)),myNorm)
924557e60a Chri*0112          ENDDO
                0113         ENDDO
                0114        ENDDO
                0115       ENDDO
12c8b75709 Jean*0116       _GLOBAL_MAX_RS( myNorm, myThid )
e6e223b277 Jean*0117       IF ( myNorm .NE. zeroRS ) THEN
78c96bee2f Alis*0118        myNorm = 1. _d 0/myNorm
924557e60a Chri*0119       ELSE
                0120        myNorm = 1. _d 0
                0121       ENDIF
                0122       DO bj=myByLo(myThid),myByHi(myThid)
                0123        DO bi=myBxLo(myThid),myBxHi(myThid)
63f3ae4253 Jean*0124         DO j=1,sNy
                0125          DO i=1,sNx
                0126           aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)*myNorm
                0127           aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)*myNorm
924557e60a Chri*0128          ENDDO
                0129         ENDDO
                0130        ENDDO
                0131       ENDDO
4606c28752 Jean*0132 
924557e60a Chri*0133 C--   Update overlap regions
                0134 CcnhDebugStarts
42bd47f06f Chri*0135 C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.1' , 1, myThid )
                0136 C     CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.1' , 1, myThid )
924557e60a Chri*0137 CcnhDebugEnds
9155c64ca4 Jean*0138       CALL EXCH_UV_XY_RS( aW2d, aS2d, .FALSE., myThid )
924557e60a Chri*0139 CcnhDebugStarts
88830be691 Alis*0140 C     CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.2' , 1, myThid )
                0141 C     CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.2' , 1, myThid )
924557e60a Chri*0142 CcnhDebugEnds
                0143 
06bb0cec77 Jean*0144       _BEGIN_MASTER(myThid)
                0145 C-- set global parameter in common block:
                0146       cg2dNorm = myNorm
aea29c8517 Alis*0147 C-- Define the solver tolerance in the appropriate Unit :
e6e223b277 Jean*0148       cg2dNormaliseRHS = cg2dTargetResWunit.LE.zeroRL
aea29c8517 Alis*0149       IF (cg2dNormaliseRHS) THEN
                0150 C-  when using a normalisation of RHS, tolerance has no unit => no conversion
                0151         cg2dTolerance = cg2dTargetResidual
                0152       ELSE
                0153 C-  convert Target-Residual (in W unit) to cg2d-solver residual unit [m^2/s^2]
                0154         cg2dTolerance = cg2dNorm * cg2dTargetResWunit
aa6b2555c8 Jean*0155      &                           * globalArea / deltaTMom
aea29c8517 Alis*0156       ENDIF
e6e223b277 Jean*0157       cg2dTolerance_sq = cg2dTolerance*cg2dTolerance
06bb0cec77 Jean*0158       _END_MASTER(myThid)
0417f72e5a Jean*0159 
                0160 CcnhDebugStarts
                0161       _BEGIN_MASTER( myThid )
                0162        WRITE(msgBuf,'(2A,1PE23.16)') 'INI_CG2D: ',
                0163      &      'CG2D normalisation factor = ', cg2dNorm
                0164        CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
                0165        IF (.NOT.cg2dNormaliseRHS) THEN
                0166         WRITE(msgBuf,'(2A,1PE22.15,A,1PE16.10,A)') 'INI_CG2D: ',
                0167      &      'cg2dTolerance =', cg2dTolerance, ' (Area=',globalArea,')'
                0168         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
                0169        ENDIF
                0170        WRITE(msgBuf,*) ' '
                0171        CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
                0172       _END_MASTER( myThid )
                0173 CcnhDebugEnds
4606c28752 Jean*0174 
924557e60a Chri*0175 C--   Initialise preconditioner
181bb5586b Chri*0176 C     Note. 20th May 1998
                0177 C           I made a weird discovery! In the model paper we argue
                0178 C           for the form of the preconditioner used here ( see
                0179 C           A Finite-volume, Incompressible Navier-Stokes Model
                0180 C           ...., Marshall et. al ). The algebra gives a simple
                0181 C           0.5 factor for the averaging of ac and aCw to get a
                0182 C           symmettric pre-conditioner. By using a factor of 0.51
                0183 C           i.e. scaling the off-diagonal terms in the
                0184 C           preconditioner down slightly I managed to get the
                0185 C           number of iterations for convergence in a test case to
                0186 C           drop form 192 -> 134! Need to investigate this further!
                0187 C           For now I have introduced a parameter cg2dpcOffDFac which
                0188 C           defaults to 0.51 but can be set at runtime.
924557e60a Chri*0189       DO bj=myByLo(myThid),myByHi(myThid)
                0190        DO bi=myBxLo(myThid),myBxHi(myThid)
63f3ae4253 Jean*0191 C-    calculate and store solver main diagonal :
aa6b2555c8 Jean*0192         DO j=0,sNy
                0193          DO i=0,sNx
                0194            ks = kSurfC(i,j,bi,bj)
63f3ae4253 Jean*0195            aC2d(i,j,bi,bj) = -(
                0196      &       aW2d(i,j,bi,bj) + aW2d(i+1,j, bi,bj)
                0197      &      +aS2d(i,j,bi,bj) + aS2d( i,j+1,bi,bj)
                0198      &      +freeSurfFac*myNorm*recip_Bo(i,j,bi,bj)*deepFac2F(ks)
aa6b2555c8 Jean*0199      &                  *rA(i,j,bi,bj)/deltaTMom/deltaTFreeSurf
63f3ae4253 Jean*0200      &                        )
                0201          ENDDO
                0202         ENDDO
                0203         DO j=1,sNy
                0204          DO i=1,sNx
                0205            aC  = aC2d( i, j, bi,bj)
                0206            aCs = aC2d( i,j-1,bi,bj)
                0207            aCw = aC2d(i-1,j, bi,bj)
e6e223b277 Jean*0208            IF ( aC .EQ. zeroRS ) THEN
63f3ae4253 Jean*0209             pC(i,j,bi,bj) = 1. _d 0
                0210            ELSE
                0211             pC(i,j,bi,bj) =  1. _d 0 / aC
                0212            ENDIF
e6e223b277 Jean*0213            IF ( aC + aCw .EQ. zeroRS ) THEN
63f3ae4253 Jean*0214             pW(i,j,bi,bj) = 0.
                0215            ELSE
                0216             pW(i,j,bi,bj) =
                0217      &         -aW2d(i,j,bi,bj)/((cg2dpcOffDFac *(aCw+aC))**2 )
                0218            ENDIF
e6e223b277 Jean*0219            IF ( aC + aCs .EQ. zeroRS ) THEN
63f3ae4253 Jean*0220             pS(i,j,bi,bj) = 0.
                0221            ELSE
                0222             pS(i,j,bi,bj) =
                0223      &         -aS2d(i,j,bi,bj)/((cg2dpcOffDFac *(aCs+aC))**2 )
                0224            ENDIF
                0225 C          pC(i,j,bi,bj) = 1.
                0226 C          pW(i,j,bi,bj) = 0.
                0227 C          pS(i,j,bi,bj) = 0.
924557e60a Chri*0228          ENDDO
                0229         ENDDO
                0230        ENDDO
                0231       ENDDO
                0232 C--   Update overlap regions
9155c64ca4 Jean*0233       CALL EXCH_XY_RS( pC, myThid )
                0234       CALL EXCH_UV_XY_RS( pW, pS, .FALSE., myThid )
255afb6daa Chri*0235 CcnhDebugStarts
88830be691 Alis*0236 C     CALL PLOT_FIELD_XYRS( pC, 'pC   INI_CG2D.2' , 1, myThid )
                0237 C     CALL PLOT_FIELD_XYRS( pW, 'pW   INI_CG2D.2' , 1, myThid )
                0238 C     CALL PLOT_FIELD_XYRS( pS, 'pS   INI_CG2D.2' , 1, myThid )
255afb6daa Chri*0239 CcnhDebugEnds
924557e60a Chri*0240 
                0241       RETURN
                0242       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!