/[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.1 - (show annotations) (download)
Thu Mar 29 15:59:21 2012 UTC (13 years, 3 months ago) by heimbach
Branch: MAIN
Initial check-in of Dan Goldberg's streamice package

1 C $Header: /u/gcmpack/MITgcm/pkg/streamice/streamice_init_varia.F,v 1.6 2011/06/29 16:24:10 dng 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 myThid,
15 I tolerance,
16 O iters )
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
205 do iter = 1, streamice_max_cg_iter
206 if (resid .gt. tolerance*resid_0) then
207
208 c to avoid using "exit"
209 iters = iters + 1
210
211 is = 1 - cg_halo
212 ie = sNx + cg_halo
213 js = 1 - cg_halo
214 je = sNy + cg_halo
215
216 DO bj = myByLo(myThid), myByHi(myThid)
217 DO bi = myBxLo(myThid), myBxHi(myThid)
218 DO j=1-OLy,sNy+OLy
219 DO i=1-OLx,sNx+OLx
220 Au_SI(i,j,bi,bj) = 0. _d 0
221 Av_SI(i,j,bi,bj) = 0. _d 0
222 ! Du_SI(i,j,bi,bj) = Real(i)
223 ! Dv_SI(i,j,bi,bj) = 0.0
224 ENDDO
225 ENDDO
226 ENDDO
227 ENDDO
228
229 IF (STREAMICE_construct_matrix) THEN
230 DO bj = myByLo(myThid), myByHi(myThid)
231 DO bi = myBxLo(myThid), myBxHi(myThid)
232 DO j=js,je
233 DO i=is,ie
234 DO colx=-1,1
235 DO coly=-1,1
236 Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) +
237 & streamice_cg_A1(i,j,bi,bj,colx,coly)*
238 & Du_SI(i+colx,j+coly,bi,bj)+
239 & streamice_cg_A2(i,j,bi,bj,colx,coly)*
240 & Dv_SI(i+colx,j+coly,bi,bj)
241 Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +
242 & streamice_cg_A3(i,j,bi,bj,colx,coly)*
243 & Du_SI(i+colx,j+coly,bi,bj)+
244 & streamice_cg_A4(i,j,bi,bj,colx,coly)*
245 & Dv_SI(i+colx,j+coly,bi,bj)
246 ENDDO
247 ENDDO
248 ENDDO
249 ENDDO
250 ENDDO
251 ENDDO
252
253 else
254
255 CALL STREAMICE_CG_ACTION( myThid,
256 O Au_SI,
257 O Av_SI,
258 I Du_SI,
259 I Dv_SI,
260 I is,ie,js,je)
261
262 ENDIF
263
264 ! ! if (iter.eq.1) then
265 ! ! CALL WRITE_FLD_XY_RL ("Au2","",Au_SI,0,myThid)
266 ! ! CALL WRITE_FLD_XY_RL ("Av2","",Av_SI,0,myThid)
267 ! ! CALL WRITE_FLD_XY_RL ("Du","",Du_SI,0,myThid)
268 ! !
269 ! ! CALL WRITE_FLD_XY_RL("Dv","",Au_SI,0,myThid)
270 ! !
271 ! ! CALL WRITE_FLD_XY_RL("DiagU1","",Diagu_SI,0,myThid)
272 ! ! CALL WRITE_FLD_XY_RL("DiagV1","",Diagv_SI,0,myThid)
273 ! ! CALL WRITE_FLD_XY_RL("DiagU2","",
274 ! ! & streamice_cg_A1(i,j,bi,bj,0,0),0,myThid)
275 ! ! CALL WRITE_FLD_XY_RL("DiagV2","",
276 ! ! & streamice_cg_A4(i,j,bi,bj,0,0),0,myThid)
277 ! ! CALL WRITE_FLD_XY_RL("DiagV2","",Diagv_SI,0,myThid)
278 ! !
279 ! ! endif
280
281
282
283
284
285 ! ! if (iter.eq.1) then
286 ! ! CALL WRITE_FLD_XY_RL ("Au1","",Au_SI,0,myThid)
287 ! ! CALL WRITE_FLD_XY_RL ("Av1","",Av_SI,0,myThid)
288 ! !
289 ! ! endif
290
291
292
293 DO bj = myByLo(myThid), myByHi(myThid)
294 DO bi = myBxLo(myThid), myBxHi(myThid)
295 dot_p1_tile(bi,bj) = 0. _d 0
296 dot_p2_tile(bi,bj) = 0. _d 0
297 ENDDO
298 ENDDO
299
300 DO bj = myByLo(myThid), myByHi(myThid)
301 DO bi = myBxLo(myThid), myBxHi(myThid)
302 DO j=1,sNy
303 DO i=1,sNx
304 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN
305 dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zu_SI(i,j,bi,bj)*
306 & Ru_SI(i,j,bi,bj)
307 dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Du_SI(i,j,bi,bj)*
308 & Au_SI(i,j,bi,bj)
309 ENDIF
310 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN
311 dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zv_SI(i,j,bi,bj)*
312 & Rv_SI(i,j,bi,bj)
313 dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Dv_SI(i,j,bi,bj)*
314 & Av_SI(i,j,bi,bj)
315 ENDIF
316 ENDDO
317 ENDDO
318 ENDDO
319 ENDDO
320
321 CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
322 CALL GLOBAL_SUM_TILE_RL( dot_p2_tile, dot_p2, myThid )
323 alpha_k = dot_p1/dot_p2
324
325 DO bj = myByLo(myThid), myByHi(myThid)
326 DO bi = myBxLo(myThid), myBxHi(myThid)
327 DO j=1-OLy,sNy+OLy
328 DO i=1-OLx,sNx+OLx
329
330 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN
331 cg_Uin(i,j,bi,bj)=cg_Uin(i,j,bi,bj)+
332 & alpha_k*Du_SI(i,j,bi,bj)
333 Ru_old_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj)
334 Zu_old_SI(i,j,bi,bj) = Zu_SI(i,j,bi,bj)
335 Ru_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj)-
336 & alpha_k*Au_SI(i,j,bi,bj)
337 Zu_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj) /
338 & DIAGu_SI(i,j,bi,bj)
339 ENDIF
340
341 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN
342 cg_Vin(i,j,bi,bj)=cg_Vin(i,j,bi,bj)+
343 & alpha_k*Dv_SI(i,j,bi,bj)
344 Rv_old_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj)
345 Zv_old_SI(i,j,bi,bj) = Zv_SI(i,j,bi,bj)
346 Rv_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj)-
347 & alpha_k*Av_SI(i,j,bi,bj)
348 Zv_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj) /
349 & DIAGv_SI(i,j,bi,bj)
350
351 ENDIF
352 ENDDO
353 ENDDO
354 ENDDO
355 ENDDO
356
357 DO bj = myByLo(myThid), myByHi(myThid)
358 DO bi = myBxLo(myThid), myBxHi(myThid)
359 dot_p1_tile(bi,bj) = 0. _d 0
360 dot_p2_tile(bi,bj) = 0. _d 0
361 ENDDO
362 ENDDO
363
364 DO bj = myByLo(myThid), myByHi(myThid)
365 DO bi = myBxLo(myThid), myBxHi(myThid)
366 DO j=1,sNy
367 DO i=1,sNx
368
369 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN
370 dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zu_SI(i,j,bi,bj)*
371 & Ru_SI(i,j,bi,bj)
372 dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Zu_old_SI(i,j,bi,bj)*
373 & Ru_old_SI(i,j,bi,bj)
374 ENDIF
375
376 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN
377 dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zv_SI(i,j,bi,bj)*
378 & Rv_SI(i,j,bi,bj)
379 dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Zv_old_SI(i,j,bi,bj)*
380 & Rv_old_SI(i,j,bi,bj)
381 ENDIF
382
383 ENDDO
384 ENDDO
385 ENDDO
386 ENDDO
387
388 CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
389 CALL GLOBAL_SUM_TILE_RL( dot_p2_tile, dot_p2, myThid )
390
391 beta_k = dot_p1/dot_p2
392
393 DO bj = myByLo(myThid), myByHi(myThid)
394 DO bi = myBxLo(myThid), myBxHi(myThid)
395 DO j=1-OLy,sNy+OLy
396 DO i=1-OLx,sNx+OLx
397 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
398 & Du_SI(i,j,bi,bj)=beta_k*Du_SI(i,j,bi,bj)+
399 & Zu_SI(i,j,bi,bj)
400 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
401 & Dv_SI(i,j,bi,bj)=beta_k*Dv_SI(i,j,bi,bj)+
402 & Zv_SI(i,j,bi,bj)
403 ENDDO
404 ENDDO
405 ENDDO
406 ENDDO
407
408 DO bj = myByLo(myThid), myByHi(myThid)
409 DO bi = myBxLo(myThid), myBxHi(myThid)
410 dot_p1_tile(bi,bj) = 0. _d 0
411 ENDDO
412 ENDDO
413
414 DO bj = myByLo(myThid), myByHi(myThid)
415 DO bi = myBxLo(myThid), myBxHi(myThid)
416 DO j=1,sNy
417 DO i=1,sNx
418 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
419 & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Ru_SI(i,j,bi,bj)**2
420 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
421 & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Rv_SI(i,j,bi,bj)**2
422 ENDDO
423 ENDDO
424 ENDDO
425 ENDDO
426
427 CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
428 resid = sqrt(dot_p1)
429
430 ! IF (iter .eq. 1) then
431 ! print *, alpha_k, beta_k, resid
432 ! ENDIF
433
434 cg_halo = cg_halo - 1
435
436 if (cg_halo .eq. 0) then
437 cg_halo = min(OLx-1,OLy-1)
438 _EXCH_XY_RL( Du_SI, myThid )
439 _EXCH_XY_RL( Dv_SI, myThid )
440 _EXCH_XY_RL( Ru_SI, myThid )
441 _EXCH_XY_RL( Rv_SI, myThid )
442 _EXCH_XY_RL( cg_Uin, myThid )
443 _EXCH_XY_RL( cg_Vin, myThid )
444 endif
445
446 endif
447 enddo ! end of CG loop
448
449 c to avoid using "exit"
450 c if iters has reached max_iters there is no convergence
451
452 IF (iters .lt. streamice_max_cg_iter) THEN
453 conv_flag = 1
454 ENDIF
455
456 DO bj = myByLo(myThid), myByHi(myThid)
457 DO bi = myBxLo(myThid), myBxHi(myThid)
458 DO j=1-OLy,sNy+OLy
459 DO i=1-OLy,sNx+OLy
460 IF (STREAMICE_umask(i,j,bi,bj).eq.3.0)
461 & cg_Uin(i,j,bi,bj)=u_bdry_values_SI(i,j,bi,bj)
462 IF (STREAMICE_vmask(i,j,bi,bj).eq.3.0)
463 & cg_Vin(i,j,bi,bj)=v_bdry_values_SI(i,j,bi,bj)
464 ENDDO
465 ENDDO
466 ENDDO
467 ENDDO
468
469 _EXCH_XY_RL( cg_Uin, myThid )
470 _EXCH_XY_RL( cg_Vin, myThid )
471
472 #endif
473 RETURN
474 END
475
476
477
478
479
480
481
482

  ViewVC Help
Powered by ViewVC 1.1.22