/[MITgcm]/MITgcm/model/src/cg3d.F
ViewVC logotype

Diff of /MITgcm/model/src/cg3d.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.7 by jmc, Thu Mar 8 20:45:53 2001 UTC revision 1.25 by jmc, Fri May 11 23:34:06 2012 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    #ifdef TARGET_NEC_SX
6    C     set a sensible default for the outer loop unrolling parameter that can
7    C     be overriden in the Makefile with the DEFINES macro or in CPP_OPTIONS.h
8    #ifndef CG3D_OUTERLOOPITERS
9    # define CG3D_OUTERLOOPITERS 10
10    #endif
11    #endif /* TARGET_NEC_SX */
12    
13  #define VERBOSE  CBOP
14    C     !ROUTINE: CG3D
15    C     !INTERFACE:
16          SUBROUTINE CG3D(
17         U                cg3d_b, cg3d_x,
18         O                firstResidual, lastResidual,
19         U                numIters,
20         I                myIter, myThid )
21    C     !DESCRIPTION: \bv
22    C     *==========================================================*
23    C     | SUBROUTINE CG3D
24    C     | o Three-dimensional grid problem conjugate-gradient
25    C     |   inverter (with preconditioner).
26    C     *==========================================================*
27    C     | Con. grad is an iterative procedure for solving Ax = b.
28    C     | It requires the A be symmetric.
29    C     | This implementation assumes A is a seven-diagonal
30    C     | matrix of the form that arises in the discrete
31    C     | representation of the del^2 operator in a
32    C     | three-dimensional space.
33    C     | Notes:
34    C     | ======
35    C     | This implementation can support shared-memory
36    C     | multi-threaded execution. In order to do this COMMON
37    C     | blocks are used for many of the arrays - even ones that
38    C     | are only used for intermedaite results. This design is
39    C     | OK if you want to all the threads to collaborate on
40    C     | solving the same problem. On the other hand if you want
41    C     | the threads to solve several different problems
42    C     | concurrently this implementation will not work.
43    C     *==========================================================*
44    C     \ev
45    
46        SUBROUTINE CG3D(    C     !USES:
      I                myThid )  
 C     /==========================================================\  
 C     | SUBROUTINE CG3D                                          |  
 C     | o Three-dimensional grid problem conjugate-gradient      |  
 C     |   inverter (with preconditioner).                        |  
 C     |==========================================================|  
 C     | Con. grad is an iterative procedure for solving Ax = b.  |  
 C     | It requires the A be symmetric.                          |  
 C     | This implementation assumes A is a five-diagonal         |  
 C     | matrix of the form that arises in the discrete           |  
 C     | representation of the del^2 operator in a                |  
 C     | two-dimensional space.                                   |  
 C     | Notes:                                                   |  
 C     | ======                                                   |  
 C     | This implementation can support shared-memory            |  
 C     | multi-threaded execution. In order to do this COMMON     |  
 C     | blocks are used for many of the arrays - even ones that  |  
 C     | are only used for intermedaite results. This design is   |  
 C     | OK if you want to all the threads to collaborate on      |  
 C     | solving the same problem. On the other hand if you want  |  
 C     | the threads to solve several different problems          |  
 C     | concurrently this implementation will not work.          |  
 C     \==========================================================/  
47        IMPLICIT NONE        IMPLICIT NONE
   
48  C     === Global data ===  C     === Global data ===
49  #include "SIZE.h"  #include "SIZE.h"
50  #include "EEPARAMS.h"  #include "EEPARAMS.h"
51  #include "PARAMS.h"  #include "PARAMS.h"
52  #include "GRID.h"  #include "GRID.h"
53    #include "SURFACE.h"
54  #include "CG3D.h"  #include "CG3D.h"
55    
56    C     !INPUT/OUTPUT PARAMETERS:
57  C     === Routine arguments ===  C     === Routine arguments ===
58  C     myThid - Thread on which I am working.  C     cg3d_b    :: The source term or "right hand side" (output: normalised RHS)
59    C     cg3d_x    :: The solution (input: first guess)
60    C     firstResidual :: the initial residual before any iterations
61    C     minResidualSq :: the lowest residual reached (squared)
62    CC    lastResidual  :: the actual residual reached
63    C     numIters  :: Inp: the maximum number of iterations allowed
64    C                  Out: the actual number of iterations used
65    CC    nIterMin  :: Inp: decide to store (if >=0) or not (if <0) lowest res. sol.
66    CC                 Out: iteration number corresponding to lowest residual
67    C     myIter    :: Current iteration number in simulation
68    C     myThid    :: my Thread Id number
69          _RL     cg3d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
70          _RL     cg3d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
71          _RL     firstResidual
72          _RL     lastResidual
73          INTEGER numIters
74          INTEGER myIter
75        INTEGER myThid        INTEGER myThid
76    
77  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
78    
79    C     !LOCAL VARIABLES:
80  C     === Local variables ====  C     === Local variables ====
81  C     actualIts      - Number of iterations taken  C     bi, bj     :: tile index in X and Y.
82  C     actualResidual - residual  C     i, j, k    :: Loop counters
83  C     bi          - Block index in X and Y.  C     it3d       :: Loop counter for CG iterations
84  C     bj  C     actualIts  :: actual CG iteration number
85  C     eta_qrN     - Used in computing search directions  C     err_sq     :: Measure of the square of the residual of Ax - b.
86  C     eta_qrNM1     suffix N and NM1 denote current and  C     eta_qrN    :: Used in computing search directions; suffix N and NM1
87  C     cgBeta        previous iterations respectively.  C     eta_qrNM1     denote current and previous iterations respectively.
88  C     alpha    C     cgBeta     :: coeff used to update conjugate direction vector "s".
89  C     sumRHS      - Sum of right-hand-side. Sometimes this is a  C     alpha      :: coeff used to update solution & residual
90  C                   useful debuggin/trouble shooting diagnostic.  C     sumRHS     :: Sum of right-hand-side. Sometimes this is a useful
91  C                   For neumann problems sumRHS needs to be ~0.  C                   debugging/trouble shooting diagnostic. For neumann problems
92  C                   or they converge at a non-zero residual.  C                   sumRHS needs to be ~0 or it converge at a non-zero residual.
93  C     err         - Measure of residual of Ax - b, usually the norm.  C     cg2d_min   :: used to store solution corresponding to lowest residual.
94  C     I, J, N     - Loop counters ( N counts CG iterations )  C     msgBuf     :: Informational/error message buffer
95          INTEGER bi, bj
96          INTEGER i, j, k, it3d
97        INTEGER actualIts        INTEGER actualIts
98        _RL    actualResidual        INTEGER km1, kp1
99        INTEGER bi, bj                      _RL     maskM1, maskP1
100        INTEGER I, J, K, it3d        _RL     cg3dTolerance_sq
101        INTEGER KM1, KP1        _RL     err_sq,  errTile(nSx,nSy)
102        _RL    err        _RL     eta_qrN, eta_qrNtile(nSx,nSy)
103        _RL    eta_qrN        _RL     eta_qrNM1
104        _RL    eta_qrNM1        _RL     cgBeta
105        _RL    cgBeta        _RL     alpha ,  alphaTile(nSx,nSy)
106        _RL    alpha        _RL     sumRHS,  sumRHStile(nSx,nSy)
107        _RL    sumRHS        _RL     rhsMax
108        _RL    rhsMax        _RL     rhsNorm
109        _RL    rhsNorm        CHARACTER*(MAX_LEN_MBUF) msgBuf
110          LOGICAL printResidual
111        INTEGER OLw        _RL     surfFac
112        INTEGER OLe  #ifdef NONLIN_FRSURF
113        INTEGER OLn        INTEGER ks
114        INTEGER OLs        _RL     surfTerm(sNx,sNy)
115        INTEGER exchWidthX  #endif /* NONLIN_FRSURF */
116        INTEGER exchWidthY  CEOP
117        INTEGER myNz  
118        _RL     topLevTerm  C--   Initialise auxiliary constant, some output variable
119          cg3dTolerance_sq = cg3dTargetResidual*cg3dTargetResidual
120          IF ( select_rStar .NE. 0 ) THEN
121            surfFac = freeSurfFac
122          ELSE
123            surfFac = 0.
124          ENDIF
125    #ifdef NONLIN_FRSURF
126          DO j=1,sNy
127           DO i=1,sNx
128             surfTerm(i,j) = 0.
129           ENDDO
130          ENDDO
131    #endif /* NONLIN_FRSURF */
132    
133  C--   Initialise inverter  C--   Initialise inverter
134        eta_qrNM1 = 1. D0        eta_qrNM1 = 1. _d 0
135    
136  C--   Normalise RHS  C--   Normalise RHS
137        rhsMax = 0. _d 0        rhsMax = 0. _d 0
138        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
139         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
140          DO K=1,Nr          DO k=1,Nr
141           DO J=1,sNy           DO j=1,sNy
142            DO I=1,sNx            DO i=1,sNx
143             cg3d_b(I,J,K,bi,bj) = cg3d_b(I,J,K,bi,bj)*cg3dNorm             cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)*cg3dNorm
144             rhsMax = MAX(ABS(cg3d_b(I,J,K,bi,bj)),rhsMax)       &                          * maskC(i,j,k,bi,bj)
145               rhsMax = MAX(ABS(cg3d_b(i,j,k,bi,bj)),rhsMax)
146            ENDDO            ENDDO
147           ENDDO           ENDDO
148          ENDDO          ENDDO
149         ENDDO         ENDDO
150        ENDDO        ENDDO
151        _GLOBAL_MAX_R8( rhsMax, myThid )        _GLOBAL_MAX_RL( rhsMax, myThid )
152        rhsNorm = 1. _d 0        rhsNorm = 1. _d 0
153        IF ( rhsMax .NE. 0. ) rhsNorm = 1. _d 0 / rhsMax        IF ( rhsMax .NE. 0. ) rhsNorm = 1. _d 0 / rhsMax
154        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
155         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
156          DO K=1,Nr          DO k=1,Nr
157           DO J=1,sNy           DO j=1,sNy
158            DO I=1,sNx            DO i=1,sNx
159             cg3d_b(I,J,K,bi,bj) = cg3d_b(I,J,K,bi,bj)*rhsNorm             cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)*rhsNorm
160             cg3d_x(I,J,K,bi,bj) = cg3d_x(I,J,K,bi,bj)*rhsNorm             cg3d_x(i,j,k,bi,bj) = cg3d_x(i,j,k,bi,bj)*rhsNorm
161            ENDDO            ENDDO
162           ENDDO           ENDDO
163          ENDDO          ENDDO
# Line 116  C--   Normalise RHS Line 165  C--   Normalise RHS
165        ENDDO        ENDDO
166    
167  C--   Update overlaps  C--   Update overlaps
168        _EXCH_XYZ_R8( cg3d_b, myThid )        _EXCH_XYZ_RL( cg3d_x, myThid )
       _EXCH_XYZ_R8( cg3d_x, myThid )  
169    
170  C--   Initial residual calculation (with free-Surface term)  C--   Initial residual calculation (with free-Surface term)
       err    = 0. _d 0  
       sumRHS = 0. _d 0  
171        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
172         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
173          DO K=1,Nr          errTile(bi,bj)    = 0. _d 0
174           KM1 = K-1          sumRHStile(bi,bj) = 0. _d 0
175           IF ( K .EQ. 1 ) KM1 = 1  #ifdef NONLIN_FRSURF
176           KP1 = K+1          IF ( select_rStar .NE. 0 ) THEN
177           IF ( K .EQ. Nr ) KP1 = 1            DO j=1,sNy
178           topLevTerm = 0.             DO i=1,sNx
179           IF ( K .EQ. 1) topLevTerm = freeSurfFac*cg3dNorm*               surfTerm(i,j) = 0.
180       &      (horiVertRatio/gravity)/deltaTMom/deltaTMom             ENDDO
181           DO J=1,sNy            ENDDO
182            DO I=1,sNx            DO k=1,Nr
183             cg3d_s(I,J,K,bi,bj) = 0.             DO j=1,sNy
184             cg3d_r(I,J,K,bi,bj) = cg3d_b(I,J,K,bi,bj) -( 0.              DO i=1,sNx
185       &     +aW3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I-1,J  ,K  ,bi,bj)               surfTerm(i,j) = surfTerm(i,j)
186       &     +aW3d(I+1,J  ,K  ,bi,bj)*cg3d_x(I+1,J  ,K  ,bi,bj)       &         +cg3d_x(i,j,k,bi,bj)*drF(k)*h0FacC(i,j,k,bi,bj)
187       &     +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J-1,K  ,bi,bj)              ENDDO
188       &     +aS3d(I  ,J+1,K  ,bi,bj)*cg3d_x(I  ,J+1,K  ,bi,bj)             ENDDO
189       &     +aV3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J  ,KM1,bi,bj)            ENDDO
190       &     +aV3d(I  ,J  ,KP1,bi,bj)*cg3d_x(I  ,J  ,KP1,bi,bj)            DO j=1,sNy
191       &     -aW3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)             DO i=1,sNx
192       &     -aW3d(I+1,J  ,K  ,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)               ks = kSurfC(i,j,bi,bj)
193       &     -aS3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)               surfTerm(i,j) = surfTerm(i,j)*cg3dNorm
194       &     -aS3d(I  ,J+1,K  ,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)       &              *recip_Rcol(i,j,bi,bj)*recip_Rcol(i,j,bi,bj)
195       &     -aV3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)       &              *rA(i,j,bi,bj)*deepFac2F(ks)
196       &     -aV3d(I  ,J  ,KP1,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)       &              *recip_Bo(i,j,bi,bj)/deltaTMom/deltaTfreesurf
197       &     -topLevTerm*_rA(I,J,bi,bj)*cg3d_x(I,J,K,bi,bj)             ENDDO
198       &     )            ENDDO
199             err = err          ENDIF
200       &     +cg3d_r(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)  #endif /* NONLIN_FRSURF */
201             sumRHS = sumRHS          DO k=1,Nr
202       &     +cg3d_b(I,J,K,bi,bj)           km1 = MAX(k-1, 1 )
203             kp1 = MIN(k+1, Nr)
204             maskM1 = 1. _d 0
205             maskP1 = 1. _d 0
206             IF ( k .EQ. 1 ) maskM1 = 0. _d 0
207             IF ( k .EQ. Nr) maskP1 = 0. _d 0
208    #ifdef TARGET_NEC_SX
209    !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
210    #endif /* TARGET_NEC_SX */
211             DO j=1,sNy
212              DO i=1,sNx
213               cg3d_r(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
214         &     -( 0.
215         &       +aW3d( i, j, k, bi,bj)*cg3d_x(i-1,j, k, bi,bj)
216         &       +aW3d(i+1,j, k, bi,bj)*cg3d_x(i+1,j, k, bi,bj)
217         &       +aS3d( i, j, k, bi,bj)*cg3d_x( i,j-1,k, bi,bj)
218         &       +aS3d( i,j+1,k, bi,bj)*cg3d_x( i,j+1,k, bi,bj)
219         &       +aV3d( i, j, k, bi,bj)*cg3d_x( i, j,km1,bi,bj)*maskM1
220         &       +aV3d( i, j,kp1,bi,bj)*cg3d_x( i, j,kp1,bi,bj)*maskP1
221         &       +aC3d( i, j, k, bi,bj)*cg3d_x( i, j, k, bi,bj)
222    #ifdef NONLIN_FRSURF
223         &       -surfFac*surfTerm(i,j)*drF(k)*h0FacC(i,j,k,bi,bj)
224    #endif /* NONLIN_FRSURF */
225         &      )
226               errTile(bi,bj) = errTile(bi,bj)
227         &                    +cg3d_r(i,j,k,bi,bj)*cg3d_r(i,j,k,bi,bj)
228               sumRHStile(bi,bj) = sumRHStile(bi,bj)+cg3d_b(i,j,k,bi,bj)
229              ENDDO
230             ENDDO
231             DO j=0,sNy+1
232              DO i=0,sNx+1
233               cg3d_s(i,j,k,bi,bj) = 0.
234            ENDDO            ENDDO
235           ENDDO           ENDDO
236          ENDDO          ENDDO
237         ENDDO         ENDDO
238        ENDDO        ENDDO
239  C     _EXCH_XYZ_R8( cg3d_r, myThid )         CALL EXCH_S3D_RL( cg3d_r, Nr, myThid )
240         OLw        = 1         CALL GLOBAL_SUM_TILE_RL( sumRHStile, sumRHS, myThid )
241         OLe        = 1         CALL GLOBAL_SUM_TILE_RL( errTile,    err_sq, myThid )
242         OLn        = 1        IF ( debugLevel.GE.debLevC .AND. diagFreq.GT.0. ) THEN
243         OLs        = 1          CALL WRITE_FLD_S3D_RL(
244         exchWidthX = 1       I        'cg3d_r_I', 'I10', 1, Nr, cg3d_r, myIter, myThid )
245         exchWidthY = 1        ENDIF
246         myNz       = Nr  
247         CALL EXCH_RL( cg3d_r,        actualIts = 0
248       I            OLw, OLe, OLs, OLn, myNz,        firstResidual = SQRT(err_sq)
249       I            exchWidthX, exchWidthY,  
250       I            FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )        printResidual = .FALSE.
251  C     _EXCH_XYZ_R8( cg3d_s, myThid )        IF ( debugLevel .GE. debLevZero ) THEN
252         OLw        = 1          _BEGIN_MASTER( myThid )
253         OLe        = 1          printResidual = printResidualFreq.GE.1
254         OLn        = 1          WRITE(standardmessageunit,'(A,1P2E22.14)')
255         OLs        = 1       &     ' cg3d: Sum(rhs),rhsMax = ',sumRHS,rhsMax
256         exchWidthX = 1          _END_MASTER( myThid )
257         exchWidthY = 1        ENDIF
258         myNz       = Nr  
259         CALL EXCH_RL( cg3d_s,        IF ( err_sq .LT. cg3dTolerance_sq ) GOTO 11
      I            OLw, OLe, OLs, OLn, myNz,  
      I            exchWidthX, exchWidthY,  
      I            FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )  
       _GLOBAL_SUM_R8( sumRHS, myThid )  
       _GLOBAL_SUM_R8( err   , myThid )  
         
       _BEGIN_MASTER( myThid )  
       write(0,'(A,1PE30.14)') ' cg3d: Sum(rhs) = ',sumRHS  
       _END_MASTER( )  
   
       actualIts      = 0  
       actualResidual = SQRT(err)  
 C     _BARRIER  
       _BEGIN_MASTER( myThid )  
        WRITE(0,'(A,I6,1PE30.14)') ' CG3D iters, err = ',  
      &  actualIts, actualResidual  
       _END_MASTER( )  
260    
261  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
262        DO 10 it3d=1, cg3dMaxIters        DO 10 it3d=1, numIters
263    
 CcnhDebugStarts  
 #ifdef VERBOSE  
        IF ( mod(it3d-1,10).EQ.0)  
      &   WRITE(0,*) ' CG3D: Iteration ',it3d-1,  
      &      ' residual = ',actualResidual  
 #endif  
 CcnhDebugEnds  
        IF ( actualResidual .LT. cg3dTargetResidual ) GOTO 11  
264  C--    Solve preconditioning equation and update  C--    Solve preconditioning equation and update
265  C--    conjugate direction vector "s".  C--    conjugate direction vector "s".
266  C      Note. On the next to loops over all tiles the inner loop ranges  C      Note. On the next two loops over all tiles the inner loop ranges
267  C            in sNx and sNy are expanded by 1 to avoid a communication  C            in sNx and sNy are expanded by 1 to avoid a communication
268  C            step. However this entails a bit of gynamastics because we only  C            step. However this entails a bit of gynamastics because we only
269  C            want eta_qrN for the interior points.  C            want eta_qrN for the interior points.
        eta_qrN = 0. _d 0  
270         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
271          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
272           DO K=1,1           eta_qrNtile(bi,bj) = 0. _d 0
273            DO J=1-1,sNy+1           DO k=1,1
274             DO I=1-1,sNx+1  #ifdef TARGET_NEC_SX
275              cg3d_q(I,J,K,bi,bj) =  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
276       &       zMC(I  ,J  ,K,bi,bj)*cg3d_r(I  ,J  ,K,bi,bj)  #endif /* TARGET_NEC_SX */
277              DO j=0,sNy+1
278               DO i=0,sNx+1
279                cg3d_q(i,j,k,bi,bj) = zMC(i,j,k,bi,bj)
280         &                        *cg3d_r(i,j,k,bi,bj)
281             ENDDO             ENDDO
282            ENDDO            ENDDO
283           ENDDO           ENDDO
284           DO K=2,Nr           DO k=2,Nr
285            DO J=1-1,sNy+1  #ifdef TARGET_NEC_SX
286             DO I=1-1,sNx+1  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
287              cg3d_q(I,J,K,bi,bj) =  #endif /* TARGET_NEC_SX */
288       &       zMC(I,J,K,bi,bj)*(cg3d_r(I,J,K  ,bi,bj)            DO j=0,sNy+1
289       &      -zML(I,J,K,bi,bj)*cg3d_q(I,J,K-1,bi,bj))             DO i=0,sNx+1
290                cg3d_q(i,j,k,bi,bj) = zMC(i,j,k,bi,bj)
291         &                      *( cg3d_r(i,j,k,bi,bj)
292         &                         -zML(i,j,k,bi,bj)*cg3d_q(i,j,k-1,bi,bj)
293         &                       )
294             ENDDO             ENDDO
295            ENDDO            ENDDO
296           ENDDO           ENDDO
297           DO K=Nr,Nr           DO k=Nr,Nr
298  caja      IF (Nr .GT. 1) THEN  #ifdef TARGET_NEC_SX
299  caja       DO J=1-1,sNy+1  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
300  caja        DO I=1-1,sNx+1  #endif /* TARGET_NEC_SX */
301  caja         cg3d_q(I,J,K,bi,bj) =            DO j=1,sNy
302  caja &        zMC(i,j,k,bi,bj)*(cg3d_r(i,j,k  ,bi,bj)             DO i=1,sNx
303  caja &       -zML(i,j,k,bi,bj)*cg3d_q(i,j,k-1,bi,bj))              eta_qrNtile(bi,bj) = eta_qrNtile(bi,bj)
304  caja        ENDDO       &                        +cg3d_q(i,j,k,bi,bj)*cg3d_r(i,j,k,bi,bj)
 caja       ENDDO  
 caja      ENDIF  
           DO J=1,sNy  
            DO I=1,sNx  
             eta_qrN = eta_qrN  
      &      +cg3d_q(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)  
305             ENDDO             ENDDO
306            ENDDO            ENDDO
307           ENDDO           ENDDO
308           DO K=Nr-1,1,-1           DO k=Nr-1,1,-1
309            DO J=1-1,sNy+1  #ifdef TARGET_NEC_SX
310             DO I=1-1,sNx+1  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
311              cg3d_q(I,J,K,bi,bj) =  #endif /* TARGET_NEC_SX */
312       &      cg3d_q(I,J,K,bi,bj)            DO j=0,sNy+1
313       &      -zMU(I,J,K,bi,bj)*cg3d_q(I,J,K+1,bi,bj)             DO i=0,sNx+1
314             ENDDO              cg3d_q(i,j,k,bi,bj) = cg3d_q(i,j,k,bi,bj)
315            ENDDO       &                         -zMU(i,j,k,bi,bj)*cg3d_q(i,j,k+1,bi,bj)
316            DO J=1,sNy             ENDDO
317             DO I=1,sNx            ENDDO
318              eta_qrN = eta_qrN  #ifdef TARGET_NEC_SX
319       &      +cg3d_q(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
320    #endif /* TARGET_NEC_SX */
321              DO j=1,sNy
322               DO i=1,sNx
323                eta_qrNtile(bi,bj) = eta_qrNtile(bi,bj)
324         &                        +cg3d_q(i,j,k,bi,bj)*cg3d_r(i,j,k,bi,bj)
325             ENDDO             ENDDO
326            ENDDO            ENDDO
327           ENDDO           ENDDO
328          ENDDO          ENDDO
329         ENDDO         ENDDO
 caja  
 caja  eta_qrN=0.  
 caja   DO bj=myByLo(myThid),myByHi(myThid)  
 caja    DO bi=myBxLo(myThid),myBxHi(myThid)  
 caja     DO K=1,Nr  
 caja      DO J=1,sNy  
 caja       DO I=1,sNx  
 caja        eta_qrN = eta_qrN  
 caja &      +cg3d_q(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)  
 caja       ENDDO  
 caja      ENDDO  
 caja     ENDDO  
 caja    ENDDO  
 caja   ENDDO  
 caja  
330    
331         _GLOBAL_SUM_R8(eta_qrN, myThid)         CALL GLOBAL_SUM_TILE_RL( eta_qrNtile,eta_qrN,myThid )
 CcnhDebugStarts  
 C      WRITE(0,*) ' CG3D: Iteration ',it3d-1,' eta_qrN = ',eta_qrN  
 CcnhDebugEnds  
332         cgBeta   = eta_qrN/eta_qrNM1         cgBeta   = eta_qrN/eta_qrNM1
333  CcnhDebugStarts  CcnhDebugStarts
334  C      WRITE(0,*) ' CG3D: Iteration ',it3d-1,' beta = ',cgBeta  c      WRITE(*,*) ' CG3D: Iteration ', it3d-1,
335    c    &            ' eta_qrN=', eta_qrN, ' beta=', cgBeta
336  CcnhDebugEnds  CcnhDebugEnds
337         eta_qrNM1 = eta_qrN         eta_qrNM1 = eta_qrN
338    
339         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
340          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
341           DO K=1,Nr           DO k=1,Nr
342            DO J=1-1,sNy+1            DO j=0,sNy+1
343             DO I=1-1,sNx+1             DO i=0,sNx+1
344              cg3d_s(I,J,K,bi,bj) = cg3d_q(I,J,K,bi,bj)              cg3d_s(i,j,k,bi,bj) = cg3d_q(i,j,k,bi,bj)
345       &                          + cgBeta*cg3d_s(I,J,K,bi,bj)       &                   + cgBeta*cg3d_s(i,j,k,bi,bj)
346             ENDDO             ENDDO
347            ENDDO            ENDDO
348           ENDDO           ENDDO
# Line 310  CcnhDebugEnds Line 351  CcnhDebugEnds
351    
352  C==    Evaluate laplace operator on conjugate gradient vector  C==    Evaluate laplace operator on conjugate gradient vector
353  C==    q = A.s  C==    q = A.s
        alpha = 0. _d 0  
        topLevTerm = freeSurfFac*cg3dNorm*  
      &      (horiVertRatio/gravity)/deltaTMom/deltaTMom  
354         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
355          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
356           IF ( Nr .GT. 1 ) THEN           alphaTile(bi,bj) = 0. _d 0
357            DO K=1,1  #ifdef NONLIN_FRSURF
358             DO J=1,sNy           IF ( select_rStar .NE. 0 ) THEN
359              DO I=1,sNx            DO j=1,sNy
360               cg3d_q(I,J,K,bi,bj) =             DO i=1,sNx
361       &       aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I-1,J  ,K  ,bi,bj)               surfTerm(i,j) = 0.
362       &      +aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I+1,J  ,K  ,bi,bj)             ENDDO
363       &      +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J-1,K  ,bi,bj)            ENDDO
364       &      +aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J+1,K  ,bi,bj)            DO k=1,Nr
365       &      +aV3d(I  ,J  ,K+1,bi,bj)*cg3d_s(I  ,J  ,K+1,bi,bj)             DO j=1,sNy
366       &      -aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)              DO i=1,sNx
367       &      -aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)               surfTerm(i,j) = surfTerm(i,j)
368       &      -aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &         +cg3d_s(i,j,k,bi,bj)*drF(k)*h0FacC(i,j,k,bi,bj)
      &      -aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)  
      &      -aV3d(I  ,J  ,K+1,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)  
      &      -topLevTerm*_rA(I,J,bi,bj)*cg3d_s(I,J,K,bi,bj)  
              alpha = alpha+cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)  
369              ENDDO              ENDDO
370             ENDDO             ENDDO
371            ENDDO            ENDDO
372              DO j=1,sNy
373               DO i=1,sNx
374                 ks = kSurfC(i,j,bi,bj)
375                 surfTerm(i,j) = surfTerm(i,j)*cg3dNorm
376         &              *recip_Rcol(i,j,bi,bj)*recip_Rcol(i,j,bi,bj)
377         &              *rA(i,j,bi,bj)*deepFac2F(ks)
378         &              *recip_Bo(i,j,bi,bj)/deltaTMom/deltaTfreesurf
379               ENDDO
380              ENDDO
381             ENDIF
382    #endif /* NONLIN_FRSURF */
383             IF ( Nr .GT. 1 ) THEN
384              k=1
385    #ifdef TARGET_NEC_SX
386    !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
387    #endif /* TARGET_NEC_SX */
388              DO j=1,sNy
389               DO i=1,sNx
390                cg3d_q(i,j,k,bi,bj) =
391         &        aW3d( i, j, k, bi,bj)*cg3d_s(i-1,j, k, bi,bj)
392         &       +aW3d(i+1,j, k, bi,bj)*cg3d_s(i+1,j, k, bi,bj)
393         &       +aS3d( i, j, k, bi,bj)*cg3d_s( i,j-1,k, bi,bj)
394         &       +aS3d( i,j+1,k, bi,bj)*cg3d_s( i,j+1,k, bi,bj)
395         &       +aV3d( i, j,k+1,bi,bj)*cg3d_s( i, j,k+1,bi,bj)
396         &       +aC3d( i, j, k, bi,bj)*cg3d_s( i, j, k, bi,bj)
397    #ifdef NONLIN_FRSURF
398         &       -surfFac*surfTerm(i,j)*drF(k)*h0FacC(i,j,k,bi,bj)
399    #endif /* NONLIN_FRSURF */
400                alphaTile(bi,bj) = alphaTile(bi,bj)
401         &             +cg3d_s(i,j,k,bi,bj)*cg3d_q(i,j,k,bi,bj)
402               ENDDO
403              ENDDO
404           ELSE           ELSE
405            DO K=1,1            k=1
406             DO J=1,sNy  #ifdef TARGET_NEC_SX
407              DO I=1,sNx  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
408               cg3d_q(I,J,K,bi,bj) =  #endif /* TARGET_NEC_SX */
409       &       aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I-1,J  ,K  ,bi,bj)            DO j=1,sNy
410       &      +aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I+1,J  ,K  ,bi,bj)             DO i=1,sNx
411       &      +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J-1,K  ,bi,bj)              cg3d_q(i,j,k,bi,bj) =
412       &      +aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J+1,K  ,bi,bj)       &        aW3d( i, j, k, bi,bj)*cg3d_s(i-1,j, k, bi,bj)
413       &      -aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &       +aW3d(i+1,j, k, bi,bj)*cg3d_s(i+1,j, k, bi,bj)
414       &      -aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &       +aS3d( i, j, k, bi,bj)*cg3d_s( i,j-1,k, bi,bj)
415       &      -aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &       +aS3d( i,j+1,k, bi,bj)*cg3d_s( i,j+1,k, bi,bj)
416       &      -aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &       +aC3d( i, j, k, bi,bj)*cg3d_s( i, j, k, bi,bj)
417       &      -topLevTerm*_rA(I,J,bi,bj)*cg3d_s(I,J,K,bi,bj)  #ifdef NONLIN_FRSURF
418               alpha = alpha+cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)       &       -surfFac*surfTerm(i,j)*drF(k)*h0FacC(i,j,k,bi,bj)
419              ENDDO  #endif /* NONLIN_FRSURF */
420                alphaTile(bi,bj) = alphaTile(bi,bj)
421         &             +cg3d_s(i,j,k,bi,bj)*cg3d_q(i,j,k,bi,bj)
422             ENDDO             ENDDO
423            ENDDO            ENDDO
424           ENDIF           ENDIF
425           DO K=2,Nr-1           DO k=2,Nr-1
426            DO J=1,sNy  #ifdef TARGET_NEC_SX
427             DO I=1,sNx  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
428              cg3d_q(I,J,K,bi,bj) =  #endif /* TARGET_NEC_SX */
429       &      aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I-1,J  ,K  ,bi,bj)            DO j=1,sNy
430       &     +aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I+1,J  ,K  ,bi,bj)             DO i=1,sNx
431       &     +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J-1,K  ,bi,bj)              cg3d_q(i,j,k,bi,bj) =
432       &     +aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J+1,K  ,bi,bj)       &        aW3d( i, j, k, bi,bj)*cg3d_s(i-1,j, k, bi,bj)
433       &     +aV3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K-1,bi,bj)       &       +aW3d(i+1,j, k, bi,bj)*cg3d_s(i+1,j, k, bi,bj)
434       &     +aV3d(I  ,J  ,K+1,bi,bj)*cg3d_s(I  ,J  ,K+1,bi,bj)       &       +aS3d( i, j, k, bi,bj)*cg3d_s( i,j-1,k, bi,bj)
435       &     -aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &       +aS3d( i,j+1,k, bi,bj)*cg3d_s( i,j+1,k, bi,bj)
436       &     -aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &       +aV3d( i, j, k, bi,bj)*cg3d_s( i, j,k-1,bi,bj)
437       &     -aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &       +aV3d( i, j,k+1,bi,bj)*cg3d_s( i, j,k+1,bi,bj)
438       &     -aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &       +aC3d( i, j, k, bi,bj)*cg3d_s( i, j, k, bi,bj)
439       &     -aV3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)  #ifdef NONLIN_FRSURF
440       &     -aV3d(I  ,J  ,K+1,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &       -surfFac*surfTerm(i,j)*drF(k)*h0FacC(i,j,k,bi,bj)
441              alpha = alpha+cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)  #endif /* NONLIN_FRSURF */
442                alphaTile(bi,bj) = alphaTile(bi,bj)
443         &             +cg3d_s(i,j,k,bi,bj)*cg3d_q(i,j,k,bi,bj)
444             ENDDO             ENDDO
445            ENDDO            ENDDO
446           ENDDO           ENDDO
447           IF ( Nr .GT. 1 ) THEN           IF ( Nr .GT. 1 ) THEN
448            DO K=Nr,Nr            k=Nr
449             DO J=1,sNy  #ifdef TARGET_NEC_SX
450              DO I=1,sNx  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
451               cg3d_q(I,J,K,bi,bj) =  #endif /* TARGET_NEC_SX */
452       &       aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I-1,J  ,K  ,bi,bj)            DO j=1,sNy
453       &      +aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I+1,J  ,K  ,bi,bj)             DO i=1,sNx
454       &      +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J-1,K  ,bi,bj)              cg3d_q(i,j,k,bi,bj) =
455       &      +aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J+1,K  ,bi,bj)       &        aW3d( i, j, k, bi,bj)*cg3d_s(i-1,j, k, bi,bj)
456       &      +aV3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K-1,bi,bj)       &       +aW3d(i+1,j, k, bi,bj)*cg3d_s(i+1,j, k, bi,bj)
457       &      -aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &       +aS3d( i, j, k, bi,bj)*cg3d_s( i,j-1,k, bi,bj)
458       &      -aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &       +aS3d( i,j+1,k, bi,bj)*cg3d_s( i,j+1,k, bi,bj)
459       &      -aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &       +aV3d( i, j, k, bi,bj)*cg3d_s( i, j,k-1,bi,bj)
460       &      -aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &       +aC3d( i, j, k, bi,bj)*cg3d_s( i, j, k, bi,bj)
461       &      -aV3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)  #ifdef NONLIN_FRSURF
462               alpha = alpha+cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)       &       -surfFac*surfTerm(i,j)*drF(k)*h0FacC(i,j,k,bi,bj)
463              ENDDO  #endif /* NONLIN_FRSURF */
464                alphaTile(bi,bj) = alphaTile(bi,bj)
465         &             +cg3d_s(i,j,k,bi,bj)*cg3d_q(i,j,k,bi,bj)
466             ENDDO             ENDDO
467            ENDDO            ENDDO
468           ENDIF           ENDIF
469          ENDDO          ENDDO
470         ENDDO         ENDDO
471         _GLOBAL_SUM_R8(alpha,myThid)         CALL GLOBAL_SUM_TILE_RL( alphaTile,  alpha,  myThid )
472  CcnhDebugStarts  CcnhDebugStarts
473  C      WRITE(0,*) ' CG3D: Iteration ',it3d-1,' SUM(s*q)= ',alpha  c      WRITE(*,*) ' CG3D: Iteration ', it3d-1,
474    c    &            ' SUM(s*q)=', alpha, ' alpha=', eta_qrN/alpha
475  CcnhDebugEnds  CcnhDebugEnds
476         alpha = eta_qrN/alpha         alpha = eta_qrN/alpha
477  CcnhDebugStarts  
478  C      WRITE(0,*) ' CG3D: Iteration ',it3d-1,' alpha= ',alpha  C==    Update simultaneously solution and residual vectors (and Iter number)
 CcnhDebugEnds  
       
 C==    Update solution and residual vectors  
479  C      Now compute "interior" points.  C      Now compute "interior" points.
        err = 0. _d 0  
480         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
481          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
482           DO K=1,Nr           errTile(bi,bj) = 0. _d 0
483            DO J=1,sNy           DO k=1,Nr
484             DO I=1,sNx  #ifdef TARGET_NEC_SX
485              cg3d_x(I,J,K,bi,bj)=cg3d_x(I,J,K,bi,bj)  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
486       &            +alpha*cg3d_s(I,J,K,bi,bj)  #endif /* TARGET_NEC_SX */
487              cg3d_r(I,J,K,bi,bj)=cg3d_r(I,J,K,bi,bj)            DO j=1,sNy
488       &            -alpha*cg3d_q(I,J,K,bi,bj)             DO i=1,sNx
489              err = err+cg3d_r(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)              cg3d_x(i,j,k,bi,bj)=cg3d_x(i,j,k,bi,bj)
490         &            +alpha*cg3d_s(i,j,k,bi,bj)
491                cg3d_r(i,j,k,bi,bj)=cg3d_r(i,j,k,bi,bj)
492         &            -alpha*cg3d_q(i,j,k,bi,bj)
493                errTile(bi,bj) = errTile(bi,bj)
494         &            +cg3d_r(i,j,k,bi,bj)*cg3d_r(i,j,k,bi,bj)
495             ENDDO             ENDDO
496            ENDDO            ENDDO
497           ENDDO           ENDDO
498          ENDDO          ENDDO
499         ENDDO         ENDDO
500           actualIts = it3d
501    
502         _GLOBAL_SUM_R8( err   , myThid )         CALL GLOBAL_SUM_TILE_RL( errTile,    err_sq, myThid )
503         err = SQRT(err)         IF ( printResidual ) THEN
504         actualIts      = it3d          IF ( MOD( it3d-1, printResidualFreq ).EQ.0 ) THEN
505         actualResidual = err           WRITE(msgBuf,'(A,I6,A,1PE21.14)')
506         IF ( actualResidual .LT. cg3dTargetResidual ) GOTO 11       &    ' cg3d: iter=', it3d, ' ; resid.= ', SQRT(err_sq)
507  C      _EXCH_XYZ_R8(cg3d_r, myThid )           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
508         OLw        = 1       &                       SQUEEZE_RIGHT, myThid )
509         OLe        = 1          ENDIF
510         OLn        = 1         ENDIF
511         OLs        = 1         IF ( err_sq .LT. cg3dTolerance_sq ) GOTO 11
512         exchWidthX = 1         CALL EXCH_S3D_RL( cg3d_r, Nr, myThid )
        exchWidthY = 1  
        myNz       = Nr  
        CALL EXCH_RL( cg3d_r,  
      I             OLw, OLe, OLs, OLn, myNz,  
      I             exchWidthX, exchWidthY,  
      I             FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )  
513    
514     10 CONTINUE     10 CONTINUE
515     11 CONTINUE     11 CONTINUE
516    
517          IF ( debugLevel.GE.debLevC .AND. diagFreq.GT.0. ) THEN
518            CALL WRITE_FLD_S3D_RL(
519         I        'cg3d_r_F', 'I10', 1, Nr, cg3d_r, myIter, myThid )
520          ENDIF
521    
522  C--   Un-normalise the answer  C--   Un-normalise the answer
523        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
524         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
525          DO K=1,Nr          DO k=1,Nr
526           DO J=1,sNy           DO j=1,sNy
527            DO I=1,sNx            DO i=1,sNx
528             cg3d_x(I,J,K,bi,bj) = cg3d_x(I,J,K,bi,bj)/rhsNorm             cg3d_x(i,j,k,bi,bj) = cg3d_x(i,j,k,bi,bj)/rhsNorm
529            ENDDO            ENDDO
530           ENDDO           ENDDO
531          ENDDO          ENDDO
532         ENDDO         ENDDO
533        ENDDO        ENDDO
534    
535  Cadj  _EXCH_XYZ_R8(cg3d_x, myThid )  C--   Return parameters to caller
536        _BEGIN_MASTER( myThid )        lastResidual = SQRT(err_sq)
537         WRITE(0,'(A,I6,1PE30.14)') ' CG3D iters, err = ',        numIters = actualIts
      &  actualIts, actualResidual  
       _END_MASTER( )  
538    
539  #endif /* ALLOW_NONHYDROSTATIC */  #endif /* ALLOW_NONHYDROSTATIC */
540    
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.25
  ViewVC Help
Powered by ViewVC 1.1.22