/[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.4 - (show annotations) (download)
Thu Jul 19 18:46:56 2012 UTC (13 years ago) by dgoldberg
Branch: MAIN
Changes since 1.3: +87 -124 lines
TAF compatibility

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

  ViewVC Help
Powered by ViewVC 1.1.22