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