/[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.24 by jmc, Wed Jun 8 01:46:34 2011 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    #ifdef TARGET_NEC_SX
6  #define VERBOSE  C     set a sensible default for the outer loop unrolling parameter that can
7    C     be overriden in the Makefile with the DEFINES macro or in CPP_OPTIONS.h
8    #ifndef CG3D_OUTERLOOPITERS
9    # define CG3D_OUTERLOOPITERS 10
10    #endif
11    #endif /* TARGET_NEC_SX */
12    
13  CBOP  CBOP
14  C     !ROUTINE: CG3D  C     !ROUTINE: CG3D
15  C     !INTERFACE:  C     !INTERFACE:
16        SUBROUTINE CG3D(          SUBROUTINE CG3D(
17       I                cg3d_b,       I                cg3d_b,
18       U                cg3d_x,       U                cg3d_x,
19       O                firstResidual,       O                firstResidual,
20       O                lastResidual,       O                lastResidual,
21       U                numIters,       U                numIters,
22       I                myThid )       I                myIter, myThid )
23  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
24  C     *==========================================================*  C     *==========================================================*
25  C     | SUBROUTINE CG3D                                            C     | SUBROUTINE CG3D
26  C     | o Three-dimensional grid problem conjugate-gradient        C     | o Three-dimensional grid problem conjugate-gradient
27  C     |   inverter (with preconditioner).                          C     |   inverter (with preconditioner).
28  C     *==========================================================*  C     *==========================================================*
29  C     | Con. grad is an iterative procedure for solving Ax = b.    C     | Con. grad is an iterative procedure for solving Ax = b.
30  C     | It requires the A be symmetric.                            C     | It requires the A be symmetric.
31  C     | This implementation assumes A is a seven-diagonal            C     | This implementation assumes A is a seven-diagonal
32  C     | matrix of the form that arises in the discrete              C     | matrix of the form that arises in the discrete
33  C     | representation of the del^2 operator in a                  C     | representation of the del^2 operator in a
34  C     | three-dimensional space.                                      C     | three-dimensional space.
35  C     | Notes:                                                      C     | Notes:
36  C     | ======                                                      C     | ======
37  C     | This implementation can support shared-memory                C     | This implementation can support shared-memory
38  C     | multi-threaded execution. In order to do this COMMON        C     | multi-threaded execution. In order to do this COMMON
39  C     | blocks are used for many of the arrays - even ones that      C     | blocks are used for many of the arrays - even ones that
40  C     | are only used for intermedaite results. This design is      C     | are only used for intermedaite results. This design is
41  C     | OK if you want to all the threads to collaborate on          C     | OK if you want to all the threads to collaborate on
42  C     | solving the same problem. On the other hand if you want      C     | solving the same problem. On the other hand if you want
43  C     | the threads to solve several different problems              C     | the threads to solve several different problems
44  C     | concurrently this implementation will not work.            C     | concurrently this implementation will not work.
45  C     *==========================================================*  C     *==========================================================*
46  C     \ev  C     \ev
47    
# Line 48  C     === Global data === Line 52  C     === Global data ===
52  #include "EEPARAMS.h"  #include "EEPARAMS.h"
53  #include "PARAMS.h"  #include "PARAMS.h"
54  #include "GRID.h"  #include "GRID.h"
55    #include "SURFACE.h"
56  #include "CG3D.h"  #include "CG3D.h"
57    
58  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
59  C     === Routine arguments ===  C     === Routine arguments ===
60  C     myThid    - Thread on which I am working.  C     cg3d_b    :: The source term or "right hand side"
61  C     cg2d_b    - The source term or "right hand side"  C     cg3d_x    :: The solution
62  C     cg2d_x    - The solution  C     firstResidual :: the initial residual before any iterations
63  C     firstResidual - the initial residual before any iterations  C     lastResidual  :: the actual residual reached
64  C     lastResidual  - the actual residual reached  C     numIters  :: Entry: the maximum number of iterations allowed
65  C     numIters  - Entry: the maximum number of iterations allowed  C               :: Exit:  the actual number of iterations used
66  C                 Exit:  the actual number of iterations used  C     myIter    :: Current iteration number in simulation
67        _RL  cg3d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)  C     myThid    :: my Thread Id number
68        _RL  cg3d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL     cg3d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
69        _RL  firstResidual        _RL     cg3d_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
70        _RL  lastResidual        _RL     firstResidual
71          _RL     lastResidual
72        INTEGER numIters        INTEGER numIters
73          INTEGER myIter
74        INTEGER myThid        INTEGER myThid
75    
   
76  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
77    
78  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
79  C     === Local variables ====  C     === Local variables ====
80  C     actualIts      - Number of iterations taken  C     actualIts      :: Number of iterations taken
81  C     actualResidual - residual  C     actualResidual :: residual
82  C     bi          - Block index in X and Y.  C     bi,bj      :: tile indices
83  C     bj  C     eta_qrN    :: Used in computing search directions
 C     eta_qrN     - Used in computing search directions  
84  C     eta_qrNM1     suffix N and NM1 denote current and  C     eta_qrNM1     suffix N and NM1 denote current and
85  C     cgBeta        previous iterations respectively.  C     cgBeta        previous iterations respectively.
86  C     alpha    C     alpha
87  C     sumRHS      - Sum of right-hand-side. Sometimes this is a  C     sumRHS     :: Sum of right-hand-side. Sometimes this is a
88  C                   useful debuggin/trouble shooting diagnostic.  C                   useful debuggin/trouble shooting diagnostic.
89  C                   For neumann problems sumRHS needs to be ~0.  C                   For neumann problems sumRHS needs to be ~0.
90  C                   or they converge at a non-zero residual.  C                   or they converge at a non-zero residual.
91  C     err         - Measure of residual of Ax - b, usually the norm.  C     err        :: Measure of residual of Ax - b, usually the norm.
92  C     I, J, N     - Loop counters ( N counts CG iterations )  C     i, j, k    :: Loop counters
93    C     it3d       :: Loop counter for CG iterations
94    C     msgBuf     :: Informational/error message buffer
95        INTEGER actualIts        INTEGER actualIts
96        _RL    actualResidual        _RL     actualResidual
97        INTEGER bi, bj                      INTEGER bi, bj
98        INTEGER I, J, K, it3d        INTEGER i, j, k, it3d
99        INTEGER KM1, KP1        INTEGER km1, kp1
100        _RL    err        _RL     maskM1, maskP1
101        _RL    eta_qrN        _RL     err,    errTile(nSx,nSy)
102        _RL    eta_qrNM1        _RL     eta_qrN,eta_qrNtile(nSx,nSy)
103        _RL    cgBeta        _RL     eta_qrNM1
104        _RL    alpha        _RL     cgBeta
105        _RL    sumRHS        _RL     alpha , alphaTile(nSx,nSy)
106        _RL    rhsMax        _RL     sumRHS, sumRHStile(nSx,nSy)
107        _RL    rhsNorm        _RL     rhsMax
108          _RL     rhsNorm
109        INTEGER OLw        CHARACTER*(MAX_LEN_MBUF) msgBuf
110        INTEGER OLe        LOGICAL printResidual
111        INTEGER OLn        _RL     surfFac
112        INTEGER OLs  #ifdef NONLIN_FRSURF
113        INTEGER exchWidthX        INTEGER ks
114        INTEGER exchWidthY        _RL     surfTerm(sNx,sNy)
115        INTEGER myNz  #endif /* NONLIN_FRSURF */
       _RL     topLevTerm  
116  CEOP  CEOP
117    
118  ceh3 needs an IF ( useNONHYDROSTATIC ) THEN        IF ( select_rStar .NE. 0 ) THEN
119            surfFac = freeSurfFac
120          ELSE
121            surfFac = 0.
122          ENDIF
123    #ifdef NONLIN_FRSURF
124          DO j=1,sNy
125           DO i=1,sNx
126             surfTerm(i,j) = 0.
127           ENDDO
128          ENDDO
129    #endif /* NONLIN_FRSURF */
130    
131  C--   Initialise inverter  C--   Initialise inverter
132        eta_qrNM1 = 1. D0        eta_qrNM1 = 1. _d 0
133    
134  C--   Normalise RHS  C--   Normalise RHS
135        rhsMax = 0. _d 0        rhsMax = 0. _d 0
136        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
137         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
138          DO K=1,Nr          DO k=1,Nr
139           DO J=1,sNy           DO j=1,sNy
140            DO I=1,sNx            DO i=1,sNx
141             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
142             rhsMax = MAX(ABS(cg3d_b(I,J,K,bi,bj)),rhsMax)       &                          * maskC(i,j,k,bi,bj)
143               rhsMax = MAX(ABS(cg3d_b(i,j,k,bi,bj)),rhsMax)
144            ENDDO            ENDDO
145           ENDDO           ENDDO
146          ENDDO          ENDDO
147         ENDDO         ENDDO
148        ENDDO        ENDDO
149        _GLOBAL_MAX_R8( rhsMax, myThid )        _GLOBAL_MAX_RL( rhsMax, myThid )
150        rhsNorm = 1. _d 0        rhsNorm = 1. _d 0
151        IF ( rhsMax .NE. 0. ) rhsNorm = 1. _d 0 / rhsMax        IF ( rhsMax .NE. 0. ) rhsNorm = 1. _d 0 / rhsMax
152        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
153         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
154          DO K=1,Nr          DO k=1,Nr
155           DO J=1,sNy           DO j=1,sNy
156            DO I=1,sNx            DO i=1,sNx
157             cg3d_b(I,J,K,bi,bj) = cg3d_b(I,J,K,bi,bj)*rhsNorm             cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)*rhsNorm
158             cg3d_x(I,J,K,bi,bj) = cg3d_x(I,J,K,bi,bj)*rhsNorm             cg3d_x(i,j,k,bi,bj) = cg3d_x(i,j,k,bi,bj)*rhsNorm
159            ENDDO            ENDDO
160           ENDDO           ENDDO
161          ENDDO          ENDDO
# Line 146  C--   Normalise RHS Line 163  C--   Normalise RHS
163        ENDDO        ENDDO
164    
165  C--   Update overlaps  C--   Update overlaps
166        _EXCH_XYZ_R8( cg3d_b, myThid )  c     _EXCH_XYZ_RL( cg3d_b, myThid )
167        _EXCH_XYZ_R8( cg3d_x, myThid )        _EXCH_XYZ_RL( cg3d_x, myThid )
168    
169  C--   Initial residual calculation (with free-Surface term)  C--   Initial residual calculation (with free-Surface term)
170        err    = 0. _d 0        err    = 0. _d 0
171        sumRHS = 0. _d 0        sumRHS = 0. _d 0
172        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
173         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
174          DO K=1,Nr          errTile(bi,bj)    = 0. _d 0
175           KM1 = K-1          sumRHStile(bi,bj) = 0. _d 0
176           IF ( K .EQ. 1 ) KM1 = 1  #ifdef NONLIN_FRSURF
177           KP1 = K+1          IF ( select_rStar .NE. 0 ) THEN
178           IF ( K .EQ. Nr ) KP1 = 1            DO j=1,sNy
179           topLevTerm = 0.             DO i=1,sNx
180           IF ( K .EQ. 1) topLevTerm = freeSurfFac*cg3dNorm*               surfTerm(i,j) = 0.
181       &      (horiVertRatio/gravity)/deltaTMom/deltaTMom             ENDDO
182           DO J=1,sNy            ENDDO
183            DO I=1,sNx            DO k=1,Nr
184             cg3d_s(I,J,K,bi,bj) = 0.             DO j=1,sNy
185             cg3d_r(I,J,K,bi,bj) = cg3d_b(I,J,K,bi,bj) -( 0.              DO i=1,sNx
186       &     +aW3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I-1,J  ,K  ,bi,bj)               surfTerm(i,j) = surfTerm(i,j)
187       &     +aW3d(I+1,J  ,K  ,bi,bj)*cg3d_x(I+1,J  ,K  ,bi,bj)       &         +cg3d_x(i,j,k,bi,bj)*drF(k)*h0FacC(i,j,k,bi,bj)
188       &     +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J-1,K  ,bi,bj)              ENDDO
189       &     +aS3d(I  ,J+1,K  ,bi,bj)*cg3d_x(I  ,J+1,K  ,bi,bj)             ENDDO
190       &     +aV3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J  ,KM1,bi,bj)            ENDDO
191       &     +aV3d(I  ,J  ,KP1,bi,bj)*cg3d_x(I  ,J  ,KP1,bi,bj)            DO j=1,sNy
192       &     -aW3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)             DO i=1,sNx
193       &     -aW3d(I+1,J  ,K  ,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)               ks = kSurfC(i,j,bi,bj)
194       &     -aS3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)               surfTerm(i,j) = surfTerm(i,j)*cg3dNorm
195       &     -aS3d(I  ,J+1,K  ,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)       &              *recip_Rcol(i,j,bi,bj)*recip_Rcol(i,j,bi,bj)
196       &     -aV3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)       &              *rA(i,j,bi,bj)*deepFac2F(ks)
197       &     -aV3d(I  ,J  ,KP1,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)       &              *recip_Bo(i,j,bi,bj)/deltaTMom/deltaTfreesurf
198       &     -topLevTerm*_rA(I,J,bi,bj)*cg3d_x(I,J,K,bi,bj)             ENDDO
199       &     )            ENDDO
200             err = err          ENDIF
201       &     +cg3d_r(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)  #endif /* NONLIN_FRSURF */
202             sumRHS = sumRHS          DO k=1,Nr
203       &     +cg3d_b(I,J,K,bi,bj)           km1 = MAX(k-1, 1 )
204             kp1 = MIN(k+1, Nr)
205             maskM1 = 1. _d 0
206             maskP1 = 1. _d 0
207             IF ( k .EQ. 1 ) maskM1 = 0. _d 0
208             IF ( k .EQ. Nr) maskP1 = 0. _d 0
209    #ifdef TARGET_NEC_SX
210    !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
211    #endif /* TARGET_NEC_SX */
212             DO j=1,sNy
213              DO i=1,sNx
214               cg3d_r(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
215         &     -( 0.
216         &       +aW3d( i, j, k, bi,bj)*cg3d_x(i-1,j, k, bi,bj)
217         &       +aW3d(i+1,j, k, bi,bj)*cg3d_x(i+1,j, k, bi,bj)
218         &       +aS3d( i, j, k, bi,bj)*cg3d_x( i,j-1,k, bi,bj)
219         &       +aS3d( i,j+1,k, bi,bj)*cg3d_x( i,j+1,k, bi,bj)
220         &       +aV3d( i, j, k, bi,bj)*cg3d_x( i, j,km1,bi,bj)*maskM1
221         &       +aV3d( i, j,kp1,bi,bj)*cg3d_x( i, j,kp1,bi,bj)*maskP1
222         &       +aC3d( i, j, k, bi,bj)*cg3d_x( i, j, k, bi,bj)
223    #ifdef NONLIN_FRSURF
224         &       -surfFac*surfTerm(i,j)*drF(k)*h0FacC(i,j,k,bi,bj)
225    #endif /* NONLIN_FRSURF */
226         &      )
227               errTile(bi,bj) = errTile(bi,bj)
228         &                    +cg3d_r(i,j,k,bi,bj)*cg3d_r(i,j,k,bi,bj)
229               sumRHStile(bi,bj) = sumRHStile(bi,bj)+cg3d_b(i,j,k,bi,bj)
230              ENDDO
231             ENDDO
232             DO J=1-1,sNy+1
233              DO I=1-1,sNx+1
234               cg3d_s(i,j,k,bi,bj) = 0.
235            ENDDO            ENDDO
236           ENDDO           ENDDO
237          ENDDO          ENDDO
238         ENDDO         ENDDO
239        ENDDO        ENDDO
240  C     _EXCH_XYZ_R8( cg3d_r, myThid )         CALL EXCH_S3D_RL( cg3d_r, Nr, myThid )
241         OLw        = 1  c      CALL EXCH_S3D_RL( cg3d_s, Nr, myThid )
242         OLe        = 1         CALL GLOBAL_SUM_TILE_RL( sumRHStile, sumRHS, myThid )
243         OLn        = 1         CALL GLOBAL_SUM_TILE_RL( errTile,    err,    myThid )
244         OLs        = 1        IF ( debugLevel.GE.debLevC .AND. diagFreq.GT.0. ) THEN
245         exchWidthX = 1          CALL WRITE_FLD_S3D_RL(
246         exchWidthY = 1       I        'cg3d_r_I', 'I10', 1, Nr, cg3d_r, myIter, myThid )
247         myNz       = Nr        ENDIF
248         CALL EXCH_RL( cg3d_r,  
249       I            OLw, OLe, OLs, OLn, myNz,        printResidual = .FALSE.
250       I            exchWidthX, exchWidthY,        IF ( debugLevel .GE. debLevZero ) THEN
251       I            FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )          _BEGIN_MASTER( myThid )
252  C     _EXCH_XYZ_R8( cg3d_s, myThid )          printResidual = printResidualFreq.GE.1
253         OLw        = 1          WRITE(standardmessageunit,'(A,1P2E22.14)')
        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 )  
       _GLOBAL_SUM_R8( sumRHS, myThid )  
       _GLOBAL_SUM_R8( err   , myThid )  
         
       _BEGIN_MASTER( myThid )  
       write(*,'(A,1P2E22.14)')  
254       &     ' cg3d: Sum(rhs),rhsMax = ',sumRHS,rhsMax       &     ' cg3d: Sum(rhs),rhsMax = ',sumRHS,rhsMax
255        _END_MASTER( )          _END_MASTER( myThid )
256          ENDIF
257    
258        actualIts      = 0        actualIts      = 0
259        actualResidual = SQRT(err)        actualResidual = SQRT(err)
 C     _BARRIER  
 c     _BEGIN_MASTER( myThid )  
 c      WRITE(*,'(A,I6,1PE30.14)') ' CG3D iters, err = ',  
 c    &  actualIts, actualResidual  
 c     _END_MASTER( )  
260        firstResidual=actualResidual        firstResidual=actualResidual
261    
262  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
263        DO 10 it3d=1, cg3dMaxIters        DO 10 it3d=1, numIters
264    
 CcnhDebugStarts  
 #ifdef VERBOSE  
 c      IF ( mod(it3d-1,10).EQ.0)  
 c    &   WRITE(*,*) ' CG3D: Iteration ',it3d-1,  
 c    &      ' residual = ',actualResidual  
 #endif  
 CcnhDebugEnds  
265         IF ( actualResidual .LT. cg3dTargetResidual ) GOTO 11         IF ( actualResidual .LT. cg3dTargetResidual ) GOTO 11
266  C--    Solve preconditioning equation and update  C--    Solve preconditioning equation and update
267  C--    conjugate direction vector "s".  C--    conjugate direction vector "s".
268  C      Note. On the next to loops over all tiles the inner loop ranges  C      Note. On the next two loops over all tiles the inner loop ranges
269  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
270  C            step. However this entails a bit of gynamastics because we only  C            step. However this entails a bit of gynamastics because we only
271  C            want eta_qrN for the interior points.  C            want eta_qrN for the interior points.
272         eta_qrN = 0. _d 0         eta_qrN = 0. _d 0
273         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
274          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
275           DO K=1,1           eta_qrNtile(bi,bj) = 0. _d 0
276            DO J=1-1,sNy+1           DO k=1,1
277             DO I=1-1,sNx+1  #ifdef TARGET_NEC_SX
278              cg3d_q(I,J,K,bi,bj) =  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
279       &       zMC(I  ,J  ,K,bi,bj)*cg3d_r(I  ,J  ,K,bi,bj)  #endif /* TARGET_NEC_SX */
280              DO j=1-1,sNy+1
281               DO i=1-1,sNx+1
282                cg3d_q(i,j,k,bi,bj) = zMC(i,j,k,bi,bj)
283         &                        *cg3d_r(i,j,k,bi,bj)
284             ENDDO             ENDDO
285            ENDDO            ENDDO
286           ENDDO           ENDDO
287           DO K=2,Nr           DO k=2,Nr
288            DO J=1-1,sNy+1  #ifdef TARGET_NEC_SX
289             DO I=1-1,sNx+1  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
290              cg3d_q(I,J,K,bi,bj) =  #endif /* TARGET_NEC_SX */
291       &       zMC(I,J,K,bi,bj)*(cg3d_r(I,J,K  ,bi,bj)            DO j=1-1,sNy+1
292       &      -zML(I,J,K,bi,bj)*cg3d_q(I,J,K-1,bi,bj))             DO i=1-1,sNx+1
293                cg3d_q(i,j,k,bi,bj) = zMC(i,j,k,bi,bj)
294         &                      *( cg3d_r(i,j,k,bi,bj)
295         &                         -zML(i,j,k,bi,bj)*cg3d_q(i,j,k-1,bi,bj)
296         &                       )
297             ENDDO             ENDDO
298            ENDDO            ENDDO
299           ENDDO           ENDDO
300           DO K=Nr,Nr           DO k=Nr,Nr
301  caja      IF (Nr .GT. 1) THEN  #ifdef TARGET_NEC_SX
302  caja       DO J=1-1,sNy+1  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
303  caja        DO I=1-1,sNx+1  #endif /* TARGET_NEC_SX */
304  caja         cg3d_q(I,J,K,bi,bj) =            DO j=1,sNy
305  caja &        zMC(i,j,k,bi,bj)*(cg3d_r(i,j,k  ,bi,bj)             DO i=1,sNx
306  caja &       -zML(i,j,k,bi,bj)*cg3d_q(i,j,k-1,bi,bj))              eta_qrNtile(bi,bj) = eta_qrNtile(bi,bj)
307  caja        ENDDO       &                        +cg3d_q(i,j,k,bi,bj)*cg3d_r(i,j,k,bi,bj)
 caja       ENDDO  
 caja      ENDIF  
           DO J=1,sNy  
            DO I=1,sNx  
             eta_qrN = eta_qrN  
      &      +cg3d_q(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)  
308             ENDDO             ENDDO
309            ENDDO            ENDDO
310           ENDDO           ENDDO
311           DO K=Nr-1,1,-1           DO k=Nr-1,1,-1
312            DO J=1-1,sNy+1  #ifdef TARGET_NEC_SX
313             DO I=1-1,sNx+1  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
314              cg3d_q(I,J,K,bi,bj) =  #endif /* TARGET_NEC_SX */
315       &      cg3d_q(I,J,K,bi,bj)            DO j=1-1,sNy+1
316       &      -zMU(I,J,K,bi,bj)*cg3d_q(I,J,K+1,bi,bj)             DO i=1-1,sNx+1
317             ENDDO              cg3d_q(i,j,k,bi,bj) = cg3d_q(i,j,k,bi,bj)
318            ENDDO       &                         -zMU(i,j,k,bi,bj)*cg3d_q(i,j,k+1,bi,bj)
319            DO J=1,sNy             ENDDO
320             DO I=1,sNx            ENDDO
321              eta_qrN = eta_qrN  #ifdef TARGET_NEC_SX
322       &      +cg3d_q(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
323    #endif /* TARGET_NEC_SX */
324              DO j=1,sNy
325               DO i=1,sNx
326                eta_qrNtile(bi,bj) = eta_qrNtile(bi,bj)
327         &                        +cg3d_q(i,j,k,bi,bj)*cg3d_r(i,j,k,bi,bj)
328             ENDDO             ENDDO
329            ENDDO            ENDDO
330           ENDDO           ENDDO
331          ENDDO          ENDDO
332         ENDDO         ENDDO
 caja  
 caja  eta_qrN=0.  
 caja   DO bj=myByLo(myThid),myByHi(myThid)  
 caja    DO bi=myBxLo(myThid),myBxHi(myThid)  
 caja     DO K=1,Nr  
 caja      DO J=1,sNy  
 caja       DO I=1,sNx  
 caja        eta_qrN = eta_qrN  
 caja &      +cg3d_q(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)  
 caja       ENDDO  
 caja      ENDDO  
 caja     ENDDO  
 caja    ENDDO  
 caja   ENDDO  
 caja  
333    
334         _GLOBAL_SUM_R8(eta_qrN, myThid)         CALL GLOBAL_SUM_TILE_RL( eta_qrNtile,eta_qrN,myThid )
335  CcnhDebugStarts  CcnhDebugStarts
336  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' eta_qrN = ',eta_qrN  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' eta_qrN = ',eta_qrN
337  CcnhDebugEnds  CcnhDebugEnds
# Line 329  CcnhDebugEnds Line 343  CcnhDebugEnds
343    
344         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
345          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
346           DO K=1,Nr           DO k=1,Nr
347            DO J=1-1,sNy+1            DO j=1-1,sNy+1
348             DO I=1-1,sNx+1             DO i=1-1,sNx+1
349              cg3d_s(I,J,K,bi,bj) = cg3d_q(I,J,K,bi,bj)              cg3d_s(i,j,k,bi,bj) = cg3d_q(i,j,k,bi,bj)
350       &                          + cgBeta*cg3d_s(I,J,K,bi,bj)       &                   + cgBeta*cg3d_s(i,j,k,bi,bj)
351             ENDDO             ENDDO
352            ENDDO            ENDDO
353           ENDDO           ENDDO
# Line 343  CcnhDebugEnds Line 357  CcnhDebugEnds
357  C==    Evaluate laplace operator on conjugate gradient vector  C==    Evaluate laplace operator on conjugate gradient vector
358  C==    q = A.s  C==    q = A.s
359         alpha = 0. _d 0         alpha = 0. _d 0
        topLevTerm = freeSurfFac*cg3dNorm*  
      &      (horiVertRatio/gravity)/deltaTMom/deltaTMom  
360         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
361          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
362           IF ( Nr .GT. 1 ) THEN           alphaTile(bi,bj) = 0. _d 0
363            DO K=1,1  #ifdef NONLIN_FRSURF
364             DO J=1,sNy           IF ( select_rStar .NE. 0 ) THEN
365              DO I=1,sNx            DO j=1,sNy
366               cg3d_q(I,J,K,bi,bj) =             DO i=1,sNx
367       &       aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I-1,J  ,K  ,bi,bj)               surfTerm(i,j) = 0.
368       &      +aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I+1,J  ,K  ,bi,bj)             ENDDO
369       &      +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J-1,K  ,bi,bj)            ENDDO
370       &      +aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J+1,K  ,bi,bj)            DO k=1,Nr
371       &      +aV3d(I  ,J  ,K+1,bi,bj)*cg3d_s(I  ,J  ,K+1,bi,bj)             DO j=1,sNy
372       &      -aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)              DO i=1,sNx
373       &      -aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)               surfTerm(i,j) = surfTerm(i,j)
374       &      -aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &         +cg3d_s(i,j,k,bi,bj)*drF(k)*h0FacC(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)  
375              ENDDO              ENDDO
376             ENDDO             ENDDO
377            ENDDO            ENDDO
378              DO j=1,sNy
379               DO i=1,sNx
380                 ks = kSurfC(i,j,bi,bj)
381                 surfTerm(i,j) = surfTerm(i,j)*cg3dNorm
382         &              *recip_Rcol(i,j,bi,bj)*recip_Rcol(i,j,bi,bj)
383         &              *rA(i,j,bi,bj)*deepFac2F(ks)
384         &              *recip_Bo(i,j,bi,bj)/deltaTMom/deltaTfreesurf
385               ENDDO
386              ENDDO
387             ENDIF
388    #endif /* NONLIN_FRSURF */
389             IF ( Nr .GT. 1 ) THEN
390              k=1
391    #ifdef TARGET_NEC_SX
392    !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
393    #endif /* TARGET_NEC_SX */
394              DO j=1,sNy
395               DO i=1,sNx
396                cg3d_q(i,j,k,bi,bj) =
397         &        aW3d( i, j, k, bi,bj)*cg3d_s(i-1,j, k, bi,bj)
398         &       +aW3d(i+1,j, k, bi,bj)*cg3d_s(i+1,j, k, bi,bj)
399         &       +aS3d( i, j, k, bi,bj)*cg3d_s( i,j-1,k, bi,bj)
400         &       +aS3d( i,j+1,k, bi,bj)*cg3d_s( i,j+1,k, bi,bj)
401         &       +aV3d( i, j,k+1,bi,bj)*cg3d_s( i, j,k+1,bi,bj)
402         &       +aC3d( i, j, k, bi,bj)*cg3d_s( i, j, k, bi,bj)
403    #ifdef NONLIN_FRSURF
404         &       -surfFac*surfTerm(i,j)*drF(k)*h0FacC(i,j,k,bi,bj)
405    #endif /* NONLIN_FRSURF */
406                alphaTile(bi,bj) = alphaTile(bi,bj)
407         &             +cg3d_s(i,j,k,bi,bj)*cg3d_q(i,j,k,bi,bj)
408               ENDDO
409              ENDDO
410           ELSE           ELSE
411            DO K=1,1            k=1
412             DO J=1,sNy  #ifdef TARGET_NEC_SX
413              DO I=1,sNx  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
414               cg3d_q(I,J,K,bi,bj) =  #endif /* TARGET_NEC_SX */
415       &       aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I-1,J  ,K  ,bi,bj)            DO j=1,sNy
416       &      +aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I+1,J  ,K  ,bi,bj)             DO i=1,sNx
417       &      +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J-1,K  ,bi,bj)              cg3d_q(i,j,k,bi,bj) =
418       &      +aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J+1,K  ,bi,bj)       &        aW3d( i, j, k, bi,bj)*cg3d_s(i-1,j, k, bi,bj)
419       &      -aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &       +aW3d(i+1,j, k, bi,bj)*cg3d_s(i+1,j, k, bi,bj)
420       &      -aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &       +aS3d( i, j, k, bi,bj)*cg3d_s( i,j-1,k, bi,bj)
421       &      -aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &       +aS3d( i,j+1,k, bi,bj)*cg3d_s( i,j+1,k, bi,bj)
422       &      -aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &       +aC3d( i, j, k, bi,bj)*cg3d_s( i, j, k, bi,bj)
423       &      -topLevTerm*_rA(I,J,bi,bj)*cg3d_s(I,J,K,bi,bj)  #ifdef NONLIN_FRSURF
424               alpha = alpha+cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)       &       -surfFac*surfTerm(i,j)*drF(k)*h0FacC(i,j,k,bi,bj)
425              ENDDO  #endif /* NONLIN_FRSURF */
426                alphaTile(bi,bj) = alphaTile(bi,bj)
427         &             +cg3d_s(i,j,k,bi,bj)*cg3d_q(i,j,k,bi,bj)
428             ENDDO             ENDDO
429            ENDDO            ENDDO
430           ENDIF           ENDIF
431           DO K=2,Nr-1           DO k=2,Nr-1
432            DO J=1,sNy  #ifdef TARGET_NEC_SX
433             DO I=1,sNx  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
434              cg3d_q(I,J,K,bi,bj) =  #endif /* TARGET_NEC_SX */
435       &      aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I-1,J  ,K  ,bi,bj)            DO j=1,sNy
436       &     +aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I+1,J  ,K  ,bi,bj)             DO i=1,sNx
437       &     +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J-1,K  ,bi,bj)              cg3d_q(i,j,k,bi,bj) =
438       &     +aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J+1,K  ,bi,bj)       &        aW3d( i, j, k, bi,bj)*cg3d_s(i-1,j, k, bi,bj)
439       &     +aV3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K-1,bi,bj)       &       +aW3d(i+1,j, k, bi,bj)*cg3d_s(i+1,j, k, bi,bj)
440       &     +aV3d(I  ,J  ,K+1,bi,bj)*cg3d_s(I  ,J  ,K+1,bi,bj)       &       +aS3d( i, j, k, bi,bj)*cg3d_s( i,j-1,k, bi,bj)
441       &     -aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &       +aS3d( i,j+1,k, bi,bj)*cg3d_s( i,j+1,k, bi,bj)
442       &     -aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &       +aV3d( i, j, k, bi,bj)*cg3d_s( i, j,k-1,bi,bj)
443       &     -aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &       +aV3d( i, j,k+1,bi,bj)*cg3d_s( i, j,k+1,bi,bj)
444       &     -aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &       +aC3d( i, j, k, bi,bj)*cg3d_s( i, j, k, bi,bj)
445       &     -aV3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)  #ifdef NONLIN_FRSURF
446       &     -aV3d(I  ,J  ,K+1,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &       -surfFac*surfTerm(i,j)*drF(k)*h0FacC(i,j,k,bi,bj)
447              alpha = alpha+cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)  #endif /* NONLIN_FRSURF */
448                alphaTile(bi,bj) = alphaTile(bi,bj)
449         &             +cg3d_s(i,j,k,bi,bj)*cg3d_q(i,j,k,bi,bj)
450             ENDDO             ENDDO
451            ENDDO            ENDDO
452           ENDDO           ENDDO
453           IF ( Nr .GT. 1 ) THEN           IF ( Nr .GT. 1 ) THEN
454            DO K=Nr,Nr            k=Nr
455             DO J=1,sNy  #ifdef TARGET_NEC_SX
456              DO I=1,sNx  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
457               cg3d_q(I,J,K,bi,bj) =  #endif /* TARGET_NEC_SX */
458       &       aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I-1,J  ,K  ,bi,bj)            DO j=1,sNy
459       &      +aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I+1,J  ,K  ,bi,bj)             DO i=1,sNx
460       &      +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J-1,K  ,bi,bj)              cg3d_q(i,j,k,bi,bj) =
461       &      +aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J+1,K  ,bi,bj)       &        aW3d( i, j, k, bi,bj)*cg3d_s(i-1,j, k, bi,bj)
462       &      +aV3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K-1,bi,bj)       &       +aW3d(i+1,j, k, bi,bj)*cg3d_s(i+1,j, k, bi,bj)
463       &      -aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &       +aS3d( i, j, k, bi,bj)*cg3d_s( i,j-1,k, bi,bj)
464       &      -aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &       +aS3d( i,j+1,k, bi,bj)*cg3d_s( i,j+1,k, bi,bj)
465       &      -aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &       +aV3d( i, j, k, bi,bj)*cg3d_s( i, j,k-1,bi,bj)
466       &      -aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)       &       +aC3d( i, j, k, bi,bj)*cg3d_s( i, j, k, bi,bj)
467       &      -aV3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)  #ifdef NONLIN_FRSURF
468               alpha = alpha+cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)       &       -surfFac*surfTerm(i,j)*drF(k)*h0FacC(i,j,k,bi,bj)
469              ENDDO  #endif /* NONLIN_FRSURF */
470                alphaTile(bi,bj) = alphaTile(bi,bj)
471         &             +cg3d_s(i,j,k,bi,bj)*cg3d_q(i,j,k,bi,bj)
472             ENDDO             ENDDO
473            ENDDO            ENDDO
474           ENDIF           ENDIF
475          ENDDO          ENDDO
476         ENDDO         ENDDO
477         _GLOBAL_SUM_R8(alpha,myThid)         CALL GLOBAL_SUM_TILE_RL( alphaTile,  alpha,  myThid )
478  CcnhDebugStarts  CcnhDebugStarts
479  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' SUM(s*q)= ',alpha  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' SUM(s*q)= ',alpha
480  CcnhDebugEnds  CcnhDebugEnds
# Line 436  CcnhDebugEnds Line 482  CcnhDebugEnds
482  CcnhDebugStarts  CcnhDebugStarts
483  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' alpha= ',alpha  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' alpha= ',alpha
484  CcnhDebugEnds  CcnhDebugEnds
485        
486  C==    Update solution and residual vectors  C==    Update solution and residual vectors
487  C      Now compute "interior" points.  C      Now compute "interior" points.
488         err = 0. _d 0         err = 0. _d 0
489         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
490          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
491           DO K=1,Nr           errTile(bi,bj) = 0. _d 0
492            DO J=1,sNy           DO k=1,Nr
493             DO I=1,sNx  #ifdef TARGET_NEC_SX
494              cg3d_x(I,J,K,bi,bj)=cg3d_x(I,J,K,bi,bj)  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
495       &            +alpha*cg3d_s(I,J,K,bi,bj)  #endif /* TARGET_NEC_SX */
496              cg3d_r(I,J,K,bi,bj)=cg3d_r(I,J,K,bi,bj)            DO j=1,sNy
497       &            -alpha*cg3d_q(I,J,K,bi,bj)             DO i=1,sNx
498              err = err+cg3d_r(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)              cg3d_x(i,j,k,bi,bj)=cg3d_x(i,j,k,bi,bj)
499         &            +alpha*cg3d_s(i,j,k,bi,bj)
500                cg3d_r(i,j,k,bi,bj)=cg3d_r(i,j,k,bi,bj)
501         &            -alpha*cg3d_q(i,j,k,bi,bj)
502                errTile(bi,bj) = errTile(bi,bj)
503         &            +cg3d_r(i,j,k,bi,bj)*cg3d_r(i,j,k,bi,bj)
504             ENDDO             ENDDO
505            ENDDO            ENDDO
506           ENDDO           ENDDO
507          ENDDO          ENDDO
508         ENDDO         ENDDO
509    
510         _GLOBAL_SUM_R8( err   , myThid )         CALL GLOBAL_SUM_TILE_RL( errTile,    err,    myThid )
511         err = SQRT(err)         err = SQRT(err)
512         actualIts      = it3d         actualIts      = it3d
513         actualResidual = err         actualResidual = err
514           IF ( printResidual ) THEN
515            IF ( MOD( it3d-1, printResidualFreq ).EQ.0 ) THEN
516             WRITE(msgBuf,'(A,I6,A,1PE21.14)')
517         &    ' cg3d: iter=', actualIts, ' ; resid.= ', actualResidual
518             CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
519         &                       SQUEEZE_RIGHT, myThid )
520            ENDIF
521           ENDIF
522         IF ( actualResidual .LT. cg3dTargetResidual ) GOTO 11         IF ( actualResidual .LT. cg3dTargetResidual ) GOTO 11
523  C      _EXCH_XYZ_R8(cg3d_r, myThid )         CALL EXCH_S3D_RL( cg3d_r, Nr, 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 )  
524    
525     10 CONTINUE     10 CONTINUE
526     11 CONTINUE     11 CONTINUE
527    
528          IF ( debugLevel.GE.debLevC .AND. diagFreq.GT.0. ) THEN
529            CALL WRITE_FLD_S3D_RL(
530         I        'cg3d_r_F', 'I10', 1, Nr, cg3d_r, myIter, myThid )
531          ENDIF
532    
533  C--   Un-normalise the answer  C--   Un-normalise the answer
534        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
535         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
536          DO K=1,Nr          DO k=1,Nr
537           DO J=1,sNy           DO j=1,sNy
538            DO I=1,sNx            DO i=1,sNx
539             cg3d_x(I,J,K,bi,bj) = cg3d_x(I,J,K,bi,bj)/rhsNorm             cg3d_x(i,j,k,bi,bj) = cg3d_x(i,j,k,bi,bj)/rhsNorm
540            ENDDO            ENDDO
541           ENDDO           ENDDO
542          ENDDO          ENDDO
543         ENDDO         ENDDO
544        ENDDO        ENDDO
545    
546  Cadj  _EXCH_XYZ_R8(cg3d_x, myThid )        lastResidual = actualResidual
547  c     _BEGIN_MASTER( myThid )        numIters = actualIts
 c      WRITE(*,'(A,I6,1PE30.14)') ' CG3D iters, err = ',  
 c    &  actualIts, actualResidual  
 c     _END_MASTER( )  
       lastResidual=actualResidual  
       numIters=actualIts  
548    
549  #endif /* ALLOW_NONHYDROSTATIC */  #endif /* ALLOW_NONHYDROSTATIC */
550    
Legend:
Removed from v.1.12.20.1  
changed lines
  Added in v.1.24
  ViewVC Help
Powered by ViewVC 1.1.22