Back to home page

MITgcm
 
 

    

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

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