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

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

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


Revision 1.1 - (show annotations) (download)
Thu Jul 19 18:56:00 2012 UTC (13 years ago) by dgoldberg
Branch: MAIN
separate subroutine for matrix-free solve; should not be used with TAF

1 C $Header: /u/gcmpack/MITgcm_contrib/dgoldberg/streamice/streamice_cg_solve.F,v 1.3 2012/05/14 16:52:29 dgoldberg 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_MATFREE (
10 U cg_Uin, ! x-velocities
11 U cg_Vin, ! y-velocities
12 I cg_Bu, ! force in x dir
13 I cg_Bv, ! force in y dir
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
45 C LOCAL VARIABLES
46 INTEGER i, j, bi, bj, cg_halo, conv_flag
47 INTEGER iter, is, js, ie, je, colx, coly, k
48 _RL dot_p1, dot_p2, alpha_k, beta_k, resid, resid_0
49 _RL dot_p1_tile (nSx,nSy)
50 _RL dot_p2_tile (nSx,nSy)
51
52 iters = streamice_max_cg_iter
53
54 #ifdef ALLOW_STREAMICE
55
56
57 conv_flag = 0
58
59 DO bj = myByLo(myThid), myByHi(myThid)
60 DO bi = myBxLo(myThid), myBxHi(myThid)
61 DO j=1,sNy
62 DO i=1,sNx
63 Zu_SI (i,j,bi,bj) = 0. _d 0
64 Zv_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 ENDDO
72 ENDDO
73 ENDDO
74 ENDDO
75
76 C FIND INITIAL RESIDUAL, and initialize r
77
78 ! #ifdef STREAMICE_CONSTRUCT_MATRIX
79
80 ! DO bj = myByLo(myThid), myByHi(myThid)
81 ! DO bi = myBxLo(myThid), myBxHi(myThid)
82 ! DO j=js,je
83 ! DO i=is,ie
84 ! DO colx=-1,1
85 ! DO coly=-1,1
86 ! Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) +
87 ! & streamice_cg_A1(i,j,bi,bj,colx,coly)*
88 ! & cg_Uin(i+colx,j+coly,bi,bj)+
89 ! & streamice_cg_A2(i,j,bi,bj,colx,coly)*
90 ! & cg_Vin(i+colx,j+coly,bi,bj)
91 ! Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +
92 ! & streamice_cg_A3(i,j,bi,bj,colx,coly)*
93 ! & cg_Uin(i+colx,j+coly,bi,bj)+
94 ! & streamice_cg_A4(i,j,bi,bj,colx,coly)*
95 ! & cg_Vin(i+colx,j+coly,bi,bj)
96 ! ENDDO
97 ! ENDDO
98 ! ENDDO
99 ! ENDDO
100 ! ENDDO
101 ! ENDDO
102
103 ! #else
104
105 CALL STREAMICE_CG_ACTION( myThid,
106 O Au_SI,
107 O Av_SI,
108 I cg_Uin,
109 I cg_Vin,
110 I 0, sNx+1, 0, sNy+1 )
111
112 ! #endif
113
114 _EXCH_XY_RL( Au_SI, myThid )
115 _EXCH_XY_RL( Av_SI, myThid )
116
117 DO bj = myByLo(myThid), myByHi(myThid)
118 DO bi = myBxLo(myThid), myBxHi(myThid)
119 DO j=1-OLy,sNy+OLy
120 DO i=1-OLx,sNx+OLx
121 Ru_SI(i,j,bi,bj)=cg_Bu(i,j,bi,bj)-
122 & Au_SI(i,j,bi,bj)
123 Rv_SI(i,j,bi,bj)=cg_Bv(i,j,bi,bj)-
124 & Av_SI(i,j,bi,bj)
125 ENDDO
126 ENDDO
127 dot_p1_tile(bi,bj) = 0. _d 0
128 dot_p2_tile(bi,bj) = 0. _d 0
129 ENDDO
130 ENDDO
131
132 DO bj = myByLo(myThid), myByHi(myThid)
133 DO bi = myBxLo(myThid), myBxHi(myThid)
134 DO j=1,sNy
135 DO i=1,sNx
136 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
137 & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Ru_SI(i,j,bi,bj)**2
138 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
139 & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Rv_SI(i,j,bi,bj)**2
140 ENDDO
141 ENDDO
142 ENDDO
143 ENDDO
144
145 CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
146 resid_0 = sqrt(dot_p1)
147
148 C CCCCCCCCCCCCCCCCCCCC
149
150 DO bj = myByLo(myThid), myByHi(myThid)
151 DO bi = myBxLo(myThid), myBxHi(myThid)
152 DO j=1-OLy,sNy+OLy
153 DO i=1-OLx,sNx+OLx
154 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
155 & Zu_SI(i,j,bi,bj)=Ru_SI(i,j,bi,bj) / DIAGu_SI(i,j,bi,bj)
156 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
157 & Zv_SI(i,j,bi,bj)=Rv_SI(i,j,bi,bj) / DIAGv_SI(i,j,bi,bj)
158 ENDDO
159 ENDDO
160 ENDDO
161 ENDDO
162
163 cg_halo = min(OLx-1,OLy-1)
164 conv_flag = 0
165
166 DO bj = myByLo(myThid), myByHi(myThid)
167 DO bi = myBxLo(myThid), myBxHi(myThid)
168 DO j=1-OLy,sNy+OLy
169 DO i=1-OLx,sNx+OLx
170 Du_SI(i,j,bi,bj)=Zu_SI(i,j,bi,bj)
171 Dv_SI(i,j,bi,bj)=Zv_SI(i,j,bi,bj)
172 ENDDO
173 ENDDO
174 ENDDO
175 ENDDO
176
177 resid = resid_0
178 iters = 0
179
180 c !!!!!!!!!!!!!!!!!!
181 c !! !!
182 c !! MAIN CG LOOP !!
183 c !! !!
184 c !!!!!!!!!!!!!!!!!!
185
186
187
188
189 c ! initially, b-grid data is valid up to 3 halo nodes out -- right? (check for MITgcm!!)
190
191 print *, "BEGINNING MAIN CG LOOP"
192
193 ! IF(STREAMICE_construct_matrix) CALL STREAMICE_CG_MAKE_A(myThid)
194
195
196 do iter = 1, streamice_max_cg_iter
197 if (resid .gt. tolerance*resid_0) then
198
199 c to avoid using "exit"
200 iters = iters + 1
201
202 is = 1 - cg_halo
203 ie = sNx + cg_halo
204 js = 1 - cg_halo
205 je = sNy + cg_halo
206
207 DO bj = myByLo(myThid), myByHi(myThid)
208 DO bi = myBxLo(myThid), myBxHi(myThid)
209 DO j=1-OLy,sNy+OLy
210 DO i=1-OLx,sNx+OLx
211 Au_SI(i,j,bi,bj) = 0. _d 0
212 Av_SI(i,j,bi,bj) = 0. _d 0
213 ENDDO
214 ENDDO
215 ENDDO
216 ENDDO
217
218 ! IF (STREAMICE_construct_matrix) THEN
219
220 ! #ifdef STREAMICE_CONSTRUCT_MATRIX
221 !
222 ! DO bj = myByLo(myThid), myByHi(myThid)
223 ! DO bi = myBxLo(myThid), myBxHi(myThid)
224 ! DO j=js,je
225 ! DO i=is,ie
226 ! DO colx=-1,1
227 ! DO coly=-1,1
228 ! Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) +
229 ! & streamice_cg_A1(i,j,bi,bj,colx,coly)*
230 ! & Du_SI(i+colx,j+coly,bi,bj)+
231 ! & streamice_cg_A2(i,j,bi,bj,colx,coly)*
232 ! & Dv_SI(i+colx,j+coly,bi,bj)
233 ! Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +
234 ! & streamice_cg_A3(i,j,bi,bj,colx,coly)*
235 ! & Du_SI(i+colx,j+coly,bi,bj)+
236 ! & streamice_cg_A4(i,j,bi,bj,colx,coly)*
237 ! & Dv_SI(i+colx,j+coly,bi,bj)
238 ! ENDDO
239 ! ENDDO
240 ! ENDDO
241 ! ENDDO
242 ! ENDDO
243 ! ENDDO
244 !
245 ! ! else
246 ! #else
247
248 CALL STREAMICE_CG_ACTION( myThid,
249 O Au_SI,
250 O Av_SI,
251 I Du_SI,
252 I Dv_SI,
253 I is,ie,js,je)
254
255 ! ENDIF
256
257 ! #endif
258
259
260 DO bj = myByLo(myThid), myByHi(myThid)
261 DO bi = myBxLo(myThid), myBxHi(myThid)
262 dot_p1_tile(bi,bj) = 0. _d 0
263 dot_p2_tile(bi,bj) = 0. _d 0
264 ENDDO
265 ENDDO
266
267 DO bj = myByLo(myThid), myByHi(myThid)
268 DO bi = myBxLo(myThid), myBxHi(myThid)
269 DO j=1,sNy
270 DO i=1,sNx
271 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN
272 dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zu_SI(i,j,bi,bj)*
273 & Ru_SI(i,j,bi,bj)
274 dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Du_SI(i,j,bi,bj)*
275 & Au_SI(i,j,bi,bj)
276 ENDIF
277 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN
278 dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zv_SI(i,j,bi,bj)*
279 & Rv_SI(i,j,bi,bj)
280 dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Dv_SI(i,j,bi,bj)*
281 & Av_SI(i,j,bi,bj)
282 ENDIF
283 ENDDO
284 ENDDO
285 ENDDO
286 ENDDO
287
288 CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
289 CALL GLOBAL_SUM_TILE_RL( dot_p2_tile, dot_p2, myThid )
290 alpha_k = dot_p1/dot_p2
291
292 DO bj = myByLo(myThid), myByHi(myThid)
293 DO bi = myBxLo(myThid), myBxHi(myThid)
294 DO j=1-OLy,sNy+OLy
295 DO i=1-OLx,sNx+OLx
296
297 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN
298 cg_Uin(i,j,bi,bj)=cg_Uin(i,j,bi,bj)+
299 & alpha_k*Du_SI(i,j,bi,bj)
300 Ru_old_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj)
301 Zu_old_SI(i,j,bi,bj) = Zu_SI(i,j,bi,bj)
302 Ru_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj)-
303 & alpha_k*Au_SI(i,j,bi,bj)
304 Zu_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj) /
305 & DIAGu_SI(i,j,bi,bj)
306 ENDIF
307
308 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN
309 cg_Vin(i,j,bi,bj)=cg_Vin(i,j,bi,bj)+
310 & alpha_k*Dv_SI(i,j,bi,bj)
311 Rv_old_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj)
312 Zv_old_SI(i,j,bi,bj) = Zv_SI(i,j,bi,bj)
313 Rv_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj)-
314 & alpha_k*Av_SI(i,j,bi,bj)
315 Zv_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj) /
316 & DIAGv_SI(i,j,bi,bj)
317
318 ENDIF
319 ENDDO
320 ENDDO
321 ENDDO
322 ENDDO
323
324 DO bj = myByLo(myThid), myByHi(myThid)
325 DO bi = myBxLo(myThid), myBxHi(myThid)
326 dot_p1_tile(bi,bj) = 0. _d 0
327 dot_p2_tile(bi,bj) = 0. _d 0
328 ENDDO
329 ENDDO
330
331 DO bj = myByLo(myThid), myByHi(myThid)
332 DO bi = myBxLo(myThid), myBxHi(myThid)
333 DO j=1,sNy
334 DO i=1,sNx
335
336 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN
337 dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zu_SI(i,j,bi,bj)*
338 & Ru_SI(i,j,bi,bj)
339 dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Zu_old_SI(i,j,bi,bj)*
340 & Ru_old_SI(i,j,bi,bj)
341 ENDIF
342
343 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN
344 dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zv_SI(i,j,bi,bj)*
345 & Rv_SI(i,j,bi,bj)
346 dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Zv_old_SI(i,j,bi,bj)*
347 & Rv_old_SI(i,j,bi,bj)
348 ENDIF
349
350 ENDDO
351 ENDDO
352 ENDDO
353 ENDDO
354
355 CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
356 CALL GLOBAL_SUM_TILE_RL( dot_p2_tile, dot_p2, myThid )
357
358 beta_k = dot_p1/dot_p2
359
360 DO bj = myByLo(myThid), myByHi(myThid)
361 DO bi = myBxLo(myThid), myBxHi(myThid)
362 DO j=1-OLy,sNy+OLy
363 DO i=1-OLx,sNx+OLx
364 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
365 & Du_SI(i,j,bi,bj)=beta_k*Du_SI(i,j,bi,bj)+
366 & Zu_SI(i,j,bi,bj)
367 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
368 & Dv_SI(i,j,bi,bj)=beta_k*Dv_SI(i,j,bi,bj)+
369 & Zv_SI(i,j,bi,bj)
370 ENDDO
371 ENDDO
372 ENDDO
373 ENDDO
374
375 DO bj = myByLo(myThid), myByHi(myThid)
376 DO bi = myBxLo(myThid), myBxHi(myThid)
377 dot_p1_tile(bi,bj) = 0. _d 0
378 ENDDO
379 ENDDO
380
381 DO bj = myByLo(myThid), myByHi(myThid)
382 DO bi = myBxLo(myThid), myBxHi(myThid)
383 DO j=1,sNy
384 DO i=1,sNx
385 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
386 & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Ru_SI(i,j,bi,bj)**2
387 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
388 & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Rv_SI(i,j,bi,bj)**2
389 ENDDO
390 ENDDO
391 ENDDO
392 ENDDO
393
394 CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
395 resid = sqrt(dot_p1)
396
397 ! IF (iter .eq. 1) then
398 ! print *, alpha_k, beta_k, resid
399 ! ENDIF
400
401 cg_halo = cg_halo - 1
402
403 if (cg_halo .eq. 0) then
404 cg_halo = min(OLx-1,OLy-1)
405 _EXCH_XY_RL( Du_SI, myThid )
406 _EXCH_XY_RL( Dv_SI, myThid )
407 _EXCH_XY_RL( Ru_SI, myThid )
408 _EXCH_XY_RL( Rv_SI, myThid )
409 _EXCH_XY_RL( cg_Uin, myThid )
410 _EXCH_XY_RL( cg_Vin, myThid )
411 endif
412
413 endif
414 enddo ! end of CG loop
415
416 c to avoid using "exit"
417 c if iters has reached max_iters there is no convergence
418
419 IF (iters .lt. streamice_max_cg_iter) THEN
420 conv_flag = 1
421 ENDIF
422
423 ! DO bj = myByLo(myThid), myByHi(myThid)
424 ! DO bi = myBxLo(myThid), myBxHi(myThid)
425 ! DO j=1-OLy,sNy+OLy
426 ! DO i=1-OLy,sNx+OLy
427 ! IF (STREAMICE_umask(i,j,bi,bj).eq.3.0)
428 ! & cg_Uin(i,j,bi,bj)=u_bdry_values_SI(i,j,bi,bj)
429 ! IF (STREAMICE_vmask(i,j,bi,bj).eq.3.0)
430 ! & cg_Vin(i,j,bi,bj)=v_bdry_values_SI(i,j,bi,bj)
431 ! ENDDO
432 ! ENDDO
433 ! ENDDO
434 ! ENDDO
435 !
436 ! _EXCH_XY_RL( cg_Uin, myThid )
437 ! _EXCH_XY_RL( cg_Vin, myThid )
438
439 #endif
440 RETURN
441 END
442
443
444
445
446
447
448
449

  ViewVC Help
Powered by ViewVC 1.1.22