model/src/ini_cg2d.F

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
79b5d5775c Jean*0149       IF ( cg2dNormaliseRHS ) THEN
aea29c8517 Alis*0150 C-  when using a normalisation of RHS, tolerance has no unit => no conversion
                0151         cg2dTolerance = cg2dTargetResidual
79b5d5775c Jean*0152       ELSEIF ( implicDiv2DFlow.GT.zeroRL ) THEN
aea29c8517 Alis*0153 C-  convert Target-Residual (in W unit) to cg2d-solver residual unit [m^2/s^2]
                0154         cg2dTolerance = cg2dNorm * cg2dTargetResWunit
79b5d5775c Jean*0155      &                           * implicDiv2DFlow
                0156      &                           * globalArea / deltaTMom
                0157       ELSE
                0158 C-  special case implicDiv2DFlow=0 :
                0159 C   not much to solve anyway, cg2dTolerance does not matter as long as it is > 0
                0160         cg2dTolerance = cg2dNorm * cg2dTargetResWunit
aa6b2555c8 Jean*0161      &                           * globalArea / deltaTMom
aea29c8517 Alis*0162       ENDIF
e6e223b277 Jean*0163       cg2dTolerance_sq = cg2dTolerance*cg2dTolerance
06bb0cec77 Jean*0164       _END_MASTER(myThid)
0417f72e5a Jean*0165
                0166 CcnhDebugStarts
                0167       _BEGIN_MASTER( myThid )
                0168        WRITE(msgBuf,'(2A,1PE23.16)') 'INI_CG2D: ',
                0169      &      'CG2D normalisation factor = ', cg2dNorm
                0170        CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
                0171        IF (.NOT.cg2dNormaliseRHS) THEN
                0172         WRITE(msgBuf,'(2A,1PE22.15,A,1PE16.10,A)') 'INI_CG2D: ',
                0173      &      'cg2dTolerance =', cg2dTolerance, ' (Area=',globalArea,')'
                0174         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
                0175        ENDIF
                0176        WRITE(msgBuf,*) ' '
                0177        CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
                0178       _END_MASTER( myThid )
                0179 CcnhDebugEnds
4606c28752 Jean*0180
924557e60a Chri*0181 C--   Initialise preconditioner
181bb5586b Chri*0182 C     Note. 20th May 1998
                0183 C           I made a weird discovery! In the model paper we argue
                0184 C           for the form of the preconditioner used here ( see
                0185 C           A Finite-volume, Incompressible Navier-Stokes Model
                0186 C           ...., Marshall et. al ). The algebra gives a simple
                0187 C           0.5 factor for the averaging of ac and aCw to get a
                0188 C           symmettric pre-conditioner. By using a factor of 0.51
                0189 C           i.e. scaling the off-diagonal terms in the
                0190 C           preconditioner down slightly I managed to get the
                0191 C           number of iterations for convergence in a test case to
                0192 C           drop form 192 -> 134! Need to investigate this further!
                0193 C           For now I have introduced a parameter cg2dpcOffDFac which
                0194 C           defaults to 0.51 but can be set at runtime.
924557e60a Chri*0195       DO bj=myByLo(myThid),myByHi(myThid)
                0196        DO bi=myBxLo(myThid),myBxHi(myThid)
63f3ae4253 Jean*0197 C-    calculate and store solver main diagonal :
aa6b2555c8 Jean*0198         DO j=0,sNy
                0199          DO i=0,sNx
                0200            ks = kSurfC(i,j,bi,bj)
63f3ae4253 Jean*0201            aC2d(i,j,bi,bj) = -(
                0202      &       aW2d(i,j,bi,bj) + aW2d(i+1,j, bi,bj)
                0203      &      +aS2d(i,j,bi,bj) + aS2d( i,j+1,bi,bj)
                0204      &      +freeSurfFac*myNorm*recip_Bo(i,j,bi,bj)*deepFac2F(ks)
aa6b2555c8 Jean*0205      &                  *rA(i,j,bi,bj)/deltaTMom/deltaTFreeSurf
63f3ae4253 Jean*0206      &                        )
                0207          ENDDO
                0208         ENDDO
                0209         DO j=1,sNy
                0210          DO i=1,sNx
                0211            aC  = aC2d( i, j, bi,bj)
                0212            aCs = aC2d( i,j-1,bi,bj)
                0213            aCw = aC2d(i-1,j, bi,bj)
e6e223b277 Jean*0214            IF ( aC .EQ. zeroRS ) THEN
63f3ae4253 Jean*0215             pC(i,j,bi,bj) = 1. _d 0
                0216            ELSE
                0217             pC(i,j,bi,bj) =  1. _d 0 / aC
                0218            ENDIF
e6e223b277 Jean*0219            IF ( aC + aCw .EQ. zeroRS ) THEN
63f3ae4253 Jean*0220             pW(i,j,bi,bj) = 0.
                0221            ELSE
                0222             pW(i,j,bi,bj) =
                0223      &         -aW2d(i,j,bi,bj)/((cg2dpcOffDFac *(aCw+aC))**2 )
                0224            ENDIF
e6e223b277 Jean*0225            IF ( aC + aCs .EQ. zeroRS ) THEN
63f3ae4253 Jean*0226             pS(i,j,bi,bj) = 0.
                0227            ELSE
                0228             pS(i,j,bi,bj) =
                0229      &         -aS2d(i,j,bi,bj)/((cg2dpcOffDFac *(aCs+aC))**2 )
                0230            ENDIF
                0231 C          pC(i,j,bi,bj) = 1.
                0232 C          pW(i,j,bi,bj) = 0.
                0233 C          pS(i,j,bi,bj) = 0.
924557e60a Chri*0234          ENDDO
                0235         ENDDO
                0236        ENDDO
                0237       ENDDO
                0238 C--   Update overlap regions
9155c64ca4 Jean*0239       CALL EXCH_XY_RS( pC, myThid )
                0240       CALL EXCH_UV_XY_RS( pW, pS, .FALSE., myThid )
255afb6daa Chri*0241 CcnhDebugStarts
88830be691 Alis*0242 C     CALL PLOT_FIELD_XYRS( pC, 'pC   INI_CG2D.2' , 1, myThid )
                0243 C     CALL PLOT_FIELD_XYRS( pW, 'pW   INI_CG2D.2' , 1, myThid )
                0244 C     CALL PLOT_FIELD_XYRS( pS, 'pS   INI_CG2D.2' , 1, myThid )
255afb6daa Chri*0245 CcnhDebugEnds
924557e60a Chri*0246
                0247       RETURN
                0248       END