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

1 C $Header: /u/gcmpack/MITgcm_contrib/dgoldberg/streamice/streamice_cg_solve.F,v 1.5 2012/07/26 16:13:18 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(
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 A_uu, ! section of matrix that multiplies u and projects on u
15 I A_uv, ! section of matrix that multiplies v and projects on u
16 I A_vu, ! section of matrix that multiplies u and projects on v
17 I A_vv, ! section of matrix that multiplies v and projects on v
18 I tolerance,
19 O iters,
20 I myThid )
21 C /============================================================\
22 C | SUBROUTINE |
23 C | o |
24 C |============================================================|
25 C | |
26 C \============================================================/
27 IMPLICIT NONE
28
29 C === Global variables ===
30 #include "SIZE.h"
31 #include "EEPARAMS.h"
32 #include "PARAMS.h"
33 #include "STREAMICE.h"
34 #include "STREAMICE_CG.h"
35
36
37 C !INPUT/OUTPUT ARGUMENTS
38 C cg_Uin, cg_Vin - input and output velocities
39 C cg_Bu, cg_Bv - driving stress
40 INTEGER myThid
41 INTEGER iters
42 _RL tolerance
43 _RL cg_Uin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
44 _RL cg_Vin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
45 _RL cg_Bu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
46 _RL cg_Bv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
47 _RL
48 & A_uu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
49 & A_vu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
50 & A_uv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
51 & A_vv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1)
52
53 C LOCAL VARIABLES
54 INTEGER i, j, bi, bj, cg_halo, conv_flag
55 INTEGER iter, is, js, ie, je, colx, coly, k
56 _RL dot_p1, dot_p2, alpha_k, beta_k, resid, resid_0
57 _RL dot_p1_tile (nSx,nSy)
58 _RL dot_p2_tile (nSx,nSy)
59 CHARACTER*(MAX_LEN_MBUF) msgBuf
60
61 iters = streamice_max_cg_iter
62
63 #ifdef ALLOW_STREAMICE
64
65
66 conv_flag = 0
67
68
69 DO bj = myByLo(myThid), myByHi(myThid)
70 DO bi = myBxLo(myThid), myBxHi(myThid)
71 DO j=1,sNy
72 DO i=1,sNx
73 Zu_SI (i,j,bi,bj) = 0. _d 0
74 Zv_SI (i,j,bi,bj) = 0. _d 0
75 Ru_SI (i,j,bi,bj) = 0. _d 0
76 Rv_SI (i,j,bi,bj) = 0. _d 0
77 Au_SI (i,j,bi,bj) = 0. _d 0
78 Av_SI (i,j,bi,bj) = 0. _d 0
79 Du_SI (i,j,bi,bj) = 0. _d 0
80 Dv_SI (i,j,bi,bj) = 0. _d 0
81 ENDDO
82 ENDDO
83 ENDDO
84 ENDDO
85
86 C FIND INITIAL RESIDUAL, and initialize r
87
88 ! #ifdef STREAMICE_CONSTRUCT_MATRIX
89
90
91
92 DO bj = myByLo(myThid), myByHi(myThid)
93 DO bi = myBxLo(myThid), myBxHi(myThid)
94 DO j=1,sNy
95 DO i=1,sNx
96 DO colx=-1,1
97 DO coly=-1,1
98 Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) +
99 & A_uu(i,j,bi,bj,colx,coly)*
100 & cg_Uin(i+colx,j+coly,bi,bj)+
101 & A_uv(i,j,bi,bj,colx,coly)*
102 & cg_Vin(i+colx,j+coly,bi,bj)
103
104
105 Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +
106 & A_vu(i,j,bi,bj,colx,coly)*
107 & cg_Uin(i+colx,j+coly,bi,bj)+
108 & A_vv(i,j,bi,bj,colx,coly)*
109 & cg_Vin(i+colx,j+coly,bi,bj)
110 ENDDO
111 ENDDO
112 ENDDO
113 ENDDO
114 ENDDO
115 ENDDO
116
117
118 _EXCH_XY_RL( Au_SI, myThid )
119 _EXCH_XY_RL( Av_SI, myThid )
120
121
122 DO bj = myByLo(myThid), myByHi(myThid)
123 DO bi = myBxLo(myThid), myBxHi(myThid)
124 DO j=1-OLy,sNy+OLy
125 DO i=1-OLx,sNx+OLx
126 Ru_SI(i,j,bi,bj)=cg_Bu(i,j,bi,bj)-
127 & Au_SI(i,j,bi,bj)
128 Rv_SI(i,j,bi,bj)=cg_Bv(i,j,bi,bj)-
129 & Av_SI(i,j,bi,bj)
130 ENDDO
131 ENDDO
132 dot_p1_tile(bi,bj) = 0. _d 0
133 dot_p2_tile(bi,bj) = 0. _d 0
134 ENDDO
135 ENDDO
136
137
138
139 DO bj = myByLo(myThid), myByHi(myThid)
140 DO bi = myBxLo(myThid), myBxHi(myThid)
141 DO j=1,sNy
142 DO i=1,sNx
143 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
144 & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Ru_SI(i,j,bi,bj)**2
145 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
146 & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Rv_SI(i,j,bi,bj)**2
147 ENDDO
148 ENDDO
149 ENDDO
150 ENDDO
151
152 CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
153 resid_0 = sqrt(dot_p1)
154
155 WRITE(msgBuf,'(A,E14.7)') 'CONJ GRAD INIT RESID, ',
156 & resid_0
157 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
158 & SQUEEZE_RIGHT , 1)
159
160 C CCCCCCCCCCCCCCCCCCCC
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 cg_halo = min(OLx-1,OLy-1)
176 conv_flag = 0
177
178 DO bj = myByLo(myThid), myByHi(myThid)
179 DO bi = myBxLo(myThid), myBxHi(myThid)
180 DO j=1-OLy,sNy+OLy
181 DO i=1-OLx,sNx+OLx
182 Du_SI(i,j,bi,bj)=Zu_SI(i,j,bi,bj)
183 Dv_SI(i,j,bi,bj)=Zv_SI(i,j,bi,bj)
184 ENDDO
185 ENDDO
186 ENDDO
187 ENDDO
188
189 resid = resid_0
190 iters = 0
191
192 c !!!!!!!!!!!!!!!!!!
193 c !! !!
194 c !! MAIN CG LOOP !!
195 c !! !!
196 c !!!!!!!!!!!!!!!!!!
197
198
199
200
201 c ! initially, b-grid data is valid up to 3 halo nodes out -- right? (check for MITgcm!!)
202
203 WRITE(msgBuf,'(A)') 'BEGINNING MAIN CG LOOP'
204 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
205 & SQUEEZE_RIGHT , 1)
206
207 ! IF(STREAMICE_construct_matrix) CALL STREAMICE_CG_MAKE_A(myThid)
208
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 ENDDO
228 ENDDO
229 ENDDO
230 ENDDO
231
232 ! IF (STREAMICE_construct_matrix) THEN
233
234 ! #ifdef STREAMICE_CONSTRUCT_MATRIX
235
236 DO bj = myByLo(myThid), myByHi(myThid)
237 DO bi = myBxLo(myThid), myBxHi(myThid)
238 DO j=js,je
239 DO i=is,ie
240 DO colx=-1,1
241 DO coly=-1,1
242 Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) +
243 & A_uu(i,j,bi,bj,colx,coly)*
244 & Du_SI(i+colx,j+coly,bi,bj)+
245 & A_uv(i,j,bi,bj,colx,coly)*
246 & Dv_SI(i+colx,j+coly,bi,bj)
247 Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +
248 & A_vu(i,j,bi,bj,colx,coly)*
249 & Du_SI(i+colx,j+coly,bi,bj)+
250 & A_vv(i,j,bi,bj,colx,coly)*
251 & Dv_SI(i+colx,j+coly,bi,bj)
252 ENDDO
253 ENDDO
254 ENDDO
255 ENDDO
256 ENDDO
257 ENDDO
258
259 ! else
260 ! #else
261 !
262 ! CALL STREAMICE_CG_ACTION( myThid,
263 ! O Au_SI,
264 ! O Av_SI,
265 ! I Du_SI,
266 ! I Dv_SI,
267 ! I is,ie,js,je)
268 !
269 ! ! ENDIF
270 !
271 ! #endif
272
273
274 DO bj = myByLo(myThid), myByHi(myThid)
275 DO bi = myBxLo(myThid), myBxHi(myThid)
276 dot_p1_tile(bi,bj) = 0. _d 0
277 dot_p2_tile(bi,bj) = 0. _d 0
278 ENDDO
279 ENDDO
280
281 DO bj = myByLo(myThid), myByHi(myThid)
282 DO bi = myBxLo(myThid), myBxHi(myThid)
283 DO j=1,sNy
284 DO i=1,sNx
285 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN
286 dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zu_SI(i,j,bi,bj)*
287 & Ru_SI(i,j,bi,bj)
288 dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Du_SI(i,j,bi,bj)*
289 & Au_SI(i,j,bi,bj)
290 ENDIF
291 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN
292 dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zv_SI(i,j,bi,bj)*
293 & Rv_SI(i,j,bi,bj)
294 dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Dv_SI(i,j,bi,bj)*
295 & Av_SI(i,j,bi,bj)
296 ENDIF
297 ENDDO
298 ENDDO
299 ENDDO
300 ENDDO
301
302 CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
303 CALL GLOBAL_SUM_TILE_RL( dot_p2_tile, dot_p2, myThid )
304 alpha_k = dot_p1/dot_p2
305
306 DO bj = myByLo(myThid), myByHi(myThid)
307 DO bi = myBxLo(myThid), myBxHi(myThid)
308 DO j=1-OLy,sNy+OLy
309 DO i=1-OLx,sNx+OLx
310
311 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN
312 cg_Uin(i,j,bi,bj)=cg_Uin(i,j,bi,bj)+
313 & alpha_k*Du_SI(i,j,bi,bj)
314 Ru_old_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj)
315 Zu_old_SI(i,j,bi,bj) = Zu_SI(i,j,bi,bj)
316 Ru_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj)-
317 & alpha_k*Au_SI(i,j,bi,bj)
318 Zu_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj) /
319 & DIAGu_SI(i,j,bi,bj)
320 ENDIF
321
322 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN
323 cg_Vin(i,j,bi,bj)=cg_Vin(i,j,bi,bj)+
324 & alpha_k*Dv_SI(i,j,bi,bj)
325 Rv_old_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj)
326 Zv_old_SI(i,j,bi,bj) = Zv_SI(i,j,bi,bj)
327 Rv_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj)-
328 & alpha_k*Av_SI(i,j,bi,bj)
329 Zv_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj) /
330 & DIAGv_SI(i,j,bi,bj)
331
332 ENDIF
333 ENDDO
334 ENDDO
335 ENDDO
336 ENDDO
337
338 DO bj = myByLo(myThid), myByHi(myThid)
339 DO bi = myBxLo(myThid), myBxHi(myThid)
340 dot_p1_tile(bi,bj) = 0. _d 0
341 dot_p2_tile(bi,bj) = 0. _d 0
342 ENDDO
343 ENDDO
344
345 DO bj = myByLo(myThid), myByHi(myThid)
346 DO bi = myBxLo(myThid), myBxHi(myThid)
347 DO j=1,sNy
348 DO i=1,sNx
349
350 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN
351 dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zu_SI(i,j,bi,bj)*
352 & Ru_SI(i,j,bi,bj)
353 dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Zu_old_SI(i,j,bi,bj)*
354 & Ru_old_SI(i,j,bi,bj)
355 ENDIF
356
357 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN
358 dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zv_SI(i,j,bi,bj)*
359 & Rv_SI(i,j,bi,bj)
360 dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Zv_old_SI(i,j,bi,bj)*
361 & Rv_old_SI(i,j,bi,bj)
362 ENDIF
363
364 ENDDO
365 ENDDO
366 ENDDO
367 ENDDO
368
369 CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
370 CALL GLOBAL_SUM_TILE_RL( dot_p2_tile, dot_p2, myThid )
371
372 beta_k = dot_p1/dot_p2
373
374 DO bj = myByLo(myThid), myByHi(myThid)
375 DO bi = myBxLo(myThid), myBxHi(myThid)
376 DO j=1-OLy,sNy+OLy
377 DO i=1-OLx,sNx+OLx
378 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
379 & Du_SI(i,j,bi,bj)=beta_k*Du_SI(i,j,bi,bj)+
380 & Zu_SI(i,j,bi,bj)
381 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
382 & Dv_SI(i,j,bi,bj)=beta_k*Dv_SI(i,j,bi,bj)+
383 & Zv_SI(i,j,bi,bj)
384 ENDDO
385 ENDDO
386 ENDDO
387 ENDDO
388
389 DO bj = myByLo(myThid), myByHi(myThid)
390 DO bi = myBxLo(myThid), myBxHi(myThid)
391 dot_p1_tile(bi,bj) = 0. _d 0
392 ENDDO
393 ENDDO
394
395 DO bj = myByLo(myThid), myByHi(myThid)
396 DO bi = myBxLo(myThid), myBxHi(myThid)
397 DO j=1,sNy
398 DO i=1,sNx
399 IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
400 & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Ru_SI(i,j,bi,bj)**2
401 IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
402 & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Rv_SI(i,j,bi,bj)**2
403 ENDDO
404 ENDDO
405 ENDDO
406 ENDDO
407
408 CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
409 resid = sqrt(dot_p1)
410
411 ! IF (iter .eq. 1) then
412 ! print *, alpha_k, beta_k, resid
413 ! ENDIF
414
415 cg_halo = cg_halo - 1
416
417 if (cg_halo .eq. 0) then
418 cg_halo = min(OLx-1,OLy-1)
419 _EXCH_XY_RL( Du_SI, myThid )
420 _EXCH_XY_RL( Dv_SI, myThid )
421 _EXCH_XY_RL( Ru_SI, myThid )
422 _EXCH_XY_RL( Rv_SI, myThid )
423 _EXCH_XY_RL( cg_Uin, myThid )
424 _EXCH_XY_RL( cg_Vin, myThid )
425 endif
426
427
428 endif
429 enddo ! end of CG loop
430
431 c to avoid using "exit"
432 c if iters has reached max_iters there is no convergence
433
434 IF (iters .lt. streamice_max_cg_iter) THEN
435 conv_flag = 1
436 ENDIF
437
438 ! DO bj = myByLo(myThid), myByHi(myThid)
439 ! DO bi = myBxLo(myThid), myBxHi(myThid)
440 ! DO j=1-OLy,sNy+OLy
441 ! DO i=1-OLy,sNx+OLy
442 ! IF (STREAMICE_umask(i,j,bi,bj).eq.3.0)
443 ! & cg_Uin(i,j,bi,bj)=u_bdry_values_SI(i,j,bi,bj)
444 ! IF (STREAMICE_vmask(i,j,bi,bj).eq.3.0)
445 ! & cg_Vin(i,j,bi,bj)=v_bdry_values_SI(i,j,bi,bj)
446 ! ENDDO
447 ! ENDDO
448 ! ENDDO
449 ! ENDDO
450 !
451 ! _EXCH_XY_RL( cg_Uin, myThid )
452 ! _EXCH_XY_RL( cg_Vin, myThid )
453
454 #endif
455 RETURN
456 END
457
458
459
460
461
462
463
464

  ViewVC Help
Powered by ViewVC 1.1.22