/[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.2 - (show annotations) (download)
Thu Jul 26 16:22:58 2012 UTC (13 years ago) by dgoldberg
Branch: MAIN
Changes since 1.1: +4 -2 lines
more replacing print statements

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

  ViewVC Help
Powered by ViewVC 1.1.22