/[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.2 - (hide annotations) (download)
Thu Jul 26 16:22:58 2012 UTC (13 years ago) by dgoldberg
Branch: MAIN
Changes since 1.1: +4 -2 lines
more replacing print statements

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

  ViewVC Help
Powered by ViewVC 1.1.22