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

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

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

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