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