/[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.4 - (show annotations) (download)
Wed Aug 27 19:29:13 2014 UTC (10 years, 10 months ago) by dgoldberg
Branch: MAIN
CVS Tags: HEAD
Changes since 1.3: +1 -1 lines
Error occurred while calculating annotation data.
updating contrib streamice repo with latest files, and separated out convergence checks; and parameterised maximum iteration counts and interface w shelfice for coupling

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

  ViewVC Help
Powered by ViewVC 1.1.22