/[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.4 - (hide annotations) (download)
Wed Aug 27 19:29:13 2014 UTC (10 years, 10 months ago) by dgoldberg
Branch: MAIN
CVS Tags: HEAD
Changes since 1.3: +1 -1 lines
updating contrib streamice repo with latest files, and separated out convergence checks; and parameterised maximum iteration counts and interface w shelfice for coupling

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

  ViewVC Help
Powered by ViewVC 1.1.22