/[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.18 by jmc, Wed Aug 23 15:22:32 2006 UTC revision 1.19 by jmc, Tue Sep 4 14:54:58 2007 UTC
# Line 15  C     !INTERFACE: Line 15  C     !INTERFACE:
15       I                myThid )       I                myThid )
16  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
17  C     *==========================================================*  C     *==========================================================*
18  C     | SUBROUTINE CG3D                                            C     | SUBROUTINE CG3D
19  C     | o Three-dimensional grid problem conjugate-gradient        C     | o Three-dimensional grid problem conjugate-gradient
20  C     |   inverter (with preconditioner).                          C     |   inverter (with preconditioner).
21  C     *==========================================================*  C     *==========================================================*
22  C     | Con. grad is an iterative procedure for solving Ax = b.    C     | Con. grad is an iterative procedure for solving Ax = b.
23  C     | It requires the A be symmetric.                            C     | It requires the A be symmetric.
24  C     | This implementation assumes A is a seven-diagonal            C     | This implementation assumes A is a seven-diagonal
25  C     | matrix of the form that arises in the discrete              C     | matrix of the form that arises in the discrete
26  C     | representation of the del^2 operator in a                  C     | representation of the del^2 operator in a
27  C     | three-dimensional space.                                      C     | three-dimensional space.
28  C     | Notes:                                                      C     | Notes:
29  C     | ======                                                      C     | ======
30  C     | This implementation can support shared-memory                C     | This implementation can support shared-memory
31  C     | multi-threaded execution. In order to do this COMMON        C     | multi-threaded execution. In order to do this COMMON
32  C     | blocks are used for many of the arrays - even ones that      C     | blocks are used for many of the arrays - even ones that
33  C     | are only used for intermedaite results. This design is      C     | are only used for intermedaite results. This design is
34  C     | OK if you want to all the threads to collaborate on          C     | OK if you want to all the threads to collaborate on
35  C     | solving the same problem. On the other hand if you want      C     | solving the same problem. On the other hand if you want
36  C     | the threads to solve several different problems              C     | the threads to solve several different problems
37  C     | concurrently this implementation will not work.            C     | concurrently this implementation will not work.
38  C     *==========================================================*  C     *==========================================================*
39  C     \ev  C     \ev
40    
# Line 75  C     bj Line 75  C     bj
75  C     eta_qrN     - Used in computing search directions  C     eta_qrN     - Used in computing search directions
76  C     eta_qrNM1     suffix N and NM1 denote current and  C     eta_qrNM1     suffix N and NM1 denote current and
77  C     cgBeta        previous iterations respectively.  C     cgBeta        previous iterations respectively.
78  C     alpha    C     alpha
79  C     sumRHS      - Sum of right-hand-side. Sometimes this is a  C     sumRHS      - Sum of right-hand-side. Sometimes this is a
80  C                   useful debuggin/trouble shooting diagnostic.  C                   useful debuggin/trouble shooting diagnostic.
81  C                   For neumann problems sumRHS needs to be ~0.  C                   For neumann problems sumRHS needs to be ~0.
# Line 84  C     err         - Measure of residual Line 84  C     err         - Measure of residual
84  C     I, J, K, N  - Loop counters ( N counts CG iterations )  C     I, J, K, N  - Loop counters ( N counts CG iterations )
85        INTEGER actualIts        INTEGER actualIts
86        _RL    actualResidual        _RL    actualResidual
87        INTEGER bi, bj                      INTEGER bi, bj
88        INTEGER I, J, K, it3d        INTEGER I, J, K, it3d
89        INTEGER Km1, Kp1        INTEGER Km1, Kp1
90        _RL    maskM1, maskP1        _RL    maskM1, maskP1
91        _RL    err, errTile        _RL    err,    errTile(nSx,nSy)
92        _RL    eta_qrN, eta_qrNtile        _RL    eta_qrN,eta_qrNtile(nSx,nSy)
93        _RL    eta_qrNM1        _RL    eta_qrNM1
94        _RL    cgBeta        _RL    cgBeta
95        _RL    alpha , alphaTile        _RL    alpha , alphaTile(nSx,nSy)
96        _RL    sumRHS, sumRHStile        _RL    sumRHS, sumRHStile(nSx,nSy)
97        _RL    rhsMax        _RL    rhsMax
98        _RL    rhsNorm        _RL    rhsNorm
99  CEOP  CEOP
# Line 142  C--   Initial residual calculation (with Line 142  C--   Initial residual calculation (with
142        sumRHS = 0. _d 0        sumRHS = 0. _d 0
143        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
144         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
145          errTile    = 0. _d 0          errTile(bi,bj)    = 0. _d 0
146          sumRHStile = 0. _d 0          sumRHStile(bi,bj) = 0. _d 0
147          DO K=1,Nr          DO K=1,Nr
148           Km1 = MAX(K-1, 1 )           Km1 = MAX(K-1, 1 )
149           Kp1 = MIN(K+1, Nr)           Kp1 = MIN(K+1, Nr)
# Line 151  C--   Initial residual calculation (with Line 151  C--   Initial residual calculation (with
151           maskP1 = 1. _d 0           maskP1 = 1. _d 0
152           IF ( K .EQ. 1 ) maskM1 = 0. _d 0           IF ( K .EQ. 1 ) maskM1 = 0. _d 0
153           IF ( K .EQ. Nr) maskP1 = 0. _d 0           IF ( K .EQ. Nr) maskP1 = 0. _d 0
154              
155           DO J=1,sNy           DO J=1,sNy
156            DO I=1,sNx            DO I=1,sNx
157             cg3d_r(I,J,K,bi,bj) = cg3d_b(I,J,K,bi,bj) -( 0.             cg3d_r(I,J,K,bi,bj) = cg3d_b(I,J,K,bi,bj) -( 0.
# Line 163  C--   Initial residual calculation (with Line 163  C--   Initial residual calculation (with
163       &     +aV3d(I  ,J  ,Kp1,bi,bj)*cg3d_x(I  ,J  ,Kp1,bi,bj)*maskP1       &     +aV3d(I  ,J  ,Kp1,bi,bj)*cg3d_x(I  ,J  ,Kp1,bi,bj)*maskP1
164       &     +aC3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)       &     +aC3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)
165       &     )       &     )
166             errTile = errTile             errTile(bi,bj) = errTile(bi,bj)
167       &     +cg3d_r(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)       &     +cg3d_r(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)
168             sumRHStile = sumRHStile             sumRHStile(bi,bj) = sumRHStile(bi,bj)
169       &     +cg3d_b(I,J,K,bi,bj)       &     +cg3d_b(I,J,K,bi,bj)
170            ENDDO            ENDDO
171           ENDDO           ENDDO
# Line 175  C--   Initial residual calculation (with Line 175  C--   Initial residual calculation (with
175            ENDDO            ENDDO
176           ENDDO           ENDDO
177          ENDDO          ENDDO
178          err    = err    + errTile  c       err    = err    + errTile(bi,bj)
179          sumRHS = sumRHS + sumRHStile  c       sumRHS = sumRHS + sumRHStile(bi,bj)
180         ENDDO         ENDDO
181        ENDDO        ENDDO
182         CALL EXCH_S3D_RL( cg3d_r, Nr, myThid )         CALL EXCH_S3D_RL( cg3d_r, Nr, myThid )
183  c      CALL EXCH_S3D_RL( cg3d_s, Nr, myThid )  c      CALL EXCH_S3D_RL( cg3d_s, Nr, myThid )
184        _GLOBAL_SUM_R8( sumRHS, myThid )  c     _GLOBAL_SUM_R8( sumRHS, myThid )
185        _GLOBAL_SUM_R8( err   , myThid )  c     _GLOBAL_SUM_R8( err   , myThid )
186           CALL GLOBAL_SUM_TILE_RL( sumRHStile, sumRHS, myThid )
187           CALL GLOBAL_SUM_TILE_RL( errTile,    err,    myThid )
188    
189        IF ( debugLevel .GE. debLevZero ) THEN        IF ( debugLevel .GE. debLevZero ) THEN
190          _BEGIN_MASTER( myThid )          _BEGIN_MASTER( myThid )
# Line 212  CcnhDebugEnds Line 214  CcnhDebugEnds
214  C--    Solve preconditioning equation and update  C--    Solve preconditioning equation and update
215  C--    conjugate direction vector "s".  C--    conjugate direction vector "s".
216  C      Note. On the next to loops over all tiles the inner loop ranges  C      Note. On the next to loops over all tiles the inner loop ranges
217  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
218  C            step. However this entails a bit of gynamastics because we only  C            step. However this entails a bit of gynamastics because we only
219  C            want eta_qrN for the interior points.  C            want eta_qrN for the interior points.
220         eta_qrN = 0. _d 0         eta_qrN = 0. _d 0
221         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
222          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
223           eta_qrNtile = 0. _d 0           eta_qrNtile(bi,bj) = 0. _d 0
224           DO K=1,1           DO K=1,1
225            DO J=1-1,sNy+1            DO J=1-1,sNy+1
226             DO I=1-1,sNx+1             DO I=1-1,sNx+1
227              cg3d_q(I,J,K,bi,bj) =              cg3d_q(I,J,K,bi,bj) =
228       &       zMC(I  ,J  ,K,bi,bj)*cg3d_r(I  ,J  ,K,bi,bj)       &       zMC(I  ,J  ,K,bi,bj)*cg3d_r(I  ,J  ,K,bi,bj)
229             ENDDO             ENDDO
230            ENDDO            ENDDO
# Line 230  C            want eta_qrN for the interi Line 232  C            want eta_qrN for the interi
232           DO K=2,Nr           DO K=2,Nr
233            DO J=1-1,sNy+1            DO J=1-1,sNy+1
234             DO I=1-1,sNx+1             DO I=1-1,sNx+1
235              cg3d_q(I,J,K,bi,bj) =              cg3d_q(I,J,K,bi,bj) =
236       &       zMC(I,J,K,bi,bj)*(cg3d_r(I,J,K  ,bi,bj)       &       zMC(I,J,K,bi,bj)*(cg3d_r(I,J,K  ,bi,bj)
237       &      -zML(I,J,K,bi,bj)*cg3d_q(I,J,K-1,bi,bj))       &      -zML(I,J,K,bi,bj)*cg3d_q(I,J,K-1,bi,bj))
238             ENDDO             ENDDO
# Line 240  C            want eta_qrN for the interi Line 242  C            want eta_qrN for the interi
242  caja      IF (Nr .GT. 1) THEN  caja      IF (Nr .GT. 1) THEN
243  caja       DO J=1-1,sNy+1  caja       DO J=1-1,sNy+1
244  caja        DO I=1-1,sNx+1  caja        DO I=1-1,sNx+1
245  caja         cg3d_q(I,J,K,bi,bj) =  caja         cg3d_q(I,J,K,bi,bj) =
246  caja &        zMC(i,j,k,bi,bj)*(cg3d_r(i,j,k  ,bi,bj)  caja &        zMC(i,j,k,bi,bj)*(cg3d_r(i,j,k  ,bi,bj)
247  caja &       -zML(i,j,k,bi,bj)*cg3d_q(i,j,k-1,bi,bj))  caja &       -zML(i,j,k,bi,bj)*cg3d_q(i,j,k-1,bi,bj))
248  caja        ENDDO  caja        ENDDO
# Line 248  caja       ENDDO Line 250  caja       ENDDO
250  caja      ENDIF  caja      ENDIF
251            DO J=1,sNy            DO J=1,sNy
252             DO I=1,sNx             DO I=1,sNx
253              eta_qrNtile = eta_qrNtile              eta_qrNtile(bi,bj) = eta_qrNtile(bi,bj)
254       &      +cg3d_q(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)       &      +cg3d_q(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)
255             ENDDO             ENDDO
256            ENDDO            ENDDO
# Line 256  caja      ENDIF Line 258  caja      ENDIF
258           DO K=Nr-1,1,-1           DO K=Nr-1,1,-1
259            DO J=1-1,sNy+1            DO J=1-1,sNy+1
260             DO I=1-1,sNx+1             DO I=1-1,sNx+1
261              cg3d_q(I,J,K,bi,bj) =              cg3d_q(I,J,K,bi,bj) =
262       &      cg3d_q(I,J,K,bi,bj)       &      cg3d_q(I,J,K,bi,bj)
263       &      -zMU(I,J,K,bi,bj)*cg3d_q(I,J,K+1,bi,bj)       &      -zMU(I,J,K,bi,bj)*cg3d_q(I,J,K+1,bi,bj)
264             ENDDO             ENDDO
265            ENDDO            ENDDO
266            DO J=1,sNy            DO J=1,sNy
267             DO I=1,sNx             DO I=1,sNx
268              eta_qrNtile = eta_qrNtile              eta_qrNtile(bi,bj) = eta_qrNtile(bi,bj)
269       &      +cg3d_q(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)       &      +cg3d_q(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)
270             ENDDO             ENDDO
271            ENDDO            ENDDO
272           ENDDO           ENDDO
273           eta_qrN = eta_qrN + eta_qrNtile  c        eta_qrN = eta_qrN + eta_qrNtile(bi,bj)
274          ENDDO          ENDDO
275         ENDDO         ENDDO
276  caja  caja
# Line 287  caja    ENDDO Line 289  caja    ENDDO
289  caja   ENDDO  caja   ENDDO
290  caja  caja
291    
292         _GLOBAL_SUM_R8(eta_qrN, myThid)  c      _GLOBAL_SUM_R8(eta_qrN, myThid)
293           CALL GLOBAL_SUM_TILE_RL( eta_qrNtile,eta_qrN,myThid )
294  CcnhDebugStarts  CcnhDebugStarts
295  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' eta_qrN = ',eta_qrN  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' eta_qrN = ',eta_qrN
296  CcnhDebugEnds  CcnhDebugEnds
# Line 302  CcnhDebugEnds Line 305  CcnhDebugEnds
305           DO K=1,Nr           DO K=1,Nr
306            DO J=1-1,sNy+1            DO J=1-1,sNy+1
307             DO I=1-1,sNx+1             DO I=1-1,sNx+1
308              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)
309       &                          + cgBeta*cg3d_s(I,J,K,bi,bj)       &                          + cgBeta*cg3d_s(I,J,K,bi,bj)
310             ENDDO             ENDDO
311            ENDDO            ENDDO
# Line 315  C==    q = A.s Line 318  C==    q = A.s
318         alpha = 0. _d 0         alpha = 0. _d 0
319         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
320          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
321           alphaTile = 0. _d 0           alphaTile(bi,bj) = 0. _d 0
322           IF ( Nr .GT. 1 ) THEN           IF ( Nr .GT. 1 ) THEN
323            DO K=1,1            DO K=1,1
324             DO J=1,sNy             DO J=1,sNy
325              DO I=1,sNx              DO I=1,sNx
326               cg3d_q(I,J,K,bi,bj) =               cg3d_q(I,J,K,bi,bj) =
327       &       aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I-1,J  ,K  ,bi,bj)       &       aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I-1,J  ,K  ,bi,bj)
328       &      +aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I+1,J  ,K  ,bi,bj)       &      +aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I+1,J  ,K  ,bi,bj)
329       &      +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J-1,K  ,bi,bj)       &      +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J-1,K  ,bi,bj)
330       &      +aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J+1,K  ,bi,bj)       &      +aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J+1,K  ,bi,bj)
331       &      +aV3d(I  ,J  ,K+1,bi,bj)*cg3d_s(I  ,J  ,K+1,bi,bj)       &      +aV3d(I  ,J  ,K+1,bi,bj)*cg3d_s(I  ,J  ,K+1,bi,bj)
332       &      +aC3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &      +aC3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)
333               alphaTile = alphaTile               alphaTile(bi,bj) = alphaTile(bi,bj)
334       &                 +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)
335              ENDDO              ENDDO
336             ENDDO             ENDDO
# Line 336  C==    q = A.s Line 339  C==    q = A.s
339            DO K=1,1            DO K=1,1
340             DO J=1,sNy             DO J=1,sNy
341              DO I=1,sNx              DO I=1,sNx
342               cg3d_q(I,J,K,bi,bj) =               cg3d_q(I,J,K,bi,bj) =
343       &       aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I-1,J  ,K  ,bi,bj)       &       aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I-1,J  ,K  ,bi,bj)
344       &      +aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I+1,J  ,K  ,bi,bj)       &      +aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I+1,J  ,K  ,bi,bj)
345       &      +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J-1,K  ,bi,bj)       &      +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J-1,K  ,bi,bj)
346       &      +aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J+1,K  ,bi,bj)       &      +aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J+1,K  ,bi,bj)
347       &      +aC3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &      +aC3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)
348               alphaTile = alphaTile               alphaTile(bi,bj) = alphaTile(bi,bj)
349       &                 +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)
350              ENDDO              ENDDO
351             ENDDO             ENDDO
# Line 351  C==    q = A.s Line 354  C==    q = A.s
354           DO K=2,Nr-1           DO K=2,Nr-1
355            DO J=1,sNy            DO J=1,sNy
356             DO I=1,sNx             DO I=1,sNx
357              cg3d_q(I,J,K,bi,bj) =              cg3d_q(I,J,K,bi,bj) =
358       &      aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I-1,J  ,K  ,bi,bj)       &      aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I-1,J  ,K  ,bi,bj)
359       &     +aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I+1,J  ,K  ,bi,bj)       &     +aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I+1,J  ,K  ,bi,bj)
360       &     +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J-1,K  ,bi,bj)       &     +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J-1,K  ,bi,bj)
# Line 359  C==    q = A.s Line 362  C==    q = A.s
362       &     +aV3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K-1,bi,bj)       &     +aV3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K-1,bi,bj)
363       &     +aV3d(I  ,J  ,K+1,bi,bj)*cg3d_s(I  ,J  ,K+1,bi,bj)       &     +aV3d(I  ,J  ,K+1,bi,bj)*cg3d_s(I  ,J  ,K+1,bi,bj)
364       &     +aC3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &     +aC3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)
365              alphaTile = alphaTile              alphaTile(bi,bj) = alphaTile(bi,bj)
366       &                +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)
367             ENDDO             ENDDO
368            ENDDO            ENDDO
# Line 368  C==    q = A.s Line 371  C==    q = A.s
371            DO K=Nr,Nr            DO K=Nr,Nr
372             DO J=1,sNy             DO J=1,sNy
373              DO I=1,sNx              DO I=1,sNx
374               cg3d_q(I,J,K,bi,bj) =               cg3d_q(I,J,K,bi,bj) =
375       &       aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I-1,J  ,K  ,bi,bj)       &       aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I-1,J  ,K  ,bi,bj)
376       &      +aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I+1,J  ,K  ,bi,bj)       &      +aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I+1,J  ,K  ,bi,bj)
377       &      +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J-1,K  ,bi,bj)       &      +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J-1,K  ,bi,bj)
378       &      +aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J+1,K  ,bi,bj)       &      +aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J+1,K  ,bi,bj)
379       &      +aV3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K-1,bi,bj)       &      +aV3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K-1,bi,bj)
380       &      +aC3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &      +aC3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)
381               alphaTile = alphaTile               alphaTile(bi,bj) = alphaTile(bi,bj)
382       &                 +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)
383              ENDDO              ENDDO
384             ENDDO             ENDDO
385            ENDDO            ENDDO
386           ENDIF           ENDIF
387           alpha = alpha + alphaTile  c        alpha = alpha + alphaTile(bi,bj)
388          ENDDO          ENDDO
389         ENDDO         ENDDO
390         _GLOBAL_SUM_R8(alpha,myThid)  c      _GLOBAL_SUM_R8(alpha,myThid)
391           CALL GLOBAL_SUM_TILE_RL( alphaTile,  alpha,  myThid )
392  CcnhDebugStarts  CcnhDebugStarts
393  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' SUM(s*q)= ',alpha  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' SUM(s*q)= ',alpha
394  CcnhDebugEnds  CcnhDebugEnds
# Line 392  CcnhDebugEnds Line 396  CcnhDebugEnds
396  CcnhDebugStarts  CcnhDebugStarts
397  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' alpha= ',alpha  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' alpha= ',alpha
398  CcnhDebugEnds  CcnhDebugEnds
399        
400  C==    Update solution and residual vectors  C==    Update solution and residual vectors
401  C      Now compute "interior" points.  C      Now compute "interior" points.
402         err = 0. _d 0         err = 0. _d 0
403         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
404          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
405          errTile    = 0. _d 0          errTile(bi,bj)    = 0. _d 0
406           DO K=1,Nr           DO K=1,Nr
407            DO J=1,sNy            DO J=1,sNy
408             DO I=1,sNx             DO I=1,sNx
# Line 406  C      Now compute "interior" points. Line 410  C      Now compute "interior" points.
410       &            +alpha*cg3d_s(I,J,K,bi,bj)       &            +alpha*cg3d_s(I,J,K,bi,bj)
411              cg3d_r(I,J,K,bi,bj)=cg3d_r(I,J,K,bi,bj)              cg3d_r(I,J,K,bi,bj)=cg3d_r(I,J,K,bi,bj)
412       &            -alpha*cg3d_q(I,J,K,bi,bj)       &            -alpha*cg3d_q(I,J,K,bi,bj)
413              errTile = errTile              errTile(bi,bj) = errTile(bi,bj)
414       &             +cg3d_r(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)       &             +cg3d_r(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)
415             ENDDO             ENDDO
416            ENDDO            ENDDO
417           ENDDO           ENDDO
418           err = err + errTile  c        err = err + errTile(bi,bj)
419          ENDDO          ENDDO
420         ENDDO         ENDDO
421    
422         _GLOBAL_SUM_R8( err   , myThid )  c      _GLOBAL_SUM_R8( err   , myThid )
423           CALL GLOBAL_SUM_TILE_RL( errTile,    err,    myThid )
424         err = SQRT(err)         err = SQRT(err)
425         actualIts      = it3d         actualIts      = it3d
426         actualResidual = err         actualResidual = err
Legend:
Removed from v.1.18  
changed lines
  Added in v.1.19
  ViewVC Help
Powered by ViewVC 1.1.22