/[MITgcm]/MITgcm_contrib/dgoldberg/streamice/streamice_cg_solve.F
ViewVC logotype

Contents of /MITgcm_contrib/dgoldberg/streamice/streamice_cg_solve.F

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


Revision 1.3 - (show annotations) (download)
Mon May 14 16:52:29 2012 UTC (13 years, 2 months ago) by dgoldberg
Branch: MAIN
Changes since 1.2: +16 -5 lines
make construct matrix option a CPP option

1 C $Header: /u/gcmpack/MITgcm_contrib/dgoldberg/streamice/streamice_cg_solve.F,v 1.2 2012/05/02 02:36:01 heimbach Exp $
2 C $Name: $
3
4 #include "STREAMICE_OPTIONS.h"
5
6 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
7
8 CBOP
9 SUBROUTINE STREAMICE_CG_SOLVE(
10 U cg_Uin,
11 U cg_Vin,
12 I cg_Bu,
13 I cg_Bv,
14 I tolerance,
15 O iters,
16 I myThid )
17 C /============================================================\
18 C | SUBROUTINE |
19 C | o |
20 C |============================================================|
21 C | |
22 C \============================================================/
23 IMPLICIT NONE
24
25 C === Global variables ===
26 #include "SIZE.h"
27 #include "EEPARAMS.h"
28 #include "PARAMS.h"
29 #include "STREAMICE.h"
30 #include "STREAMICE_CG.h"
31
32
33 C !INPUT/OUTPUT ARGUMENTS
34 C cg_Uin, cg_Vin - input and output velocities
35 C cg_Bu, cg_Bv - driving stress
36 INTEGER myThid
37 INTEGER iters
38 _RL tolerance
39 _RL cg_Uin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
40 _RL cg_Vin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
41 _RL cg_Bu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
42 _RL cg_Bv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
43
44 C LOCAL VARIABLES
45 INTEGER i, j, bi, bj, cg_halo, conv_flag
46 INTEGER iter, is, js, ie, je, colx, coly, k
47 _RL dot_p1, dot_p2, resid_0, alpha_k, beta_k, resid
48 _RL dot_p1_tile (nSx,nSy)
49 _RL dot_p2_tile (nSx,nSy)
50
51 iters = streamice_max_cg_iter
52
53 #ifdef ALLOW_STREAMICE
54
55 conv_flag = 0
56
57 DO bj = myByLo(myThid), myByHi(myThid)
58 DO bi = myBxLo(myThid), myBxHi(myThid)
59 DO j=1,sNy
60 DO i=1,sNx
61 Zu_SI (i,j,bi,bj) = 0. _d 0
62 Zv_SI (i,j,bi,bj) = 0. _d 0
63 DIAGu_SI (i,j,bi,bj) = 0. _d 0
64 DIAGv_SI (i,j,bi,bj) = 0. _d 0
65 Ru_SI (i,j,bi,bj) = 0. _d 0
66 Rv_SI (i,j,bi,bj) = 0. _d 0
67 Au_SI (i,j,bi,bj) = 0. _d 0
68 Av_SI (i,j,bi,bj) = 0. _d 0
69 Du_SI (i,j,bi,bj) = 0. _d 0
70 Dv_SI (i,j,bi,bj) = 0. _d 0
71 ubd_SI (i,j,bi,bj) = 0. _d 0
72 vbd_SI (i,j,bi,bj) = 0. _d 0
73 ENDDO
74 ENDDO
75 ENDDO
76 ENDDO
77
78 C PREAMBLE: get bdry contribution, find matrix diagonal,
79 C initialize iterates R (residual), Z, and D
80
81 CALL STREAMICE_CG_BOUND_VALS( myThid,
82 O ubd_SI,
83 O vbd_SI)
84
85 DO bj = myByLo(myThid), myByHi(myThid)
86 DO bi = myBxLo(myThid), myBxHi(myThid)
87 DO j=1-OLy,sNy+OLy
88 DO i=1-OLx,sNx+OLx
89 RHSu_SI (i,j,bi,bj) = cg_Bu(i,j,bi,bj)
90 & - ubd_SI(i,j,bi,bj)
91 RHSv_SI (i,j,bi,bj) = cg_Bv(i,j,bi,bj)
92 & - vbd_SI(i,j,bi,bj)
93 ENDDO
94 ENDDO
95 ENDDO
96 ENDDO
97
98 _EXCH_XY_RL( RHSu_SI, myThid )
99 _EXCH_XY_RL( RHSv_SI, myThid )
100
101 CALL STREAMICE_CG_ADIAG( myThid,
102 O DIAGu_SI,
103 O DIAGv_SI)
104
105 _EXCH_XY_RL( DIAGu_SI, myThid )
106 _EXCH_XY_RL( DIAGv_SI, myThid )
107
108 C FIX PROBLEM WITH PRECOND LATER
109
110 ! DO bj = myByLo(myThid), myByHi(myThid)
111 ! DO bi = myBxLo(myThid), myBxHi(myThid)
112 ! DO j=1-OLy,sNy+OLy
113 ! DO i=1-OLx,sNx+OLx
114 ! DIAGu_SI(i,j,bi,bj)=1.0
115 ! DIAGv_SI(i,j,bi,bj)=1.0
116 ! ENDDO
117 ! ENDDO
118 ! ENDDO
119 ! ENDDO
120
121 CALL STREAMICE_CG_ACTION( myThid,
122 O Au_SI,
123 O Av_SI,
124 I cg_Uin,
125 I cg_Vin,
126 I 0, sNx+1, 0, sNy+1 )
127
128 _EXCH_XY_RL( Au_SI, myThid )
129 _EXCH_XY_RL( Av_SI, myThid )
130
131 DO bj = myByLo(myThid), myByHi(myThid)
132 DO bi = myBxLo(myThid), myBxHi(myThid)
133 DO j=1-OLy,sNy+OLy
134 DO i=1-OLx,sNx+OLx
135 Ru_SI(i,j,bi,bj)=RHSu_SI(i,j,bi,bj)-
136 & Au_SI(i,j,bi,bj)
137 Rv_SI(i,j,bi,bj)=RHSv_SI(i,j,bi,bj)-
138 & Av_SI(i,j,bi,bj)
139 ENDDO
140 ENDDO
141 dot_p1_tile(bi,bj) = 0. _d 0
142 dot_p2_tile(bi,bj) = 0. _d 0
143 ENDDO
144 ENDDO
145
146 DO bj = myByLo(myThid), myByHi(myThid)
147 DO bi = myBxLo(myThid), myBxHi(myThid)
148 DO j=1,sNy
149 DO i=1,sNx
150 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
151 & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Ru_SI(i,j,bi,bj)**2
152 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
153 & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Rv_SI(i,j,bi,bj)**2
154 ENDDO
155 ENDDO
156 ENDDO
157 ENDDO
158
159 CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
160 resid_0 = sqrt(dot_p1)
161
162 DO bj = myByLo(myThid), myByHi(myThid)
163 DO bi = myBxLo(myThid), myBxHi(myThid)
164 DO j=1-OLy,sNy+OLy
165 DO i=1-OLx,sNx+OLx
166 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
167 & Zu_SI(i,j,bi,bj)=Ru_SI(i,j,bi,bj) / DIAGu_SI(i,j,bi,bj)
168 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
169 & Zv_SI(i,j,bi,bj)=Rv_SI(i,j,bi,bj) / DIAGv_SI(i,j,bi,bj)
170 ENDDO
171 ENDDO
172 ENDDO
173 ENDDO
174
175
176 cg_halo = min(OLx-1,OLy-1)
177 conv_flag = 0
178
179 DO bj = myByLo(myThid), myByHi(myThid)
180 DO bi = myBxLo(myThid), myBxHi(myThid)
181 DO j=1-OLy,sNy+OLy
182 DO i=1-OLx,sNx+OLx
183 Du_SI(i,j,bi,bj)=Zu_SI(i,j,bi,bj)
184 Dv_SI(i,j,bi,bj)=Zv_SI(i,j,bi,bj)
185 ENDDO
186 ENDDO
187 ENDDO
188 ENDDO
189
190 resid = resid_0
191 iters = 0
192
193 c !!!!!!!!!!!!!!!!!!
194 c !! !!
195 c !! MAIN CG LOOP !!
196 c !! !!
197 c !!!!!!!!!!!!!!!!!!
198
199 c ! initially, b-grid data is valid up to 3 halo nodes out -- right? (check for MITgcm!!)
200
201 print *, "BEGINNING MAIN CG LOOP"
202
203 ! IF(STREAMICE_construct_matrix) CALL STREAMICE_CG_MAKE_A(myThid)
204 #ifdef STREAMICE_CONSTRUCT_MATRIX
205
206 CALL STREAMICE_CG_MAKE_A(myThid)
207
208 #endif
209
210 do iter = 1, streamice_max_cg_iter
211 if (resid .gt. tolerance*resid_0) then
212
213 c to avoid using "exit"
214 iters = iters + 1
215
216 is = 1 - cg_halo
217 ie = sNx + cg_halo
218 js = 1 - cg_halo
219 je = sNy + cg_halo
220
221 DO bj = myByLo(myThid), myByHi(myThid)
222 DO bi = myBxLo(myThid), myBxHi(myThid)
223 DO j=1-OLy,sNy+OLy
224 DO i=1-OLx,sNx+OLx
225 Au_SI(i,j,bi,bj) = 0. _d 0
226 Av_SI(i,j,bi,bj) = 0. _d 0
227 ! Du_SI(i,j,bi,bj) = Real(i)
228 ! Dv_SI(i,j,bi,bj) = 0.0
229 ENDDO
230 ENDDO
231 ENDDO
232 ENDDO
233
234 ! IF (STREAMICE_construct_matrix) THEN
235
236 #ifdef STREAMICE_CONSTRUCT_MATRIX
237
238 DO bj = myByLo(myThid), myByHi(myThid)
239 DO bi = myBxLo(myThid), myBxHi(myThid)
240 DO j=js,je
241 DO i=is,ie
242 DO colx=-1,1
243 DO coly=-1,1
244 Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) +
245 & streamice_cg_A1(i,j,bi,bj,colx,coly)*
246 & Du_SI(i+colx,j+coly,bi,bj)+
247 & streamice_cg_A2(i,j,bi,bj,colx,coly)*
248 & Dv_SI(i+colx,j+coly,bi,bj)
249 Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +
250 & streamice_cg_A3(i,j,bi,bj,colx,coly)*
251 & Du_SI(i+colx,j+coly,bi,bj)+
252 & streamice_cg_A4(i,j,bi,bj,colx,coly)*
253 & Dv_SI(i+colx,j+coly,bi,bj)
254 ENDDO
255 ENDDO
256 ENDDO
257 ENDDO
258 ENDDO
259 ENDDO
260
261 ! else
262 #else
263
264 CALL STREAMICE_CG_ACTION( myThid,
265 O Au_SI,
266 O Av_SI,
267 I Du_SI,
268 I Dv_SI,
269 I is,ie,js,je)
270
271 ! ENDIF
272
273 #endif
274
275 ! ! if (iter.eq.1) then
276 ! ! CALL WRITE_FLD_XY_RL ("Au2","",Au_SI,0,myThid)
277 ! ! CALL WRITE_FLD_XY_RL ("Av2","",Av_SI,0,myThid)
278 ! ! CALL WRITE_FLD_XY_RL ("Du","",Du_SI,0,myThid)
279 ! !
280 ! ! CALL WRITE_FLD_XY_RL("Dv","",Au_SI,0,myThid)
281 ! !
282 ! ! CALL WRITE_FLD_XY_RL("DiagU1","",Diagu_SI,0,myThid)
283 ! ! CALL WRITE_FLD_XY_RL("DiagV1","",Diagv_SI,0,myThid)
284 ! ! CALL WRITE_FLD_XY_RL("DiagU2","",
285 ! ! & streamice_cg_A1(i,j,bi,bj,0,0),0,myThid)
286 ! ! CALL WRITE_FLD_XY_RL("DiagV2","",
287 ! ! & streamice_cg_A4(i,j,bi,bj,0,0),0,myThid)
288 ! ! CALL WRITE_FLD_XY_RL("DiagV2","",Diagv_SI,0,myThid)
289 ! !
290 ! ! endif
291
292
293
294
295
296 ! ! if (iter.eq.1) then
297 ! ! CALL WRITE_FLD_XY_RL ("Au1","",Au_SI,0,myThid)
298 ! ! CALL WRITE_FLD_XY_RL ("Av1","",Av_SI,0,myThid)
299 ! !
300 ! ! endif
301
302
303
304 DO bj = myByLo(myThid), myByHi(myThid)
305 DO bi = myBxLo(myThid), myBxHi(myThid)
306 dot_p1_tile(bi,bj) = 0. _d 0
307 dot_p2_tile(bi,bj) = 0. _d 0
308 ENDDO
309 ENDDO
310
311 DO bj = myByLo(myThid), myByHi(myThid)
312 DO bi = myBxLo(myThid), myBxHi(myThid)
313 DO j=1,sNy
314 DO i=1,sNx
315 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN
316 dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zu_SI(i,j,bi,bj)*
317 & Ru_SI(i,j,bi,bj)
318 dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Du_SI(i,j,bi,bj)*
319 & Au_SI(i,j,bi,bj)
320 ENDIF
321 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN
322 dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zv_SI(i,j,bi,bj)*
323 & Rv_SI(i,j,bi,bj)
324 dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Dv_SI(i,j,bi,bj)*
325 & Av_SI(i,j,bi,bj)
326 ENDIF
327 ENDDO
328 ENDDO
329 ENDDO
330 ENDDO
331
332 CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
333 CALL GLOBAL_SUM_TILE_RL( dot_p2_tile, dot_p2, myThid )
334 alpha_k = dot_p1/dot_p2
335
336 DO bj = myByLo(myThid), myByHi(myThid)
337 DO bi = myBxLo(myThid), myBxHi(myThid)
338 DO j=1-OLy,sNy+OLy
339 DO i=1-OLx,sNx+OLx
340
341 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN
342 cg_Uin(i,j,bi,bj)=cg_Uin(i,j,bi,bj)+
343 & alpha_k*Du_SI(i,j,bi,bj)
344 Ru_old_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj)
345 Zu_old_SI(i,j,bi,bj) = Zu_SI(i,j,bi,bj)
346 Ru_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj)-
347 & alpha_k*Au_SI(i,j,bi,bj)
348 Zu_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj) /
349 & DIAGu_SI(i,j,bi,bj)
350 ENDIF
351
352 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN
353 cg_Vin(i,j,bi,bj)=cg_Vin(i,j,bi,bj)+
354 & alpha_k*Dv_SI(i,j,bi,bj)
355 Rv_old_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj)
356 Zv_old_SI(i,j,bi,bj) = Zv_SI(i,j,bi,bj)
357 Rv_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj)-
358 & alpha_k*Av_SI(i,j,bi,bj)
359 Zv_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj) /
360 & DIAGv_SI(i,j,bi,bj)
361
362 ENDIF
363 ENDDO
364 ENDDO
365 ENDDO
366 ENDDO
367
368 DO bj = myByLo(myThid), myByHi(myThid)
369 DO bi = myBxLo(myThid), myBxHi(myThid)
370 dot_p1_tile(bi,bj) = 0. _d 0
371 dot_p2_tile(bi,bj) = 0. _d 0
372 ENDDO
373 ENDDO
374
375 DO bj = myByLo(myThid), myByHi(myThid)
376 DO bi = myBxLo(myThid), myBxHi(myThid)
377 DO j=1,sNy
378 DO i=1,sNx
379
380 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN
381 dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zu_SI(i,j,bi,bj)*
382 & Ru_SI(i,j,bi,bj)
383 dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Zu_old_SI(i,j,bi,bj)*
384 & Ru_old_SI(i,j,bi,bj)
385 ENDIF
386
387 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN
388 dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zv_SI(i,j,bi,bj)*
389 & Rv_SI(i,j,bi,bj)
390 dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Zv_old_SI(i,j,bi,bj)*
391 & Rv_old_SI(i,j,bi,bj)
392 ENDIF
393
394 ENDDO
395 ENDDO
396 ENDDO
397 ENDDO
398
399 CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
400 CALL GLOBAL_SUM_TILE_RL( dot_p2_tile, dot_p2, myThid )
401
402 beta_k = dot_p1/dot_p2
403
404 DO bj = myByLo(myThid), myByHi(myThid)
405 DO bi = myBxLo(myThid), myBxHi(myThid)
406 DO j=1-OLy,sNy+OLy
407 DO i=1-OLx,sNx+OLx
408 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
409 & Du_SI(i,j,bi,bj)=beta_k*Du_SI(i,j,bi,bj)+
410 & Zu_SI(i,j,bi,bj)
411 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
412 & Dv_SI(i,j,bi,bj)=beta_k*Dv_SI(i,j,bi,bj)+
413 & Zv_SI(i,j,bi,bj)
414 ENDDO
415 ENDDO
416 ENDDO
417 ENDDO
418
419 DO bj = myByLo(myThid), myByHi(myThid)
420 DO bi = myBxLo(myThid), myBxHi(myThid)
421 dot_p1_tile(bi,bj) = 0. _d 0
422 ENDDO
423 ENDDO
424
425 DO bj = myByLo(myThid), myByHi(myThid)
426 DO bi = myBxLo(myThid), myBxHi(myThid)
427 DO j=1,sNy
428 DO i=1,sNx
429 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
430 & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Ru_SI(i,j,bi,bj)**2
431 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
432 & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Rv_SI(i,j,bi,bj)**2
433 ENDDO
434 ENDDO
435 ENDDO
436 ENDDO
437
438 CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
439 resid = sqrt(dot_p1)
440
441 ! IF (iter .eq. 1) then
442 ! print *, alpha_k, beta_k, resid
443 ! ENDIF
444
445 cg_halo = cg_halo - 1
446
447 if (cg_halo .eq. 0) then
448 cg_halo = min(OLx-1,OLy-1)
449 _EXCH_XY_RL( Du_SI, myThid )
450 _EXCH_XY_RL( Dv_SI, myThid )
451 _EXCH_XY_RL( Ru_SI, myThid )
452 _EXCH_XY_RL( Rv_SI, myThid )
453 _EXCH_XY_RL( cg_Uin, myThid )
454 _EXCH_XY_RL( cg_Vin, myThid )
455 endif
456
457 endif
458 enddo ! end of CG loop
459
460 c to avoid using "exit"
461 c if iters has reached max_iters there is no convergence
462
463 IF (iters .lt. streamice_max_cg_iter) THEN
464 conv_flag = 1
465 ENDIF
466
467 DO bj = myByLo(myThid), myByHi(myThid)
468 DO bi = myBxLo(myThid), myBxHi(myThid)
469 DO j=1-OLy,sNy+OLy
470 DO i=1-OLy,sNx+OLy
471 IF (STREAMICE_umask(i,j,bi,bj).eq.3.0)
472 & cg_Uin(i,j,bi,bj)=u_bdry_values_SI(i,j,bi,bj)
473 IF (STREAMICE_vmask(i,j,bi,bj).eq.3.0)
474 & cg_Vin(i,j,bi,bj)=v_bdry_values_SI(i,j,bi,bj)
475 ENDDO
476 ENDDO
477 ENDDO
478 ENDDO
479
480 _EXCH_XY_RL( cg_Uin, myThid )
481 _EXCH_XY_RL( cg_Vin, myThid )
482
483 #endif
484 RETURN
485 END
486
487
488
489
490
491
492
493

  ViewVC Help
Powered by ViewVC 1.1.22