/[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.18 by jmc, Wed Aug 23 15:22:32 2006 UTC revision 1.23 by jmc, Sun Jan 17 21:55:48 2010 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    #ifdef TARGET_NEC_SX
6    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
# Line 12  C     !INTERFACE: Line 19  C     !INTERFACE:
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 45  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     cg3d_b    - The source term or "right hand side"  C     cg3d_x    :: The solution
62  C     cg3d_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, K, 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    maskM1, maskP1        _RL     maskM1, maskP1
101        _RL    err, errTile        _RL     err,    errTile(nSx,nSy)
102        _RL    eta_qrN, eta_qrNtile        _RL     eta_qrN,eta_qrNtile(nSx,nSy)
103        _RL    eta_qrNM1        _RL     eta_qrNM1
104        _RL    cgBeta        _RL     cgBeta
105        _RL    alpha , alphaTile        _RL     alpha , alphaTile(nSx,nSy)
106        _RL    sumRHS, sumRHStile        _RL     sumRHS, sumRHStile(nSx,nSy)
107        _RL    rhsMax        _RL     rhsMax
108        _RL    rhsNorm        _RL     rhsNorm
109          CHARACTER*(MAX_LEN_MBUF) msgBuf
110          _RL     surfFac
111    #ifdef NONLIN_FRSURF
112          INTEGER ks
113          _RL     surfTerm(sNx,sNy)
114    #endif /* NONLIN_FRSURF */
115  CEOP  CEOP
116    
117          IF ( select_rStar .NE. 0 ) THEN
118            surfFac = freeSurfFac
119          ELSE
120            surfFac = 0.
121          ENDIF
122    #ifdef NONLIN_FRSURF
123          DO j=1,sNy
124           DO i=1,sNx
125             surfTerm(i,j) = 0.
126           ENDDO
127          ENDDO
128    #endif /* NONLIN_FRSURF */
129    
130  C--   Initialise inverter  C--   Initialise inverter
131        eta_qrNM1 = 1. D0        eta_qrNM1 = 1. _d 0
132    
133  C--   Normalise RHS  C--   Normalise RHS
134        rhsMax = 0. _d 0        rhsMax = 0. _d 0
135        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
136         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
137          DO K=1,Nr          DO k=1,Nr
138           DO J=1,sNy           DO j=1,sNy
139            DO I=1,sNx            DO i=1,sNx
140             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
141       &                          * maskC(I,J,K,bi,bj)       &                          * maskC(i,j,k,bi,bj)
142             rhsMax = MAX(ABS(cg3d_b(I,J,K,bi,bj)),rhsMax)             rhsMax = MAX(ABS(cg3d_b(i,j,k,bi,bj)),rhsMax)
143            ENDDO            ENDDO
144           ENDDO           ENDDO
145          ENDDO          ENDDO
146         ENDDO         ENDDO
147        ENDDO        ENDDO
148        _GLOBAL_MAX_R8( rhsMax, myThid )        _GLOBAL_MAX_RL( rhsMax, myThid )
149        rhsNorm = 1. _d 0        rhsNorm = 1. _d 0
150        IF ( rhsMax .NE. 0. ) rhsNorm = 1. _d 0 / rhsMax        IF ( rhsMax .NE. 0. ) rhsNorm = 1. _d 0 / rhsMax
151        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
152         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
153          DO K=1,Nr          DO k=1,Nr
154           DO J=1,sNy           DO j=1,sNy
155            DO I=1,sNx            DO i=1,sNx
156             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
157             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
158            ENDDO            ENDDO
159           ENDDO           ENDDO
160          ENDDO          ENDDO
# Line 134  C--   Normalise RHS Line 162  C--   Normalise RHS
162        ENDDO        ENDDO
163    
164  C--   Update overlaps  C--   Update overlaps
165  c     _EXCH_XYZ_R8( cg3d_b, myThid )  c     _EXCH_XYZ_RL( cg3d_b, myThid )
166        _EXCH_XYZ_R8( cg3d_x, myThid )        _EXCH_XYZ_RL( cg3d_x, myThid )
167    
168  C--   Initial residual calculation (with free-Surface term)  C--   Initial residual calculation (with free-Surface term)
169        err    = 0. _d 0        err    = 0. _d 0
170        sumRHS = 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 = MAX(K-1, 1 )          IF ( select_rStar .NE. 0 ) THEN
177           Kp1 = MIN(K+1, Nr)            DO j=1,sNy
178               DO i=1,sNx
179                 surfTerm(i,j) = 0.
180               ENDDO
181              ENDDO
182              DO k=1,Nr
183               DO j=1,sNy
184                DO i=1,sNx
185                 surfTerm(i,j) = surfTerm(i,j)
186         &         +cg3d_x(i,j,k,bi,bj)*drF(k)*h0FacC(i,j,k,bi,bj)
187                ENDDO
188               ENDDO
189              ENDDO
190              DO j=1,sNy
191               DO i=1,sNx
192                 ks = ksurfC(i,j,bi,bj)
193                 surfTerm(i,j) = surfTerm(i,j)*cg3dNorm
194         &              *recip_Rcol(i,j,bi,bj)*recip_Rcol(i,j,bi,bj)
195         &              *rA(i,j,bi,bj)*deepFac2F(ks)
196         &              *recip_Bo(i,j,bi,bj)/deltaTMom/deltaTfreesurf
197               ENDDO
198              ENDDO
199            ENDIF
200    #endif /* NONLIN_FRSURF */
201            DO k=1,Nr
202             km1 = MAX(k-1, 1 )
203             kp1 = MIN(k+1, Nr)
204           maskM1 = 1. _d 0           maskM1 = 1. _d 0
205           maskP1 = 1. _d 0           maskP1 = 1. _d 0
206           IF ( K .EQ. 1 ) maskM1 = 0. _d 0           IF ( k .EQ. 1 ) maskM1 = 0. _d 0
207           IF ( K .EQ. Nr) maskP1 = 0. _d 0           IF ( k .EQ. Nr) maskP1 = 0. _d 0
208              #ifdef TARGET_NEC_SX
209           DO J=1,sNy  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
210            DO I=1,sNx  #endif /* TARGET_NEC_SX */
211             cg3d_r(I,J,K,bi,bj) = cg3d_b(I,J,K,bi,bj) -( 0.           DO j=1,sNy
212       &     +aW3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I-1,J  ,K  ,bi,bj)            DO i=1,sNx
213       &     +aW3d(I+1,J  ,K  ,bi,bj)*cg3d_x(I+1,J  ,K  ,bi,bj)             cg3d_r(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
214       &     +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J-1,K  ,bi,bj)       &     -( 0.
215       &     +aS3d(I  ,J+1,K  ,bi,bj)*cg3d_x(I  ,J+1,K  ,bi,bj)       &       +aW3d( i, j, k, bi,bj)*cg3d_x(i-1,j, k, bi,bj)
216       &     +aV3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J  ,Km1,bi,bj)*maskM1       &       +aW3d(i+1,j, k, bi,bj)*cg3d_x(i+1,j, k, bi,bj)
217       &     +aV3d(I  ,J  ,Kp1,bi,bj)*cg3d_x(I  ,J  ,Kp1,bi,bj)*maskP1       &       +aS3d( i, j, k, bi,bj)*cg3d_x( i,j-1,k, bi,bj)
218       &     +aC3d(I  ,J  ,K  ,bi,bj)*cg3d_x(I  ,J  ,K  ,bi,bj)       &       +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             errTile = errTile       &       +aV3d( i, j,kp1,bi,bj)*cg3d_x( i, j,kp1,bi,bj)*maskP1
221       &     +cg3d_r(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)       &       +aC3d( i, j, k, bi,bj)*cg3d_x( i, j, k, bi,bj)
222             sumRHStile = sumRHStile  #ifdef NONLIN_FRSURF
223       &     +cg3d_b(I,J,K,bi,bj)       &       -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            ENDDO
230           ENDDO           ENDDO
231           DO J=1-1,sNy+1           DO J=1-1,sNy+1
232            DO I=1-1,sNx+1            DO I=1-1,sNx+1
233             cg3d_s(I,J,K,bi,bj) = 0.             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         CALL EXCH_S3D_RL( cg3d_r, Nr, myThid )         CALL EXCH_S3D_RL( cg3d_r, Nr, myThid )
240  c      CALL EXCH_S3D_RL( cg3d_s, Nr, myThid )  c      CALL EXCH_S3D_RL( cg3d_s, Nr, myThid )
241        _GLOBAL_SUM_R8( sumRHS, myThid )         CALL GLOBAL_SUM_TILE_RL( sumRHStile, sumRHS, myThid )
242        _GLOBAL_SUM_R8( err   , myThid )         CALL GLOBAL_SUM_TILE_RL( errTile,    err,    myThid )
243          IF ( debugLevel.GT.debLevB .AND. diagFreq.GT.0. ) THEN
244            CALL WRITE_FLD_S3D_RL(
245         I        'cg3d_r_I', 'I10', 1, Nr, cg3d_r, myIter, myThid )
246          ENDIF
247    
248        IF ( debugLevel .GE. debLevZero ) THEN        IF ( debugLevel .GE. debLevZero ) THEN
249          _BEGIN_MASTER( myThid )          _BEGIN_MASTER( myThid )
250          write(standardmessageunit,'(A,1P2E22.14)')          WRITE(standardmessageunit,'(A,1P2E22.14)')
251       &     ' cg3d: Sum(rhs),rhsMax = ',sumRHS,rhsMax       &     ' cg3d: Sum(rhs),rhsMax = ',sumRHS,rhsMax
252          _END_MASTER( myThid )          _END_MASTER( myThid )
253        ENDIF        ENDIF
254    
255        actualIts      = 0        actualIts      = 0
256        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( myThid )  
257        firstResidual=actualResidual        firstResidual=actualResidual
258    
259  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  C     >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
# Line 211  CcnhDebugEnds Line 267  CcnhDebugEnds
267         IF ( actualResidual .LT. cg3dTargetResidual ) GOTO 11         IF ( actualResidual .LT. cg3dTargetResidual ) GOTO 11
268  C--    Solve preconditioning equation and update  C--    Solve preconditioning equation and update
269  C--    conjugate direction vector "s".  C--    conjugate direction vector "s".
270  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
271  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
272  C            step. However this entails a bit of gynamastics because we only  C            step. However this entails a bit of gynamastics because we only
273  C            want eta_qrN for the interior points.  C            want eta_qrN for the interior points.
274         eta_qrN = 0. _d 0         eta_qrN = 0. _d 0
275         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
276          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
277           eta_qrNtile = 0. _d 0           eta_qrNtile(bi,bj) = 0. _d 0
278           DO K=1,1           DO k=1,1
279            DO J=1-1,sNy+1  #ifdef TARGET_NEC_SX
280             DO I=1-1,sNx+1  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
281              cg3d_q(I,J,K,bi,bj) =  #endif /* TARGET_NEC_SX */
282       &       zMC(I  ,J  ,K,bi,bj)*cg3d_r(I  ,J  ,K,bi,bj)            DO j=1-1,sNy+1
283               DO i=1-1,sNx+1
284                cg3d_q(i,j,k,bi,bj) = zMC(i,j,k,bi,bj)
285         &                        *cg3d_r(i,j,k,bi,bj)
286             ENDDO             ENDDO
287            ENDDO            ENDDO
288           ENDDO           ENDDO
289           DO K=2,Nr           DO k=2,Nr
290            DO J=1-1,sNy+1  #ifdef TARGET_NEC_SX
291             DO I=1-1,sNx+1  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
292              cg3d_q(I,J,K,bi,bj) =  #endif /* TARGET_NEC_SX */
293       &       zMC(I,J,K,bi,bj)*(cg3d_r(I,J,K  ,bi,bj)            DO j=1-1,sNy+1
294       &      -zML(I,J,K,bi,bj)*cg3d_q(I,J,K-1,bi,bj))             DO i=1-1,sNx+1
295                cg3d_q(i,j,k,bi,bj) = zMC(i,j,k,bi,bj)
296         &                      *( cg3d_r(i,j,k,bi,bj)
297         &                         -zML(i,j,k,bi,bj)*cg3d_q(i,j,k-1,bi,bj)
298         &                       )
299             ENDDO             ENDDO
300            ENDDO            ENDDO
301           ENDDO           ENDDO
302           DO K=Nr,Nr           DO k=Nr,Nr
303  caja      IF (Nr .GT. 1) THEN  #ifdef TARGET_NEC_SX
304  caja       DO J=1-1,sNy+1  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
305  caja        DO I=1-1,sNx+1  #endif /* TARGET_NEC_SX */
306  caja         cg3d_q(I,J,K,bi,bj) =            DO j=1,sNy
307  caja &        zMC(i,j,k,bi,bj)*(cg3d_r(i,j,k  ,bi,bj)             DO i=1,sNx
308  caja &       -zML(i,j,k,bi,bj)*cg3d_q(i,j,k-1,bi,bj))              eta_qrNtile(bi,bj) = eta_qrNtile(bi,bj)
309  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)  
310             ENDDO             ENDDO
311            ENDDO            ENDDO
312           ENDDO           ENDDO
313           DO K=Nr-1,1,-1           DO k=Nr-1,1,-1
314            DO J=1-1,sNy+1  #ifdef TARGET_NEC_SX
315             DO I=1-1,sNx+1  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
316              cg3d_q(I,J,K,bi,bj) =  #endif /* TARGET_NEC_SX */
317       &      cg3d_q(I,J,K,bi,bj)            DO j=1-1,sNy+1
318       &      -zMU(I,J,K,bi,bj)*cg3d_q(I,J,K+1,bi,bj)             DO i=1-1,sNx+1
319             ENDDO              cg3d_q(i,j,k,bi,bj) = cg3d_q(i,j,k,bi,bj)
320            ENDDO       &                         -zMU(i,j,k,bi,bj)*cg3d_q(i,j,k+1,bi,bj)
321            DO J=1,sNy             ENDDO
322             DO I=1,sNx            ENDDO
323              eta_qrNtile = eta_qrNtile  #ifdef TARGET_NEC_SX
324       &      +cg3d_q(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
325    #endif /* TARGET_NEC_SX */
326              DO j=1,sNy
327               DO i=1,sNx
328                eta_qrNtile(bi,bj) = eta_qrNtile(bi,bj)
329         &                        +cg3d_q(i,j,k,bi,bj)*cg3d_r(i,j,k,bi,bj)
330             ENDDO             ENDDO
331            ENDDO            ENDDO
332           ENDDO           ENDDO
          eta_qrN = eta_qrN + eta_qrNtile  
333          ENDDO          ENDDO
334         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  
335    
336         _GLOBAL_SUM_R8(eta_qrN, myThid)         CALL GLOBAL_SUM_TILE_RL( eta_qrNtile,eta_qrN,myThid )
337  CcnhDebugStarts  CcnhDebugStarts
338  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' eta_qrN = ',eta_qrN  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' eta_qrN = ',eta_qrN
339  CcnhDebugEnds  CcnhDebugEnds
# Line 299  CcnhDebugEnds Line 345  CcnhDebugEnds
345    
346         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
347          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
348           DO K=1,Nr           DO k=1,Nr
349            DO J=1-1,sNy+1            DO j=1-1,sNy+1
350             DO I=1-1,sNx+1             DO i=1-1,sNx+1
351              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)
352       &                          + cgBeta*cg3d_s(I,J,K,bi,bj)       &                   + cgBeta*cg3d_s(i,j,k,bi,bj)
353             ENDDO             ENDDO
354            ENDDO            ENDDO
355           ENDDO           ENDDO
# Line 315  C==    q = A.s Line 361  C==    q = A.s
361         alpha = 0. _d 0         alpha = 0. _d 0
362         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
363          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
364           alphaTile = 0. _d 0           alphaTile(bi,bj) = 0. _d 0
365           IF ( Nr .GT. 1 ) THEN  #ifdef NONLIN_FRSURF
366            DO K=1,1           IF ( select_rStar .NE. 0 ) THEN
367             DO J=1,sNy            DO j=1,sNy
368              DO I=1,sNx             DO i=1,sNx
369               cg3d_q(I,J,K,bi,bj) =               surfTerm(i,j) = 0.
370       &       aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I-1,J  ,K  ,bi,bj)             ENDDO
371       &      +aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I+1,J  ,K  ,bi,bj)            ENDDO
372       &      +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J-1,K  ,bi,bj)            DO k=1,Nr
373       &      +aS3d(I  ,J+1,K  ,bi,bj)*cg3d_s(I  ,J+1,K  ,bi,bj)             DO j=1,sNy
374       &      +aV3d(I  ,J  ,K+1,bi,bj)*cg3d_s(I  ,J  ,K+1,bi,bj)              DO i=1,sNx
375       &      +aC3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J  ,K  ,bi,bj)               surfTerm(i,j) = surfTerm(i,j)
376               alphaTile = alphaTile       &         +cg3d_s(i,j,k,bi,bj)*drF(k)*h0FacC(i,j,k,bi,bj)
      &                 +cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)  
377              ENDDO              ENDDO
378             ENDDO             ENDDO
379            ENDDO            ENDDO
380              DO j=1,sNy
381               DO i=1,sNx
382                 ks = ksurfC(i,j,bi,bj)
383                 surfTerm(i,j) = surfTerm(i,j)*cg3dNorm
384         &              *recip_Rcol(i,j,bi,bj)*recip_Rcol(i,j,bi,bj)
385         &              *rA(i,j,bi,bj)*deepFac2F(ks)
386         &              *recip_Bo(i,j,bi,bj)/deltaTMom/deltaTfreesurf
387               ENDDO
388              ENDDO
389             ENDIF
390    #endif /* NONLIN_FRSURF */
391             IF ( Nr .GT. 1 ) THEN
392              k=1
393    #ifdef TARGET_NEC_SX
394    !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
395    #endif /* TARGET_NEC_SX */
396              DO j=1,sNy
397               DO i=1,sNx
398                cg3d_q(i,j,k,bi,bj) =
399         &        aW3d( i, j, k, bi,bj)*cg3d_s(i-1,j, k, bi,bj)
400         &       +aW3d(i+1,j, k, bi,bj)*cg3d_s(i+1,j, k, bi,bj)
401         &       +aS3d( i, j, k, bi,bj)*cg3d_s( i,j-1,k, bi,bj)
402         &       +aS3d( i,j+1,k, bi,bj)*cg3d_s( i,j+1,k, bi,bj)
403         &       +aV3d( i, j,k+1,bi,bj)*cg3d_s( i, j,k+1,bi,bj)
404         &       +aC3d( i, j, k, bi,bj)*cg3d_s( i, j, k, bi,bj)
405    #ifdef NONLIN_FRSURF
406         &       -surfFac*surfTerm(i,j)*drF(k)*h0FacC(i,j,k,bi,bj)
407    #endif /* NONLIN_FRSURF */
408                alphaTile(bi,bj) = alphaTile(bi,bj)
409         &             +cg3d_s(i,j,k,bi,bj)*cg3d_q(i,j,k,bi,bj)
410               ENDDO
411              ENDDO
412           ELSE           ELSE
413            DO K=1,1            k=1
414             DO J=1,sNy  #ifdef TARGET_NEC_SX
415              DO I=1,sNx  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
416               cg3d_q(I,J,K,bi,bj) =  #endif /* TARGET_NEC_SX */
417       &       aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I-1,J  ,K  ,bi,bj)            DO j=1,sNy
418       &      +aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I+1,J  ,K  ,bi,bj)             DO i=1,sNx
419       &      +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J-1,K  ,bi,bj)              cg3d_q(i,j,k,bi,bj) =
420       &      +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)
421       &      +aC3d(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)
422               alphaTile = alphaTile       &       +aS3d( i, j, k, bi,bj)*cg3d_s( i,j-1,k, bi,bj)
423       &                 +cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)       &       +aS3d( i,j+1,k, bi,bj)*cg3d_s( i,j+1,k, bi,bj)
424              ENDDO       &       +aC3d( i, j, k, bi,bj)*cg3d_s( i, j, k, bi,bj)
425    #ifdef NONLIN_FRSURF
426         &       -surfFac*surfTerm(i,j)*drF(k)*h0FacC(i,j,k,bi,bj)
427    #endif /* NONLIN_FRSURF */
428                alphaTile(bi,bj) = alphaTile(bi,bj)
429         &             +cg3d_s(i,j,k,bi,bj)*cg3d_q(i,j,k,bi,bj)
430             ENDDO             ENDDO
431            ENDDO            ENDDO
432           ENDIF           ENDIF
433           DO K=2,Nr-1           DO k=2,Nr-1
434            DO J=1,sNy  #ifdef TARGET_NEC_SX
435             DO I=1,sNx  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
436              cg3d_q(I,J,K,bi,bj) =  #endif /* TARGET_NEC_SX */
437       &      aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I-1,J  ,K  ,bi,bj)            DO j=1,sNy
438       &     +aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I+1,J  ,K  ,bi,bj)             DO i=1,sNx
439       &     +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J-1,K  ,bi,bj)              cg3d_q(i,j,k,bi,bj) =
440       &     +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)
441       &     +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)
442       &     +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)
443       &     +aC3d(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)
444              alphaTile = alphaTile       &       +aV3d( i, j, k, bi,bj)*cg3d_s( i, j,k-1,bi,bj)
445       &                +cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)       &       +aV3d( i, j,k+1,bi,bj)*cg3d_s( i, j,k+1,bi,bj)
446         &       +aC3d( i, j, k, bi,bj)*cg3d_s( i, j, k, bi,bj)
447    #ifdef NONLIN_FRSURF
448         &       -surfFac*surfTerm(i,j)*drF(k)*h0FacC(i,j,k,bi,bj)
449    #endif /* NONLIN_FRSURF */
450                alphaTile(bi,bj) = alphaTile(bi,bj)
451         &             +cg3d_s(i,j,k,bi,bj)*cg3d_q(i,j,k,bi,bj)
452             ENDDO             ENDDO
453            ENDDO            ENDDO
454           ENDDO           ENDDO
455           IF ( Nr .GT. 1 ) THEN           IF ( Nr .GT. 1 ) THEN
456            DO K=Nr,Nr            k=Nr
457             DO J=1,sNy  #ifdef TARGET_NEC_SX
458              DO I=1,sNx  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
459               cg3d_q(I,J,K,bi,bj) =  #endif /* TARGET_NEC_SX */
460       &       aW3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I-1,J  ,K  ,bi,bj)            DO j=1,sNy
461       &      +aW3d(I+1,J  ,K  ,bi,bj)*cg3d_s(I+1,J  ,K  ,bi,bj)             DO i=1,sNx
462       &      +aS3d(I  ,J  ,K  ,bi,bj)*cg3d_s(I  ,J-1,K  ,bi,bj)              cg3d_q(i,j,k,bi,bj) =
463       &      +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)
464       &      +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)
465       &      +aC3d(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)
466               alphaTile = alphaTile       &       +aS3d( i,j+1,k, bi,bj)*cg3d_s( i,j+1,k, bi,bj)
467       &                 +cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)       &       +aV3d( i, j, k, bi,bj)*cg3d_s( i, j,k-1,bi,bj)
468              ENDDO       &       +aC3d( i, j, k, bi,bj)*cg3d_s( i, j, k, bi,bj)
469    #ifdef NONLIN_FRSURF
470         &       -surfFac*surfTerm(i,j)*drF(k)*h0FacC(i,j,k,bi,bj)
471    #endif /* NONLIN_FRSURF */
472                alphaTile(bi,bj) = alphaTile(bi,bj)
473         &             +cg3d_s(i,j,k,bi,bj)*cg3d_q(i,j,k,bi,bj)
474             ENDDO             ENDDO
475            ENDDO            ENDDO
476           ENDIF           ENDIF
          alpha = alpha + alphaTile  
477          ENDDO          ENDDO
478         ENDDO         ENDDO
479         _GLOBAL_SUM_R8(alpha,myThid)         CALL GLOBAL_SUM_TILE_RL( alphaTile,  alpha,  myThid )
480  CcnhDebugStarts  CcnhDebugStarts
481  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' SUM(s*q)= ',alpha  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' SUM(s*q)= ',alpha
482  CcnhDebugEnds  CcnhDebugEnds
# Line 392  CcnhDebugEnds Line 484  CcnhDebugEnds
484  CcnhDebugStarts  CcnhDebugStarts
485  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' alpha= ',alpha  C      WRITE(*,*) ' CG3D: Iteration ',it3d-1,' alpha= ',alpha
486  CcnhDebugEnds  CcnhDebugEnds
487        
488  C==    Update solution and residual vectors  C==    Update solution and residual vectors
489  C      Now compute "interior" points.  C      Now compute "interior" points.
490         err = 0. _d 0         err = 0. _d 0
491         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
492          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
493          errTile    = 0. _d 0           errTile(bi,bj) = 0. _d 0
494           DO K=1,Nr           DO k=1,Nr
495            DO J=1,sNy  #ifdef TARGET_NEC_SX
496             DO I=1,sNx  !CDIR OUTERUNROLL=CG3D_OUTERLOOPITERS
497              cg3d_x(I,J,K,bi,bj)=cg3d_x(I,J,K,bi,bj)  #endif /* TARGET_NEC_SX */
498       &            +alpha*cg3d_s(I,J,K,bi,bj)            DO j=1,sNy
499              cg3d_r(I,J,K,bi,bj)=cg3d_r(I,J,K,bi,bj)             DO i=1,sNx
500       &            -alpha*cg3d_q(I,J,K,bi,bj)              cg3d_x(i,j,k,bi,bj)=cg3d_x(i,j,k,bi,bj)
501              errTile = errTile       &            +alpha*cg3d_s(i,j,k,bi,bj)
502       &             +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)
503         &            -alpha*cg3d_q(i,j,k,bi,bj)
504                errTile(bi,bj) = errTile(bi,bj)
505         &            +cg3d_r(i,j,k,bi,bj)*cg3d_r(i,j,k,bi,bj)
506             ENDDO             ENDDO
507            ENDDO            ENDDO
508           ENDDO           ENDDO
          err = err + errTile  
509          ENDDO          ENDDO
510         ENDDO         ENDDO
511    
512         _GLOBAL_SUM_R8( err   , myThid )         CALL GLOBAL_SUM_TILE_RL( errTile,    err,    myThid )
513         err = SQRT(err)         err = SQRT(err)
514         actualIts      = it3d         actualIts      = it3d
515         actualResidual = err         actualResidual = err
516           IF ( debugLevel.GT.debLevB ) THEN
517    c       IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
518    c    &     ) THEN
519            _BEGIN_MASTER( myThid )
520             WRITE(msgBuf,'(A,I6,A,1PE21.14)')
521         &    ' cg3d: iter=', actualIts, ' ; resid.= ', actualResidual
522             CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
523         &                       SQUEEZE_RIGHT, myThid )
524            _END_MASTER( myThid )
525    c       ENDIF
526           ENDIF
527         IF ( actualResidual .LT. cg3dTargetResidual ) GOTO 11         IF ( actualResidual .LT. cg3dTargetResidual ) GOTO 11
528         CALL EXCH_S3D_RL( cg3d_r, Nr, myThid )         CALL EXCH_S3D_RL( cg3d_r, Nr, myThid )
529    
530     10 CONTINUE     10 CONTINUE
531     11 CONTINUE     11 CONTINUE
532    
533          IF ( debugLevel.GT.debLevB .AND. diagFreq.GT.0. ) THEN
534            CALL WRITE_FLD_S3D_RL(
535         I        'cg3d_r_F', 'I10', 1, Nr, cg3d_r, myIter, myThid )
536          ENDIF
537    
538  C--   Un-normalise the answer  C--   Un-normalise the answer
539        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
540         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
541          DO K=1,Nr          DO k=1,Nr
542           DO J=1,sNy           DO j=1,sNy
543            DO I=1,sNx            DO i=1,sNx
544             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
545            ENDDO            ENDDO
546           ENDDO           ENDDO
547          ENDDO          ENDDO
548         ENDDO         ENDDO
549        ENDDO        ENDDO
550    
551  Cadj  _EXCH_XYZ_R8(cg3d_x, myThid )        lastResidual = actualResidual
552  c     _BEGIN_MASTER( myThid )        numIters = actualIts
 c      WRITE(*,'(A,I6,1PE30.14)') ' CG3D iters, err = ',  
 c    &  actualIts, actualResidual  
 c     _END_MASTER( myThid )  
       lastResidual=actualResidual  
       numIters=actualIts  
553    
554  #endif /* ALLOW_NONHYDROSTATIC */  #endif /* ALLOW_NONHYDROSTATIC */
555    
Legend:
Removed from v.1.18  
changed lines
  Added in v.1.23
  ViewVC Help
Powered by ViewVC 1.1.22