/[MITgcm]/MITgcm/model/src/cg2d.F
ViewVC logotype

Diff of /MITgcm/model/src/cg2d.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.10 by cnh, Sun Sep 6 17:35:19 1998 UTC revision 1.34 by cnh, Wed Sep 26 18:09:14 2001 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_EEOPTIONS.h"  #include "CPP_OPTIONS.h"
5    
6    CBOP
7    C     !ROUTINE: CG2D
8    C     !INTERFACE:
9        SUBROUTINE CG2D(          SUBROUTINE CG2D(  
10         I                cg2d_b,
11         U                cg2d_x,
12         O                firstResidual,
13         O                lastResidual,
14         U                numIters,
15       I                myThid )       I                myThid )
16  C     /==========================================================\  C     !DESCRIPTION: \bv
17  C     | SUBROUTINE CG2D                                          |  C     *==========================================================*
18  C     | o Two-dimensional grid problem conjugate-gradient        |  C     | SUBROUTINE CG2D                                          
19  C     |   inverter (with preconditioner).                        |  C     | o Two-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 five-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     | two-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    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"
46  #include "PARAMS.h"  #include "PARAMS.h"
47  #include "GRID.h"  #include "GRID.h"
48  #include "CG2D.h"  #include "CG2D.h"
49    #include "SURFACE.h"
50    
51    C     !INPUT/OUTPUT PARAMETERS:
52  C     === Routine arguments ===  C     === Routine arguments ===
53  C     myThid - Thread on which I am working.  C     myThid    - Thread on which I am working.
54    C     cg2d_b    - The source term or "right hand side"
55    C     cg2d_x    - The solution
56    C     firstResidual - the initial residual before any iterations
57    C     lastResidual  - the actual residual reached
58    C     numIters  - Entry: the maximum number of iterations allowed
59    C                 Exit:  the actual number of iterations used
60          _RL  cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
61          _RL  cg2d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
62          _RL  firstResidual
63          _RL  lastResidual
64          INTEGER numIters
65        INTEGER myThid        INTEGER myThid
66    
67    C     !LOCAL VARIABLES:
68  C     === Local variables ====  C     === Local variables ====
69  C     actualIts      - Number of iterations taken  C     actualIts      - Number of iterations taken
70  C     actualResidual - residual  C     actualResidual - residual
71  C     bi          - Block index in X and Y.  C     bi          - Block index in X and Y.
72  C     bj  C     bj
73  C     etaN        - Used in computing search directions  C     eta_qrN     - Used in computing search directions
74  C     etaNM1        suffix N and NM1 denote current and  C     eta_qrNM1     suffix N and NM1 denote current and
75  C     cgBeta        previous iterations respectively.  C     cgBeta        previous iterations respectively.
76  C     alpha    C     alpha  
77  C     sumRHS      - Sum of right-hand-side. Sometimes this is a  C     sumRHS      - Sum of right-hand-side. Sometimes this is a
# Line 54  C                   or they converge at Line 81  C                   or they converge at
81  C     err         - Measure of residual of Ax - b, usually the norm.  C     err         - Measure of residual of Ax - b, usually the norm.
82  C     I, J, N     - Loop counters ( N counts CG iterations )  C     I, J, N     - Loop counters ( N counts CG iterations )
83        INTEGER actualIts        INTEGER actualIts
84        REAL    actualResidual        _RL    actualResidual
85        INTEGER bi, bj                      INTEGER bi, bj              
86        INTEGER I, J, it2d        INTEGER I, J, it2d
87        REAL    err        _RL    err
88        REAL    etaN        _RL    eta_qrN
89        REAL    etaNM1        _RL    eta_qrNM1
90        REAL    cgBeta        _RL    cgBeta
91        REAL    alpha        _RL    alpha
92        REAL    sumRHS        _RL    sumRHS
93        REAL    rhsMax        _RL    rhsMax
94        REAL    rhsNorm        _RL    rhsNorm
95    
96          INTEGER OLw
97          INTEGER OLe
98          INTEGER OLn
99          INTEGER OLs
100          INTEGER exchWidthX
101          INTEGER exchWidthY
102          INTEGER myNz
103    CEOP
104    
105    
106    CcnhDebugStarts
107    C     CHARACTER*(MAX_LEN_FNAM) suff
108    CcnhDebugEnds
109    
110    
111  C--   Initialise inverter  C--   Initialise inverter
112        etaNBuf(1,myThid)   = 0. D0        eta_qrNM1 = 1. _d 0
       errBuf(1,myThid)    = 0. D0  
       sumRHSBuf(1,myThid) = 0. D0  
       etaNM1              = 1. D0  
113    
114  CcnhDebugStarts  CcnhDebugStarts
115        _EXCH_XY_R8( cg2d_b, myThid )  C     _EXCH_XY_R8( cg2d_b, myThid )
116        CALL PLOT_FIELD_XYRL( cg2d_b, 'CG2D.0 CG2D_B' , 1, myThid )  C     CALL PLOT_FIELD_XYRL( cg2d_b, 'CG2D.0 CG2D_B' , 1, myThid )
117    C     suff = 'unnormalised'
118    C     CALL WRITE_FLD_XY_RL (  'cg2d_b.',suff,    cg2d_b, 1, myThid)
119    C     STOP
120  CcnhDebugEnds  CcnhDebugEnds
121    
122  C--   Normalise RHS  C--   Normalise RHS
# Line 89  C--   Normalise RHS Line 131  C--   Normalise RHS
131          ENDDO          ENDDO
132         ENDDO         ENDDO
133        ENDDO        ENDDO
134        rhsMaxBuf(1,myThid) = rhsMax  
135        _GLOBAL_MAX_R8( rhsMaxbuf, rhsMax, myThid )        IF (cg2dNormaliseRHS) THEN
136        rhsMax =  rhsMaxBuf(1,1)  C-  Normalise RHS :
137    #ifdef LETS_MAKE_JAM
138    C     _GLOBAL_MAX_R8( rhsMax, myThid )
139          rhsMax=1.
140    #else
141          _GLOBAL_MAX_R8( rhsMax, myThid )
142    Catm  rhsMax=1.
143    #endif
144        rhsNorm = 1. _d 0        rhsNorm = 1. _d 0
145        IF ( rhsMax .NE. 0. ) rhsNorm = 1. _d 0 / rhsMax        IF ( rhsMax .NE. 0. ) rhsNorm = 1. _d 0 / rhsMax
146        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
# Line 104  C--   Normalise RHS Line 153  C--   Normalise RHS
153          ENDDO          ENDDO
154         ENDDO         ENDDO
155        ENDDO        ENDDO
156    C- end Normalise RHS
157          ENDIF
158    
159  C--   Update overlaps  C--   Update overlaps
160        _EXCH_XY_R8( cg2d_b, myThid )        _EXCH_XY_R8( cg2d_b, myThid )
161        _EXCH_XY_R8( cg2d_x, myThid )        _EXCH_XY_R8( cg2d_x, myThid )
162  CcnhDebugStarts  CcnhDebugStarts
163        CALL PLOT_FIELD_XYRL( cg2d_b, 'CG2D.1 CG2D_B' , 1, myThid )  C     CALL PLOT_FIELD_XYRL( cg2d_b, 'CG2D.1 CG2D_B' , 1, myThid )
164    C     suff = 'normalised'
165    C     CALL WRITE_FLD_XY_RL (  'cg2d_b.',suff,    cg2d_b, 1, myThid)
166  CcnhDebugEnds  CcnhDebugEnds
167    
168  C--   Initial residual calculation  C--   Initial residual calculation
# Line 129  C--   Initial residual calculation Line 182  C--   Initial residual calculation
182       &    -aW2d(I+1,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)       &    -aW2d(I+1,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)
183       &    -aS2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)       &    -aS2d(I  ,J  ,bi,bj)*cg2d_x(I  ,J  ,bi,bj)
184       &    -aS2d(I  ,J+1,bi,bj)*cg2d_x(I  ,J  ,bi,bj)       &    -aS2d(I  ,J+1,bi,bj)*cg2d_x(I  ,J  ,bi,bj)
185       &    -freeSurfFac*_rA(i,j,bi,bj)* horiVertRatio*       &    -freeSurfFac*_rA(i,j,bi,bj)*recip_Bo(i,j,bi,bj)*
186       &     cg2d_x(I  ,J  ,bi,bj)/deltaTMom/deltaTMom*cg2dNorm       &     cg2d_x(I  ,J  ,bi,bj)/deltaTMom/deltaTMom*cg2dNorm
187       &    )       &    )
188            err            = err            +            err            = err            +
# Line 140  C--   Initial residual calculation Line 193  C--   Initial residual calculation
193          ENDDO          ENDDO
194         ENDDO         ENDDO
195        ENDDO        ENDDO
196        _EXCH_XY_R8( cg2d_r, myThid )  C     _EXCH_XY_R8( cg2d_r, myThid )
197        _EXCH_XY_R8( cg2d_s, myThid )  #ifdef LETS_MAKE_JAM
198        sumRHSBuf(1,myThid) = sumRHS        CALL EXCH_XY_O1_R8_JAM( cg2d_r )
199        _GLOBAL_SUM_R8( sumRHSBuf , sumRHS, myThid )  #else
200        sumRHS = sumRHSBuf(1,1)         OLw        = 1
201        errBuf(1,myThid)    = err         OLe        = 1
202  C     WRITE(6,*) ' mythid, err = ', mythid, SQRT(err)         OLn        = 1
203        _GLOBAL_SUM_R8( errBuf    , err   , myThid )         OLs        = 1
204        err    = errBuf(1,1)         exchWidthX = 1
205        write(0,*) 'cg2d: Sum(rhs) = ',sumRHS         exchWidthY = 1
206           myNz       = 1
207        actualIts      = 0        IF (useCubedSphereExchange) THEN
208        actualResidual = SQRT(err)         CALL EXCH_RL_CUBE( cg2d_r,
209  C     _BARRIER       I            OLw, OLe, OLs, OLn, myNz,
210         I            exchWidthX, exchWidthY,
211         I            FORWARD_SIMULATION, EXCH_UPDATE_CORNERS, myThid )
212          ELSE
213           CALL EXCH_RL( cg2d_r,
214         I            OLw, OLe, OLs, OLn, myNz,
215         I            exchWidthX, exchWidthY,
216         I            FORWARD_SIMULATION, EXCH_UPDATE_CORNERS, myThid )
217          ENDIF
218    #endif
219    C     _EXCH_XY_R8( cg2d_s, myThid )
220    #ifdef LETS_MAKE_JAM
221          CALL EXCH_XY_O1_R8_JAM( cg2d_s )
222    #else
223           OLw        = 1
224           OLe        = 1
225           OLn        = 1
226           OLs        = 1
227           exchWidthX = 1
228           exchWidthY = 1
229           myNz       = 1
230          IF (useCubedSphereExchange) THEN
231           CALL EXCH_RL_CUBE( cg2d_s,
232         I            OLw, OLe, OLs, OLn, myNz,
233         I            exchWidthX, exchWidthY,
234         I            FORWARD_SIMULATION, EXCH_UPDATE_CORNERS, myThid )
235          ELSE
236           CALL EXCH_RL( cg2d_s,
237         I            OLw, OLe, OLs, OLn, myNz,
238         I            exchWidthX, exchWidthY,
239         I            FORWARD_SIMULATION, EXCH_UPDATE_CORNERS, myThid )
240          ENDIF
241    #endif
242           _GLOBAL_SUM_R8( sumRHS, myThid )
243           _GLOBAL_SUM_R8( err   , myThid )
244           err = SQRT(err)
245           actualIts      = 0
246           actualResidual = err
247          
248        _BEGIN_MASTER( myThid )        _BEGIN_MASTER( myThid )
249         WRITE(0,*) ' CG2D iters, err = ', actualIts, actualResidual        write(*,'(A,1P2E22.14)')' cg2d: Sum(rhs),rhsMax = ',
250         &                                  sumRHS,rhsMax
251        _END_MASTER( )        _END_MASTER( )
252    C     _BARRIER
253    c     _BEGIN_MASTER( myThid )
254    c      WRITE(*,'(A,I6,1PE30.14)') ' CG2D iters, err = ',
255    c    & actualIts, actualResidual
256    c     _END_MASTER( )
257          firstResidual=actualResidual
258    
259          IF ( err .LT. cg2dTolerance ) GOTO 11
260    
261  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
262        DO 10 it2d=1, cg2dMaxIters        DO 10 it2d=1, numIters
263    
264  CcnhDebugStarts  CcnhDebugStarts
265  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' residual = ',actualResidual  C      WRITE(*,*) ' CG2D: Iteration ',it2d-1,' residual = ',
266    C    &  actualResidual
267  CcnhDebugEnds  CcnhDebugEnds
        IF ( err .LT. cg2dTargetResidual ) GOTO 11  
268  C--    Solve preconditioning equation and update  C--    Solve preconditioning equation and update
269  C--    conjugate direction vector "s".  C--    conjugate direction vector "s".
270         etaN = 0. _d 0         eta_qrN = 0. _d 0
271         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
272          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
273           DO J=1,sNy           DO J=1,sNy
# Line 181  C--    conjugate direction vector "s". Line 281  C--    conjugate direction vector "s".
281  CcnhDebugStarts  CcnhDebugStarts
282  C          cg2d_q(I,J,bi,bj) = cg2d_r(I  ,J  ,bi,bj)  C          cg2d_q(I,J,bi,bj) = cg2d_r(I  ,J  ,bi,bj)
283  CcnhDebugEnds  CcnhDebugEnds
284             etaN = etaN             eta_qrN = eta_qrN
285       &     +cg2d_q(I,J,bi,bj)*cg2d_r(I,J,bi,bj)       &     +cg2d_q(I,J,bi,bj)*cg2d_r(I,J,bi,bj)
286            ENDDO            ENDDO
287           ENDDO           ENDDO
288          ENDDO          ENDDO
289         ENDDO         ENDDO
290    
291         etanBuf(1,myThid) = etaN         _GLOBAL_SUM_R8(eta_qrN, myThid)
        _GLOBAL_SUM_R8(etaNbuf,etaN, myThid)  
        etaN = etaNBuf(1,1)  
292  CcnhDebugStarts  CcnhDebugStarts
293  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' etaN = ',etaN  C      WRITE(*,*) ' CG2D: Iteration ',it2d-1,' eta_qrN = ',eta_qrN
294  CcnhDebugEnds  CcnhDebugEnds
295         cgBeta   = etaN/etaNM1         cgBeta   = eta_qrN/eta_qrNM1
296  CcnhDebugStarts  CcnhDebugStarts
297  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' beta = ',cgBeta  C      WRITE(*,*) ' CG2D: Iteration ',it2d-1,' beta = ',cgBeta
298  CcnhDebugEnds  CcnhDebugEnds
299         etaNM1 = etaN         eta_qrNM1 = eta_qrN
300    
301         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
302          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
303           DO J=1,sNy           DO J=1,sNy
304            DO I=1,sNx            DO I=1,sNx
305             cg2d_s(I,J,bi,bj) = cg2d_q(I,J,bi,bj) + cgBeta*cg2d_s(I,J,bi,bj)             cg2d_s(I,J,bi,bj) = cg2d_q(I,J,bi,bj)
306         &                       + cgBeta*cg2d_s(I,J,bi,bj)
307            ENDDO            ENDDO
308           ENDDO           ENDDO
309          ENDDO          ENDDO
# Line 212  CcnhDebugEnds Line 311  CcnhDebugEnds
311    
312  C--    Do exchanges that require messages i.e. between  C--    Do exchanges that require messages i.e. between
313  C--    processes.  C--    processes.
314         _EXCH_XY_R8( cg2d_s, myThid )  C      _EXCH_XY_R8( cg2d_s, myThid )
315    #ifdef LETS_MAKE_JAM
316          CALL EXCH_XY_O1_R8_JAM( cg2d_s )
317    #else
318           OLw        = 1
319           OLe        = 1
320           OLn        = 1
321           OLs        = 1
322           exchWidthX = 1
323           exchWidthY = 1
324           myNz       = 1
325          IF (useCubedSphereExchange) THEN
326           CALL EXCH_RL_CUBE( cg2d_s,
327         I            OLw, OLe, OLs, OLn, myNz,
328         I            exchWidthX, exchWidthY,
329         I            FORWARD_SIMULATION, EXCH_UPDATE_CORNERS, myThid )
330          ELSE
331           CALL EXCH_RL( cg2d_s,
332         I            OLw, OLe, OLs, OLn, myNz,
333         I            exchWidthX, exchWidthY,
334         I            FORWARD_SIMULATION, EXCH_UPDATE_CORNERS, myThid )
335          ENDIF
336    #endif
337    
338  C==    Evaluate laplace operator on conjugate gradient vector  C==    Evaluate laplace operator on conjugate gradient vector
339  C==    q = A.s  C==    q = A.s
# Line 230  C==    q = A.s Line 351  C==    q = A.s
351       &    -aW2d(I+1,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)       &    -aW2d(I+1,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)
352       &    -aS2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)       &    -aS2d(I  ,J  ,bi,bj)*cg2d_s(I  ,J  ,bi,bj)
353       &    -aS2d(I  ,J+1,bi,bj)*cg2d_s(I  ,J  ,bi,bj)       &    -aS2d(I  ,J+1,bi,bj)*cg2d_s(I  ,J  ,bi,bj)
354       &    -freeSurfFac*_rA(i,j,bi,bj)* horiVertRatio*       &    -freeSurfFac*_rA(i,j,bi,bj)*recip_Bo(i,j,bi,bj)*
355       &     cg2d_s(I  ,J  ,bi,bj)/deltaTMom/deltaTMom*cg2dNorm       &     cg2d_s(I  ,J  ,bi,bj)/deltaTMom/deltaTMom*cg2dNorm
356            alpha = alpha+cg2d_s(I,J,bi,bj)*cg2d_q(I,J,bi,bj)            alpha = alpha+cg2d_s(I,J,bi,bj)*cg2d_q(I,J,bi,bj)
357            ENDDO            ENDDO
358           ENDDO           ENDDO
359          ENDDO          ENDDO
360         ENDDO         ENDDO
361         alphaBuf(1,myThid) = alpha         _GLOBAL_SUM_R8(alpha,myThid)
        _GLOBAL_SUM_R8(alphaBuf,alpha,myThid)  
        alpha = alphaBuf(1,1)  
362  CcnhDebugStarts  CcnhDebugStarts
363  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' SUM(s*q)= ',alpha  C      WRITE(*,*) ' CG2D: Iteration ',it2d-1,' SUM(s*q)= ',alpha
364  CcnhDebugEnds  CcnhDebugEnds
365         alpha = etaN/alpha         alpha = eta_qrN/alpha
366  CcnhDebugStarts  CcnhDebugStarts
367  C      WRITE(0,*) ' CG2D: Iteration ',it2d-1,' alpha= ',alpha  C      WRITE(*,*) ' CG2D: Iteration ',it2d-1,' alpha= ',alpha
368  CcnhDebugEnds  CcnhDebugEnds
369            
370  C==    Update solution and residual vectors  C==    Update solution and residual vectors
# Line 263  C      Now compute "interior" points. Line 382  C      Now compute "interior" points.
382          ENDDO          ENDDO
383         ENDDO         ENDDO
384    
385         errBuf(1,myThid) = err         _GLOBAL_SUM_R8( err   , myThid )
        _GLOBAL_SUM_R8( errBuf    , err   , myThid )  
        err = errBuf(1,1)  
386         err = SQRT(err)         err = SQRT(err)
387         actualIts      = it2d         actualIts      = it2d
388         actualResidual = err         actualResidual = err
389         IF ( err .LT. cg2dTargetResidual ) GOTO 11         IF ( err .LT. cg2dTolerance ) GOTO 11
390         _EXCH_XY_R8(cg2d_r, myThid )  C      _EXCH_XY_R8(cg2d_r, myThid )
391    #ifdef LETS_MAKE_JAM
392          CALL EXCH_XY_O1_R8_JAM( cg2d_r )
393    #else
394           OLw        = 1
395           OLe        = 1
396           OLn        = 1
397           OLs        = 1
398           exchWidthX = 1
399           exchWidthY = 1
400           myNz       = 1
401          IF (useCubedSphereExchange) THEN
402           CALL EXCH_RL_CUBE( cg2d_r,
403         I            OLw, OLe, OLs, OLn, myNz,
404         I            exchWidthX, exchWidthY,
405         I            FORWARD_SIMULATION, EXCH_UPDATE_CORNERS, myThid )
406          ELSE
407           CALL EXCH_RL( cg2d_r,
408         I            OLw, OLe, OLs, OLn, myNz,
409         I            exchWidthX, exchWidthY,
410         I            FORWARD_SIMULATION, EXCH_UPDATE_CORNERS, myThid )
411          ENDIF
412    #endif
413    
414     10 CONTINUE     10 CONTINUE
415     11 CONTINUE     11 CONTINUE
416    
417          IF (cg2dNormaliseRHS) THEN
418  C--   Un-normalise the answer  C--   Un-normalise the answer
419        DO bj=myByLo(myThid),myByHi(myThid)          DO bj=myByLo(myThid),myByHi(myThid)
420         DO bi=myBxLo(myThid),myBxHi(myThid)           DO bi=myBxLo(myThid),myBxHi(myThid)
421          DO J=1,sNy            DO J=1,sNy
422           DO I=1,sNx             DO I=1,sNx
423            cg2d_x(I  ,J  ,bi,bj) = cg2d_x(I  ,J  ,bi,bj)/rhsNorm              cg2d_x(I  ,J  ,bi,bj) = cg2d_x(I  ,J  ,bi,bj)/rhsNorm
424               ENDDO
425              ENDDO
426           ENDDO           ENDDO
427          ENDDO          ENDDO
428         ENDDO        ENDIF
       ENDDO  
429    
430        _EXCH_XY_R8(cg2d_x, myThid )  C     The following exchange was moved up to solve_for_pressure
431        _BEGIN_MASTER( myThid )  C     for compatibility with TAMC.
432         WRITE(0,*) ' CG2D iters, err = ', actualIts, actualResidual  C     _EXCH_XY_R8(cg2d_x, myThid )
433        _END_MASTER( )  c     _BEGIN_MASTER( myThid )
434    c      WRITE(*,'(A,I6,1PE30.14)') ' CG2D iters, err = ',
435    c    & actualIts, actualResidual
436    c     _END_MASTER( )
437    
438    C--   Return parameters to caller
439          lastResidual=actualResidual
440          numIters=actualIts
441    
442  CcnhDebugStarts  CcnhDebugStarts
443  C     CALL PLOT_FIELD_XYRL( cg2d_x, 'CALC_MOM_RHS CG2D_X' , 1, myThid )  C     CALL PLOT_FIELD_XYRL( cg2d_x, 'CALC_MOM_RHS CG2D_X' , 1, myThid )
# Line 312  C        ENDDO Line 461  C        ENDDO
461  C       ENDDO  C       ENDDO
462  C      ENDDO  C      ENDDO
463  C     ENDDO  C     ENDDO
464  C     errBuf(1,myThid)    = err  C     _GLOBAL_SUM_R8( err   , myThid )
465  C     _GLOBAL_SUM_R8( errBuf    , err   , myThid )  C     write(*,*) 'cg2d: Ax - b = ',SQRT(err)
 C     err    = errBuf(1,1)  
 C     write(0,*) 'cg2d: Ax - b = ',SQRT(err)  
466  CcnhDebugEnds  CcnhDebugEnds
467    
468          RETURN
469        END        END
Legend:
Removed from v.1.10  
changed lines
  Added in v.1.34
  ViewVC Help
Powered by ViewVC 1.1.22