/[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.4 by adcroft, Mon Mar 27 22:25:43 2000 UTC revision 1.12.20.1 by edhill, Thu Oct 2 18:10:45 2003 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4    #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    
7  #define VERBOSE  #define VERBOSE
8    
9    CBOP
10    C     !ROUTINE: CG3D
11    C     !INTERFACE:
12        SUBROUTINE CG3D(          SUBROUTINE CG3D(  
13         I                cg3d_b,
14         U                cg3d_x,
15         O                firstResidual,
16         O                lastResidual,
17         U                numIters,
18       I                myThid )       I                myThid )
19  C     /==========================================================\  C     !DESCRIPTION: \bv
20  C     | SUBROUTINE CG3D                                          |  C     *==========================================================*
21  C     | o Three-dimensional grid problem conjugate-gradient      |  C     | SUBROUTINE CG3D                                          
22  C     |   inverter (with preconditioner).                        |  C     | o Three-dimensional grid problem conjugate-gradient      
23  C     |==========================================================|  C     |   inverter (with preconditioner).                        
24  C     | Con. grad is an iterative procedure for solving Ax = b.  |  C     *==========================================================*
25  C     | It requires the A be symmetric.                          |  C     | Con. grad is an iterative procedure for solving Ax = b.  
26  C     | This implementation assumes A is a five-diagonal         |  C     | It requires the A be symmetric.                          
27  C     | matrix of the form that arises in the discrete           |  C     | This implementation assumes A is a seven-diagonal          
28  C     | representation of the del^2 operator in a                |  C     | matrix of the form that arises in the discrete            
29  C     | two-dimensional space.                                   |  C     | representation of the del^2 operator in a                
30  C     | Notes:                                                   |  C     | three-dimensional space.                                    
31  C     | ======                                                   |  C     | Notes:                                                    
32  C     | This implementation can support shared-memory            |  C     | ======                                                    
33  C     | multi-threaded execution. In order to do this COMMON     |  C     | This implementation can support shared-memory              
34  C     | blocks are used for many of the arrays - even ones that  |  C     | multi-threaded execution. In order to do this COMMON      
35  C     | are only used for intermedaite results. This design is   |  C     | blocks are used for many of the arrays - even ones that    
36  C     | OK if you want to all the threads to collaborate on      |  C     | are only used for intermedaite results. This design is    
37  C     | solving the same problem. On the other hand if you want  |  C     | OK if you want to all the threads to collaborate on        
38  C     | the threads to solve several different problems          |  C     | solving the same problem. On the other hand if you want    
39  C     | concurrently this implementation will not work.          |  C     | the threads to solve several different problems            
40  C     \==========================================================/  C     | concurrently this implementation will not work.          
41        IMPLICIT NONE  C     *==========================================================*
42    C     \ev
43    
44    C     !USES:
45          IMPLICIT NONE
46  C     === Global data ===  C     === Global data ===
47  #include "SIZE.h"  #include "SIZE.h"
48  #include "EEPARAMS.h"  #include "EEPARAMS.h"
# Line 37  C     === Global data === Line 50  C     === Global data ===
50  #include "GRID.h"  #include "GRID.h"
51  #include "CG3D.h"  #include "CG3D.h"
52    
53    C     !INPUT/OUTPUT PARAMETERS:
54  C     === Routine arguments ===  C     === Routine arguments ===
55  C     myThid - Thread on which I am working.  C     myThid    - Thread on which I am working.
56    C     cg2d_b    - The source term or "right hand side"
57    C     cg2d_x    - The solution
58    C     firstResidual - the initial residual before any iterations
59    C     lastResidual  - the actual residual reached
60    C     numIters  - Entry: the maximum number of iterations allowed
61    C                 Exit:  the actual number of iterations used
62          _RL  cg3d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
63          _RL  cg3d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
64          _RL  firstResidual
65          _RL  lastResidual
66          INTEGER numIters
67        INTEGER myThid        INTEGER myThid
68    
69    
70  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
71    
72    C     !LOCAL VARIABLES:
73  C     === Local variables ====  C     === Local variables ====
74  C     actualIts      - Number of iterations taken  C     actualIts      - Number of iterations taken
75  C     actualResidual - residual  C     actualResidual - residual
76  C     bi          - Block index in X and Y.  C     bi          - Block index in X and Y.
77  C     bj  C     bj
78  C     etaN        - Used in computing search directions  C     eta_qrN     - Used in computing search directions
79  C     etaNM1        suffix N and NM1 denote current and  C     eta_qrNM1     suffix N and NM1 denote current and
80  C     cgBeta        previous iterations respectively.  C     cgBeta        previous iterations respectively.
81  C     alpha    C     alpha  
82  C     sumRHS      - Sum of right-hand-side. Sometimes this is a  C     sumRHS      - Sum of right-hand-side. Sometimes this is a
# Line 64  C     I, J, N     - Loop counters ( N co Line 91  C     I, J, N     - Loop counters ( N co
91        INTEGER I, J, K, it3d        INTEGER I, J, K, it3d
92        INTEGER KM1, KP1        INTEGER KM1, KP1
93        _RL    err        _RL    err
94        _RL    etaN        _RL    eta_qrN
95        _RL    etaNM1        _RL    eta_qrNM1
96        _RL    cgBeta        _RL    cgBeta
97        _RL    alpha        _RL    alpha
98        _RL    sumRHS        _RL    sumRHS
# Line 79  C     I, J, N     - Loop counters ( N co Line 106  C     I, J, N     - Loop counters ( N co
106        INTEGER exchWidthX        INTEGER exchWidthX
107        INTEGER exchWidthY        INTEGER exchWidthY
108        INTEGER myNz        INTEGER myNz
109        _RL     topLevFac        _RL     topLevTerm
110    CEOP
111    
112    ceh3 needs an IF ( useNONHYDROSTATIC ) THEN
113    
114    
115  C--   Initialise inverter  C--   Initialise inverter
116        etaNM1              = 1. D0        eta_qrNM1 = 1. D0
117    
118  C--   Normalise RHS  C--   Normalise RHS
119        rhsMax = 0. _d 0        rhsMax = 0. _d 0
# Line 118  C--   Update overlaps Line 149  C--   Update overlaps
149        _EXCH_XYZ_R8( cg3d_b, myThid )        _EXCH_XYZ_R8( cg3d_b, myThid )
150        _EXCH_XYZ_R8( cg3d_x, myThid )        _EXCH_XYZ_R8( cg3d_x, myThid )
151    
152  #ifdef NONO  C--   Initial residual calculation (with free-Surface term)
 CcnhDebugStarts  
 C--   Initial residual calculation  
       err    = 0. _d 0  
       sumRHS = 0. _d 0  
       DO bj=myByLo(myThid),myByHi(myThid)  
        DO bi=myBxLo(myThid),myBxHi(myThid)  
          DO J=1,sNy  
           DO I=1,sNx  
           alpha = 0. _d 0  
           DO K=1,Nr  
            KM1 = K-1  
            IF ( KM1 .EQ. 0 )    KM1 = 1  
            KP1 = K+1  
            IF ( KP1 .EQ. Nr+1 ) KP1 = 1  
             cg3d_r(I,J,K,bi,bj) = cg3d_b(I,J,K,bi,bj) -( 0.  
      &     +aW3d(I  ,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)  
      &     +aS3d(I  ,J  ,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)  
      &     +aV3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J  ,KM1,bi,bj)  
      &     +aV3d(I  ,J  ,KP1,bi,bj)*cg3d_x(I  ,J  ,KP1,bi,bj)  
      &     -aW3d(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)  
      &     )  
            alpha = alpha  
      &     +cg3d_r(I,J,K,bi,bj)  
            sumRHS = sumRHS  
      &     +cg3d_b(I,J,K,bi,bj)  
           ENDDO  
           err = err + alpha*alpha  
          ENDDO  
         ENDDO  
        ENDDO  
       ENDDO  
       WRITE(6,*) 'DEBUG  mythid, err = ', mythid, SQRT(err)  
       _GLOBAL_SUM_R8( err   , myThid )  
       _GLOBAL_SUM_R8( sumRHS   , myThid )  
       _BEGIN_MASTER( myThid )  
       write(0,*) 'DEBUG cg3d: Sum(rhs) = ',sumRHS  
       _END_MASTER( )  
       actualIts      = 0  
       actualResidual = SQRT(err)  
 C     _BARRIER  
       _BEGIN_MASTER( myThid )  
        WRITE(0,'(A,I6,1PE30.14)') 'DEBUG  CG3D iters, err = ',  
      &  actualIts, actualResidual  
       _END_MASTER( )  
 CcnhDebugEnds  
 #endif  
   
 C--   Initial residual calculation  
153        err    = 0. _d 0        err    = 0. _d 0
154        sumRHS = 0. _d 0        sumRHS = 0. _d 0
155        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
# Line 183  C--   Initial residual calculation Line 159  C--   Initial residual calculation
159           IF ( K .EQ. 1 ) KM1 = 1           IF ( K .EQ. 1 ) KM1 = 1
160           KP1 = K+1           KP1 = K+1
161           IF ( K .EQ. Nr ) KP1 = 1           IF ( K .EQ. Nr ) KP1 = 1
162           topLevFac = 0.           topLevTerm = 0.
163           IF ( K .EQ. 1) topLevFac = 1.           IF ( K .EQ. 1) topLevTerm = freeSurfFac*cg3dNorm*
164         &      (horiVertRatio/gravity)/deltaTMom/deltaTMom
165           DO J=1,sNy           DO J=1,sNy
166            DO I=1,sNx            DO I=1,sNx
167             cg3d_s(I,J,K,bi,bj) = 0.             cg3d_s(I,J,K,bi,bj) = 0.
# Line 201  C--   Initial residual calculation Line 178  C--   Initial residual calculation
178       &     -aS3d(I  ,J+1,K  ,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)       &     -aS3d(I  ,J+1,K  ,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)
179       &     -aV3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)       &     -aV3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)
180       &     -aV3d(I  ,J  ,KP1,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)       &     -aV3d(I  ,J  ,KP1,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)
181       &     -freeSurfFac*_rA(i,j,bi,bj)* horiVertRatio*       &     -topLevTerm*_rA(I,J,bi,bj)*cg3d_x(I,J,K,bi,bj)
      &      cg3d_x(I  ,J  ,K,bi,bj)/deltaTMom/deltaTMom*cg3dNorm  
      &      *topLevFac  
182       &     )       &     )
183             err = err             err = err
184       &     +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)
# Line 242  C     _EXCH_XYZ_R8( cg3d_s, myThid ) Line 217  C     _EXCH_XYZ_R8( cg3d_s, myThid )
217        _GLOBAL_SUM_R8( err   , myThid )        _GLOBAL_SUM_R8( err   , myThid )
218                
219        _BEGIN_MASTER( myThid )        _BEGIN_MASTER( myThid )
220        write(0,'(A,1PE30.14)') ' cg3d: Sum(rhs) = ',sumRHS        write(*,'(A,1P2E22.14)')
221         &     ' cg3d: Sum(rhs),rhsMax = ',sumRHS,rhsMax
222        _END_MASTER( )        _END_MASTER( )
223    
224        actualIts      = 0        actualIts      = 0
225        actualResidual = SQRT(err)        actualResidual = SQRT(err)
226  C     _BARRIER  C     _BARRIER
227        _BEGIN_MASTER( myThid )  c     _BEGIN_MASTER( myThid )
228         WRITE(0,'(A,I6,1PE30.14)') ' CG3D iters, err = ',  c      WRITE(*,'(A,I6,1PE30.14)') ' CG3D iters, err = ',
229       &  actualIts, actualResidual  c    &  actualIts, actualResidual
230        _END_MASTER( )  c     _END_MASTER( )
231          firstResidual=actualResidual
232    
233  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
234        DO 10 it3d=1, cg3dMaxIters        DO 10 it3d=1, cg3dMaxIters
235    
236  CcnhDebugStarts  CcnhDebugStarts
237  #ifdef VERBOSE  #ifdef VERBOSE
238         IF ( mod(it3d-1,10).EQ.0)  c      IF ( mod(it3d-1,10).EQ.0)
239       &   WRITE(0,*) ' CG3D: Iteration ',it3d-1,  c    &   WRITE(*,*) ' CG3D: Iteration ',it3d-1,
240       &      ' residual = ',actualResidual  c    &      ' residual = ',actualResidual
241  #endif  #endif
242  CcnhDebugEnds  CcnhDebugEnds
243         IF ( actualResidual .LT. cg3dTargetResidual ) GOTO 11         IF ( actualResidual .LT. cg3dTargetResidual ) GOTO 11
# Line 269  C--    conjugate direction vector "s". Line 246  C--    conjugate direction vector "s".
246  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
247  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
248  C            step. However this entails a bit of gynamastics because we only  C            step. However this entails a bit of gynamastics because we only
249  C            want etaN for the interior points.  C            want eta_qrN for the interior points.
250         etaN = 0. _d 0         eta_qrN = 0. _d 0
251         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
252          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
253           DO K=1,1           DO K=1,1
# Line 302  caja       ENDDO Line 279  caja       ENDDO
279  caja      ENDIF  caja      ENDIF
280            DO J=1,sNy            DO J=1,sNy
281             DO I=1,sNx             DO I=1,sNx
282              etaN = etaN              eta_qrN = eta_qrN
283       &      +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)
284             ENDDO             ENDDO
285            ENDDO            ENDDO
# Line 317  caja      ENDIF Line 294  caja      ENDIF
294            ENDDO            ENDDO
295            DO J=1,sNy            DO J=1,sNy
296             DO I=1,sNx             DO I=1,sNx
297              etaN = etaN              eta_qrN = eta_qrN
298       &      +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)
299             ENDDO             ENDDO
300            ENDDO            ENDDO
# Line 325  caja      ENDIF Line 302  caja      ENDIF
302          ENDDO          ENDDO
303         ENDDO         ENDDO
304  caja  caja
305  caja  etaN=0.  caja  eta_qrN=0.
306  caja   DO bj=myByLo(myThid),myByHi(myThid)  caja   DO bj=myByLo(myThid),myByHi(myThid)
307  caja    DO bi=myBxLo(myThid),myBxHi(myThid)  caja    DO bi=myBxLo(myThid),myBxHi(myThid)
308  caja     DO K=1,Nr  caja     DO K=1,Nr
309  caja      DO J=1,sNy  caja      DO J=1,sNy
310  caja       DO I=1,sNx  caja       DO I=1,sNx
311  caja        etaN = etaN  caja        eta_qrN = eta_qrN
312  caja &      +cg3d_q(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)  caja &      +cg3d_q(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)
313  caja       ENDDO  caja       ENDDO
314  caja      ENDDO  caja      ENDDO
# Line 340  caja    ENDDO Line 317  caja    ENDDO
317  caja   ENDDO  caja   ENDDO
318  caja  caja
319    
320         _GLOBAL_SUM_R8(etaN, myThid)         _GLOBAL_SUM_R8(eta_qrN, myThid)
321  CcnhDebugStarts  CcnhDebugStarts
322  C      WRITE(0,*) ' CG3D: Iteration ',it3d-1,' etaN = ',etaN  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' eta_qrN = ',eta_qrN
323  CcnhDebugEnds  CcnhDebugEnds
324         cgBeta   = etaN/etaNM1         cgBeta   = eta_qrN/eta_qrNM1
325  CcnhDebugStarts  CcnhDebugStarts
326  C      WRITE(0,*) ' CG3D: Iteration ',it3d-1,' beta = ',cgBeta  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' beta = ',cgBeta
327  CcnhDebugEnds  CcnhDebugEnds
328         etaNM1 = etaN         eta_qrNM1 = eta_qrN
329    
330         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
331          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
# Line 366  CcnhDebugEnds Line 343  CcnhDebugEnds
343  C==    Evaluate laplace operator on conjugate gradient vector  C==    Evaluate laplace operator on conjugate gradient vector
344  C==    q = A.s  C==    q = A.s
345         alpha = 0. _d 0         alpha = 0. _d 0
346           topLevTerm = freeSurfFac*cg3dNorm*
347         &      (horiVertRatio/gravity)/deltaTMom/deltaTMom
348         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
349          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
350           IF ( Nr .GT. 1 ) THEN           IF ( Nr .GT. 1 ) THEN
# Line 383  C==    q = A.s Line 362  C==    q = A.s
362       &      -aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &      -aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)
363       &      -aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &      -aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)
364       &      -aV3d(I  ,J  ,K+1,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &      -aV3d(I  ,J  ,K+1,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)
365       &      -freeSurfFac*_rA(i,j,bi,bj)* horiVertRatio*       &      -topLevTerm*_rA(I,J,bi,bj)*cg3d_s(I,J,K,bi,bj)
      &       cg3d_s(I  ,J  ,K,bi,bj)/deltaTMom/deltaTMom*cg3dNorm  
366               alpha = alpha+cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)               alpha = alpha+cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)
367              ENDDO              ENDDO
368             ENDDO             ENDDO
# Line 402  C==    q = A.s Line 380  C==    q = A.s
380       &      -aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &      -aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)
381       &      -aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &      -aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)
382       &      -aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &      -aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)
383       &      -freeSurfFac*_rA(i,j,bi,bj)* horiVertRatio*       &      -topLevTerm*_rA(I,J,bi,bj)*cg3d_s(I,J,K,bi,bj)
      &       cg3d_s(I  ,J  ,K,bi,bj)/deltaTMom/deltaTMom*cg3dNorm  
384               alpha = alpha+cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)               alpha = alpha+cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)
385              ENDDO              ENDDO
386             ENDDO             ENDDO
# Line 453  C==    q = A.s Line 430  C==    q = A.s
430         ENDDO         ENDDO
431         _GLOBAL_SUM_R8(alpha,myThid)         _GLOBAL_SUM_R8(alpha,myThid)
432  CcnhDebugStarts  CcnhDebugStarts
433  C      WRITE(0,*) ' CG3D: Iteration ',it3d-1,' SUM(s*q)= ',alpha  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' SUM(s*q)= ',alpha
434  CcnhDebugEnds  CcnhDebugEnds
435         alpha = etaN/alpha         alpha = eta_qrN/alpha
436  CcnhDebugStarts  CcnhDebugStarts
437  C      WRITE(0,*) ' CG3D: Iteration ',it3d-1,' alpha= ',alpha  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' alpha= ',alpha
438  CcnhDebugEnds  CcnhDebugEnds
439            
440  C==    Update solution and residual vectors  C==    Update solution and residual vectors
# Line 514  C--   Un-normalise the answer Line 491  C--   Un-normalise the answer
491        ENDDO        ENDDO
492    
493  Cadj  _EXCH_XYZ_R8(cg3d_x, myThid )  Cadj  _EXCH_XYZ_R8(cg3d_x, myThid )
494        _BEGIN_MASTER( myThid )  c     _BEGIN_MASTER( myThid )
495         WRITE(0,'(A,I6,1PE30.14)') ' CG3D iters, err = ',  c      WRITE(*,'(A,I6,1PE30.14)') ' CG3D iters, err = ',
496       &  actualIts, actualResidual  c    &  actualIts, actualResidual
497        _END_MASTER( )  c     _END_MASTER( )
498          lastResidual=actualResidual
499  CcnhDebugStarts        numIters=actualIts
 C     CALL PLOT_FIELD_XYRL( cg2d_x, 'CALC_MOM_RHS CG2D_X' , 1, myThid )  
 C     err    = 0. _d 0  
 C     DO bj=myByLo(myThid),myByHi(myThid)  
 C      DO bi=myBxLo(myThid),myBxHi(myThid)  
 C       DO J=1,sNy  
 C        DO I=1,sNx  
 C         cg2d_r(I,J,bi,bj) = cg2d_b(I,J,bi,bj) -  
 C    &    (aW2d(I  ,J  ,bi,bj)*cg2d_x(I-1,J  ,bi,bj)  
 C    &    +aW2d(I+1,J  ,bi,bj)*cg2d_x(I+1,J  ,bi,bj)  
 C    &    +aS2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J-1,bi,bj)  
 C    &    +aS2d(I  ,J+1,bi,bj)*cg2d_x(I  ,J+1,bi,bj)  
 C    &    -aW2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)  
 C    &    -aW2d(I+1,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)  
 C    &    -aS2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)  
 C    &    -aS2d(I  ,J+1,bi,bj)*cg2d_x(I  ,J  ,bi,bj))  
 C         err            = err            +  
 C    &     cg2d_r(I,J,bi,bj)*cg2d_r(I,J,bi,bj)  
 C        ENDDO  
 C       ENDDO  
 C      ENDDO  
 C     ENDDO  
 C     _GLOBAL_SUM_R8( err   , myThid )  
 C     write(0,*) 'cg2d: Ax - b = ',SQRT(err)  
 CcnhDebugEnds  
500    
501  #endif /* ALLOW_NONHYDROSTATIC */  #endif /* ALLOW_NONHYDROSTATIC */
502    
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.12.20.1
  ViewVC Help
Powered by ViewVC 1.1.22