model/src/cg2d.F

1dbaea09ee Chri*0001 #include "CPP_OPTIONS.h"
2b7cba8d4b Mart*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 CG2D_OUTERLOOPITERS
                0006 # define CG2D_OUTERLOOPITERS 10
                0007 #endif
                0008 #endif /* TARGET_NEC_SX */
924557e60a Chri*0009
9366854e02 Chri*0010 CBOP
                0011 C     !ROUTINE: CG2D
                0012 C     !INTERFACE:
4606c28752 Jean*0013       SUBROUTINE CG2D(
2739a7f265 Jean*0014      U                cg2d_b, cg2d_x,
                0015      O                firstResidual, minResidualSq, lastResidual,
                0016      U                numIters, nIterMin,
924557e60a Chri*0017      I                myThid )
9366854e02 Chri*0018 C     !DESCRIPTION: \bv
aecc8b0f47 Mart*0019 C     *================================================================*
4606c28752 Jean*0020 C     | SUBROUTINE CG2D
aecc8b0f47 Mart*0021 C     | o Two-dimensional grid problem conjugate-gradient inverter
                0022 C     |   (with preconditioner).
                0023 C     *================================================================*
4606c28752 Jean*0024 C     | Con. grad is an iterative procedure for solving Ax = b.
                0025 C     | It requires the A be symmetric.
aecc8b0f47 Mart*0026 C     | This implementation assumes A is a five-diagonal matrix
                0027 C     | of the form that arises in the discrete representation of
                0028 C     | the del^2 operator in a two-dimensional space.
                0029 C     *================================================================*
9366854e02 Chri*0030 C     \ev
924557e60a Chri*0031
9366854e02 Chri*0032 C     !USES:
                0033       IMPLICIT NONE
924557e60a Chri*0034 C     === Global data ===
                0035 #include "SIZE.h"
                0036 #include "EEPARAMS.h"
                0037 #include "PARAMS.h"
aea29c8517 Alis*0038 #include "CG2D.h"
924557e60a Chri*0039
9366854e02 Chri*0040 C     !INPUT/OUTPUT PARAMETERS:
2739a7f265 Jean*0041 C     cg2d_b    :: The source term or "right hand side" (output: normalised RHS)
                0042 C     cg2d_x    :: The solution (input: first guess)
4606c28752 Jean*0043 C     firstResidual :: the initial residual before any iterations
2739a7f265 Jean*0044 C     minResidualSq :: the lowest residual reached (squared)
4606c28752 Jean*0045 C     lastResidual  :: the actual residual reached
2739a7f265 Jean*0046 C     numIters  :: Inp: the maximum number of iterations allowed
                0047 C                  Out: the actual number of iterations used
                0048 C     nIterMin  :: Inp: decide to store (if >=0) or not (if <0) lowest res. sol.
                0049 C                  Out: iteration number corresponding to lowest residual
9155c64ca4 Jean*0050 C     myThid    :: Thread on which I am working.
46dc4f419b Chri*0051       _RL  cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
030bea3287 Alis*0052       _RL  cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
aea29c8517 Alis*0053       _RL  firstResidual
2739a7f265 Jean*0054       _RL  minResidualSq
aea29c8517 Alis*0055       _RL  lastResidual
030bea3287 Alis*0056       INTEGER numIters
2739a7f265 Jean*0057       INTEGER nIterMin
030bea3287 Alis*0058       INTEGER myThid
924557e60a Chri*0059
9366854e02 Chri*0060 C     !LOCAL VARIABLES:
2739a7f265 Jean*0061 C     bi, bj     :: tile index in X and Y.
                0062 C     i, j, it2d :: Loop counters ( it2d counts CG iterations )
                0063 C     actualIts  :: actual CG iteration number
                0064 C     err_sq     :: Measure of the square of the residual of Ax - b.
                0065 C     eta_qrN    :: Used in computing search directions; suffix N and NM1
                0066 C     eta_qrNM1     denote current and previous iterations respectively.
                0067 C     cgBeta     :: coeff used to update conjugate direction vector "s".
                0068 C     alpha      :: coeff used to update solution & residual
                0069 C     sumRHS     :: Sum of right-hand-side. Sometimes this is a useful
                0070 C                   debugging/trouble shooting diagnostic. For neumann problems
                0071 C                   sumRHS needs to be ~0 or it converge at a non-zero residual.
                0072 C     cg2d_min   :: used to store solution corresponding to lowest residual.
                0073 C     msgBuf     :: Informational/error message buffer
aecc8b0f47 Mart*0074 C--   local working array (used to be in CG2D.h common block:
                0075 C     cg2d_q     :: Intermediate matrix-vector product term
                0076 C     cg2d_r     ::   *same*
                0077 C     cg2d_s     ::   *same*
4606c28752 Jean*0078       INTEGER bi, bj
2739a7f265 Jean*0079       INTEGER i, j, it2d
                0080       INTEGER actualIts
                0081       _RL    err_sq,  errTile(nSx,nSy)
                0082       _RL    eta_qrN, eta_qrNtile(nSx,nSy)
980a7d9757 Jean*0083       _RL    eta_qrNM1
46dc4f419b Chri*0084       _RL    cgBeta
2739a7f265 Jean*0085       _RL    alpha,   alphaTile(nSx,nSy)
                0086       _RL    sumRHS,  sumRHStile(nSx,nSy)
46dc4f419b Chri*0087       _RL    rhsMax
                0088       _RL    rhsNorm
2739a7f265 Jean*0089       _RL    cg2d_min(1:sNx,1:sNy,nSx,nSy)
d68b36a485 Mart*0090       _RL    cg2d_q  (1:sNx,1:sNy,nSx,nSy)
aecc8b0f47 Mart*0091       _RL    cg2d_r(1-1:sNx+1,1-1:sNy+1,nSx,nSy)
                0092       _RL    cg2d_s(1-1:sNx+1,1-1:sNy+1,nSx,nSy)
d439db03c9 Oliv*0093 #ifdef CG2D_SINGLECPU_SUM
                0094       _RL    localBuf(1:sNx,1:sNy,nSx,nSy)
                0095 #endif
1f36709e5c Jean*0096       CHARACTER*(MAX_LEN_MBUF) msgBuf
764aea3440 Jean*0097       LOGICAL printResidual
9366854e02 Chri*0098 CEOP
a85d6ab24e Chri*0099
2739a7f265 Jean*0100 C--   Initialise auxiliary constant, some output variable and inverter
                0101       minResidualSq = -1. _d 0
                0102       eta_qrNM1     =  1. _d 0
dbef787bda Chri*0103
924557e60a Chri*0104 C--   Normalise RHS
                0105       rhsMax = 0. _d 0
                0106       DO bj=myByLo(myThid),myByHi(myThid)
                0107        DO bi=myBxLo(myThid),myBxHi(myThid)
2739a7f265 Jean*0108         DO j=1,sNy
                0109          DO i=1,sNx
                0110           cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)*cg2dNorm
                0111           rhsMax = MAX(ABS(cg2d_b(i,j,bi,bj)),rhsMax)
924557e60a Chri*0112          ENDDO
                0113         ENDDO
                0114        ENDDO
                0115       ENDDO
aea29c8517 Alis*0116
                0117       IF (cg2dNormaliseRHS) THEN
                0118 C-  Normalise RHS :
e6e223b277 Jean*0119        _GLOBAL_MAX_RL( rhsMax, myThid )
                0120        rhsNorm = 1. _d 0
                0121        IF ( rhsMax .NE. 0. ) rhsNorm = 1. _d 0 / rhsMax
                0122        DO bj=myByLo(myThid),myByHi(myThid)
                0123         DO bi=myBxLo(myThid),myBxHi(myThid)
                0124          DO j=1,sNy
                0125           DO i=1,sNx
                0126            cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)*rhsNorm
                0127            cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)*rhsNorm
                0128           ENDDO
924557e60a Chri*0129          ENDDO
                0130         ENDDO
                0131        ENDDO
aea29c8517 Alis*0132 C- end Normalise RHS
                0133       ENDIF
924557e60a Chri*0134
                0135 C--   Update overlaps
9155c64ca4 Jean*0136       CALL EXCH_XY_RL( cg2d_x, myThid )
924557e60a Chri*0137
                0138 C--   Initial residual calculation
                0139       DO bj=myByLo(myThid),myByHi(myThid)
                0140        DO bi=myBxLo(myThid),myBxHi(myThid)
ff51fd05c6 Jean*0141 C-    Initialise local working arrays:
                0142         DO j=0,sNy+1
                0143          DO i=0,sNx+1
                0144           cg2d_r(i,j,bi,bj) = 0. _d 0
                0145           cg2d_s(i,j,bi,bj) = 0. _d 0
                0146          ENDDO
                0147         ENDDO
2739a7f265 Jean*0148         IF ( nIterMin.GE.0 ) THEN
ff51fd05c6 Jean*0149 C-    Initialise saved solution
2739a7f265 Jean*0150          DO j=1,sNy
                0151           DO i=1,sNx
ff51fd05c6 Jean*0152            cg2d_min(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
2739a7f265 Jean*0153           ENDDO
                0154          ENDDO
                0155         ENDIF
a619b19bc9 Jean*0156         sumRHStile(bi,bj) = 0. _d 0
                0157         errTile(bi,bj)    = 0. _d 0
2b7cba8d4b Mart*0158 #ifdef TARGET_NEC_SX
                0159 !CDIR OUTERUNROLL=CG2D_OUTERLOOPITERS
                0160 #endif /* TARGET_NEC_SX */
2739a7f265 Jean*0161         DO j=1,sNy
                0162          DO i=1,sNx
                0163           cg2d_r(i,j,bi,bj) = cg2d_b(i,j,bi,bj) -
                0164      &    (aW2d(i  ,j  ,bi,bj)*cg2d_x(i-1,j  ,bi,bj)
                0165      &    +aW2d(i+1,j  ,bi,bj)*cg2d_x(i+1,j  ,bi,bj)
                0166      &    +aS2d(i  ,j  ,bi,bj)*cg2d_x(i  ,j-1,bi,bj)
                0167      &    +aS2d(i  ,j+1,bi,bj)*cg2d_x(i  ,j+1,bi,bj)
                0168      &    +aC2d(i  ,j  ,bi,bj)*cg2d_x(i  ,j  ,bi,bj)
c0a4efc370 Chri*0169      &    )
d439db03c9 Oliv*0170 #ifdef CG2D_SINGLECPU_SUM
2739a7f265 Jean*0171           localBuf(i,j,bi,bj) = cg2d_r(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
d439db03c9 Oliv*0172 #else
a619b19bc9 Jean*0173           errTile(bi,bj)    = errTile(bi,bj)
2739a7f265 Jean*0174      &                      + cg2d_r(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
d439db03c9 Oliv*0175 #endif
d68b36a485 Mart*0176           sumRHStile(bi,bj) = sumRHStile(bi,bj) + cg2d_b(i,j,bi,bj)
924557e60a Chri*0177          ENDDO
                0178         ENDDO
                0179        ENDDO
                0180       ENDDO
9155c64ca4 Jean*0181       CALL EXCH_S3D_RL( cg2d_r, 1, myThid )
d439db03c9 Oliv*0182 #ifdef CG2D_SINGLECPU_SUM
2739a7f265 Jean*0183       CALL GLOBAL_SUM_SINGLECPU_RL(localBuf, err_sq, 0, 0, myThid)
d439db03c9 Oliv*0184 #else
2739a7f265 Jean*0185       CALL GLOBAL_SUM_TILE_RL( errTile,    err_sq, myThid )
d439db03c9 Oliv*0186 #endif
d68b36a485 Mart*0187       CALL GLOBAL_SUM_TILE_RL( sumRHStile, sumRHS, myThid )
2739a7f265 Jean*0188       actualIts = 0
                0189       firstResidual = SQRT(err_sq)
                0190       IF ( nIterMin.GE.0 ) THEN
                0191         nIterMin = 0
                0192         minResidualSq = err_sq
                0193       ENDIF
d607aa1d50 Dimi*0194
764aea3440 Jean*0195       printResidual = .FALSE.
9155c64ca4 Jean*0196       IF ( debugLevel .GE. debLevZero ) THEN
494ad43bae Patr*0197         _BEGIN_MASTER( myThid )
764aea3440 Jean*0198         printResidual = printResidualFreq.GE.1
e6e223b277 Jean*0199         WRITE(standardMessageUnit,'(A,1P2E22.14)')
4606c28752 Jean*0200      &  ' cg2d: Sum(rhs),rhsMax = ', sumRHS,rhsMax
ea6e02f692 Ed H*0201         _END_MASTER( myThid )
9155c64ca4 Jean*0202       ENDIF
aea29c8517 Alis*0203
2739a7f265 Jean*0204       IF ( err_sq .LT. cg2dTolerance_sq ) GOTO 11
924557e60a Chri*0205
                0206 C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
030bea3287 Alis*0207       DO 10 it2d=1, numIters
924557e60a Chri*0208
                0209 C--    Solve preconditioning equation and update
                0210 C--    conjugate direction vector "s".
                0211        DO bj=myByLo(myThid),myByHi(myThid)
                0212         DO bi=myBxLo(myThid),myBxHi(myThid)
a619b19bc9 Jean*0213          eta_qrNtile(bi,bj) = 0. _d 0
2b7cba8d4b Mart*0214 #ifdef TARGET_NEC_SX
                0215 !CDIR OUTERUNROLL=CG2D_OUTERLOOPITERS
                0216 #endif /* TARGET_NEC_SX */
2739a7f265 Jean*0217          DO j=1,sNy
                0218           DO i=1,sNx
                0219            cg2d_q(i,j,bi,bj) =
                0220      &      pC(i  ,j  ,bi,bj)*cg2d_r(i  ,j  ,bi,bj)
                0221      &     +pW(i  ,j  ,bi,bj)*cg2d_r(i-1,j  ,bi,bj)
                0222      &     +pW(i+1,j  ,bi,bj)*cg2d_r(i+1,j  ,bi,bj)
                0223      &     +pS(i  ,j  ,bi,bj)*cg2d_r(i  ,j-1,bi,bj)
                0224      &     +pS(i  ,j+1,bi,bj)*cg2d_r(i  ,j+1,bi,bj)
c0a4efc370 Chri*0225 CcnhDebugStarts
2739a7f265 Jean*0226 c          cg2d_q(i,j,bi,bj) = cg2d_r(j  ,j  ,bi,bj)
c0a4efc370 Chri*0227 CcnhDebugEnds
d439db03c9 Oliv*0228 #ifdef CG2D_SINGLECPU_SUM
2739a7f265 Jean*0229           localBuf(i,j,bi,bj) =
                0230      &      cg2d_q(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
d439db03c9 Oliv*0231 #else
a619b19bc9 Jean*0232            eta_qrNtile(bi,bj) = eta_qrNtile(bi,bj)
2739a7f265 Jean*0233      &     +cg2d_q(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
d439db03c9 Oliv*0234 #endif
924557e60a Chri*0235           ENDDO
                0236          ENDDO
                0237         ENDDO
                0238        ENDDO
                0239
d439db03c9 Oliv*0240 #ifdef CG2D_SINGLECPU_SUM
71d7d0e7f1 Jean*0241        CALL GLOBAL_SUM_SINGLECPU_RL( localBuf,eta_qrN,0,0,myThid )
d439db03c9 Oliv*0242 #else
a619b19bc9 Jean*0243        CALL GLOBAL_SUM_TILE_RL( eta_qrNtile,eta_qrN,myThid )
d439db03c9 Oliv*0244 #endif
980a7d9757 Jean*0245        cgBeta   = eta_qrN/eta_qrNM1
924557e60a Chri*0246 CcnhDebugStarts
2739a7f265 Jean*0247 c      WRITE(*,*) ' CG2D: Iteration ', it2d-1,
                0248 c    &            ' eta_qrN=', eta_qrN, ' beta=', cgBeta
924557e60a Chri*0249 CcnhDebugEnds
980a7d9757 Jean*0250        eta_qrNM1 = eta_qrN
924557e60a Chri*0251
                0252        DO bj=myByLo(myThid),myByHi(myThid)
                0253         DO bi=myBxLo(myThid),myBxHi(myThid)
2739a7f265 Jean*0254          DO j=1,sNy
                0255           DO i=1,sNx
                0256            cg2d_s(i,j,bi,bj) = cg2d_q(i,j,bi,bj)
                0257      &                       + cgBeta*cg2d_s(i,j,bi,bj)
924557e60a Chri*0258           ENDDO
                0259          ENDDO
                0260         ENDDO
                0261        ENDDO
                0262
764aea3440 Jean*0263 C--    Do exchanges that require messages i.e. between processes.
9155c64ca4 Jean*0264        CALL EXCH_S3D_RL( cg2d_s, 1, myThid )
924557e60a Chri*0265
                0266 C==    Evaluate laplace operator on conjugate gradient vector
                0267 C==    q = A.s
                0268        DO bj=myByLo(myThid),myByHi(myThid)
                0269         DO bi=myBxLo(myThid),myBxHi(myThid)
a619b19bc9 Jean*0270          alphaTile(bi,bj) = 0. _d 0
2b7cba8d4b Mart*0271 #ifdef TARGET_NEC_SX
                0272 !CDIR OUTERUNROLL=CG2D_OUTERLOOPITERS
                0273 #endif /* TARGET_NEC_SX */
2739a7f265 Jean*0274          DO j=1,sNy
                0275           DO i=1,sNx
                0276            cg2d_q(i,j,bi,bj) =
                0277      &     aW2d(i  ,j  ,bi,bj)*cg2d_s(i-1,j  ,bi,bj)
                0278      &    +aW2d(i+1,j  ,bi,bj)*cg2d_s(i+1,j  ,bi,bj)
                0279      &    +aS2d(i  ,j  ,bi,bj)*cg2d_s(i  ,j-1,bi,bj)
                0280      &    +aS2d(i  ,j+1,bi,bj)*cg2d_s(i  ,j+1,bi,bj)
                0281      &    +aC2d(i  ,j  ,bi,bj)*cg2d_s(i  ,j  ,bi,bj)
d439db03c9 Oliv*0282 #ifdef CG2D_SINGLECPU_SUM
2739a7f265 Jean*0283           localBuf(i,j,bi,bj) = cg2d_s(i,j,bi,bj)*cg2d_q(i,j,bi,bj)
d439db03c9 Oliv*0284 #else
a619b19bc9 Jean*0285           alphaTile(bi,bj) = alphaTile(bi,bj)
2739a7f265 Jean*0286      &                     + cg2d_s(i,j,bi,bj)*cg2d_q(i,j,bi,bj)
d439db03c9 Oliv*0287 #endif
924557e60a Chri*0288           ENDDO
                0289          ENDDO
                0290         ENDDO
                0291        ENDDO
d439db03c9 Oliv*0292 #ifdef CG2D_SINGLECPU_SUM
71d7d0e7f1 Jean*0293        CALL GLOBAL_SUM_SINGLECPU_RL(localBuf, alpha, 0, 0, myThid)
d439db03c9 Oliv*0294 #else
a619b19bc9 Jean*0295        CALL GLOBAL_SUM_TILE_RL( alphaTile,  alpha,  myThid )
d439db03c9 Oliv*0296 #endif
924557e60a Chri*0297 CcnhDebugStarts
2739a7f265 Jean*0298 c      WRITE(*,*) ' CG2D: Iteration ', it2d-1,
                0299 c    &            ' SUM(s*q)=', alpha, ' alpha=', eta_qrN/alpha
924557e60a Chri*0300 CcnhDebugEnds
980a7d9757 Jean*0301        alpha = eta_qrN/alpha
4606c28752 Jean*0302
2739a7f265 Jean*0303 C==    Update simultaneously solution and residual vectors (and Iter number)
924557e60a Chri*0304 C      Now compute "interior" points.
                0305        DO bj=myByLo(myThid),myByHi(myThid)
                0306         DO bi=myBxLo(myThid),myBxHi(myThid)
a619b19bc9 Jean*0307          errTile(bi,bj) = 0. _d 0
2739a7f265 Jean*0308          DO j=1,sNy
                0309           DO i=1,sNx
                0310            cg2d_x(i,j,bi,bj)=cg2d_x(i,j,bi,bj)+alpha*cg2d_s(i,j,bi,bj)
                0311            cg2d_r(i,j,bi,bj)=cg2d_r(i,j,bi,bj)-alpha*cg2d_q(i,j,bi,bj)
d439db03c9 Oliv*0312 #ifdef CG2D_SINGLECPU_SUM
2739a7f265 Jean*0313            localBuf(i,j,bi,bj) = cg2d_r(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
d439db03c9 Oliv*0314 #else
a619b19bc9 Jean*0315            errTile(bi,bj) = errTile(bi,bj)
2739a7f265 Jean*0316      &                    + cg2d_r(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
d439db03c9 Oliv*0317 #endif
924557e60a Chri*0318           ENDDO
                0319          ENDDO
                0320         ENDDO
                0321        ENDDO
2739a7f265 Jean*0322        actualIts = it2d
924557e60a Chri*0323
d439db03c9 Oliv*0324 #ifdef CG2D_SINGLECPU_SUM
2739a7f265 Jean*0325        CALL GLOBAL_SUM_SINGLECPU_RL(localBuf, err_sq, 0, 0, myThid)
d439db03c9 Oliv*0326 #else
2739a7f265 Jean*0327        CALL GLOBAL_SUM_TILE_RL( errTile,    err_sq,    myThid )
d439db03c9 Oliv*0328 #endif
764aea3440 Jean*0329        IF ( printResidual ) THEN
                0330         IF ( MOD( it2d-1, printResidualFreq ).EQ.0 ) THEN
                0331          WRITE(msgBuf,'(A,I6,A,1PE21.14)')
2739a7f265 Jean*0332      &    ' cg2d: iter=', it2d, ' ; resid.= ', SQRT(err_sq)
764aea3440 Jean*0333          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0334      &                       SQUEEZE_RIGHT, myThid )
                0335         ENDIF
1f36709e5c Jean*0336        ENDIF
2739a7f265 Jean*0337        IF ( err_sq .LT. cg2dTolerance_sq ) GOTO 11
                0338        IF ( err_sq .LT. minResidualSq ) THEN
                0339 C-     Store lowest residual solution
                0340          minResidualSq = err_sq
                0341          nIterMin = it2d
                0342          DO bj=myByLo(myThid),myByHi(myThid)
                0343           DO bi=myBxLo(myThid),myBxHi(myThid)
                0344            DO j=1,sNy
                0345             DO i=1,sNx
                0346              cg2d_min(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
                0347             ENDDO
                0348            ENDDO
                0349           ENDDO
                0350          ENDDO
                0351        ENDIF
1f36709e5c Jean*0352
9155c64ca4 Jean*0353        CALL EXCH_S3D_RL( cg2d_r, 1, myThid )
a85d6ab24e Chri*0354
924557e60a Chri*0355    10 CONTINUE
                0356    11 CONTINUE
                0357
2739a7f265 Jean*0358       IF ( nIterMin.GE.0 .AND. err_sq .GT. minResidualSq ) THEN
                0359 C-    use the lowest residual solution (instead of current one = last residual)
                0360         DO bj=myByLo(myThid),myByHi(myThid)
                0361          DO bi=myBxLo(myThid),myBxHi(myThid)
                0362           DO j=1,sNy
                0363            DO i=1,sNx
                0364              cg2d_x(i,j,bi,bj) = cg2d_min(i,j,bi,bj)
                0365            ENDDO
                0366           ENDDO
                0367          ENDDO
                0368         ENDDO
                0369       ENDIF
                0370
aea29c8517 Alis*0371       IF (cg2dNormaliseRHS) THEN
924557e60a Chri*0372 C--   Un-normalise the answer
e6e223b277 Jean*0373 c       rhsMax = 1. _d 0 / rhsNorm
aea29c8517 Alis*0374         DO bj=myByLo(myThid),myByHi(myThid)
                0375          DO bi=myBxLo(myThid),myBxHi(myThid)
2739a7f265 Jean*0376           DO j=1,sNy
                0377            DO i=1,sNx
                0378             cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)/rhsNorm
e6e223b277 Jean*0379 c           cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)*rhsMax
aea29c8517 Alis*0380            ENDDO
                0381           ENDDO
924557e60a Chri*0382          ENDDO
                0383         ENDDO
aea29c8517 Alis*0384       ENDIF
924557e60a Chri*0385
030bea3287 Alis*0386 C--   Return parameters to caller
2739a7f265 Jean*0387       lastResidual = SQRT(err_sq)
764aea3440 Jean*0388       numIters = actualIts
924557e60a Chri*0389
                0390 CcnhDebugStarts
2739a7f265 Jean*0391 c     _EXCH_XY_RL(cg2d_x, myThid )
                0392 c     CALL PLOT_FIELD_XYRL( cg2d_x, 'CALC_MOM_RHS CG2D_X' , 1, myThid )
                0393 c     err_sq = 0. _d 0
                0394 c     DO bj=myByLo(myThid),myByHi(myThid)
                0395 c      DO bi=myBxLo(myThid),myBxHi(myThid)
                0396 c       DO j=1,sNy
                0397 c        DO i=1,sNx
                0398 c         cg2d_r(i,j,bi,bj) = cg2d_b(i,j,bi,bj) -
                0399 c    &    (aW2d(i  ,j  ,bi,bj)*cg2d_x(i-1,j  ,bi,bj)
                0400 c    &    +aW2d(i+1,j  ,bi,bj)*cg2d_x(i+1,j  ,bi,bj)
                0401 c    &    +aS2d(i  ,j  ,bi,bj)*cg2d_x(i  ,j-1,bi,bj)
                0402 c    &    +aS2d(i  ,j+1,bi,bj)*cg2d_x(i  ,j+1,bi,bj)
                0403 c    &    +aC2d(i  ,j  ,bi,bj)*cg2d_x(i  ,j  ,bi,bj)
                0404 c    &    )
                0405 c         err_sq = err_sq + cg2d_r(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
                0406 c        ENDDO
                0407 c       ENDDO
                0408 c      ENDDO
                0409 c     ENDDO
                0410 c     _GLOBAL_SUM_RL( err_sq, myThid )
                0411 c     write(*,*) 'cg2d: Ax - b = ',SQRT(err_sq)
924557e60a Chri*0412 CcnhDebugEnds
                0413
88830be691 Alis*0414       RETURN
924557e60a Chri*0415       END