/[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.12.20.1 by edhill, Thu Oct 2 18:10:45 2003 UTC revision 1.17 by jmc, Tue Dec 20 20:31:28 2005 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
 #include "PACKAGES_CONFIG.h"  
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    
 #define VERBOSE  
   
6  CBOP  CBOP
7  C     !ROUTINE: CG3D  C     !ROUTINE: CG3D
8  C     !INTERFACE:  C     !INTERFACE:
# Line 53  C     === Global data === Line 50  C     === Global data ===
50  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
51  C     === Routine arguments ===  C     === Routine arguments ===
52  C     myThid    - Thread on which I am working.  C     myThid    - Thread on which I am working.
53  C     cg2d_b    - The source term or "right hand side"  C     cg3d_b    - The source term or "right hand side"
54  C     cg2d_x    - The solution  C     cg3d_x    - The solution
55  C     firstResidual - the initial residual before any iterations  C     firstResidual - the initial residual before any iterations
56  C     lastResidual  - the actual residual reached  C     lastResidual  - the actual residual reached
57  C     numIters  - Entry: the maximum number of iterations allowed  C     numIters  - Entry: the maximum number of iterations allowed
# Line 84  C                   useful debuggin/trou Line 81  C                   useful debuggin/trou
81  C                   For neumann problems sumRHS needs to be ~0.  C                   For neumann problems sumRHS needs to be ~0.
82  C                   or they converge at a non-zero residual.  C                   or they converge at a non-zero residual.
83  C     err         - Measure of residual of Ax - b, usually the norm.  C     err         - Measure of residual of Ax - b, usually the norm.
84  C     I, J, 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    err        _RL    maskM1, maskP1
91        _RL    eta_qrN        _RL    err, errTile
92          _RL    eta_qrN, eta_qrNtile
93        _RL    eta_qrNM1        _RL    eta_qrNM1
94        _RL    cgBeta        _RL    cgBeta
95        _RL    alpha        _RL    alpha , alphaTile
96        _RL    sumRHS        _RL    sumRHS, sumRHStile
97        _RL    rhsMax        _RL    rhsMax
98        _RL    rhsNorm        _RL    rhsNorm
   
       INTEGER OLw  
       INTEGER OLe  
       INTEGER OLn  
       INTEGER OLs  
       INTEGER exchWidthX  
       INTEGER exchWidthY  
       INTEGER myNz  
       _RL     topLevTerm  
99  CEOP  CEOP
100    
 ceh3 needs an IF ( useNONHYDROSTATIC ) THEN  
   
101    
102  C--   Initialise inverter  C--   Initialise inverter
103        eta_qrNM1 = 1. D0        eta_qrNM1 = 1. D0
# Line 123  C--   Normalise RHS Line 110  C--   Normalise RHS
110           DO J=1,sNy           DO J=1,sNy
111            DO I=1,sNx            DO I=1,sNx
112             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
113         &                          * maskC(I,J,K,bi,bj)
114             rhsMax = MAX(ABS(cg3d_b(I,J,K,bi,bj)),rhsMax)             rhsMax = MAX(ABS(cg3d_b(I,J,K,bi,bj)),rhsMax)
115            ENDDO            ENDDO
116           ENDDO           ENDDO
# Line 146  C--   Normalise RHS Line 134  C--   Normalise RHS
134        ENDDO        ENDDO
135    
136  C--   Update overlaps  C--   Update overlaps
137        _EXCH_XYZ_R8( cg3d_b, myThid )  c     _EXCH_XYZ_R8( cg3d_b, myThid )
138        _EXCH_XYZ_R8( cg3d_x, myThid )        _EXCH_XYZ_R8( cg3d_x, myThid )
139    
140  C--   Initial residual calculation (with free-Surface term)  C--   Initial residual calculation (with free-Surface term)
# Line 154  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
146            sumRHStile = 0. _d 0
147          DO K=1,Nr          DO K=1,Nr
148           KM1 = K-1           Km1 = MAX(K-1, 1 )
149           IF ( K .EQ. 1 ) KM1 = 1           Kp1 = MIN(K+1, Nr)
150           KP1 = K+1           maskM1 = 1. _d 0
151           IF ( K .EQ. Nr ) KP1 = 1           maskP1 = 1. _d 0
152           topLevTerm = 0.           IF ( K .EQ. 1 ) maskM1 = 0. _d 0
153           IF ( K .EQ. 1) topLevTerm = freeSurfFac*cg3dNorm*           IF ( K .EQ. Nr) maskP1 = 0. _d 0
154       &      (horiVertRatio/gravity)/deltaTMom/deltaTMom            
155           DO J=1,sNy           DO J=1,sNy
156            DO I=1,sNx            DO I=1,sNx
            cg3d_s(I,J,K,bi,bj) = 0.  
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.
158       &     +aW3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I-1,J  ,K  ,bi,bj)       &     +aW3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I-1,J  ,K  ,bi,bj)
159       &     +aW3d(I+1,J  ,K  ,bi,bj)*cg3d_x(I+1,J  ,K  ,bi,bj)       &     +aW3d(I+1,J  ,K  ,bi,bj)*cg3d_x(I+1,J  ,K  ,bi,bj)
160       &     +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J-1,K  ,bi,bj)       &     +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J-1,K  ,bi,bj)
161       &     +aS3d(I  ,J+1,K  ,bi,bj)*cg3d_x(I  ,J+1,K  ,bi,bj)       &     +aS3d(I  ,J+1,K  ,bi,bj)*cg3d_x(I  ,J+1,K  ,bi,bj)
162       &     +aV3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J  ,KM1,bi,bj)       &     +aV3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J  ,Km1,bi,bj)*maskM1
163       &     +aV3d(I  ,J  ,KP1,bi,bj)*cg3d_x(I  ,J  ,KP1,bi,bj)       &     +aV3d(I  ,J  ,Kp1,bi,bj)*cg3d_x(I  ,J  ,Kp1,bi,bj)*maskP1
164       &     -aW3d(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)
      &     -aW3d(I+1,J  ,K  ,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)  
      &     -aS3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)  
      &     -aS3d(I  ,J+1,K  ,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)  
      &     -aV3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)  
      &     -aV3d(I  ,J  ,KP1,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)  
      &     -topLevTerm*_rA(I,J,bi,bj)*cg3d_x(I,J,K,bi,bj)  
165       &     )       &     )
166             err = err             errTile = errTile
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             sumRHS = sumRHS             sumRHStile = sumRHStile
169       &     +cg3d_b(I,J,K,bi,bj)       &     +cg3d_b(I,J,K,bi,bj)
170            ENDDO            ENDDO
171           ENDDO           ENDDO
172             DO J=1-1,sNy+1
173              DO I=1-1,sNx+1
174               cg3d_s(I,J,K,bi,bj) = 0.
175              ENDDO
176             ENDDO
177          ENDDO          ENDDO
178            err    = err    + errTile
179            sumRHS = sumRHS + sumRHStile
180         ENDDO         ENDDO
181        ENDDO        ENDDO
182  C     _EXCH_XYZ_R8( cg3d_r, myThid )         CALL EXCH_S3D_RL( cg3d_r, myThid )
183         OLw        = 1  c      CALL EXCH_S3D_RL( cg3d_s, myThid )
        OLe        = 1  
        OLn        = 1  
        OLs        = 1  
        exchWidthX = 1  
        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 )  
 C     _EXCH_XYZ_R8( cg3d_s, myThid )  
        OLw        = 1  
        OLe        = 1  
        OLn        = 1  
        OLs        = 1  
        exchWidthX = 1  
        exchWidthY = 1  
        myNz       = Nr  
        CALL EXCH_RL( cg3d_s,  
      I            OLw, OLe, OLs, OLn, myNz,  
      I            exchWidthX, exchWidthY,  
      I            FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )  
184        _GLOBAL_SUM_R8( sumRHS, myThid )        _GLOBAL_SUM_R8( sumRHS, myThid )
185        _GLOBAL_SUM_R8( err   , myThid )        _GLOBAL_SUM_R8( err   , myThid )
186                
187        _BEGIN_MASTER( myThid )        IF ( debugLevel .GE. debLevZero ) THEN
188        write(*,'(A,1P2E22.14)')          _BEGIN_MASTER( myThid )
189            write(standardmessageunit,'(A,1P2E22.14)')
190       &     ' cg3d: Sum(rhs),rhsMax = ',sumRHS,rhsMax       &     ' cg3d: Sum(rhs),rhsMax = ',sumRHS,rhsMax
191        _END_MASTER( )          _END_MASTER( myThid )
192          ENDIF
193    
194        actualIts      = 0        actualIts      = 0
195        actualResidual = SQRT(err)        actualResidual = SQRT(err)
# Line 227  C     _BARRIER Line 197  C     _BARRIER
197  c     _BEGIN_MASTER( myThid )  c     _BEGIN_MASTER( myThid )
198  c      WRITE(*,'(A,I6,1PE30.14)') ' CG3D iters, err = ',  c      WRITE(*,'(A,I6,1PE30.14)') ' CG3D iters, err = ',
199  c    &  actualIts, actualResidual  c    &  actualIts, actualResidual
200  c     _END_MASTER( )  c     _END_MASTER( myThid )
201        firstResidual=actualResidual        firstResidual=actualResidual
202    
203  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
204        DO 10 it3d=1, cg3dMaxIters        DO 10 it3d=1, numIters
205    
206  CcnhDebugStarts  CcnhDebugStarts
 #ifdef VERBOSE  
207  c      IF ( mod(it3d-1,10).EQ.0)  c      IF ( mod(it3d-1,10).EQ.0)
208  c    &   WRITE(*,*) ' CG3D: Iteration ',it3d-1,  c    &   WRITE(*,*) ' CG3D: Iteration ',it3d-1,
209  c    &      ' residual = ',actualResidual  c    &      ' residual = ',actualResidual
 #endif  
210  CcnhDebugEnds  CcnhDebugEnds
211         IF ( actualResidual .LT. cg3dTargetResidual ) GOTO 11         IF ( actualResidual .LT. cg3dTargetResidual ) GOTO 11
212  C--    Solve preconditioning equation and update  C--    Solve preconditioning equation and update
# Line 250  C            want eta_qrN for the interi Line 218  C            want eta_qrN for the interi
218         eta_qrN = 0. _d 0         eta_qrN = 0. _d 0
219         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
220          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
221             eta_qrNtile = 0. _d 0
222           DO K=1,1           DO K=1,1
223            DO J=1-1,sNy+1            DO J=1-1,sNy+1
224             DO I=1-1,sNx+1             DO I=1-1,sNx+1
# Line 279  caja       ENDDO Line 248  caja       ENDDO
248  caja      ENDIF  caja      ENDIF
249            DO J=1,sNy            DO J=1,sNy
250             DO I=1,sNx             DO I=1,sNx
251              eta_qrN = eta_qrN              eta_qrNtile = eta_qrNtile
252       &      +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)
253             ENDDO             ENDDO
254            ENDDO            ENDDO
# Line 294  caja      ENDIF Line 263  caja      ENDIF
263            ENDDO            ENDDO
264            DO J=1,sNy            DO J=1,sNy
265             DO I=1,sNx             DO I=1,sNx
266              eta_qrN = eta_qrN              eta_qrNtile = eta_qrNtile
267       &      +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)
268             ENDDO             ENDDO
269            ENDDO            ENDDO
270           ENDDO           ENDDO
271             eta_qrN = eta_qrN + eta_qrNtile
272          ENDDO          ENDDO
273         ENDDO         ENDDO
274  caja  caja
# Line 343  CcnhDebugEnds Line 313  CcnhDebugEnds
313  C==    Evaluate laplace operator on conjugate gradient vector  C==    Evaluate laplace operator on conjugate gradient vector
314  C==    q = A.s  C==    q = A.s
315         alpha = 0. _d 0         alpha = 0. _d 0
        topLevTerm = freeSurfFac*cg3dNorm*  
      &      (horiVertRatio/gravity)/deltaTMom/deltaTMom  
316         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
317          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
318             alphaTile = 0. _d 0
319           IF ( Nr .GT. 1 ) THEN           IF ( Nr .GT. 1 ) THEN
320            DO K=1,1            DO K=1,1
321             DO J=1,sNy             DO J=1,sNy
# Line 357  C==    q = A.s Line 326  C==    q = A.s
326       &      +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)
327       &      +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)
328       &      +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)
329       &      -aW3d(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)
330       &      -aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)               alphaTile = alphaTile
331       &      -aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &                 +cg3d_s(I,J,K,bi,bj)*cg3d_q(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)  
332              ENDDO              ENDDO
333             ENDDO             ENDDO
334            ENDDO            ENDDO
# Line 376  C==    q = A.s Line 341  C==    q = A.s
341       &      +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)
342       &      +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)
343       &      +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)
344       &      -aW3d(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)
345       &      -aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)               alphaTile = alphaTile
346       &      -aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &                 +cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)
      &      -aS3d(I  ,J+1,K  ,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)  
347              ENDDO              ENDDO
348             ENDDO             ENDDO
349            ENDDO            ENDDO
# Line 396  C==    q = A.s Line 358  C==    q = A.s
358       &     +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)
359       &     +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)
360       &     +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)
361       &     -aW3d(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)
362       &     -aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)              alphaTile = alphaTile
363       &     -aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &                +cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)
      &     -aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)  
      &     -aV3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)  
      &     -aV3d(I  ,J  ,K+1,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)  
364             ENDDO             ENDDO
365            ENDDO            ENDDO
366           ENDDO           ENDDO
# Line 416  C==    q = A.s Line 374  C==    q = A.s
374       &      +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)
375       &      +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)
376       &      +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)
377       &      -aW3d(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)
378       &      -aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)               alphaTile = alphaTile
379       &      -aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &                 +cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)
      &      -aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)  
      &      -aV3d(I  ,J  ,K  ,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)  
380              ENDDO              ENDDO
381             ENDDO             ENDDO
382            ENDDO            ENDDO
383           ENDIF           ENDIF
384             alpha = alpha + alphaTile
385          ENDDO          ENDDO
386         ENDDO         ENDDO
387         _GLOBAL_SUM_R8(alpha,myThid)         _GLOBAL_SUM_R8(alpha,myThid)
# Line 442  C      Now compute "interior" points. Line 398  C      Now compute "interior" points.
398         err = 0. _d 0         err = 0. _d 0
399         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
400          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
401            errTile    = 0. _d 0
402           DO K=1,Nr           DO K=1,Nr
403            DO J=1,sNy            DO J=1,sNy
404             DO I=1,sNx             DO I=1,sNx
# Line 449  C      Now compute "interior" points. Line 406  C      Now compute "interior" points.
406       &            +alpha*cg3d_s(I,J,K,bi,bj)       &            +alpha*cg3d_s(I,J,K,bi,bj)
407              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)
408       &            -alpha*cg3d_q(I,J,K,bi,bj)       &            -alpha*cg3d_q(I,J,K,bi,bj)
409              err = err+cg3d_r(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)              errTile = errTile
410         &             +cg3d_r(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)
411             ENDDO             ENDDO
412            ENDDO            ENDDO
413           ENDDO           ENDDO
414             err = err + errTile
415          ENDDO          ENDDO
416         ENDDO         ENDDO
417    
# Line 461  C      Now compute "interior" points. Line 420  C      Now compute "interior" points.
420         actualIts      = it3d         actualIts      = it3d
421         actualResidual = err         actualResidual = err
422         IF ( actualResidual .LT. cg3dTargetResidual ) GOTO 11         IF ( actualResidual .LT. cg3dTargetResidual ) GOTO 11
423  C      _EXCH_XYZ_R8(cg3d_r, myThid )         CALL EXCH_S3D_RL( cg3d_r, 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, myNz,  
      I             exchWidthX, exchWidthY,  
      I             FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )  
424    
425     10 CONTINUE     10 CONTINUE
426     11 CONTINUE     11 CONTINUE
# Line 494  Cadj  _EXCH_XYZ_R8(cg3d_x, myThid ) Line 442  Cadj  _EXCH_XYZ_R8(cg3d_x, myThid )
442  c     _BEGIN_MASTER( myThid )  c     _BEGIN_MASTER( myThid )
443  c      WRITE(*,'(A,I6,1PE30.14)') ' CG3D iters, err = ',  c      WRITE(*,'(A,I6,1PE30.14)') ' CG3D iters, err = ',
444  c    &  actualIts, actualResidual  c    &  actualIts, actualResidual
445  c     _END_MASTER( )  c     _END_MASTER( myThid )
446        lastResidual=actualResidual        lastResidual=actualResidual
447        numIters=actualIts        numIters=actualIts
448    
Legend:
Removed from v.1.12.20.1  
changed lines
  Added in v.1.17
  ViewVC Help
Powered by ViewVC 1.1.22