/[MITgcm]/MITgcm_contrib/dgoldberg/streamice/streamice_cg_solve.F
ViewVC logotype

Annotation 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 - (hide 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 dgoldberg 1.4 C $Header: /u/gcmpack/MITgcm_contrib/dgoldberg/streamice/streamice_cg_solve.F,v 1.3 2012/05/14 16:52:29 dgoldberg Exp $
2 heimbach 1.1 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 dgoldberg 1.4 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 heimbach 1.1 I tolerance,
19 heimbach 1.2 O iters,
20     I myThid )
21 heimbach 1.1 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 dgoldberg 1.4 _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 heimbach 1.1
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 dgoldberg 1.4 _RL dot_p1, dot_p2, alpha_k, beta_k, resid, resid_0
57 heimbach 1.1 _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 dgoldberg 1.4
65 heimbach 1.1 conv_flag = 0
66    
67 dgoldberg 1.4
68 heimbach 1.1 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 dgoldberg 1.4
85     C FIND INITIAL RESIDUAL, and initialize r
86 heimbach 1.1
87 dgoldberg 1.4 ! #ifdef STREAMICE_CONSTRUCT_MATRIX
88 heimbach 1.1
89 dgoldberg 1.4
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 heimbach 1.1
103 dgoldberg 1.4
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 heimbach 1.1 ENDDO
114     ENDDO
115    
116 dgoldberg 1.4
117 heimbach 1.1 _EXCH_XY_RL( Au_SI, myThid )
118     _EXCH_XY_RL( Av_SI, myThid )
119    
120 dgoldberg 1.4
121 heimbach 1.1 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 dgoldberg 1.4 Ru_SI(i,j,bi,bj)=cg_Bu(i,j,bi,bj)-
126 heimbach 1.1 & Au_SI(i,j,bi,bj)
127 dgoldberg 1.4 Rv_SI(i,j,bi,bj)=cg_Bv(i,j,bi,bj)-
128 heimbach 1.1 & 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 dgoldberg 1.4
137    
138 heimbach 1.1 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 dgoldberg 1.4 C CCCCCCCCCCCCCCCCCCCC
155    
156 heimbach 1.1 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 dgoldberg 1.4
192    
193    
194 heimbach 1.1
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 dgoldberg 1.3 ! IF(STREAMICE_construct_matrix) CALL STREAMICE_CG_MAKE_A(myThid)
200    
201 heimbach 1.1
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 dgoldberg 1.3 ! IF (STREAMICE_construct_matrix) THEN
225    
226 dgoldberg 1.4 ! #ifdef STREAMICE_CONSTRUCT_MATRIX
227 dgoldberg 1.3
228 heimbach 1.1 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 dgoldberg 1.4 & A_uu(i,j,bi,bj,colx,coly)*
236 heimbach 1.1 & Du_SI(i+colx,j+coly,bi,bj)+
237 dgoldberg 1.4 & A_uv(i,j,bi,bj,colx,coly)*
238 heimbach 1.1 & Dv_SI(i+colx,j+coly,bi,bj)
239     Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +
240 dgoldberg 1.4 & A_vu(i,j,bi,bj,colx,coly)*
241 heimbach 1.1 & Du_SI(i+colx,j+coly,bi,bj)+
242 dgoldberg 1.4 & A_vv(i,j,bi,bj,colx,coly)*
243 heimbach 1.1 & Dv_SI(i+colx,j+coly,bi,bj)
244     ENDDO
245     ENDDO
246     ENDDO
247     ENDDO
248     ENDDO
249     ENDDO
250    
251 dgoldberg 1.3 ! else
252 dgoldberg 1.4 ! #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 heimbach 1.1
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 dgoldberg 1.4
420 heimbach 1.1 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 dgoldberg 1.4 ! 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 heimbach 1.1
446     #endif
447     RETURN
448     END
449    
450    
451    
452    
453    
454    
455    
456    

  ViewVC Help
Powered by ViewVC 1.1.22