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

Annotation 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.1 - (hide annotations) (download)
Thu Jul 19 18:56:00 2012 UTC (13 years ago) by dgoldberg
Branch: MAIN
separate subroutine for matrix-free solve; should not be used with TAF

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

  ViewVC Help
Powered by ViewVC 1.1.22