Back to home page

MITgcm
 
 

    

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

view on githubraw file Latest commit e6e223b2 on 2024-09-19 22:01:15 UTC
015c19607c Jean*0001 #include "CPP_OPTIONS.h"
                0002 #ifdef TARGET_NEC_SX
                0003 C     set a sensible default for the outer loop unrolling parameter that can
                0004 C     be overriden in the Makefile with the DEFINES macro or in CPP_OPTIONS.h
                0005 #ifndef CG3D_OUTERLOOPITERS
                0006 # define CG3D_OUTERLOOPITERS 10
                0007 #endif
                0008 #endif /* TARGET_NEC_SX */
                0009 
                0010 CBOP
                0011 C     !ROUTINE: CG3D_EX0
                0012 C     !INTERFACE:
                0013       SUBROUTINE CG3D_EX0(
                0014      U                cg3d_b, cg3d_x,
                0015      O                firstResidual, lastResidual,
                0016      U                numIters,
                0017      I                myIter, myThid )
                0018 C     !DESCRIPTION: \bv
                0019 C     *==========================================================*
                0020 C     | SUBROUTINE CG3D_EX0
                0021 C     | o Three-dimensional grid problem conjugate-gradient
                0022 C     |   inverter (with preconditioner).
                0023 C     | This is the disconnected-tile version (each tile treated
                0024 C     |   independently as a small domain, with locally periodic
                0025 C     |   BC at the edges.
                0026 C     *==========================================================*
                0027 C     | Con. grad is an iterative procedure for solving Ax = b.
                0028 C     | It requires the A be symmetric.
                0029 C     | This implementation assumes A is a seven-diagonal matrix
                0030 C     | of the form that arises in the discrete representation of
                0031 C     | the del^2 operator in a three-dimensional space.
                0032 C     *==========================================================*
                0033 C     \ev
                0034 
                0035 C     !USES:
                0036       IMPLICIT NONE
                0037 C     === Global data ===
                0038 #include "SIZE.h"
                0039 #include "EEPARAMS.h"
e6e223b277 Jean*0040 #ifdef ALLOW_NONHYDROSTATIC
                0041 # include "PARAMS.h"
                0042 # include "GRID.h"
                0043 # include "SURFACE.h"
                0044 # include "CG3D.h"
                0045 #endif /* ALLOW_NONHYDROSTATIC */
015c19607c Jean*0046 
                0047 C     !INPUT/OUTPUT PARAMETERS:
                0048 C     cg3d_b    :: The source term or "right hand side" (output: normalised RHS)
                0049 C     cg3d_x    :: The solution (input: first guess)
                0050 C     firstResidual :: the initial residual before any iterations
                0051 C     minResidualSq :: the lowest residual reached (squared)
                0052 CC    lastResidual  :: the actual residual reached
                0053 C     numIters  :: Inp: the maximum number of iterations allowed
                0054 C                  Out: the actual number of iterations used
                0055 CC    nIterMin  :: Inp: decide to store (if >=0) or not (if <0) lowest res. sol.
                0056 CC                 Out: iteration number corresponding to lowest residual
                0057 C     myIter    :: Current iteration number in simulation
                0058 C     myThid    :: my Thread Id number
                0059       _RL     cg3d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0060       _RL     cg3d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0061       _RL     firstResidual
                0062       _RL     lastResidual
                0063       INTEGER numIters
                0064       INTEGER myIter
                0065       INTEGER myThid
                0066 
                0067 #ifdef ALLOW_NONHYDROSTATIC
                0068 C     !LOCAL VARIABLES:
                0069 C     bi, bj     :: tile index in X and Y.
                0070 C     i, j, k    :: Loop counters
                0071 C     it3d       :: Loop counter for CG iterations
                0072 C     actualIts  :: actual CG iteration number
                0073 C     err_sq     :: Measure of the square of the residual of Ax - b.
                0074 C     eta_qrN    :: Used in computing search directions; suffix N and NM1
                0075 C     eta_qrNM1     denote current and previous iterations respectively.
                0076 C     cgBeta     :: coeff used to update conjugate direction vector "s".
                0077 C     alpha      :: coeff used to update solution & residual
                0078 C     sumRHS     :: Sum of right-hand-side. Sometimes this is a useful
                0079 C                   debugging/trouble shooting diagnostic. For neumann problems
                0080 C                   sumRHS needs to be ~0 or it converge at a non-zero residual.
                0081 C     cg2d_min   :: used to store solution corresponding to lowest residual.
                0082 C     msgBuf     :: Informational/error message buffer
                0083       INTEGER bi, bj
                0084       INTEGER i, j, k, it3d
                0085       INTEGER actualIts(nSx,nSy)
                0086       INTEGER km1, kp1
                0087       _RL     maskM1, maskP1
                0088       _RL     err_sq,  errTile(nSx,nSy)
                0089       _RL     eta_qrNtile(nSx,nSy)
                0090       _RL     eta_qrNM1(nSx,nSy)
                0091       _RL     cgBeta
                0092       _RL     alpha ,  alphaTile(nSx,nSy)
                0093       _RL     sumRHS,  sumRHStile(nSx,nSy)
                0094       _RL     rhsMax, rhsMaxLoc
                0095       _RL     rhsNorm(nSx,nSy)
                0096       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0097       LOGICAL printResidual
                0098       _RL     surfFac
                0099 #ifdef NONLIN_FRSURF
                0100       INTEGER ks
                0101       _RL     surfTerm(sNx,sNy)
                0102 #endif /* NONLIN_FRSURF */
                0103 CEOP
                0104 
                0105 C--   Initialise auxiliary constant, some output variable
                0106       IF ( select_rStar .NE. 0 ) THEN
                0107         surfFac = freeSurfFac
                0108       ELSE
                0109         surfFac = 0.
                0110       ENDIF
                0111 #ifdef NONLIN_FRSURF
                0112       DO j=1,sNy
                0113        DO i=1,sNx
                0114          surfTerm(i,j) = 0.
                0115        ENDDO
                0116       ENDDO
                0117 #endif /* NONLIN_FRSURF */
                0118 
                0119 C--   Initialise inverter and   Normalise RHS
                0120       rhsMax = 0. _d 0
                0121       DO bj=myByLo(myThid),myByHi(myThid)
                0122        DO bi=myBxLo(myThid),myBxHi(myThid)
                0123 
                0124         actualIts(bi,bj) = 0
                0125         eta_qrNM1(bi,bj) =  1. _d 0
                0126         rhsMaxLoc = 0. _d 0
                0127         DO k=1,Nr
                0128          DO j=1,sNy
                0129           DO i=1,sNx
                0130            cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)*cg3dNorm
                0131      &                          * maskC(i,j,k,bi,bj)
                0132            rhsMaxLoc = MAX(ABS(cg3d_b(i,j,k,bi,bj)),rhsMaxLoc)
                0133           ENDDO
                0134          ENDDO
                0135         ENDDO
                0136         rhsNorm(bi,bj) = 1. _d 0
                0137         IF ( rhsMaxLoc .NE. 0. ) rhsNorm(bi,bj) = 1. _d 0 / rhsMaxLoc
                0138         DO k=1,Nr
                0139          DO j=1,sNy
                0140           DO i=1,sNx
                0141            cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)*rhsNorm(bi,bj)
                0142            cg3d_x(i,j,k,bi,bj) = cg3d_x(i,j,k,bi,bj)*rhsNorm(bi,bj)
                0143           ENDDO
                0144          ENDDO
                0145         ENDDO
                0146         rhsMax = MAX( rhsMaxLoc, rhsMax )
                0147 
                0148        ENDDO
                0149       ENDDO
                0150       _GLOBAL_MAX_RL( rhsMax, myThid )
                0151 
                0152 C--   Update overlaps
                0153       _EXCH_XYZ_RL( cg3d_x, myThid )
                0154 
                0155 C--   Initial residual calculation (with free-Surface term)
                0156       err_sq = 0.
                0157       sumRHS = 0.
                0158       DO bj=myByLo(myThid),myByHi(myThid)
                0159        DO bi=myBxLo(myThid),myBxHi(myThid)
                0160         errTile(bi,bj)    = 0. _d 0
                0161         sumRHStile(bi,bj) = 0. _d 0
                0162 #ifdef NONLIN_FRSURF
                0163         IF ( select_rStar .NE. 0 ) THEN
                0164           DO j=1,sNy
                0165            DO i=1,sNx
                0166              surfTerm(i,j) = 0.
                0167            ENDDO
                0168           ENDDO
                0169           DO k=1,Nr
                0170            DO j=1,sNy
                0171             DO i=1,sNx
                0172              surfTerm(i,j) = surfTerm(i,j)
                0173      &         +cg3d_x(i,j,k,bi,bj)*drF(k)*h0FacC(i,j,k,bi,bj)
                0174             ENDDO
                0175            ENDDO
                0176           ENDDO
                0177           DO j=1,sNy
                0178            DO i=1,sNx
                0179              ks = kSurfC(i,j,bi,bj)
                0180              surfTerm(i,j) = surfTerm(i,j)*cg3dNorm
                0181      &              *recip_Rcol(i,j,bi,bj)*recip_Rcol(i,j,bi,bj)
                0182      &              *rA(i,j,bi,bj)*deepFac2F(ks)
e6e223b277 Jean*0183      &              *recip_Bo(i,j,bi,bj)/deltaTMom/deltaTFreeSurf
015c19607c Jean*0184            ENDDO
                0185           ENDDO
                0186         ENDIF
                0187 #endif /* NONLIN_FRSURF */
                0188         DO k=1,Nr
                0189          km1 = MAX(k-1, 1 )
                0190          kp1 = MIN(k+1, Nr)
                0191          maskM1 = 1. _d 0
                0192          maskP1 = 1. _d 0
                0193          IF ( k .EQ. 1 ) maskM1 = 0. _d 0
                0194          IF ( k .EQ. Nr) maskP1 = 0. _d 0
                0195 #ifdef TARGET_NEC_SX
                0196 !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
                0197 #endif /* TARGET_NEC_SX */
                0198          DO j=1,sNy
                0199           DO i=1,sNx
                0200            cg3d_r(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
                0201      &     -( 0.
                0202      &       +aW3d( i, j, k, bi,bj)*cg3d_x(i-1,j, k, bi,bj)
                0203      &       +aW3d(i+1,j, k, bi,bj)*cg3d_x(i+1,j, k, bi,bj)
                0204      &       +aS3d( i, j, k, bi,bj)*cg3d_x( i,j-1,k, bi,bj)
                0205      &       +aS3d( i,j+1,k, bi,bj)*cg3d_x( i,j+1,k, bi,bj)
                0206      &       +aV3d( i, j, k, bi,bj)*cg3d_x( i, j,km1,bi,bj)*maskM1
                0207      &       +aV3d( i, j,kp1,bi,bj)*cg3d_x( i, j,kp1,bi,bj)*maskP1
                0208      &       +aC3d( i, j, k, bi,bj)*cg3d_x( i, j, k, bi,bj)
                0209 #ifdef NONLIN_FRSURF
                0210      &       -surfFac*surfTerm(i,j)*drF(k)*h0FacC(i,j,k,bi,bj)
                0211 #endif /* NONLIN_FRSURF */
                0212      &      )
                0213            errTile(bi,bj) = errTile(bi,bj)
                0214      &                    +cg3d_r(i,j,k,bi,bj)*cg3d_r(i,j,k,bi,bj)
                0215            sumRHStile(bi,bj) = sumRHStile(bi,bj)+cg3d_b(i,j,k,bi,bj)
                0216           ENDDO
                0217          ENDDO
                0218          DO j=0,sNy+1
                0219           DO i=0,sNx+1
                0220            cg3d_s(i,j,k,bi,bj) = 0.
                0221           ENDDO
                0222          ENDDO
                0223         ENDDO
                0224         err_sq = MAX( errTile(bi,bj), err_sq )
                0225         sumRHS = MAX( ABS(sumRHStile(bi,bj)), sumRHS )
                0226        ENDDO
                0227       ENDDO
                0228        CALL EXCH_S3D_RL( cg3d_r, Nr, myThid )
                0229       _GLOBAL_MAX_RL( err_sq, myThid )
                0230       _GLOBAL_MAX_RL( sumRHS, myThid )
                0231       IF ( debugLevel.GE.debLevC .AND. diagFreq.GT.0. ) THEN
                0232         CALL WRITE_FLD_S3D_RL(
                0233      I        'cg3d_r_I', 'I10', 1, Nr, cg3d_r, myIter, myThid )
                0234       ENDIF
                0235       firstResidual = SQRT(err_sq)
                0236 
                0237       printResidual = .FALSE.
                0238       IF ( debugLevel .GE. debLevZero ) THEN
                0239         _BEGIN_MASTER( myThid )
                0240         printResidual = printResidualFreq.GE.1
e6e223b277 Jean*0241         WRITE(standardMessageUnit,'(A,1P2E22.14)')
015c19607c Jean*0242      &     ' cg3d: Sum(rhs),rhsMax = ',sumRHS,rhsMax
                0243         _END_MASTER( myThid )
                0244       ENDIF
                0245 
                0246 c     IF ( err_sq .LT. cg3dTolerance_sq ) GOTO 11
                0247 
                0248 C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
                0249       DO it3d=1, numIters
                0250       IF ( err_sq.GE.cg3dTolerance_sq ) THEN
                0251        err_sq = 0. _d 0
                0252 
                0253        DO bj=myByLo(myThid),myByHi(myThid)
                0254        DO bi=myBxLo(myThid),myBxHi(myThid)
                0255         IF ( errTile(bi,bj).GE.cg3dTolerance_sq ) THEN
                0256 
                0257 C--    Solve preconditioning equation and update
                0258 C--    conjugate direction vector "s".
                0259 C      Note. On the next two loops over all tiles the inner loop ranges
                0260 C            in sNx and sNy are expanded by 1 to avoid a communication
                0261 C            step. However this entails a bit of gynamastics because we only
                0262 C            want eta_qrN for the interior points.
                0263          eta_qrNtile(bi,bj) = 0. _d 0
                0264          DO k=1,1
                0265 #ifdef TARGET_NEC_SX
                0266 !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
                0267 #endif /* TARGET_NEC_SX */
                0268           DO j=0,sNy+1
                0269            DO i=0,sNx+1
                0270             cg3d_q(i,j,k,bi,bj) = zMC(i,j,k,bi,bj)
                0271      &                        *cg3d_r(i,j,k,bi,bj)
                0272            ENDDO
                0273           ENDDO
                0274          ENDDO
                0275          DO k=2,Nr
                0276 #ifdef TARGET_NEC_SX
                0277 !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
                0278 #endif /* TARGET_NEC_SX */
                0279           DO j=0,sNy+1
                0280            DO i=0,sNx+1
                0281             cg3d_q(i,j,k,bi,bj) = zMC(i,j,k,bi,bj)
                0282      &                      *( cg3d_r(i,j,k,bi,bj)
                0283      &                         -zML(i,j,k,bi,bj)*cg3d_q(i,j,k-1,bi,bj)
                0284      &                       )
                0285            ENDDO
                0286           ENDDO
                0287          ENDDO
                0288          DO k=Nr,Nr
                0289 #ifdef TARGET_NEC_SX
                0290 !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
                0291 #endif /* TARGET_NEC_SX */
                0292           DO j=1,sNy
                0293            DO i=1,sNx
                0294             eta_qrNtile(bi,bj) = eta_qrNtile(bi,bj)
                0295      &                        +cg3d_q(i,j,k,bi,bj)*cg3d_r(i,j,k,bi,bj)
                0296            ENDDO
                0297           ENDDO
                0298          ENDDO
                0299          DO k=Nr-1,1,-1
                0300 #ifdef TARGET_NEC_SX
                0301 !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
                0302 #endif /* TARGET_NEC_SX */
                0303           DO j=0,sNy+1
                0304            DO i=0,sNx+1
                0305             cg3d_q(i,j,k,bi,bj) = cg3d_q(i,j,k,bi,bj)
                0306      &                         -zMU(i,j,k,bi,bj)*cg3d_q(i,j,k+1,bi,bj)
                0307            ENDDO
                0308           ENDDO
                0309 #ifdef TARGET_NEC_SX
                0310 !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
                0311 #endif /* TARGET_NEC_SX */
                0312           DO j=1,sNy
                0313            DO i=1,sNx
                0314             eta_qrNtile(bi,bj) = eta_qrNtile(bi,bj)
                0315      &                        +cg3d_q(i,j,k,bi,bj)*cg3d_r(i,j,k,bi,bj)
                0316            ENDDO
                0317           ENDDO
                0318          ENDDO
                0319 
                0320          cgBeta   = eta_qrNtile(bi,bj)/eta_qrNM1(bi,bj)
                0321          eta_qrNM1(bi,bj) = eta_qrNtile(bi,bj)
                0322 
                0323          DO k=1,Nr
                0324           DO j=0,sNy+1
                0325            DO i=0,sNx+1
                0326             cg3d_s(i,j,k,bi,bj) = cg3d_q(i,j,k,bi,bj)
                0327      &                   + cgBeta*cg3d_s(i,j,k,bi,bj)
                0328            ENDDO
                0329           ENDDO
                0330          ENDDO
                0331 
                0332 C==    Evaluate laplace operator on conjugate gradient vector
                0333 C==    q = A.s
                0334          alphaTile(bi,bj) = 0. _d 0
                0335 #ifdef NONLIN_FRSURF
                0336          IF ( select_rStar .NE. 0 ) THEN
                0337           DO j=1,sNy
                0338            DO i=1,sNx
                0339              surfTerm(i,j) = 0.
                0340            ENDDO
                0341           ENDDO
                0342           DO k=1,Nr
                0343            DO j=1,sNy
                0344             DO i=1,sNx
                0345              surfTerm(i,j) = surfTerm(i,j)
                0346      &         +cg3d_s(i,j,k,bi,bj)*drF(k)*h0FacC(i,j,k,bi,bj)
                0347             ENDDO
                0348            ENDDO
                0349           ENDDO
                0350           DO j=1,sNy
                0351            DO i=1,sNx
                0352              ks = kSurfC(i,j,bi,bj)
                0353              surfTerm(i,j) = surfTerm(i,j)*cg3dNorm
                0354      &              *recip_Rcol(i,j,bi,bj)*recip_Rcol(i,j,bi,bj)
                0355      &              *rA(i,j,bi,bj)*deepFac2F(ks)
e6e223b277 Jean*0356      &              *recip_Bo(i,j,bi,bj)/deltaTMom/deltaTFreeSurf
015c19607c Jean*0357            ENDDO
                0358           ENDDO
                0359          ENDIF
                0360 #endif /* NONLIN_FRSURF */
                0361          IF ( Nr .GT. 1 ) THEN
                0362           k=1
                0363 #ifdef TARGET_NEC_SX
                0364 !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
                0365 #endif /* TARGET_NEC_SX */
                0366           DO j=1,sNy
                0367            DO i=1,sNx
                0368             cg3d_q(i,j,k,bi,bj) =
                0369      &        aW3d( i, j, k, bi,bj)*cg3d_s(i-1,j, k, bi,bj)
                0370      &       +aW3d(i+1,j, k, bi,bj)*cg3d_s(i+1,j, k, bi,bj)
                0371      &       +aS3d( i, j, k, bi,bj)*cg3d_s( i,j-1,k, bi,bj)
                0372      &       +aS3d( i,j+1,k, bi,bj)*cg3d_s( i,j+1,k, bi,bj)
                0373      &       +aV3d( i, j,k+1,bi,bj)*cg3d_s( i, j,k+1,bi,bj)
                0374      &       +aC3d( i, j, k, bi,bj)*cg3d_s( i, j, k, bi,bj)
                0375 #ifdef NONLIN_FRSURF
                0376      &       -surfFac*surfTerm(i,j)*drF(k)*h0FacC(i,j,k,bi,bj)
                0377 #endif /* NONLIN_FRSURF */
                0378             alphaTile(bi,bj) = alphaTile(bi,bj)
                0379      &             +cg3d_s(i,j,k,bi,bj)*cg3d_q(i,j,k,bi,bj)
                0380            ENDDO
                0381           ENDDO
                0382          ELSE
                0383           k=1
                0384 #ifdef TARGET_NEC_SX
                0385 !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
                0386 #endif /* TARGET_NEC_SX */
                0387           DO j=1,sNy
                0388            DO i=1,sNx
                0389             cg3d_q(i,j,k,bi,bj) =
                0390      &        aW3d( i, j, k, bi,bj)*cg3d_s(i-1,j, k, bi,bj)
                0391      &       +aW3d(i+1,j, k, bi,bj)*cg3d_s(i+1,j, k, bi,bj)
                0392      &       +aS3d( i, j, k, bi,bj)*cg3d_s( i,j-1,k, bi,bj)
                0393      &       +aS3d( i,j+1,k, bi,bj)*cg3d_s( i,j+1,k, bi,bj)
                0394      &       +aC3d( i, j, k, bi,bj)*cg3d_s( i, j, k, bi,bj)
                0395 #ifdef NONLIN_FRSURF
                0396      &       -surfFac*surfTerm(i,j)*drF(k)*h0FacC(i,j,k,bi,bj)
                0397 #endif /* NONLIN_FRSURF */
                0398             alphaTile(bi,bj) = alphaTile(bi,bj)
                0399      &             +cg3d_s(i,j,k,bi,bj)*cg3d_q(i,j,k,bi,bj)
                0400            ENDDO
                0401           ENDDO
                0402          ENDIF
                0403          DO k=2,Nr-1
                0404 #ifdef TARGET_NEC_SX
                0405 !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
                0406 #endif /* TARGET_NEC_SX */
                0407           DO j=1,sNy
                0408            DO i=1,sNx
                0409             cg3d_q(i,j,k,bi,bj) =
                0410      &        aW3d( i, j, k, bi,bj)*cg3d_s(i-1,j, k, bi,bj)
                0411      &       +aW3d(i+1,j, k, bi,bj)*cg3d_s(i+1,j, k, bi,bj)
                0412      &       +aS3d( i, j, k, bi,bj)*cg3d_s( i,j-1,k, bi,bj)
                0413      &       +aS3d( i,j+1,k, bi,bj)*cg3d_s( i,j+1,k, bi,bj)
                0414      &       +aV3d( i, j, k, bi,bj)*cg3d_s( i, j,k-1,bi,bj)
                0415      &       +aV3d( i, j,k+1,bi,bj)*cg3d_s( i, j,k+1,bi,bj)
                0416      &       +aC3d( i, j, k, bi,bj)*cg3d_s( i, j, k, bi,bj)
                0417 #ifdef NONLIN_FRSURF
                0418      &       -surfFac*surfTerm(i,j)*drF(k)*h0FacC(i,j,k,bi,bj)
                0419 #endif /* NONLIN_FRSURF */
                0420             alphaTile(bi,bj) = alphaTile(bi,bj)
                0421      &             +cg3d_s(i,j,k,bi,bj)*cg3d_q(i,j,k,bi,bj)
                0422            ENDDO
                0423           ENDDO
                0424          ENDDO
                0425          IF ( Nr .GT. 1 ) THEN
                0426           k=Nr
                0427 #ifdef TARGET_NEC_SX
                0428 !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
                0429 #endif /* TARGET_NEC_SX */
                0430           DO j=1,sNy
                0431            DO i=1,sNx
                0432             cg3d_q(i,j,k,bi,bj) =
                0433      &        aW3d( i, j, k, bi,bj)*cg3d_s(i-1,j, k, bi,bj)
                0434      &       +aW3d(i+1,j, k, bi,bj)*cg3d_s(i+1,j, k, bi,bj)
                0435      &       +aS3d( i, j, k, bi,bj)*cg3d_s( i,j-1,k, bi,bj)
                0436      &       +aS3d( i,j+1,k, bi,bj)*cg3d_s( i,j+1,k, bi,bj)
                0437      &       +aV3d( i, j, k, bi,bj)*cg3d_s( i, j,k-1,bi,bj)
                0438      &       +aC3d( i, j, k, bi,bj)*cg3d_s( i, j, k, bi,bj)
                0439 #ifdef NONLIN_FRSURF
                0440      &       -surfFac*surfTerm(i,j)*drF(k)*h0FacC(i,j,k,bi,bj)
                0441 #endif /* NONLIN_FRSURF */
                0442             alphaTile(bi,bj) = alphaTile(bi,bj)
                0443      &             +cg3d_s(i,j,k,bi,bj)*cg3d_q(i,j,k,bi,bj)
                0444            ENDDO
                0445           ENDDO
                0446          ENDIF
                0447          alpha = eta_qrNtile(bi,bj)/alphaTile(bi,bj)
                0448 
                0449 C==    Update simultaneously solution and residual vectors (and Iter number)
                0450 C      Now compute "interior" points.
                0451          errTile(bi,bj) = 0. _d 0
                0452          DO k=1,Nr
                0453 #ifdef TARGET_NEC_SX
                0454 !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
                0455 #endif /* TARGET_NEC_SX */
                0456           DO j=1,sNy
                0457            DO i=1,sNx
                0458             cg3d_x(i,j,k,bi,bj)=cg3d_x(i,j,k,bi,bj)
                0459      &            +alpha*cg3d_s(i,j,k,bi,bj)
                0460             cg3d_r(i,j,k,bi,bj)=cg3d_r(i,j,k,bi,bj)
                0461      &            -alpha*cg3d_q(i,j,k,bi,bj)
                0462             errTile(bi,bj) = errTile(bi,bj)
                0463      &            +cg3d_r(i,j,k,bi,bj)*cg3d_r(i,j,k,bi,bj)
                0464            ENDDO
                0465           ENDDO
                0466          ENDDO
                0467          actualIts(bi,bj) = it3d
                0468 
                0469          IF ( printResidual ) THEN
                0470           IF ( MOD( it3d-1, printResidualFreq ).EQ.0 ) THEN
                0471            WRITE(msgBuf,'(A,2I4,A,I6,A,1PE21.14)') ' cg3d(bi,bj=', bi,
                0472      &      bj, '): iter=', it3d, ' ; resid.= ', SQRT(errTile(bi,bj))
                0473            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0474      &                         SQUEEZE_RIGHT, myThid )
                0475           ENDIF
                0476          ENDIF
                0477 
                0478          CALL FILL_HALO_LOCAL_RL(
                0479      U                  cg3d_r(0,0,1,bi,bj),
                0480      I                  1, 1, 1, 1, Nr,
                0481      I                  EXCH_IGNORE_CORNERS, bi, bj, myThid )
                0482 
                0483         ENDIF
                0484         err_sq = MAX( errTile(bi,bj), err_sq )
                0485 C-    end bi,bj loops
                0486        ENDDO
                0487        ENDDO
                0488 C-    end cg-2d iteration loop
                0489       ENDIF
                0490       ENDDO
                0491 
                0492 c  11 CONTINUE
                0493 
                0494       IF ( debugLevel.GE.debLevC .AND. diagFreq.GT.0. ) THEN
                0495         CALL WRITE_FLD_S3D_RL(
                0496      I        'cg3d_r_F', 'I10', 1, Nr, cg3d_r, myIter, myThid )
                0497       ENDIF
                0498 
                0499 C--   Un-normalise the answer
                0500       numIters = 0
                0501       DO bj=myByLo(myThid),myByHi(myThid)
                0502        DO bi=myBxLo(myThid),myBxHi(myThid)
                0503         DO k=1,Nr
                0504          DO j=1,sNy
                0505           DO i=1,sNx
                0506            cg3d_x(i,j,k,bi,bj) = cg3d_x(i,j,k,bi,bj)/rhsNorm(bi,bj)
                0507           ENDDO
                0508          ENDDO
                0509         ENDDO
                0510         numIters = MAX( actualIts(bi,bj), numIters )
                0511        ENDDO
                0512       ENDDO
                0513 
                0514 C--   Return parameters to caller
                0515 C     return largest Iter # and Max residual     in numIters and lastResidual
                0516       _GLOBAL_MAX_RL( err_sq, myThid )
                0517       lastResidual = SQRT(err_sq)
                0518       alpha = numIters
                0519       _GLOBAL_MAX_RL( alpha, myThid )
                0520       numIters = NINT( alpha )
                0521 
                0522 #endif /* ALLOW_NONHYDROSTATIC */
                0523 
                0524       RETURN
                0525       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!