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