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

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

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

Revision 1.12 - (hide annotations) (download)
Wed Sep 26 18:09:14 2001 UTC (24 years, 1 month ago) by cnh
Branch: MAIN
CVS Tags: checkpoint46n_post, checkpoint47e_post, checkpoint44e_post, checkpoint46l_post, checkpoint46g_pre, checkpoint47c_post, release1_p13_pre, checkpoint50c_post, checkpoint46f_post, checkpoint48e_post, checkpoint50c_pre, checkpoint44f_post, checkpoint46b_post, checkpoint43a-release1mods, ecco_c50_e32, ecco_c50_e33, ecco_c50_e30, ecco_c50_e31, release1_p13, checkpoint48i_post, checkpoint46l_pre, chkpt44d_post, checkpoint51, checkpoint50, release1_p8, release1_p9, checkpoint50d_post, release1_p1, release1_p2, release1_p3, release1_p4, release1_p5, release1_p6, release1_p7, checkpoint50b_pre, checkpoint44e_pre, checkpoint51f_post, release1_b1, ecco_c51_e34d, ecco_c51_e34e, ecco_c51_e34f, ecco_c51_e34g, ecco_c51_e34a, ecco_c51_e34b, ecco_c51_e34c, checkpoint48b_post, checkpoint43, checkpoint51d_post, checkpoint48c_pre, checkpoint47d_pre, release1_chkpt44d_post, checkpoint47a_post, checkpoint48d_pre, checkpoint51j_post, checkpoint47i_post, release1_p11, checkpoint47d_post, icebear5, icebear4, icebear3, icebear2, checkpoint46d_pre, checkpoint48d_post, release1-branch_tutorials, checkpoint48f_post, checkpoint45d_post, checkpoint46j_pre, chkpt44a_post, checkpoint44h_pre, checkpoint48h_post, ecco_c50_e29, checkpoint51b_pre, checkpoint46a_post, checkpoint47g_post, checkpoint46j_post, checkpoint51h_pre, checkpoint46k_post, ecco_c50_e28, chkpt44c_pre, checkpoint48a_post, checkpoint45a_post, checkpoint50f_post, checkpoint50a_post, checkpoint50f_pre, ecco_c44_e19, ecco_c44_e18, ecco_c44_e17, ecco_c44_e16, release1_p12, release1_p10, release1_p16, release1_p17, release1_p14, release1_p15, checkpoint47j_post, ecco_c50_e33a, branch-exfmods-tag, checkpoint44g_post, branchpoint-genmake2, checkpoint46e_pre, checkpoint48c_post, checkpoint45b_post, checkpoint46b_pre, release1-branch-end, release1_final_v1, checkpoint51b_post, checkpoint51c_post, checkpoint46c_pre, checkpoint46, checkpoint47b_post, checkpoint44b_post, ecco_c51_e34, checkpoint46h_pre, checkpoint46m_post, checkpoint46a_pre, checkpoint50g_post, checkpoint45c_post, ecco_ice2, ecco_ice1, checkpoint44h_post, checkpoint46g_post, release1_p12_pre, ecco_c44_e22, checkpoint50h_post, checkpoint50e_pre, checkpoint50i_post, ecco_c44_e25, checkpoint51i_pre, checkpoint47f_post, checkpoint50e_post, chkpt44a_pre, checkpoint46i_post, ecco_c44_e23, ecco_c44_e20, ecco_c44_e21, ecco_c44_e26, ecco_c44_e27, ecco_c44_e24, checkpoint46c_post, ecco-branch-mod1, ecco-branch-mod2, ecco-branch-mod3, ecco-branch-mod4, ecco-branch-mod5, checkpoint50d_pre, checkpoint46e_post, release1_beta1, checkpoint51e_post, checkpoint44b_pre, checkpoint42, checkpoint41, checkpoint47, checkpoint44, checkpoint45, checkpoint48, checkpoint49, checkpoint46h_post, checkpoint51f_pre, chkpt44c_post, checkpoint48g_post, checkpoint47h_post, checkpoint44f_pre, checkpoint51g_post, checkpoint46d_post, checkpoint50b_post, release1-branch_branchpoint, checkpoint51a_post
Branch point for: c24_e25_ice, branch-exfmods-curt, release1_final, release1-branch, branch-genmake2, release1, ecco-branch, release1_50yr, icebear, release1_coupled
Changes since 1.11: +33 -24 lines 
Bringing comments up to data and formatting for document extraction.
1 cnh 1.12 C $Header: /u/gcmpack/models/MITgcmUV/model/src/cg3d.F,v 1.11 2001/06/29 17:14:49 adcroft Exp $
2 adcroft 1.10 C $Name: $
3 adcroft 1.1
4     #include "CPP_OPTIONS.h"
5    
6     #define VERBOSE
7    
8 cnh 1.12 CBOP
9     C !ROUTINE: CG3D
10     C !INTERFACE:
11 adcroft 1.1 SUBROUTINE CG3D(
12 adcroft 1.11 I cg3d_b,
13     U cg3d_x,
14     O firstResidual,
15     O lastResidual,
16     U numIters,
17 adcroft 1.1 I myThid )
18 cnh 1.12 C !DESCRIPTION: \bv
19     C *==========================================================*
20     C | SUBROUTINE CG3D
21     C | o Three-dimensional grid problem conjugate-gradient
22     C | inverter (with preconditioner).
23     C *==========================================================*
24     C | Con. grad is an iterative procedure for solving Ax = b.
25     C | It requires the A be symmetric.
26     C | This implementation assumes A is a seven-diagonal
27     C | matrix of the form that arises in the discrete
28     C | representation of the del^2 operator in a
29     C | three-dimensional space.
30     C | Notes:
31     C | ======
32     C | This implementation can support shared-memory
33     C | multi-threaded execution. In order to do this COMMON
34     C | blocks are used for many of the arrays - even ones that
35     C | are only used for intermedaite results. This design is
36     C | OK if you want to all the threads to collaborate on
37     C | solving the same problem. On the other hand if you want
38     C | the threads to solve several different problems
39     C | concurrently this implementation will not work.
40     C *==========================================================*
41     C \ev
42    
43     C !USES:
44 adcroft 1.1 IMPLICIT NONE
45     C === Global data ===
46     #include "SIZE.h"
47     #include "EEPARAMS.h"
48     #include "PARAMS.h"
49     #include "GRID.h"
50     #include "CG3D.h"
51    
52 cnh 1.12 C !INPUT/OUTPUT PARAMETERS:
53 adcroft 1.1 C === Routine arguments ===
54 adcroft 1.11 C myThid - Thread on which I am working.
55     C cg2d_b - The source term or "right hand side"
56     C cg2d_x - The solution
57     C firstResidual - the initial residual before any iterations
58     C lastResidual - the actual residual reached
59     C numIters - Entry: the maximum number of iterations allowed
60     C Exit: the actual number of iterations used
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 adcroft 1.1 INTEGER myThid
67    
68 adcroft 1.11
69 adcroft 1.4 #ifdef ALLOW_NONHYDROSTATIC
70    
71 cnh 1.12 C !LOCAL VARIABLES:
72 adcroft 1.1 C === Local variables ====
73     C actualIts - Number of iterations taken
74     C actualResidual - residual
75     C bi - Block index in X and Y.
76     C bj
77 jmc 1.6 C eta_qrN - Used in computing search directions
78     C eta_qrNM1 suffix N and NM1 denote current and
79 adcroft 1.1 C cgBeta previous iterations respectively.
80     C alpha
81     C sumRHS - Sum of right-hand-side. Sometimes this is a
82     C useful debuggin/trouble shooting diagnostic.
83     C For neumann problems sumRHS needs to be ~0.
84     C or they converge at a non-zero residual.
85     C err - Measure of residual of Ax - b, usually the norm.
86     C I, J, N - Loop counters ( N counts CG iterations )
87     INTEGER actualIts
88     _RL actualResidual
89     INTEGER bi, bj
90     INTEGER I, J, K, it3d
91     INTEGER KM1, KP1
92     _RL err
93 jmc 1.6 _RL eta_qrN
94     _RL eta_qrNM1
95 adcroft 1.1 _RL cgBeta
96     _RL alpha
97     _RL sumRHS
98     _RL rhsMax
99     _RL rhsNorm
100    
101     INTEGER OLw
102     INTEGER OLe
103     INTEGER OLn
104     INTEGER OLs
105     INTEGER exchWidthX
106     INTEGER exchWidthY
107     INTEGER myNz
108 jmc 1.7 _RL topLevTerm
109 cnh 1.12 CEOP
110 adcroft 1.1
111     C-- Initialise inverter
112 jmc 1.6 eta_qrNM1 = 1. D0
113 adcroft 1.1
114     C-- Normalise RHS
115     rhsMax = 0. _d 0
116     DO bj=myByLo(myThid),myByHi(myThid)
117     DO bi=myBxLo(myThid),myBxHi(myThid)
118     DO K=1,Nr
119     DO J=1,sNy
120     DO I=1,sNx
121     cg3d_b(I,J,K,bi,bj) = cg3d_b(I,J,K,bi,bj)*cg3dNorm
122     rhsMax = MAX(ABS(cg3d_b(I,J,K,bi,bj)),rhsMax)
123     ENDDO
124     ENDDO
125     ENDDO
126     ENDDO
127     ENDDO
128 adcroft 1.2 _GLOBAL_MAX_R8( rhsMax, myThid )
129 adcroft 1.1 rhsNorm = 1. _d 0
130     IF ( rhsMax .NE. 0. ) rhsNorm = 1. _d 0 / rhsMax
131     DO bj=myByLo(myThid),myByHi(myThid)
132     DO bi=myBxLo(myThid),myBxHi(myThid)
133     DO K=1,Nr
134     DO J=1,sNy
135     DO I=1,sNx
136     cg3d_b(I,J,K,bi,bj) = cg3d_b(I,J,K,bi,bj)*rhsNorm
137     cg3d_x(I,J,K,bi,bj) = cg3d_x(I,J,K,bi,bj)*rhsNorm
138     ENDDO
139     ENDDO
140     ENDDO
141     ENDDO
142     ENDDO
143    
144     C-- Update overlaps
145     _EXCH_XYZ_R8( cg3d_b, myThid )
146     _EXCH_XYZ_R8( cg3d_x, myThid )
147    
148 jmc 1.7 C-- Initial residual calculation (with free-Surface term)
149 adcroft 1.1 err = 0. _d 0
150     sumRHS = 0. _d 0
151     DO bj=myByLo(myThid),myByHi(myThid)
152     DO bi=myBxLo(myThid),myBxHi(myThid)
153     DO K=1,Nr
154     KM1 = K-1
155     IF ( K .EQ. 1 ) KM1 = 1
156     KP1 = K+1
157     IF ( K .EQ. Nr ) KP1 = 1
158 jmc 1.7 topLevTerm = 0.
159     IF ( K .EQ. 1) topLevTerm = freeSurfFac*cg3dNorm*
160     & (horiVertRatio/gravity)/deltaTMom/deltaTMom
161 adcroft 1.1 DO J=1,sNy
162     DO I=1,sNx
163     cg3d_s(I,J,K,bi,bj) = 0.
164     cg3d_r(I,J,K,bi,bj) = cg3d_b(I,J,K,bi,bj) -( 0.
165     & +aW3d(I ,J ,K ,bi,bj)*cg3d_x(I-1,J ,K ,bi,bj)
166     & +aW3d(I+1,J ,K ,bi,bj)*cg3d_x(I+1,J ,K ,bi,bj)
167     & +aS3d(I ,J ,K ,bi,bj)*cg3d_x(I ,J-1,K ,bi,bj)
168     & +aS3d(I ,J+1,K ,bi,bj)*cg3d_x(I ,J+1,K ,bi,bj)
169     & +aV3d(I ,J ,K ,bi,bj)*cg3d_x(I ,J ,KM1,bi,bj)
170     & +aV3d(I ,J ,KP1,bi,bj)*cg3d_x(I ,J ,KP1,bi,bj)
171     & -aW3d(I ,J ,K ,bi,bj)*cg3d_x(I ,J ,K ,bi,bj)
172     & -aW3d(I+1,J ,K ,bi,bj)*cg3d_x(I ,J ,K ,bi,bj)
173     & -aS3d(I ,J ,K ,bi,bj)*cg3d_x(I ,J ,K ,bi,bj)
174     & -aS3d(I ,J+1,K ,bi,bj)*cg3d_x(I ,J ,K ,bi,bj)
175     & -aV3d(I ,J ,K ,bi,bj)*cg3d_x(I ,J ,K ,bi,bj)
176     & -aV3d(I ,J ,KP1,bi,bj)*cg3d_x(I ,J ,K ,bi,bj)
177 jmc 1.7 & -topLevTerm*_rA(I,J,bi,bj)*cg3d_x(I,J,K,bi,bj)
178 adcroft 1.1 & )
179     err = err
180     & +cg3d_r(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)
181     sumRHS = sumRHS
182     & +cg3d_b(I,J,K,bi,bj)
183     ENDDO
184     ENDDO
185     ENDDO
186     ENDDO
187     ENDDO
188     C _EXCH_XYZ_R8( cg3d_r, myThid )
189     OLw = 1
190     OLe = 1
191     OLn = 1
192     OLs = 1
193     exchWidthX = 1
194     exchWidthY = 1
195     myNz = Nr
196     CALL EXCH_RL( cg3d_r,
197 adcroft 1.10 I OLw, OLe, OLs, OLn, myNz,
198 adcroft 1.1 I exchWidthX, exchWidthY,
199     I FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )
200     C _EXCH_XYZ_R8( cg3d_s, myThid )
201     OLw = 1
202     OLe = 1
203     OLn = 1
204     OLs = 1
205     exchWidthX = 1
206     exchWidthY = 1
207     myNz = Nr
208     CALL EXCH_RL( cg3d_s,
209 adcroft 1.10 I OLw, OLe, OLs, OLn, myNz,
210 adcroft 1.1 I exchWidthX, exchWidthY,
211     I FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )
212 adcroft 1.2 _GLOBAL_SUM_R8( sumRHS, myThid )
213     _GLOBAL_SUM_R8( err , myThid )
214 adcroft 1.1
215     _BEGIN_MASTER( myThid )
216 adcroft 1.11 write(*,'(A,1P2E22.14)')
217     & ' cg3d: Sum(rhs),rhsMax = ',sumRHS,rhsMax
218 adcroft 1.1 _END_MASTER( )
219    
220     actualIts = 0
221     actualResidual = SQRT(err)
222     C _BARRIER
223 adcroft 1.11 c _BEGIN_MASTER( myThid )
224     c WRITE(*,'(A,I6,1PE30.14)') ' CG3D iters, err = ',
225     c & actualIts, actualResidual
226     c _END_MASTER( )
227     firstResidual=actualResidual
228 adcroft 1.1
229     C >>>>>>>>>>>>>>> BEGIN SOLVER <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
230     DO 10 it3d=1, cg3dMaxIters
231    
232     CcnhDebugStarts
233     #ifdef VERBOSE
234 adcroft 1.11 c IF ( mod(it3d-1,10).EQ.0)
235     c & WRITE(*,*) ' CG3D: Iteration ',it3d-1,
236     c & ' residual = ',actualResidual
237 adcroft 1.1 #endif
238     CcnhDebugEnds
239     IF ( actualResidual .LT. cg3dTargetResidual ) GOTO 11
240     C-- Solve preconditioning equation and update
241     C-- conjugate direction vector "s".
242     C Note. On the next to loops over all tiles the inner loop ranges
243     C in sNx and sNy are expanded by 1 to avoid a communication
244     C step. However this entails a bit of gynamastics because we only
245 jmc 1.6 C want eta_qrN for the interior points.
246     eta_qrN = 0. _d 0
247 adcroft 1.1 DO bj=myByLo(myThid),myByHi(myThid)
248     DO bi=myBxLo(myThid),myBxHi(myThid)
249     DO K=1,1
250     DO J=1-1,sNy+1
251     DO I=1-1,sNx+1
252     cg3d_q(I,J,K,bi,bj) =
253     & zMC(I ,J ,K,bi,bj)*cg3d_r(I ,J ,K,bi,bj)
254     ENDDO
255     ENDDO
256     ENDDO
257     DO K=2,Nr
258     DO J=1-1,sNy+1
259     DO I=1-1,sNx+1
260     cg3d_q(I,J,K,bi,bj) =
261     & zMC(I,J,K,bi,bj)*(cg3d_r(I,J,K ,bi,bj)
262     & -zML(I,J,K,bi,bj)*cg3d_q(I,J,K-1,bi,bj))
263     ENDDO
264     ENDDO
265     ENDDO
266     DO K=Nr,Nr
267     caja IF (Nr .GT. 1) THEN
268     caja DO J=1-1,sNy+1
269     caja DO I=1-1,sNx+1
270     caja cg3d_q(I,J,K,bi,bj) =
271     caja & zMC(i,j,k,bi,bj)*(cg3d_r(i,j,k ,bi,bj)
272     caja & -zML(i,j,k,bi,bj)*cg3d_q(i,j,k-1,bi,bj))
273     caja ENDDO
274     caja ENDDO
275     caja ENDIF
276     DO J=1,sNy
277     DO I=1,sNx
278 jmc 1.6 eta_qrN = eta_qrN
279 adcroft 1.1 & +cg3d_q(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)
280     ENDDO
281     ENDDO
282     ENDDO
283     DO K=Nr-1,1,-1
284     DO J=1-1,sNy+1
285     DO I=1-1,sNx+1
286     cg3d_q(I,J,K,bi,bj) =
287     & cg3d_q(I,J,K,bi,bj)
288     & -zMU(I,J,K,bi,bj)*cg3d_q(I,J,K+1,bi,bj)
289     ENDDO
290     ENDDO
291     DO J=1,sNy
292     DO I=1,sNx
293 jmc 1.6 eta_qrN = eta_qrN
294 adcroft 1.1 & +cg3d_q(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)
295     ENDDO
296     ENDDO
297     ENDDO
298     ENDDO
299     ENDDO
300     caja
301 jmc 1.6 caja eta_qrN=0.
302 adcroft 1.1 caja DO bj=myByLo(myThid),myByHi(myThid)
303     caja DO bi=myBxLo(myThid),myBxHi(myThid)
304     caja DO K=1,Nr
305     caja DO J=1,sNy
306     caja DO I=1,sNx
307 jmc 1.6 caja eta_qrN = eta_qrN
308 adcroft 1.1 caja & +cg3d_q(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)
309     caja ENDDO
310     caja ENDDO
311     caja ENDDO
312     caja ENDDO
313     caja ENDDO
314     caja
315    
316 jmc 1.6 _GLOBAL_SUM_R8(eta_qrN, myThid)
317 adcroft 1.1 CcnhDebugStarts
318 heimbach 1.8 C WRITE(*,*) ' CG3D: Iteration ',it3d-1,' eta_qrN = ',eta_qrN
319 adcroft 1.1 CcnhDebugEnds
320 jmc 1.6 cgBeta = eta_qrN/eta_qrNM1
321 adcroft 1.1 CcnhDebugStarts
322 heimbach 1.8 C WRITE(*,*) ' CG3D: Iteration ',it3d-1,' beta = ',cgBeta
323 adcroft 1.1 CcnhDebugEnds
324 jmc 1.6 eta_qrNM1 = eta_qrN
325 adcroft 1.1
326     DO bj=myByLo(myThid),myByHi(myThid)
327     DO bi=myBxLo(myThid),myBxHi(myThid)
328     DO K=1,Nr
329     DO J=1-1,sNy+1
330     DO I=1-1,sNx+1
331     cg3d_s(I,J,K,bi,bj) = cg3d_q(I,J,K,bi,bj)
332     & + cgBeta*cg3d_s(I,J,K,bi,bj)
333     ENDDO
334     ENDDO
335     ENDDO
336     ENDDO
337     ENDDO
338    
339     C== Evaluate laplace operator on conjugate gradient vector
340     C== q = A.s
341     alpha = 0. _d 0
342 jmc 1.7 topLevTerm = freeSurfFac*cg3dNorm*
343     & (horiVertRatio/gravity)/deltaTMom/deltaTMom
344 adcroft 1.1 DO bj=myByLo(myThid),myByHi(myThid)
345     DO bi=myBxLo(myThid),myBxHi(myThid)
346     IF ( Nr .GT. 1 ) THEN
347     DO K=1,1
348     DO J=1,sNy
349     DO I=1,sNx
350     cg3d_q(I,J,K,bi,bj) =
351     & aW3d(I ,J ,K ,bi,bj)*cg3d_s(I-1,J ,K ,bi,bj)
352     & +aW3d(I+1,J ,K ,bi,bj)*cg3d_s(I+1,J ,K ,bi,bj)
353     & +aS3d(I ,J ,K ,bi,bj)*cg3d_s(I ,J-1,K ,bi,bj)
354     & +aS3d(I ,J+1,K ,bi,bj)*cg3d_s(I ,J+1,K ,bi,bj)
355     & +aV3d(I ,J ,K+1,bi,bj)*cg3d_s(I ,J ,K+1,bi,bj)
356     & -aW3d(I ,J ,K ,bi,bj)*cg3d_s(I ,J ,K ,bi,bj)
357     & -aW3d(I+1,J ,K ,bi,bj)*cg3d_s(I ,J ,K ,bi,bj)
358     & -aS3d(I ,J ,K ,bi,bj)*cg3d_s(I ,J ,K ,bi,bj)
359     & -aS3d(I ,J+1,K ,bi,bj)*cg3d_s(I ,J ,K ,bi,bj)
360     & -aV3d(I ,J ,K+1,bi,bj)*cg3d_s(I ,J ,K ,bi,bj)
361 jmc 1.7 & -topLevTerm*_rA(I,J,bi,bj)*cg3d_s(I,J,K,bi,bj)
362 adcroft 1.1 alpha = alpha+cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)
363     ENDDO
364     ENDDO
365     ENDDO
366     ELSE
367     DO K=1,1
368     DO J=1,sNy
369     DO I=1,sNx
370     cg3d_q(I,J,K,bi,bj) =
371     & aW3d(I ,J ,K ,bi,bj)*cg3d_s(I-1,J ,K ,bi,bj)
372     & +aW3d(I+1,J ,K ,bi,bj)*cg3d_s(I+1,J ,K ,bi,bj)
373     & +aS3d(I ,J ,K ,bi,bj)*cg3d_s(I ,J-1,K ,bi,bj)
374     & +aS3d(I ,J+1,K ,bi,bj)*cg3d_s(I ,J+1,K ,bi,bj)
375     & -aW3d(I ,J ,K ,bi,bj)*cg3d_s(I ,J ,K ,bi,bj)
376     & -aW3d(I+1,J ,K ,bi,bj)*cg3d_s(I ,J ,K ,bi,bj)
377     & -aS3d(I ,J ,K ,bi,bj)*cg3d_s(I ,J ,K ,bi,bj)
378     & -aS3d(I ,J+1,K ,bi,bj)*cg3d_s(I ,J ,K ,bi,bj)
379 jmc 1.7 & -topLevTerm*_rA(I,J,bi,bj)*cg3d_s(I,J,K,bi,bj)
380 adcroft 1.1 alpha = alpha+cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)
381     ENDDO
382     ENDDO
383     ENDDO
384     ENDIF
385     DO K=2,Nr-1
386     DO J=1,sNy
387     DO I=1,sNx
388     cg3d_q(I,J,K,bi,bj) =
389     & aW3d(I ,J ,K ,bi,bj)*cg3d_s(I-1,J ,K ,bi,bj)
390     & +aW3d(I+1,J ,K ,bi,bj)*cg3d_s(I+1,J ,K ,bi,bj)
391     & +aS3d(I ,J ,K ,bi,bj)*cg3d_s(I ,J-1,K ,bi,bj)
392     & +aS3d(I ,J+1,K ,bi,bj)*cg3d_s(I ,J+1,K ,bi,bj)
393     & +aV3d(I ,J ,K ,bi,bj)*cg3d_s(I ,J ,K-1,bi,bj)
394     & +aV3d(I ,J ,K+1,bi,bj)*cg3d_s(I ,J ,K+1,bi,bj)
395     & -aW3d(I ,J ,K ,bi,bj)*cg3d_s(I ,J ,K ,bi,bj)
396     & -aW3d(I+1,J ,K ,bi,bj)*cg3d_s(I ,J ,K ,bi,bj)
397     & -aS3d(I ,J ,K ,bi,bj)*cg3d_s(I ,J ,K ,bi,bj)
398     & -aS3d(I ,J+1,K ,bi,bj)*cg3d_s(I ,J ,K ,bi,bj)
399     & -aV3d(I ,J ,K ,bi,bj)*cg3d_s(I ,J ,K ,bi,bj)
400     & -aV3d(I ,J ,K+1,bi,bj)*cg3d_s(I ,J ,K ,bi,bj)
401     alpha = alpha+cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)
402     ENDDO
403     ENDDO
404     ENDDO
405     IF ( Nr .GT. 1 ) THEN
406     DO K=Nr,Nr
407     DO J=1,sNy
408     DO I=1,sNx
409     cg3d_q(I,J,K,bi,bj) =
410     & aW3d(I ,J ,K ,bi,bj)*cg3d_s(I-1,J ,K ,bi,bj)
411     & +aW3d(I+1,J ,K ,bi,bj)*cg3d_s(I+1,J ,K ,bi,bj)
412     & +aS3d(I ,J ,K ,bi,bj)*cg3d_s(I ,J-1,K ,bi,bj)
413     & +aS3d(I ,J+1,K ,bi,bj)*cg3d_s(I ,J+1,K ,bi,bj)
414     & +aV3d(I ,J ,K ,bi,bj)*cg3d_s(I ,J ,K-1,bi,bj)
415     & -aW3d(I ,J ,K ,bi,bj)*cg3d_s(I ,J ,K ,bi,bj)
416     & -aW3d(I+1,J ,K ,bi,bj)*cg3d_s(I ,J ,K ,bi,bj)
417     & -aS3d(I ,J ,K ,bi,bj)*cg3d_s(I ,J ,K ,bi,bj)
418     & -aS3d(I ,J+1,K ,bi,bj)*cg3d_s(I ,J ,K ,bi,bj)
419     & -aV3d(I ,J ,K ,bi,bj)*cg3d_s(I ,J ,K ,bi,bj)
420     alpha = alpha+cg3d_s(I,J,K,bi,bj)*cg3d_q(I,J,K,bi,bj)
421     ENDDO
422     ENDDO
423     ENDDO
424     ENDIF
425     ENDDO
426     ENDDO
427 adcroft 1.2 _GLOBAL_SUM_R8(alpha,myThid)
428 adcroft 1.1 CcnhDebugStarts
429 heimbach 1.8 C WRITE(*,*) ' CG3D: Iteration ',it3d-1,' SUM(s*q)= ',alpha
430 adcroft 1.1 CcnhDebugEnds
431 jmc 1.6 alpha = eta_qrN/alpha
432 adcroft 1.1 CcnhDebugStarts
433 heimbach 1.8 C WRITE(*,*) ' CG3D: Iteration ',it3d-1,' alpha= ',alpha
434 adcroft 1.1 CcnhDebugEnds
435    
436     C== Update solution and residual vectors
437     C Now compute "interior" points.
438     err = 0. _d 0
439     DO bj=myByLo(myThid),myByHi(myThid)
440     DO bi=myBxLo(myThid),myBxHi(myThid)
441     DO K=1,Nr
442     DO J=1,sNy
443     DO I=1,sNx
444     cg3d_x(I,J,K,bi,bj)=cg3d_x(I,J,K,bi,bj)
445     & +alpha*cg3d_s(I,J,K,bi,bj)
446     cg3d_r(I,J,K,bi,bj)=cg3d_r(I,J,K,bi,bj)
447     & -alpha*cg3d_q(I,J,K,bi,bj)
448     err = err+cg3d_r(I,J,K,bi,bj)*cg3d_r(I,J,K,bi,bj)
449     ENDDO
450     ENDDO
451     ENDDO
452     ENDDO
453     ENDDO
454    
455 adcroft 1.2 _GLOBAL_SUM_R8( err , myThid )
456 adcroft 1.1 err = SQRT(err)
457     actualIts = it3d
458     actualResidual = err
459     IF ( actualResidual .LT. cg3dTargetResidual ) GOTO 11
460     C _EXCH_XYZ_R8(cg3d_r, myThid )
461     OLw = 1
462     OLe = 1
463     OLn = 1
464     OLs = 1
465     exchWidthX = 1
466     exchWidthY = 1
467     myNz = Nr
468     CALL EXCH_RL( cg3d_r,
469 adcroft 1.10 I OLw, OLe, OLs, OLn, myNz,
470 adcroft 1.1 I exchWidthX, exchWidthY,
471     I FORWARD_SIMULATION, EXCH_IGNORE_CORNERS, myThid )
472    
473     10 CONTINUE
474     11 CONTINUE
475    
476     C-- Un-normalise the answer
477     DO bj=myByLo(myThid),myByHi(myThid)
478     DO bi=myBxLo(myThid),myBxHi(myThid)
479     DO K=1,Nr
480     DO J=1,sNy
481     DO I=1,sNx
482     cg3d_x(I,J,K,bi,bj) = cg3d_x(I,J,K,bi,bj)/rhsNorm
483     ENDDO
484     ENDDO
485     ENDDO
486     ENDDO
487     ENDDO
488    
489 adcroft 1.3 Cadj _EXCH_XYZ_R8(cg3d_x, myThid )
490 adcroft 1.11 c _BEGIN_MASTER( myThid )
491     c WRITE(*,'(A,I6,1PE30.14)') ' CG3D iters, err = ',
492     c & actualIts, actualResidual
493     c _END_MASTER( )
494     lastResidual=actualResidual
495     numIters=actualIts
496 adcroft 1.1
497     #endif /* ALLOW_NONHYDROSTATIC */
498    
499     RETURN
500     END
  ViewVC Help
Powered by ViewVC 1.1.22