/[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.1 - (hide annotations) (download)
Thu Mar 29 15:59:21 2012 UTC (13 years, 3 months ago) by heimbach
Branch: MAIN
Initial check-in of Dan Goldberg's streamice package

1 heimbach 1.1 C $Header: /u/gcmpack/MITgcm/pkg/streamice/streamice_init_varia.F,v 1.6 2011/06/29 16:24:10 dng 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,
11     U cg_Vin,
12     I cg_Bu,
13     I cg_Bv,
14     I myThid,
15     I tolerance,
16     O iters )
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     C LOCAL VARIABLES
45     INTEGER i, j, bi, bj, cg_halo, conv_flag
46     INTEGER iter, is, js, ie, je, colx, coly, k
47     _RL dot_p1, dot_p2, resid_0, alpha_k, beta_k, resid
48     _RL dot_p1_tile (nSx,nSy)
49     _RL dot_p2_tile (nSx,nSy)
50    
51     iters = streamice_max_cg_iter
52    
53     #ifdef ALLOW_STREAMICE
54    
55     conv_flag = 0
56    
57     DO bj = myByLo(myThid), myByHi(myThid)
58     DO bi = myBxLo(myThid), myBxHi(myThid)
59     DO j=1,sNy
60     DO i=1,sNx
61     Zu_SI (i,j,bi,bj) = 0. _d 0
62     Zv_SI (i,j,bi,bj) = 0. _d 0
63     DIAGu_SI (i,j,bi,bj) = 0. _d 0
64     DIAGv_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     ubd_SI (i,j,bi,bj) = 0. _d 0
72     vbd_SI (i,j,bi,bj) = 0. _d 0
73     ENDDO
74     ENDDO
75     ENDDO
76     ENDDO
77    
78     C PREAMBLE: get bdry contribution, find matrix diagonal,
79     C initialize iterates R (residual), Z, and D
80    
81     CALL STREAMICE_CG_BOUND_VALS( myThid,
82     O ubd_SI,
83     O vbd_SI)
84    
85     DO bj = myByLo(myThid), myByHi(myThid)
86     DO bi = myBxLo(myThid), myBxHi(myThid)
87     DO j=1-OLy,sNy+OLy
88     DO i=1-OLx,sNx+OLx
89     RHSu_SI (i,j,bi,bj) = cg_Bu(i,j,bi,bj)
90     & - ubd_SI(i,j,bi,bj)
91     RHSv_SI (i,j,bi,bj) = cg_Bv(i,j,bi,bj)
92     & - vbd_SI(i,j,bi,bj)
93     ENDDO
94     ENDDO
95     ENDDO
96     ENDDO
97    
98     _EXCH_XY_RL( RHSu_SI, myThid )
99     _EXCH_XY_RL( RHSv_SI, myThid )
100    
101     CALL STREAMICE_CG_ADIAG( myThid,
102     O DIAGu_SI,
103     O DIAGv_SI)
104    
105     _EXCH_XY_RL( DIAGu_SI, myThid )
106     _EXCH_XY_RL( DIAGv_SI, myThid )
107    
108     C FIX PROBLEM WITH PRECOND LATER
109    
110     ! DO bj = myByLo(myThid), myByHi(myThid)
111     ! DO bi = myBxLo(myThid), myBxHi(myThid)
112     ! DO j=1-OLy,sNy+OLy
113     ! DO i=1-OLx,sNx+OLx
114     ! DIAGu_SI(i,j,bi,bj)=1.0
115     ! DIAGv_SI(i,j,bi,bj)=1.0
116     ! ENDDO
117     ! ENDDO
118     ! ENDDO
119     ! ENDDO
120    
121     CALL STREAMICE_CG_ACTION( myThid,
122     O Au_SI,
123     O Av_SI,
124     I cg_Uin,
125     I cg_Vin,
126     I 0, sNx+1, 0, sNy+1 )
127    
128     _EXCH_XY_RL( Au_SI, myThid )
129     _EXCH_XY_RL( Av_SI, myThid )
130    
131     DO bj = myByLo(myThid), myByHi(myThid)
132     DO bi = myBxLo(myThid), myBxHi(myThid)
133     DO j=1-OLy,sNy+OLy
134     DO i=1-OLx,sNx+OLx
135     Ru_SI(i,j,bi,bj)=RHSu_SI(i,j,bi,bj)-
136     & Au_SI(i,j,bi,bj)
137     Rv_SI(i,j,bi,bj)=RHSv_SI(i,j,bi,bj)-
138     & Av_SI(i,j,bi,bj)
139     ENDDO
140     ENDDO
141     dot_p1_tile(bi,bj) = 0. _d 0
142     dot_p2_tile(bi,bj) = 0. _d 0
143     ENDDO
144     ENDDO
145    
146     DO bj = myByLo(myThid), myByHi(myThid)
147     DO bi = myBxLo(myThid), myBxHi(myThid)
148     DO j=1,sNy
149     DO i=1,sNx
150     IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
151     & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Ru_SI(i,j,bi,bj)**2
152     IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
153     & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Rv_SI(i,j,bi,bj)**2
154     ENDDO
155     ENDDO
156     ENDDO
157     ENDDO
158    
159     CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
160     resid_0 = sqrt(dot_p1)
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    
176     cg_halo = min(OLx-1,OLy-1)
177     conv_flag = 0
178    
179     DO bj = myByLo(myThid), myByHi(myThid)
180     DO bi = myBxLo(myThid), myBxHi(myThid)
181     DO j=1-OLy,sNy+OLy
182     DO i=1-OLx,sNx+OLx
183     Du_SI(i,j,bi,bj)=Zu_SI(i,j,bi,bj)
184     Dv_SI(i,j,bi,bj)=Zv_SI(i,j,bi,bj)
185     ENDDO
186     ENDDO
187     ENDDO
188     ENDDO
189    
190     resid = resid_0
191     iters = 0
192    
193     c !!!!!!!!!!!!!!!!!!
194     c !! !!
195     c !! MAIN CG LOOP !!
196     c !! !!
197     c !!!!!!!!!!!!!!!!!!
198    
199     c ! initially, b-grid data is valid up to 3 halo nodes out -- right? (check for MITgcm!!)
200    
201     print *, "BEGINNING MAIN CG LOOP"
202    
203     IF(STREAMICE_construct_matrix) CALL STREAMICE_CG_MAKE_A(myThid)
204    
205     do iter = 1, streamice_max_cg_iter
206     if (resid .gt. tolerance*resid_0) then
207    
208     c to avoid using "exit"
209     iters = iters + 1
210    
211     is = 1 - cg_halo
212     ie = sNx + cg_halo
213     js = 1 - cg_halo
214     je = sNy + cg_halo
215    
216     DO bj = myByLo(myThid), myByHi(myThid)
217     DO bi = myBxLo(myThid), myBxHi(myThid)
218     DO j=1-OLy,sNy+OLy
219     DO i=1-OLx,sNx+OLx
220     Au_SI(i,j,bi,bj) = 0. _d 0
221     Av_SI(i,j,bi,bj) = 0. _d 0
222     ! Du_SI(i,j,bi,bj) = Real(i)
223     ! Dv_SI(i,j,bi,bj) = 0.0
224     ENDDO
225     ENDDO
226     ENDDO
227     ENDDO
228    
229     IF (STREAMICE_construct_matrix) THEN
230     DO bj = myByLo(myThid), myByHi(myThid)
231     DO bi = myBxLo(myThid), myBxHi(myThid)
232     DO j=js,je
233     DO i=is,ie
234     DO colx=-1,1
235     DO coly=-1,1
236     Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) +
237     & streamice_cg_A1(i,j,bi,bj,colx,coly)*
238     & Du_SI(i+colx,j+coly,bi,bj)+
239     & streamice_cg_A2(i,j,bi,bj,colx,coly)*
240     & Dv_SI(i+colx,j+coly,bi,bj)
241     Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +
242     & streamice_cg_A3(i,j,bi,bj,colx,coly)*
243     & Du_SI(i+colx,j+coly,bi,bj)+
244     & streamice_cg_A4(i,j,bi,bj,colx,coly)*
245     & Dv_SI(i+colx,j+coly,bi,bj)
246     ENDDO
247     ENDDO
248     ENDDO
249     ENDDO
250     ENDDO
251     ENDDO
252    
253     else
254    
255     CALL STREAMICE_CG_ACTION( myThid,
256     O Au_SI,
257     O Av_SI,
258     I Du_SI,
259     I Dv_SI,
260     I is,ie,js,je)
261    
262     ENDIF
263    
264     ! ! if (iter.eq.1) then
265     ! ! CALL WRITE_FLD_XY_RL ("Au2","",Au_SI,0,myThid)
266     ! ! CALL WRITE_FLD_XY_RL ("Av2","",Av_SI,0,myThid)
267     ! ! CALL WRITE_FLD_XY_RL ("Du","",Du_SI,0,myThid)
268     ! !
269     ! ! CALL WRITE_FLD_XY_RL("Dv","",Au_SI,0,myThid)
270     ! !
271     ! ! CALL WRITE_FLD_XY_RL("DiagU1","",Diagu_SI,0,myThid)
272     ! ! CALL WRITE_FLD_XY_RL("DiagV1","",Diagv_SI,0,myThid)
273     ! ! CALL WRITE_FLD_XY_RL("DiagU2","",
274     ! ! & streamice_cg_A1(i,j,bi,bj,0,0),0,myThid)
275     ! ! CALL WRITE_FLD_XY_RL("DiagV2","",
276     ! ! & streamice_cg_A4(i,j,bi,bj,0,0),0,myThid)
277     ! ! CALL WRITE_FLD_XY_RL("DiagV2","",Diagv_SI,0,myThid)
278     ! !
279     ! ! endif
280    
281    
282    
283    
284    
285     ! ! if (iter.eq.1) then
286     ! ! CALL WRITE_FLD_XY_RL ("Au1","",Au_SI,0,myThid)
287     ! ! CALL WRITE_FLD_XY_RL ("Av1","",Av_SI,0,myThid)
288     ! !
289     ! ! endif
290    
291    
292    
293     DO bj = myByLo(myThid), myByHi(myThid)
294     DO bi = myBxLo(myThid), myBxHi(myThid)
295     dot_p1_tile(bi,bj) = 0. _d 0
296     dot_p2_tile(bi,bj) = 0. _d 0
297     ENDDO
298     ENDDO
299    
300     DO bj = myByLo(myThid), myByHi(myThid)
301     DO bi = myBxLo(myThid), myBxHi(myThid)
302     DO j=1,sNy
303     DO i=1,sNx
304     IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN
305     dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zu_SI(i,j,bi,bj)*
306     & Ru_SI(i,j,bi,bj)
307     dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Du_SI(i,j,bi,bj)*
308     & Au_SI(i,j,bi,bj)
309     ENDIF
310     IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN
311     dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zv_SI(i,j,bi,bj)*
312     & Rv_SI(i,j,bi,bj)
313     dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Dv_SI(i,j,bi,bj)*
314     & Av_SI(i,j,bi,bj)
315     ENDIF
316     ENDDO
317     ENDDO
318     ENDDO
319     ENDDO
320    
321     CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
322     CALL GLOBAL_SUM_TILE_RL( dot_p2_tile, dot_p2, myThid )
323     alpha_k = dot_p1/dot_p2
324    
325     DO bj = myByLo(myThid), myByHi(myThid)
326     DO bi = myBxLo(myThid), myBxHi(myThid)
327     DO j=1-OLy,sNy+OLy
328     DO i=1-OLx,sNx+OLx
329    
330     IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN
331     cg_Uin(i,j,bi,bj)=cg_Uin(i,j,bi,bj)+
332     & alpha_k*Du_SI(i,j,bi,bj)
333     Ru_old_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj)
334     Zu_old_SI(i,j,bi,bj) = Zu_SI(i,j,bi,bj)
335     Ru_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj)-
336     & alpha_k*Au_SI(i,j,bi,bj)
337     Zu_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj) /
338     & DIAGu_SI(i,j,bi,bj)
339     ENDIF
340    
341     IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN
342     cg_Vin(i,j,bi,bj)=cg_Vin(i,j,bi,bj)+
343     & alpha_k*Dv_SI(i,j,bi,bj)
344     Rv_old_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj)
345     Zv_old_SI(i,j,bi,bj) = Zv_SI(i,j,bi,bj)
346     Rv_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj)-
347     & alpha_k*Av_SI(i,j,bi,bj)
348     Zv_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj) /
349     & DIAGv_SI(i,j,bi,bj)
350    
351     ENDIF
352     ENDDO
353     ENDDO
354     ENDDO
355     ENDDO
356    
357     DO bj = myByLo(myThid), myByHi(myThid)
358     DO bi = myBxLo(myThid), myBxHi(myThid)
359     dot_p1_tile(bi,bj) = 0. _d 0
360     dot_p2_tile(bi,bj) = 0. _d 0
361     ENDDO
362     ENDDO
363    
364     DO bj = myByLo(myThid), myByHi(myThid)
365     DO bi = myBxLo(myThid), myBxHi(myThid)
366     DO j=1,sNy
367     DO i=1,sNx
368    
369     IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN
370     dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zu_SI(i,j,bi,bj)*
371     & Ru_SI(i,j,bi,bj)
372     dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Zu_old_SI(i,j,bi,bj)*
373     & Ru_old_SI(i,j,bi,bj)
374     ENDIF
375    
376     IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN
377     dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zv_SI(i,j,bi,bj)*
378     & Rv_SI(i,j,bi,bj)
379     dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Zv_old_SI(i,j,bi,bj)*
380     & Rv_old_SI(i,j,bi,bj)
381     ENDIF
382    
383     ENDDO
384     ENDDO
385     ENDDO
386     ENDDO
387    
388     CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
389     CALL GLOBAL_SUM_TILE_RL( dot_p2_tile, dot_p2, myThid )
390    
391     beta_k = dot_p1/dot_p2
392    
393     DO bj = myByLo(myThid), myByHi(myThid)
394     DO bi = myBxLo(myThid), myBxHi(myThid)
395     DO j=1-OLy,sNy+OLy
396     DO i=1-OLx,sNx+OLx
397     IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
398     & Du_SI(i,j,bi,bj)=beta_k*Du_SI(i,j,bi,bj)+
399     & Zu_SI(i,j,bi,bj)
400     IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
401     & Dv_SI(i,j,bi,bj)=beta_k*Dv_SI(i,j,bi,bj)+
402     & Zv_SI(i,j,bi,bj)
403     ENDDO
404     ENDDO
405     ENDDO
406     ENDDO
407    
408     DO bj = myByLo(myThid), myByHi(myThid)
409     DO bi = myBxLo(myThid), myBxHi(myThid)
410     dot_p1_tile(bi,bj) = 0. _d 0
411     ENDDO
412     ENDDO
413    
414     DO bj = myByLo(myThid), myByHi(myThid)
415     DO bi = myBxLo(myThid), myBxHi(myThid)
416     DO j=1,sNy
417     DO i=1,sNx
418     IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
419     & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Ru_SI(i,j,bi,bj)**2
420     IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
421     & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Rv_SI(i,j,bi,bj)**2
422     ENDDO
423     ENDDO
424     ENDDO
425     ENDDO
426    
427     CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
428     resid = sqrt(dot_p1)
429    
430     ! IF (iter .eq. 1) then
431     ! print *, alpha_k, beta_k, resid
432     ! ENDIF
433    
434     cg_halo = cg_halo - 1
435    
436     if (cg_halo .eq. 0) then
437     cg_halo = min(OLx-1,OLy-1)
438     _EXCH_XY_RL( Du_SI, myThid )
439     _EXCH_XY_RL( Dv_SI, myThid )
440     _EXCH_XY_RL( Ru_SI, myThid )
441     _EXCH_XY_RL( Rv_SI, myThid )
442     _EXCH_XY_RL( cg_Uin, myThid )
443     _EXCH_XY_RL( cg_Vin, myThid )
444     endif
445    
446     endif
447     enddo ! end of CG loop
448    
449     c to avoid using "exit"
450     c if iters has reached max_iters there is no convergence
451    
452     IF (iters .lt. streamice_max_cg_iter) THEN
453     conv_flag = 1
454     ENDIF
455    
456     DO bj = myByLo(myThid), myByHi(myThid)
457     DO bi = myBxLo(myThid), myBxHi(myThid)
458     DO j=1-OLy,sNy+OLy
459     DO i=1-OLy,sNx+OLy
460     IF (STREAMICE_umask(i,j,bi,bj).eq.3.0)
461     & cg_Uin(i,j,bi,bj)=u_bdry_values_SI(i,j,bi,bj)
462     IF (STREAMICE_vmask(i,j,bi,bj).eq.3.0)
463     & cg_Vin(i,j,bi,bj)=v_bdry_values_SI(i,j,bi,bj)
464     ENDDO
465     ENDDO
466     ENDDO
467     ENDDO
468    
469     _EXCH_XY_RL( cg_Uin, myThid )
470     _EXCH_XY_RL( cg_Vin, myThid )
471    
472     #endif
473     RETURN
474     END
475    
476    
477    
478    
479    
480    
481    
482    

  ViewVC Help
Powered by ViewVC 1.1.22