/[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.13 - (hide annotations) (download)
Thu Oct 9 04:19:18 2003 UTC (20 years, 7 months ago) by edhill
Branch: MAIN
CVS Tags: checkpoint51k_post, checkpoint52l_pre, hrcube4, hrcube5, checkpoint52d_pre, checkpoint56b_post, checkpoint52j_pre, checkpoint51o_pre, checkpoint54d_post, checkpoint54e_post, checkpoint51l_post, checkpoint52l_post, checkpoint52k_post, checkpoint55, checkpoint54, checkpoint57, checkpoint56, checkpoint53, checkpoint52, checkpoint52f_post, checkpoint54f_post, checkpoint51t_post, checkpoint51n_post, checkpoint55i_post, checkpoint52i_pre, hrcube_1, hrcube_2, hrcube_3, checkpoint51s_post, checkpoint55c_post, checkpoint52e_pre, checkpoint52e_post, checkpoint51n_pre, checkpoint53d_post, checkpoint57a_post, checkpoint52b_pre, checkpoint54b_post, checkpoint51l_pre, checkpoint52m_post, checkpoint55g_post, checkpoint51q_post, checkpoint52b_post, checkpoint52c_post, checkpoint52f_pre, checkpoint55d_post, checkpoint54a_pre, checkpoint53c_post, checkpoint55d_pre, checkpoint55j_post, checkpoint54a_post, checkpoint55h_post, checkpoint51r_post, checkpoint51i_post, checkpoint55b_post, checkpoint53a_post, checkpoint55f_post, checkpoint52d_post, checkpoint53g_post, checkpoint52a_pre, checkpoint52i_post, checkpoint52h_pre, checkpoint56a_post, checkpoint53f_post, checkpoint52j_post, branch-netcdf, checkpoint52n_post, checkpoint53b_pre, checkpoint56c_post, checkpoint57a_pre, checkpoint55a_post, checkpoint51o_post, checkpoint53b_post, checkpoint52a_post, ecco_c52_e35, checkpoint51m_post, checkpoint53d_pre, checkpoint55e_post, checkpoint54c_post, checkpoint51p_post, checkpoint51u_post
Branch point for: branch-nonh, tg2-branch, netcdf-sm0, checkpoint51n_branch
Changes since 1.12: +5 -1 lines 
 o first check-in for the "branch-genmake2" merge
 o verification suite as run on shelley (gcc 3.2.2):

Wed Oct  8 23:42:29 EDT 2003
                T           S           U           V
G D M    c        m  s        m  s        m  s        m  s
E p a R  g  m  m  e  .  m  m  e  .  m  m  e  .  m  m  e  .
N n k u  2  i  a  a  d  i  a  a  d  i  a  a  d  i  a  a  d
2 d e n  d  n  x  n  .  n  x  n  .  n  x  n  .  n  x  n  .

OPTFILE=NONE

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