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