/[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.2 by adcroft, Tue May 18 17:36:55 1999 UTC revision 1.16 by jmc, Tue Nov 8 06:18:10 2005 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    
6  #define VERBOSE  CBOP
7    C     !ROUTINE: CG3D
8    C     !INTERFACE:
9        SUBROUTINE CG3D(          SUBROUTINE CG3D(  
10         I                cg3d_b,
11         U                cg3d_x,
12         O                firstResidual,
13         O                lastResidual,
14         U                numIters,
15       I                myThid )       I                myThid )
16  C     /==========================================================\  C     !DESCRIPTION: \bv
17  C     | SUBROUTINE CG3D                                          |  C     *==========================================================*
18  C     | o Three-dimensional grid problem conjugate-gradient      |  C     | SUBROUTINE CG3D                                          
19  C     |   inverter (with preconditioner).                        |  C     | o Three-dimensional grid problem conjugate-gradient      
20  C     |==========================================================|  C     |   inverter (with preconditioner).                        
21  C     | Con. grad is an iterative procedure for solving Ax = b.  |  C     *==========================================================*
22  C     | It requires the A be symmetric.                          |  C     | Con. grad is an iterative procedure for solving Ax = b.  
23  C     | This implementation assumes A is a five-diagonal         |  C     | It requires the A be symmetric.                          
24  C     | matrix of the form that arises in the discrete           |  C     | This implementation assumes A is a seven-diagonal          
25  C     | representation of the del^2 operator in a                |  C     | matrix of the form that arises in the discrete            
26  C     | two-dimensional space.                                   |  C     | representation of the del^2 operator in a                
27  C     | Notes:                                                   |  C     | three-dimensional space.                                    
28  C     | ======                                                   |  C     | Notes:                                                    
29  C     | This implementation can support shared-memory            |  C     | ======                                                    
30  C     | multi-threaded execution. In order to do this COMMON     |  C     | This implementation can support shared-memory              
31  C     | blocks are used for many of the arrays - even ones that  |  C     | multi-threaded execution. In order to do this COMMON      
32  C     | are only used for intermedaite results. This design is   |  C     | blocks are used for many of the arrays - even ones that    
33  C     | OK if you want to all the threads to collaborate on      |  C     | are only used for intermedaite results. This design is    
34  C     | solving the same problem. On the other hand if you want  |  C     | OK if you want to all the threads to collaborate on        
35  C     | the threads to solve several different problems          |  C     | solving the same problem. On the other hand if you want    
36  C     | concurrently this implementation will not work.          |  C     | the threads to solve several different problems            
37  C     \==========================================================/  C     | concurrently this implementation will not work.          
38        IMPLICIT NONE  C     *==========================================================*
39    C     \ev
40    
41    C     !USES:
42          IMPLICIT NONE
43  C     === Global data ===  C     === Global data ===
44  #include "SIZE.h"  #include "SIZE.h"
45  #include "EEPARAMS.h"  #include "EEPARAMS.h"
# Line 37  C     === Global data === Line 47  C     === Global data ===
47  #include "GRID.h"  #include "GRID.h"
48  #include "CG3D.h"  #include "CG3D.h"
49    
50    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"
54    C     cg2d_x    - The solution
55    C     firstResidual - the initial residual before any iterations
56    C     lastResidual  - the actual residual reached
57    C     numIters  - Entry: the maximum number of iterations allowed
58    C                 Exit:  the actual number of iterations used
59          _RL  cg3d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
60          _RL  cg3d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
61          _RL  firstResidual
62          _RL  lastResidual
63          INTEGER numIters
64        INTEGER myThid        INTEGER myThid
65    
66    
67    #ifdef ALLOW_NONHYDROSTATIC
68    
69    C     !LOCAL VARIABLES:
70  C     === Local variables ====  C     === Local variables ====
71  C     actualIts      - Number of iterations taken  C     actualIts      - Number of iterations taken
72  C     actualResidual - residual  C     actualResidual - residual
73  C     bi          - Block index in X and Y.  C     bi          - Block index in X and Y.
74  C     bj  C     bj
75  C     etaN        - Used in computing search directions  C     eta_qrN     - Used in computing search directions
76  C     etaNM1        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
# Line 61  C     I, J, N     - Loop counters ( N co Line 87  C     I, J, N     - Loop counters ( N co
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    err, errTile
91        _RL    etaN        _RL    eta_qrN, eta_qrNtile
92        _RL    etaNM1        _RL    eta_qrNM1
93        _RL    cgBeta        _RL    cgBeta
94        _RL    alpha        _RL    alpha , alphaTile
95        _RL    sumRHS        _RL    sumRHS, sumRHStile
96        _RL    rhsMax        _RL    rhsMax
97        _RL    rhsNorm        _RL    rhsNorm
98          _RL    topLevTerm
99        INTEGER OLw  CEOP
       INTEGER OLe  
       INTEGER OLn  
       INTEGER OLs  
       INTEGER exchWidthX  
       INTEGER exchWidthY  
       INTEGER myNz  
       _RL     topLevFac  
100    
101    
 CcnhDebugStarts  
       CHARACTER*(MAX_LEN_FNAM) suff  
 CcnhDebugEnds  
   
 #ifdef ALLOW_NONHYDROSTATIC  
   
102  C--   Initialise inverter  C--   Initialise inverter
103        etaNM1              = 1. D0        eta_qrNM1 = 1. D0
104    
105  C--   Normalise RHS  C--   Normalise RHS
106        rhsMax = 0. _d 0        rhsMax = 0. _d 0
# Line 120  C--   Normalise RHS Line 133  C--   Normalise RHS
133        ENDDO        ENDDO
134    
135  C--   Update overlaps  C--   Update overlaps
136        _EXCH_XYZ_R8( cg3d_b, myThid )  c     _EXCH_XYZ_R8( cg3d_b, myThid )
137        _EXCH_XYZ_R8( cg3d_x, myThid )        _EXCH_XYZ_R8( cg3d_x, myThid )
138    
139  #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  
140        err    = 0. _d 0        err    = 0. _d 0
141        sumRHS = 0. _d 0        sumRHS = 0. _d 0
142        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
143         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
144            errTile    = 0. _d 0
145            sumRHStile = 0. _d 0
146          DO K=1,Nr          DO K=1,Nr
147           KM1 = K-1           KM1 = K-1
148           IF ( K .EQ. 1 ) KM1 = 1           IF ( K .EQ. 1 ) KM1 = 1
149           KP1 = K+1           KP1 = K+1
150           IF ( K .EQ. Nr ) KP1 = 1           IF ( K .EQ. Nr ) KP1 = 1
151           topLevFac = 0.           topLevTerm = 0.
152           IF ( K .EQ. 1) topLevFac = 1.           IF ( K .EQ. 1) topLevTerm = freeSurfFac*cg3dNorm*
153         &      (horiVertRatio/gravity)/deltaTMom/deltaTMom
154           DO J=1,sNy           DO J=1,sNy
155            DO I=1,sNx            DO I=1,sNx
            cg3d_s(I,J,K,bi,bj) = 0.  
156             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.
157       &     +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)
158       &     +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)
# Line 206  C--   Initial residual calculation Line 166  C--   Initial residual calculation
166       &     -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)
167       &     -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)
168       &     -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)
169       &     -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  
170       &     )       &     )
171             err = err             errTile = errTile
172       &     +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)
173             sumRHS = sumRHS             sumRHStile = sumRHStile
174       &     +cg3d_b(I,J,K,bi,bj)       &     +cg3d_b(I,J,K,bi,bj)
175            ENDDO            ENDDO
176           ENDDO           ENDDO
177             DO J=1-1,sNy+1
178              DO I=1-1,sNx+1
179               cg3d_s(I,J,K,bi,bj) = 0.
180              ENDDO
181             ENDDO
182          ENDDO          ENDDO
183            err    = err    + errTile
184            sumRHS = sumRHS + sumRHStile
185         ENDDO         ENDDO
186        ENDDO        ENDDO
187  C     _EXCH_XYZ_R8( cg3d_r, myThid )  C     _EXCH_XYZ_R8( cg3d_r, myThid )
188         OLw        = 1         CALL EXCH_S3D_RL( cg3d_r, 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 )  
189  C     _EXCH_XYZ_R8( cg3d_s, myThid )  C     _EXCH_XYZ_R8( cg3d_s, myThid )
190         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_s,  
      I            OLw, OLe, OLs, OLn, myNz,  
      I            exchWidthX, exchWidthY,  
      I            FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )  
191        _GLOBAL_SUM_R8( sumRHS, myThid )        _GLOBAL_SUM_R8( sumRHS, myThid )
192        _GLOBAL_SUM_R8( err   , myThid )        _GLOBAL_SUM_R8( err   , myThid )
193                
194        _BEGIN_MASTER( myThid )        IF ( debugLevel .GE. debLevZero ) THEN
195        write(0,'(A,1PE30.14)') ' cg3d: Sum(rhs) = ',sumRHS          _BEGIN_MASTER( myThid )
196        _END_MASTER( )          write(standardmessageunit,'(A,1P2E22.14)')
197         &     ' cg3d: Sum(rhs),rhsMax = ',sumRHS,rhsMax
198            _END_MASTER( myThid )
199          ENDIF
200    
201        actualIts      = 0        actualIts      = 0
202        actualResidual = SQRT(err)        actualResidual = SQRT(err)
203  C     _BARRIER  C     _BARRIER
204        _BEGIN_MASTER( myThid )  c     _BEGIN_MASTER( myThid )
205         WRITE(0,'(A,I6,1PE30.14)') ' CG3D iters, err = ',  c      WRITE(*,'(A,I6,1PE30.14)') ' CG3D iters, err = ',
206       &  actualIts, actualResidual  c    &  actualIts, actualResidual
207        _END_MASTER( )  c     _END_MASTER( myThid )
208          firstResidual=actualResidual
209    
210  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
211        DO 10 it3d=1, cg3dMaxIters        DO 10 it3d=1, cg3dMaxIters
212    
213  CcnhDebugStarts  CcnhDebugStarts
214  #ifdef VERBOSE  c      IF ( mod(it3d-1,10).EQ.0)
215         IF ( mod(it3d-1,10).EQ.0)  c    &   WRITE(*,*) ' CG3D: Iteration ',it3d-1,
216       &   WRITE(0,*) ' CG3D: Iteration ',it3d-1,  c    &      ' residual = ',actualResidual
      &      ' residual = ',actualResidual  
 #endif  
217  CcnhDebugEnds  CcnhDebugEnds
218         IF ( actualResidual .LT. cg3dTargetResidual ) GOTO 11         IF ( actualResidual .LT. cg3dTargetResidual ) GOTO 11
219  C--    Solve preconditioning equation and update  C--    Solve preconditioning equation and update
# Line 274  C--    conjugate direction vector "s". Line 221  C--    conjugate direction vector "s".
221  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
222  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
223  C            step. However this entails a bit of gynamastics because we only  C            step. However this entails a bit of gynamastics because we only
224  C            want etaN for the interior points.  C            want eta_qrN for the interior points.
225         etaN = 0. _d 0         eta_qrN = 0. _d 0
226         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
227          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
228             eta_qrNtile = 0. _d 0
229           DO K=1,1           DO K=1,1
230            DO J=1-1,sNy+1            DO J=1-1,sNy+1
231             DO I=1-1,sNx+1             DO I=1-1,sNx+1
# Line 307  caja       ENDDO Line 255  caja       ENDDO
255  caja      ENDIF  caja      ENDIF
256            DO J=1,sNy            DO J=1,sNy
257             DO I=1,sNx             DO I=1,sNx
258              etaN = etaN              eta_qrNtile = eta_qrNtile
259       &      +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)
260             ENDDO             ENDDO
261            ENDDO            ENDDO
# Line 322  caja      ENDIF Line 270  caja      ENDIF
270            ENDDO            ENDDO
271            DO J=1,sNy            DO J=1,sNy
272             DO I=1,sNx             DO I=1,sNx
273              etaN = etaN              eta_qrNtile = eta_qrNtile
274       &      +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)
275             ENDDO             ENDDO
276            ENDDO            ENDDO
277           ENDDO           ENDDO
278             eta_qrN = eta_qrN + eta_qrNtile
279          ENDDO          ENDDO
280         ENDDO         ENDDO
281  caja  caja
282  caja  etaN=0.  caja  eta_qrN=0.
283  caja   DO bj=myByLo(myThid),myByHi(myThid)  caja   DO bj=myByLo(myThid),myByHi(myThid)
284  caja    DO bi=myBxLo(myThid),myBxHi(myThid)  caja    DO bi=myBxLo(myThid),myBxHi(myThid)
285  caja     DO K=1,Nr  caja     DO K=1,Nr
286  caja      DO J=1,sNy  caja      DO J=1,sNy
287  caja       DO I=1,sNx  caja       DO I=1,sNx
288  caja        etaN = etaN  caja        eta_qrN = eta_qrN
289  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)
290  caja       ENDDO  caja       ENDDO
291  caja      ENDDO  caja      ENDDO
# Line 345  caja    ENDDO Line 294  caja    ENDDO
294  caja   ENDDO  caja   ENDDO
295  caja  caja
296    
297         _GLOBAL_SUM_R8(etaN, myThid)         _GLOBAL_SUM_R8(eta_qrN, myThid)
298  CcnhDebugStarts  CcnhDebugStarts
299  C      WRITE(0,*) ' CG3D: Iteration ',it3d-1,' etaN = ',etaN  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' eta_qrN = ',eta_qrN
300  CcnhDebugEnds  CcnhDebugEnds
301         cgBeta   = etaN/etaNM1         cgBeta   = eta_qrN/eta_qrNM1
302  CcnhDebugStarts  CcnhDebugStarts
303  C      WRITE(0,*) ' CG3D: Iteration ',it3d-1,' beta = ',cgBeta  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' beta = ',cgBeta
304  CcnhDebugEnds  CcnhDebugEnds
305         etaNM1 = etaN         eta_qrNM1 = eta_qrN
306    
307         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
308          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
# Line 371  CcnhDebugEnds Line 320  CcnhDebugEnds
320  C==    Evaluate laplace operator on conjugate gradient vector  C==    Evaluate laplace operator on conjugate gradient vector
321  C==    q = A.s  C==    q = A.s
322         alpha = 0. _d 0         alpha = 0. _d 0
323           topLevTerm = freeSurfFac*cg3dNorm*
324         &      (horiVertRatio/gravity)/deltaTMom/deltaTMom
325         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
326          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
327             alphaTile = 0. _d 0
328           IF ( Nr .GT. 1 ) THEN           IF ( Nr .GT. 1 ) THEN
329            DO K=1,1            DO K=1,1
330             DO J=1,sNy             DO J=1,sNy
# Line 388  C==    q = A.s Line 340  C==    q = A.s
340       &      -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)
341       &      -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)
342       &      -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)
343       &      -freeSurfFac*_rA(i,j,bi,bj)* horiVertRatio*       &      -topLevTerm*_rA(I,J,bi,bj)*cg3d_s(I,J,K,bi,bj)
344       &       cg3d_s(I  ,J  ,K,bi,bj)/deltaTMom/deltaTMom*cg3dNorm               alphaTile = alphaTile
345               alpha = alpha+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)
346              ENDDO              ENDDO
347             ENDDO             ENDDO
348            ENDDO            ENDDO
# Line 407  C==    q = A.s Line 359  C==    q = A.s
359       &      -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)
360       &      -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)
361       &      -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)
362       &      -freeSurfFac*_rA(i,j,bi,bj)* horiVertRatio*       &      -topLevTerm*_rA(I,J,bi,bj)*cg3d_s(I,J,K,bi,bj)
363       &       cg3d_s(I  ,J  ,K,bi,bj)/deltaTMom/deltaTMom*cg3dNorm               alphaTile = alphaTile
364               alpha = alpha+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)
365              ENDDO              ENDDO
366             ENDDO             ENDDO
367            ENDDO            ENDDO
# Line 430  C==    q = A.s Line 382  C==    q = A.s
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       &     -aV3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &     -aV3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)
384       &     -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)
385              alpha = alpha+cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)              alphaTile = alphaTile
386         &                +cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)
387             ENDDO             ENDDO
388            ENDDO            ENDDO
389           ENDDO           ENDDO
# Line 449  C==    q = A.s Line 402  C==    q = A.s
402       &      -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)
403       &      -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)
404       &      -aV3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &      -aV3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)
405               alpha = alpha+cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)               alphaTile = alphaTile
406         &                 +cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)
407              ENDDO              ENDDO
408             ENDDO             ENDDO
409            ENDDO            ENDDO
410           ENDIF           ENDIF
411             alpha = alpha + alphaTile
412          ENDDO          ENDDO
413         ENDDO         ENDDO
414         _GLOBAL_SUM_R8(alpha,myThid)         _GLOBAL_SUM_R8(alpha,myThid)
415  CcnhDebugStarts  CcnhDebugStarts
416  C      WRITE(0,*) ' CG3D: Iteration ',it3d-1,' SUM(s*q)= ',alpha  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' SUM(s*q)= ',alpha
417  CcnhDebugEnds  CcnhDebugEnds
418         alpha = etaN/alpha         alpha = eta_qrN/alpha
419  CcnhDebugStarts  CcnhDebugStarts
420  C      WRITE(0,*) ' CG3D: Iteration ',it3d-1,' alpha= ',alpha  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' alpha= ',alpha
421  CcnhDebugEnds  CcnhDebugEnds
422            
423  C==    Update solution and residual vectors  C==    Update solution and residual vectors
# Line 470  C      Now compute "interior" points. Line 425  C      Now compute "interior" points.
425         err = 0. _d 0         err = 0. _d 0
426         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
427          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
428            errTile    = 0. _d 0
429           DO K=1,Nr           DO K=1,Nr
430            DO J=1,sNy            DO J=1,sNy
431             DO I=1,sNx             DO I=1,sNx
# Line 477  C      Now compute "interior" points. Line 433  C      Now compute "interior" points.
433       &            +alpha*cg3d_s(I,J,K,bi,bj)       &            +alpha*cg3d_s(I,J,K,bi,bj)
434              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)
435       &            -alpha*cg3d_q(I,J,K,bi,bj)       &            -alpha*cg3d_q(I,J,K,bi,bj)
436              err = err+cg3d_r(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)             errTile = errTile
437         &             +cg3d_r(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)
438             ENDDO             ENDDO
439            ENDDO            ENDDO
440           ENDDO           ENDDO
441             err = err + errTile
442          ENDDO          ENDDO
443         ENDDO         ENDDO
444    
# Line 490  C      Now compute "interior" points. Line 448  C      Now compute "interior" points.
448         actualResidual = err         actualResidual = err
449         IF ( actualResidual .LT. cg3dTargetResidual ) GOTO 11         IF ( actualResidual .LT. cg3dTargetResidual ) GOTO 11
450  C      _EXCH_XYZ_R8(cg3d_r, myThid )  C      _EXCH_XYZ_R8(cg3d_r, myThid )
451         OLw        = 1         CALL EXCH_S3D_RL( cg3d_r, 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 )  
452    
453     10 CONTINUE     10 CONTINUE
454     11 CONTINUE     11 CONTINUE
# Line 518  C--   Un-normalise the answer Line 466  C--   Un-normalise the answer
466         ENDDO         ENDDO
467        ENDDO        ENDDO
468    
469        _EXCH_XYZ_R8(cg3d_x, myThid )  Cadj  _EXCH_XYZ_R8(cg3d_x, myThid )
470        _BEGIN_MASTER( myThid )  c     _BEGIN_MASTER( myThid )
471         WRITE(0,'(A,I6,1PE30.14)') ' CG3D iters, err = ',  c      WRITE(*,'(A,I6,1PE30.14)') ' CG3D iters, err = ',
472       &  actualIts, actualResidual  c    &  actualIts, actualResidual
473        _END_MASTER( )  c     _END_MASTER( myThid )
474          lastResidual=actualResidual
475  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  
476    
477  #endif /* ALLOW_NONHYDROSTATIC */  #endif /* ALLOW_NONHYDROSTATIC */
478    
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.16
  ViewVC Help
Powered by ViewVC 1.1.22