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    cg2dTolerance_sq
                0082       _RL    err_sq,  errTile(nSx,nSy)
                0083       _RL    eta_qrN, eta_qrNtile(nSx,nSy)
980a7d9757 Jean*0084       _RL    eta_qrNM1
46dc4f419b Chri*0085       _RL    cgBeta
2739a7f265 Jean*0086       _RL    alpha,   alphaTile(nSx,nSy)
                0087       _RL    sumRHS,  sumRHStile(nSx,nSy)
46dc4f419b Chri*0088       _RL    rhsMax
                0089       _RL    rhsNorm
2739a7f265 Jean*0090       _RL    cg2d_min(1:sNx,1:sNy,nSx,nSy)
d68b36a485 Mart*0091       _RL    cg2d_q  (1:sNx,1:sNy,nSx,nSy)
aecc8b0f47 Mart*0092       _RL    cg2d_r(1-1:sNx+1,1-1:sNy+1,nSx,nSy)
                0093       _RL    cg2d_s(1-1:sNx+1,1-1:sNy+1,nSx,nSy)
d439db03c9 Oliv*0094 #ifdef CG2D_SINGLECPU_SUM
                0095       _RL    localBuf(1:sNx,1:sNy,nSx,nSy)
                0096 #endif
1f36709e5c Jean*0097       CHARACTER*(MAX_LEN_MBUF) msgBuf
764aea3440 Jean*0098       LOGICAL printResidual
9366854e02 Chri*0099 CEOP
a85d6ab24e Chri*0100
2739a7f265 Jean*0101 C--   Initialise auxiliary constant, some output variable and inverter
                0102       cg2dTolerance_sq = cg2dTolerance*cg2dTolerance
                0103       minResidualSq = -1. _d 0
                0104       eta_qrNM1     =  1. _d 0
dbef787bda Chri*0105
924557e60a Chri*0106 C--   Normalise RHS
                0107       rhsMax = 0. _d 0
                0108       DO bj=myByLo(myThid),myByHi(myThid)
                0109        DO bi=myBxLo(myThid),myBxHi(myThid)
2739a7f265 Jean*0110         DO j=1,sNy
                0111          DO i=1,sNx
                0112           cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)*cg2dNorm
                0113           rhsMax = MAX(ABS(cg2d_b(i,j,bi,bj)),rhsMax)
924557e60a Chri*0114          ENDDO
                0115         ENDDO
                0116        ENDDO
                0117       ENDDO
aea29c8517 Alis*0118
                0119       IF (cg2dNormaliseRHS) THEN
                0120 C-  Normalise RHS :
12c8b75709 Jean*0121       _GLOBAL_MAX_RL( rhsMax, myThid )
924557e60a Chri*0122       rhsNorm = 1. _d 0
                0123       IF ( rhsMax .NE. 0. ) rhsNorm = 1. _d 0 / rhsMax
                0124       DO bj=myByLo(myThid),myByHi(myThid)
                0125        DO bi=myBxLo(myThid),myBxHi(myThid)
2739a7f265 Jean*0126         DO j=1,sNy
                0127          DO i=1,sNx
                0128           cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)*rhsNorm
                0129           cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)*rhsNorm
924557e60a Chri*0130          ENDDO
                0131         ENDDO
                0132        ENDDO
                0133       ENDDO
aea29c8517 Alis*0134 C- end Normalise RHS
                0135       ENDIF
924557e60a Chri*0136
                0137 C--   Update overlaps
9155c64ca4 Jean*0138       CALL EXCH_XY_RL( cg2d_x, myThid )
924557e60a Chri*0139
                0140 C--   Initial residual calculation
                0141       DO bj=myByLo(myThid),myByHi(myThid)
                0142        DO bi=myBxLo(myThid),myBxHi(myThid)
ff51fd05c6 Jean*0143 C-    Initialise local working arrays:
                0144         DO j=0,sNy+1
                0145          DO i=0,sNx+1
                0146           cg2d_r(i,j,bi,bj) = 0. _d 0
                0147           cg2d_s(i,j,bi,bj) = 0. _d 0
                0148          ENDDO
                0149         ENDDO
2739a7f265 Jean*0150         IF ( nIterMin.GE.0 ) THEN
ff51fd05c6 Jean*0151 C-    Initialise saved solution
2739a7f265 Jean*0152          DO j=1,sNy
                0153           DO i=1,sNx
ff51fd05c6 Jean*0154            cg2d_min(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
2739a7f265 Jean*0155           ENDDO
                0156          ENDDO
                0157         ENDIF
a619b19bc9 Jean*0158         sumRHStile(bi,bj) = 0. _d 0
                0159         errTile(bi,bj)    = 0. _d 0
2b7cba8d4b Mart*0160 #ifdef TARGET_NEC_SX
                0161 !CDIR OUTERUNROLL=CG2D_OUTERLOOPITERS
                0162 #endif /* TARGET_NEC_SX */
2739a7f265 Jean*0163         DO j=1,sNy
                0164          DO i=1,sNx
                0165           cg2d_r(i,j,bi,bj) = cg2d_b(i,j,bi,bj) -
                0166      &    (aW2d(i  ,j  ,bi,bj)*cg2d_x(i-1,j  ,bi,bj)
                0167      &    +aW2d(i+1,j  ,bi,bj)*cg2d_x(i+1,j  ,bi,bj)
                0168      &    +aS2d(i  ,j  ,bi,bj)*cg2d_x(i  ,j-1,bi,bj)
                0169      &    +aS2d(i  ,j+1,bi,bj)*cg2d_x(i  ,j+1,bi,bj)
                0170      &    +aC2d(i  ,j  ,bi,bj)*cg2d_x(i  ,j  ,bi,bj)
c0a4efc370 Chri*0171      &    )
d439db03c9 Oliv*0172 #ifdef CG2D_SINGLECPU_SUM
2739a7f265 Jean*0173           localBuf(i,j,bi,bj) = cg2d_r(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
d439db03c9 Oliv*0174 #else
a619b19bc9 Jean*0175           errTile(bi,bj)    = errTile(bi,bj)
2739a7f265 Jean*0176      &                      + cg2d_r(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
d439db03c9 Oliv*0177 #endif
d68b36a485 Mart*0178           sumRHStile(bi,bj) = sumRHStile(bi,bj) + cg2d_b(i,j,bi,bj)
924557e60a Chri*0179          ENDDO
                0180         ENDDO
                0181        ENDDO
                0182       ENDDO
9155c64ca4 Jean*0183       CALL EXCH_S3D_RL( cg2d_r, 1, myThid )
d439db03c9 Oliv*0184 #ifdef CG2D_SINGLECPU_SUM
2739a7f265 Jean*0185       CALL GLOBAL_SUM_SINGLECPU_RL(localBuf, err_sq, 0, 0, myThid)
d439db03c9 Oliv*0186 #else
2739a7f265 Jean*0187       CALL GLOBAL_SUM_TILE_RL( errTile,    err_sq, myThid )
d439db03c9 Oliv*0188 #endif
d68b36a485 Mart*0189       CALL GLOBAL_SUM_TILE_RL( sumRHStile, sumRHS, myThid )
2739a7f265 Jean*0190       actualIts = 0
                0191       firstResidual = SQRT(err_sq)
                0192       IF ( nIterMin.GE.0 ) THEN
                0193         nIterMin = 0
                0194         minResidualSq = err_sq
                0195       ENDIF
d607aa1d50 Dimi*0196
764aea3440 Jean*0197       printResidual = .FALSE.
9155c64ca4 Jean*0198       IF ( debugLevel .GE. debLevZero ) THEN
494ad43bae Patr*0199         _BEGIN_MASTER( myThid )
764aea3440 Jean*0200         printResidual = printResidualFreq.GE.1
4606c28752 Jean*0201         WRITE(standardmessageunit,'(A,1P2E22.14)')
                0202      &  ' cg2d: Sum(rhs),rhsMax = ', sumRHS,rhsMax
ea6e02f692 Ed H*0203         _END_MASTER( myThid )
9155c64ca4 Jean*0204       ENDIF
aea29c8517 Alis*0205
2739a7f265 Jean*0206       IF ( err_sq .LT. cg2dTolerance_sq ) GOTO 11
924557e60a Chri*0207
                0208 C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
030bea3287 Alis*0209       DO 10 it2d=1, numIters
924557e60a Chri*0210
                0211 C--    Solve preconditioning equation and update
                0212 C--    conjugate direction vector "s".
                0213        DO bj=myByLo(myThid),myByHi(myThid)
                0214         DO bi=myBxLo(myThid),myBxHi(myThid)
a619b19bc9 Jean*0215          eta_qrNtile(bi,bj) = 0. _d 0
2b7cba8d4b Mart*0216 #ifdef TARGET_NEC_SX
                0217 !CDIR OUTERUNROLL=CG2D_OUTERLOOPITERS
                0218 #endif /* TARGET_NEC_SX */
2739a7f265 Jean*0219          DO j=1,sNy
                0220           DO i=1,sNx
                0221            cg2d_q(i,j,bi,bj) =
                0222      &      pC(i  ,j  ,bi,bj)*cg2d_r(i  ,j  ,bi,bj)
                0223      &     +pW(i  ,j  ,bi,bj)*cg2d_r(i-1,j  ,bi,bj)
                0224      &     +pW(i+1,j  ,bi,bj)*cg2d_r(i+1,j  ,bi,bj)
                0225      &     +pS(i  ,j  ,bi,bj)*cg2d_r(i  ,j-1,bi,bj)
                0226      &     +pS(i  ,j+1,bi,bj)*cg2d_r(i  ,j+1,bi,bj)
c0a4efc370 Chri*0227 CcnhDebugStarts
2739a7f265 Jean*0228 c          cg2d_q(i,j,bi,bj) = cg2d_r(j  ,j  ,bi,bj)
c0a4efc370 Chri*0229 CcnhDebugEnds
d439db03c9 Oliv*0230 #ifdef CG2D_SINGLECPU_SUM
2739a7f265 Jean*0231           localBuf(i,j,bi,bj) =
                0232      &      cg2d_q(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
d439db03c9 Oliv*0233 #else
a619b19bc9 Jean*0234            eta_qrNtile(bi,bj) = eta_qrNtile(bi,bj)
2739a7f265 Jean*0235      &     +cg2d_q(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
d439db03c9 Oliv*0236 #endif
924557e60a Chri*0237           ENDDO
                0238          ENDDO
                0239         ENDDO
                0240        ENDDO
                0241
d439db03c9 Oliv*0242 #ifdef CG2D_SINGLECPU_SUM
71d7d0e7f1 Jean*0243        CALL GLOBAL_SUM_SINGLECPU_RL( localBuf,eta_qrN,0,0,myThid )
d439db03c9 Oliv*0244 #else
a619b19bc9 Jean*0245        CALL GLOBAL_SUM_TILE_RL( eta_qrNtile,eta_qrN,myThid )
d439db03c9 Oliv*0246 #endif
980a7d9757 Jean*0247        cgBeta   = eta_qrN/eta_qrNM1
924557e60a Chri*0248 CcnhDebugStarts
2739a7f265 Jean*0249 c      WRITE(*,*) ' CG2D: Iteration ', it2d-1,
                0250 c    &            ' eta_qrN=', eta_qrN, ' beta=', cgBeta
924557e60a Chri*0251 CcnhDebugEnds
980a7d9757 Jean*0252        eta_qrNM1 = eta_qrN
924557e60a Chri*0253
                0254        DO bj=myByLo(myThid),myByHi(myThid)
                0255         DO bi=myBxLo(myThid),myBxHi(myThid)
2739a7f265 Jean*0256          DO j=1,sNy
                0257           DO i=1,sNx
                0258            cg2d_s(i,j,bi,bj) = cg2d_q(i,j,bi,bj)
                0259      &                       + cgBeta*cg2d_s(i,j,bi,bj)
924557e60a Chri*0260           ENDDO
                0261          ENDDO
                0262         ENDDO
                0263        ENDDO
                0264
764aea3440 Jean*0265 C--    Do exchanges that require messages i.e. between processes.
9155c64ca4 Jean*0266        CALL EXCH_S3D_RL( cg2d_s, 1, myThid )
924557e60a Chri*0267
                0268 C==    Evaluate laplace operator on conjugate gradient vector
                0269 C==    q = A.s
                0270        DO bj=myByLo(myThid),myByHi(myThid)
                0271         DO bi=myBxLo(myThid),myBxHi(myThid)
a619b19bc9 Jean*0272          alphaTile(bi,bj) = 0. _d 0
2b7cba8d4b Mart*0273 #ifdef TARGET_NEC_SX
                0274 !CDIR OUTERUNROLL=CG2D_OUTERLOOPITERS
                0275 #endif /* TARGET_NEC_SX */
2739a7f265 Jean*0276          DO j=1,sNy
                0277           DO i=1,sNx
                0278            cg2d_q(i,j,bi,bj) =
                0279      &     aW2d(i  ,j  ,bi,bj)*cg2d_s(i-1,j  ,bi,bj)
                0280      &    +aW2d(i+1,j  ,bi,bj)*cg2d_s(i+1,j  ,bi,bj)
                0281      &    +aS2d(i  ,j  ,bi,bj)*cg2d_s(i  ,j-1,bi,bj)
                0282      &    +aS2d(i  ,j+1,bi,bj)*cg2d_s(i  ,j+1,bi,bj)
                0283      &    +aC2d(i  ,j  ,bi,bj)*cg2d_s(i  ,j  ,bi,bj)
d439db03c9 Oliv*0284 #ifdef CG2D_SINGLECPU_SUM
2739a7f265 Jean*0285           localBuf(i,j,bi,bj) = cg2d_s(i,j,bi,bj)*cg2d_q(i,j,bi,bj)
d439db03c9 Oliv*0286 #else
a619b19bc9 Jean*0287           alphaTile(bi,bj) = alphaTile(bi,bj)
2739a7f265 Jean*0288      &                     + cg2d_s(i,j,bi,bj)*cg2d_q(i,j,bi,bj)
d439db03c9 Oliv*0289 #endif
924557e60a Chri*0290           ENDDO
                0291          ENDDO
                0292         ENDDO
                0293        ENDDO
d439db03c9 Oliv*0294 #ifdef CG2D_SINGLECPU_SUM
71d7d0e7f1 Jean*0295        CALL GLOBAL_SUM_SINGLECPU_RL(localBuf, alpha, 0, 0, myThid)
d439db03c9 Oliv*0296 #else
a619b19bc9 Jean*0297        CALL GLOBAL_SUM_TILE_RL( alphaTile,  alpha,  myThid )
d439db03c9 Oliv*0298 #endif
924557e60a Chri*0299 CcnhDebugStarts
2739a7f265 Jean*0300 c      WRITE(*,*) ' CG2D: Iteration ', it2d-1,
                0301 c    &            ' SUM(s*q)=', alpha, ' alpha=', eta_qrN/alpha
924557e60a Chri*0302 CcnhDebugEnds
980a7d9757 Jean*0303        alpha = eta_qrN/alpha
4606c28752 Jean*0304
2739a7f265 Jean*0305 C==    Update simultaneously solution and residual vectors (and Iter number)
924557e60a Chri*0306 C      Now compute "interior" points.
                0307        DO bj=myByLo(myThid),myByHi(myThid)
                0308         DO bi=myBxLo(myThid),myBxHi(myThid)
a619b19bc9 Jean*0309          errTile(bi,bj) = 0. _d 0
2739a7f265 Jean*0310          DO j=1,sNy
                0311           DO i=1,sNx
                0312            cg2d_x(i,j,bi,bj)=cg2d_x(i,j,bi,bj)+alpha*cg2d_s(i,j,bi,bj)
                0313            cg2d_r(i,j,bi,bj)=cg2d_r(i,j,bi,bj)-alpha*cg2d_q(i,j,bi,bj)
d439db03c9 Oliv*0314 #ifdef CG2D_SINGLECPU_SUM
2739a7f265 Jean*0315            localBuf(i,j,bi,bj) = cg2d_r(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
d439db03c9 Oliv*0316 #else
a619b19bc9 Jean*0317            errTile(bi,bj) = errTile(bi,bj)
2739a7f265 Jean*0318      &                    + cg2d_r(i,j,bi,bj)*cg2d_r(i,j,bi,bj)
d439db03c9 Oliv*0319 #endif
924557e60a Chri*0320           ENDDO
                0321          ENDDO
                0322         ENDDO
                0323        ENDDO
2739a7f265 Jean*0324        actualIts = it2d
924557e60a Chri*0325
d439db03c9 Oliv*0326 #ifdef CG2D_SINGLECPU_SUM
2739a7f265 Jean*0327        CALL GLOBAL_SUM_SINGLECPU_RL(localBuf, err_sq, 0, 0, myThid)
d439db03c9 Oliv*0328 #else
2739a7f265 Jean*0329        CALL GLOBAL_SUM_TILE_RL( errTile,    err_sq,    myThid )
d439db03c9 Oliv*0330 #endif
764aea3440 Jean*0331        IF ( printResidual ) THEN
                0332         IF ( MOD( it2d-1, printResidualFreq ).EQ.0 ) THEN
                0333          WRITE(msgBuf,'(A,I6,A,1PE21.14)')
2739a7f265 Jean*0334      &    ' cg2d: iter=', it2d, ' ; resid.= ', SQRT(err_sq)
764aea3440 Jean*0335          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0336      &                       SQUEEZE_RIGHT, myThid )
                0337         ENDIF
1f36709e5c Jean*0338        ENDIF
2739a7f265 Jean*0339        IF ( err_sq .LT. cg2dTolerance_sq ) GOTO 11
                0340        IF ( err_sq .LT. minResidualSq ) THEN
                0341 C-     Store lowest residual solution
                0342          minResidualSq = err_sq
                0343          nIterMin = it2d
                0344          DO bj=myByLo(myThid),myByHi(myThid)
                0345           DO bi=myBxLo(myThid),myBxHi(myThid)
                0346            DO j=1,sNy
                0347             DO i=1,sNx
                0348              cg2d_min(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
                0349             ENDDO
                0350            ENDDO
                0351           ENDDO
                0352          ENDDO
                0353        ENDIF
1f36709e5c Jean*0354
9155c64ca4 Jean*0355        CALL EXCH_S3D_RL( cg2d_r, 1, myThid )
a85d6ab24e Chri*0356
924557e60a Chri*0357    10 CONTINUE
                0358    11 CONTINUE
                0359
2739a7f265 Jean*0360       IF ( nIterMin.GE.0 .AND. err_sq .GT. minResidualSq ) THEN
                0361 C-    use the lowest residual solution (instead of current one = last residual)
                0362         DO bj=myByLo(myThid),myByHi(myThid)
                0363          DO bi=myBxLo(myThid),myBxHi(myThid)
                0364           DO j=1,sNy
                0365            DO i=1,sNx
                0366              cg2d_x(i,j,bi,bj) = cg2d_min(i,j,bi,bj)
                0367            ENDDO
                0368           ENDDO
                0369          ENDDO
                0370         ENDDO
                0371       ENDIF
                0372
aea29c8517 Alis*0373       IF (cg2dNormaliseRHS) THEN
924557e60a Chri*0374 C--   Un-normalise the answer
aea29c8517 Alis*0375         DO bj=myByLo(myThid),myByHi(myThid)
                0376          DO bi=myBxLo(myThid),myBxHi(myThid)
2739a7f265 Jean*0377           DO j=1,sNy
                0378            DO i=1,sNx
                0379             cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)/rhsNorm
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