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 |
|