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

Annotation of /MITgcm_contrib/dgoldberg/streamice/adstreamice_cg_solve.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.4 - (hide annotations) (download)
Thu Mar 7 15:23:19 2013 UTC (12 years, 4 months ago) by dgoldberg
Branch: MAIN
Changes since 1.3: +1 -2 lines
bug fixes, GL smoothing, changes for controlling bathym with const surf elev

1 dgoldberg 1.4 C $Header: /u/gcmpack/MITgcm_contrib/dgoldberg/streamice/adstreamice_cg_solve.F,v 1.3 2012/07/26 16:13:18 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 ADSTREAMICE_CG_SOLVE(
10     U U_state, ! velocities - need to be recalc'ed
11     I cg_Bu, ! adjoint of vel (input)
12     U V_state, ! velocities - need to be recalc'ed
13     I cg_Bv, ! adjoint of vel (input)
14     I Bu_state, ! to recalc velocities
15     U cg_Uin, ! adjoint of RHS (output)
16     I Bv_state, ! to recalc velocities
17     U cg_Vin, ! adjoint of RHS (output)
18     I A_uu, ! section of matrix that multiplies u and projects on u
19     U adA_uu, ! adjoint of matrix coeffs (output)
20     I A_uv, ! section of matrix that multiplies v and projects on u
21     U adA_uv, ! adjoint of matrix coeffs (output)
22     I A_vu, ! section of matrix that multiplies u and projects on v
23     U adA_vu, ! adjoint of matrix coeffs (output)
24     I A_vv, ! section of matrix that multiplies v and projects on v
25     U adA_vv, ! adjoint of matrix coeffs (output)
26     I tolerance,
27     I iters,
28     I myThid )
29     C /============================================================\
30     C | SUBROUTINE |
31     C | o |
32     C |============================================================|
33     C | |
34     C \============================================================/
35     IMPLICIT NONE
36    
37     C === Global variables ===
38     #include "SIZE.h"
39     #include "EEPARAMS.h"
40     #include "PARAMS.h"
41     #include "STREAMICE.h"
42     #include "STREAMICE_CG.h"
43    
44    
45     C !INPUT/OUTPUT ARGUMENTS
46     C cg_Uin, cg_Vin - input and output velocities
47     C cg_Bu, cg_Bv - driving stress
48     _RL cg_Uin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
49     _RL cg_Vin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
50     _RL cg_Bu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
51     _RL cg_Bv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
52     _RL U_state (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
53     _RL V_state (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
54     _RL Bu_state (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
55     _RL Bv_state (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
56     _RL
57     & A_uu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
58     & A_vu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
59     & A_uv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
60     & A_vv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
61     & adA_uu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
62     & adA_vu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
63     & adA_uv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
64     & adA_vv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1)
65     _RL tolerance
66     INTEGER iters
67     INTEGER myThid
68    
69     C LOCAL VARIABLES
70     INTEGER i, j, bi, bj, cg_halo, conv_flag, tmpiter
71     INTEGER iter, is, js, ie, je, colx, coly, k
72     _RL Utemp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
73     _RL Vtemp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
74     _RL UtempSt (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
75     _RL VtempSt (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
76     _RL dot_p1, dot_p2, alpha_k, beta_k, resid, resid_0
77     _RL dot_p1_tile (nSx,nSy)
78     _RL dot_p2_tile (nSx,nSy)
79     _RL ad_tolerance
80 dgoldberg 1.3 CHARACTER*(MAX_LEN_MBUF) msgBuf
81 dgoldberg 1.1
82     ! iters = streamice_max_cg_iter
83    
84     #ifdef ALLOW_STREAMICE
85    
86 dgoldberg 1.3 WRITE(msgBuf,'(A)') 'CALLING MANUAL CG ADJOINT'
87     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
88     & SQUEEZE_RIGHT , 1)
89    
90 dgoldberg 1.1
91     conv_flag = 0
92     ad_tolerance = 1.e-14
93    
94     DO bj = myByLo(myThid), myByHi(myThid)
95     DO bi = myBxLo(myThid), myBxHi(myThid)
96     DO j=1-Oly,sNy+Oly
97     DO i=1-Olx,sNx+Olx
98     Utemp (i,j,bi,bj) =
99     & cg_Uin (i,j,bi,bj)
100     Vtemp (i,j,bi,bj) =
101     & cg_Vin (i,j,bi,bj)
102     UtempSt (i,j,bi,bj) =
103     & U_state (i,j,bi,bj)
104     VtempSt (i,j,bi,bj) =
105     & V_state (i,j,bi,bj)
106     ENDDO
107     ENDDO
108     ENDDO
109     ENDDO
110    
111 dgoldberg 1.2
112    
113 dgoldberg 1.1 CALL STREAMICE_CG_SOLVE(
114     & U_state,
115     & V_state,
116     & Bu_state,
117     & Bv_state,
118     & A_uu,
119     & A_uv,
120     & A_vu,
121     & A_vv,
122     & tolerance,
123     & tmpiter,
124     & myThid )
125    
126     tmpiter = 0
127    
128 dgoldberg 1.2 _EXCH_XY_RL( cg_Bu, myThid )
129     _EXCH_XY_RL( cg_Bv, myThid )
130    
131 dgoldberg 1.1 CALL STREAMICE_CG_SOLVE(
132     & cg_Uin,
133     & cg_Vin,
134     & cg_Bu,
135     & cg_Bv,
136     & A_uu,
137     & A_uv,
138     & A_vu,
139     & A_vv,
140     & ad_tolerance,
141     & tmpiter,
142     & myThid )
143    
144 dgoldberg 1.2 _EXCH_XY_RL( cg_Uin, myThid )
145     _EXCH_XY_RL( cg_Vin, myThid )
146     _EXCH_XY_RL( U_state, myThid )
147     _EXCH_XY_RL( V_state, myThid )
148    
149 dgoldberg 1.1 ! DO bj = myByLo(myThid), myByHi(myThid)
150     ! DO bi = myBxLo(myThid), myBxHi(myThid)
151     ! DO j=1,sNy
152     ! DO i=1,sNx
153     ! Zu_SI (i,j,bi,bj) = 0. _d 0
154     ! Zv_SI (i,j,bi,bj) = 0. _d 0
155     ! Ru_SI (i,j,bi,bj) = 0. _d 0
156     ! Rv_SI (i,j,bi,bj) = 0. _d 0
157     ! Au_SI (i,j,bi,bj) = 0. _d 0
158     ! Av_SI (i,j,bi,bj) = 0. _d 0
159     ! Du_SI (i,j,bi,bj) = 0. _d 0
160     ! Dv_SI (i,j,bi,bj) = 0. _d 0
161     ! ENDDO
162     ! ENDDO
163     ! ENDDO
164     ! ENDDO
165     !
166     ! C FIND INITIAL RESIDUAL, and initialize r
167     !
168     ! ! #ifdef STREAMICE_CONSTRUCT_MATRIX
169     !
170     !
171     !
172     ! DO bj = myByLo(myThid), myByHi(myThid)
173     ! DO bi = myBxLo(myThid), myBxHi(myThid)
174     ! DO j=js,je
175     ! DO i=is,ie
176     ! DO colx=-1,1
177     ! DO coly=-1,1
178     ! Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) +
179     ! & A_uu(i,j,bi,bj,colx,coly)*
180     ! & cg_Uin(i+colx,j+coly,bi,bj)+
181     ! & A_uv(i,j,bi,bj,colx,coly)*
182     ! & cg_Vin(i+colx,j+coly,bi,bj)
183     ! Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +
184     ! & A_vu(i,j,bi,bj,colx,coly)*
185     ! & cg_Uin(i+colx,j+coly,bi,bj)+
186     ! & A_vv(i,j,bi,bj,colx,coly)*
187     ! & cg_Vin(i+colx,j+coly,bi,bj)
188     ! ENDDO
189     ! ENDDO
190     ! ENDDO
191     ! ENDDO
192     ! ENDDO
193     ! ENDDO
194     !
195     !
196     !
197     ! ! #else
198     ! !
199     ! ! CALL STREAMICE_CG_ACTION( myThid,
200     ! ! O Au_SI,
201     ! ! O Av_SI,
202     ! ! I cg_Uin,
203     ! ! I cg_Vin,
204     ! ! I 0, sNx+1, 0, sNy+1 )
205     ! !
206     ! ! #endif
207     !
208     ! _EXCH_XY_RL( Au_SI, myThid )
209     ! _EXCH_XY_RL( Av_SI, myThid )
210     !
211     ! DO bj = myByLo(myThid), myByHi(myThid)
212     ! DO bi = myBxLo(myThid), myBxHi(myThid)
213     ! DO j=1-OLy,sNy+OLy
214     ! DO i=1-OLx,sNx+OLx
215     ! Ru_SI(i,j,bi,bj)=cg_Bu(i,j,bi,bj)-
216     ! & Au_SI(i,j,bi,bj)
217     ! Rv_SI(i,j,bi,bj)=cg_Bv(i,j,bi,bj)-
218     ! & Av_SI(i,j,bi,bj)
219     ! ENDDO
220     ! ENDDO
221     ! dot_p1_tile(bi,bj) = 0. _d 0
222     ! dot_p2_tile(bi,bj) = 0. _d 0
223     ! ENDDO
224     ! ENDDO
225     !
226     ! DO bj = myByLo(myThid), myByHi(myThid)
227     ! DO bi = myBxLo(myThid), myBxHi(myThid)
228     ! DO j=1,sNy
229     ! DO i=1,sNx
230     ! IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
231     ! & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Ru_SI(i,j,bi,bj)**2
232     ! IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
233     ! & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Rv_SI(i,j,bi,bj)**2
234     ! ENDDO
235     ! ENDDO
236     ! ENDDO
237     ! ENDDO
238     !
239     ! CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
240     ! resid_0 = sqrt(dot_p1)
241     !
242     ! C CCCCCCCCCCCCCCCCCCCC
243     !
244     ! DO bj = myByLo(myThid), myByHi(myThid)
245     ! DO bi = myBxLo(myThid), myBxHi(myThid)
246     ! DO j=1-OLy,sNy+OLy
247     ! DO i=1-OLx,sNx+OLx
248     ! IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
249     ! & Zu_SI(i,j,bi,bj)=Ru_SI(i,j,bi,bj) / DIAGu_SI(i,j,bi,bj)
250     ! IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
251     ! & Zv_SI(i,j,bi,bj)=Rv_SI(i,j,bi,bj) / DIAGv_SI(i,j,bi,bj)
252     ! ENDDO
253     ! ENDDO
254     ! ENDDO
255     ! ENDDO
256     !
257     ! cg_halo = min(OLx-1,OLy-1)
258     ! conv_flag = 0
259     !
260     ! DO bj = myByLo(myThid), myByHi(myThid)
261     ! DO bi = myBxLo(myThid), myBxHi(myThid)
262     ! DO j=1-OLy,sNy+OLy
263     ! DO i=1-OLx,sNx+OLx
264     ! Du_SI(i,j,bi,bj)=Zu_SI(i,j,bi,bj)
265     ! Dv_SI(i,j,bi,bj)=Zv_SI(i,j,bi,bj)
266     ! ENDDO
267     ! ENDDO
268     ! ENDDO
269     ! ENDDO
270     !
271     ! resid = resid_0
272     ! iters = 0
273     !
274     ! c !!!!!!!!!!!!!!!!!!
275     ! c !! !!
276     ! c !! MAIN CG LOOP !!
277     ! c !! !!
278     ! c !!!!!!!!!!!!!!!!!!
279     !
280     !
281     !
282     !
283     ! c ! initially, b-grid data is valid up to 3 halo nodes out -- right? (check for MITgcm!!)
284     !
285     !
286     !
287     ! ! IF(STREAMICE_construct_matrix) CALL STREAMICE_CG_MAKE_A(myThid)
288     !
289     !
290     ! do iter = 1, streamice_max_cg_iter
291     ! if (resid .gt. tolerance*resid_0) then
292     !
293     ! c to avoid using "exit"
294     ! iters = iters + 1
295     !
296     ! is = 1 - cg_halo
297     ! ie = sNx + cg_halo
298     ! js = 1 - cg_halo
299     ! je = sNy + cg_halo
300     !
301     ! DO bj = myByLo(myThid), myByHi(myThid)
302     ! DO bi = myBxLo(myThid), myBxHi(myThid)
303     ! DO j=1-OLy,sNy+OLy
304     ! DO i=1-OLx,sNx+OLx
305     ! Au_SI(i,j,bi,bj) = 0. _d 0
306     ! Av_SI(i,j,bi,bj) = 0. _d 0
307     ! ENDDO
308     ! ENDDO
309     ! ENDDO
310     ! ENDDO
311     !
312     ! ! IF (STREAMICE_construct_matrix) THEN
313     !
314     ! ! #ifdef STREAMICE_CONSTRUCT_MATRIX
315     !
316     ! DO bj = myByLo(myThid), myByHi(myThid)
317     ! DO bi = myBxLo(myThid), myBxHi(myThid)
318     ! DO j=js,je
319     ! DO i=is,ie
320     ! DO colx=-1,1
321     ! DO coly=-1,1
322     ! Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) +
323     ! & A_uu(i,j,bi,bj,colx,coly)*
324     ! & Du_SI(i+colx,j+coly,bi,bj)+
325     ! & A_uv(i,j,bi,bj,colx,coly)*
326     ! & Dv_SI(i+colx,j+coly,bi,bj)
327     ! Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) +
328     ! & A_vu(i,j,bi,bj,colx,coly)*
329     ! & Du_SI(i+colx,j+coly,bi,bj)+
330     ! & A_vv(i,j,bi,bj,colx,coly)*
331     ! & Dv_SI(i+colx,j+coly,bi,bj)
332     ! ENDDO
333     ! ENDDO
334     ! ENDDO
335     ! ENDDO
336     ! ENDDO
337     ! ENDDO
338     !
339     ! ! else
340     ! ! #else
341     ! !
342     ! ! CALL STREAMICE_CG_ACTION( myThid,
343     ! ! O Au_SI,
344     ! ! O Av_SI,
345     ! ! I Du_SI,
346     ! ! I Dv_SI,
347     ! ! I is,ie,js,je)
348     ! !
349     ! ! ! ENDIF
350     ! !
351     ! ! #endif
352     !
353     !
354     ! DO bj = myByLo(myThid), myByHi(myThid)
355     ! DO bi = myBxLo(myThid), myBxHi(myThid)
356     ! dot_p1_tile(bi,bj) = 0. _d 0
357     ! dot_p2_tile(bi,bj) = 0. _d 0
358     ! ENDDO
359     ! ENDDO
360     !
361     ! DO bj = myByLo(myThid), myByHi(myThid)
362     ! DO bi = myBxLo(myThid), myBxHi(myThid)
363     ! DO j=1,sNy
364     ! DO i=1,sNx
365     ! IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN
366     ! dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zu_SI(i,j,bi,bj)*
367     ! & Ru_SI(i,j,bi,bj)
368     ! dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Du_SI(i,j,bi,bj)*
369     ! & Au_SI(i,j,bi,bj)
370     ! ENDIF
371     ! IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN
372     ! dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zv_SI(i,j,bi,bj)*
373     ! & Rv_SI(i,j,bi,bj)
374     ! dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Dv_SI(i,j,bi,bj)*
375     ! & Av_SI(i,j,bi,bj)
376     ! ENDIF
377     ! ENDDO
378     ! ENDDO
379     ! ENDDO
380     ! ENDDO
381     !
382     ! CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
383     ! CALL GLOBAL_SUM_TILE_RL( dot_p2_tile, dot_p2, myThid )
384     ! alpha_k = dot_p1/dot_p2
385     !
386     ! DO bj = myByLo(myThid), myByHi(myThid)
387     ! DO bi = myBxLo(myThid), myBxHi(myThid)
388     ! DO j=1-OLy,sNy+OLy
389     ! DO i=1-OLx,sNx+OLx
390     !
391     ! IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN
392     ! cg_Uin(i,j,bi,bj)=cg_Uin(i,j,bi,bj)+
393     ! & alpha_k*Du_SI(i,j,bi,bj)
394     ! Ru_old_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj)
395     ! Zu_old_SI(i,j,bi,bj) = Zu_SI(i,j,bi,bj)
396     ! Ru_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj)-
397     ! & alpha_k*Au_SI(i,j,bi,bj)
398     ! Zu_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj) /
399     ! & DIAGu_SI(i,j,bi,bj)
400     ! ENDIF
401     !
402     ! IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN
403     ! cg_Vin(i,j,bi,bj)=cg_Vin(i,j,bi,bj)+
404     ! & alpha_k*Dv_SI(i,j,bi,bj)
405     ! Rv_old_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj)
406     ! Zv_old_SI(i,j,bi,bj) = Zv_SI(i,j,bi,bj)
407     ! Rv_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj)-
408     ! & alpha_k*Av_SI(i,j,bi,bj)
409     ! Zv_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj) /
410     ! & DIAGv_SI(i,j,bi,bj)
411     !
412     ! ENDIF
413     ! ENDDO
414     ! ENDDO
415     ! ENDDO
416     ! ENDDO
417     !
418     ! DO bj = myByLo(myThid), myByHi(myThid)
419     ! DO bi = myBxLo(myThid), myBxHi(myThid)
420     ! dot_p1_tile(bi,bj) = 0. _d 0
421     ! dot_p2_tile(bi,bj) = 0. _d 0
422     ! ENDDO
423     ! ENDDO
424     !
425     ! DO bj = myByLo(myThid), myByHi(myThid)
426     ! DO bi = myBxLo(myThid), myBxHi(myThid)
427     ! DO j=1,sNy
428     ! DO i=1,sNx
429     !
430     ! IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN
431     ! dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zu_SI(i,j,bi,bj)*
432     ! & Ru_SI(i,j,bi,bj)
433     ! dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Zu_old_SI(i,j,bi,bj)*
434     ! & Ru_old_SI(i,j,bi,bj)
435     ! ENDIF
436     !
437     ! IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN
438     ! dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zv_SI(i,j,bi,bj)*
439     ! & Rv_SI(i,j,bi,bj)
440     ! dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Zv_old_SI(i,j,bi,bj)*
441     ! & Rv_old_SI(i,j,bi,bj)
442     ! ENDIF
443     !
444     ! ENDDO
445     ! ENDDO
446     ! ENDDO
447     ! ENDDO
448     !
449     ! CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
450     ! CALL GLOBAL_SUM_TILE_RL( dot_p2_tile, dot_p2, myThid )
451     !
452     ! beta_k = dot_p1/dot_p2
453     !
454     ! DO bj = myByLo(myThid), myByHi(myThid)
455     ! DO bi = myBxLo(myThid), myBxHi(myThid)
456     ! DO j=1-OLy,sNy+OLy
457     ! DO i=1-OLx,sNx+OLx
458     ! IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
459     ! & Du_SI(i,j,bi,bj)=beta_k*Du_SI(i,j,bi,bj)+
460     ! & Zu_SI(i,j,bi,bj)
461     ! IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
462     ! & Dv_SI(i,j,bi,bj)=beta_k*Dv_SI(i,j,bi,bj)+
463     ! & Zv_SI(i,j,bi,bj)
464     ! ENDDO
465     ! ENDDO
466     ! ENDDO
467     ! ENDDO
468     !
469     ! DO bj = myByLo(myThid), myByHi(myThid)
470     ! DO bi = myBxLo(myThid), myBxHi(myThid)
471     ! dot_p1_tile(bi,bj) = 0. _d 0
472     ! ENDDO
473     ! ENDDO
474     !
475     ! DO bj = myByLo(myThid), myByHi(myThid)
476     ! DO bi = myBxLo(myThid), myBxHi(myThid)
477     ! DO j=1,sNy
478     ! DO i=1,sNx
479     ! IF (STREAMICE_umask(i,j,bi,bj).eq.1.0)
480     ! & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Ru_SI(i,j,bi,bj)**2
481     ! IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0)
482     ! & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Rv_SI(i,j,bi,bj)**2
483     ! ENDDO
484     ! ENDDO
485     ! ENDDO
486     ! ENDDO
487     !
488     ! CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid )
489     ! resid = sqrt(dot_p1)
490     !
491     ! ! IF (iter .eq. 1) then
492     ! ! print *, alpha_k, beta_k, resid
493     ! ! ENDIF
494     !
495     ! cg_halo = cg_halo - 1
496     !
497     ! if (cg_halo .eq. 0) then
498     ! cg_halo = min(OLx-1,OLy-1)
499     ! _EXCH_XY_RL( Du_SI, myThid )
500     ! _EXCH_XY_RL( Dv_SI, myThid )
501     ! _EXCH_XY_RL( Ru_SI, myThid )
502     ! _EXCH_XY_RL( Rv_SI, myThid )
503     ! _EXCH_XY_RL( cg_Uin, myThid )
504     ! _EXCH_XY_RL( cg_Vin, myThid )
505     ! endif
506     !
507     ! endif
508     ! enddo ! end of CG loop
509    
510    
511    
512     c to avoid using "exit"
513     c if iters has reached max_iters there is no convergence
514    
515     ! IF (iters .lt. streamice_max_cg_iter) THEN
516     ! conv_flag = 1
517     ! ENDIF
518    
519     DO bj = myByLo(myThid), myByHi(myThid)
520     DO bi = myBxLo(myThid), myBxHi(myThid)
521     DO j=1-Oly,sNy+Oly
522     DO i=1-Olx,sNx+Olx
523     DO colx=-1,1
524     DO coly=-1,1
525    
526     if (STREAMICE_umask(i,j,bi,bj).eq.1) then
527     if (STREAMICE_umask(i+colx,j+coly,bi,bj).eq.1) then
528     adA_uu(i,j,bi,bj,colx,coly) =
529     & adA_uu(i,j,bi,bj,colx,coly) -
530     & cg_Uin(i,j,bi,bj) *
531     & U_state(i+colx,j+coly,bi,bj)
532     ! print *,"ADA",cg_Uin(i,j,bi,bj)
533     endif
534     if (STREAMICE_vmask(i+colx,j+coly,bi,bj).eq.1) then
535     adA_uv(i,j,bi,bj,colx,coly) =
536     & adA_uv(i,j,bi,bj,colx,coly) -
537     & cg_Uin(i,j,bi,bj) *
538     & V_state(i+colx,j+coly,bi,bj)
539     endif
540     endif
541    
542     if (STREAMICE_vmask(i,j,bi,bj).eq.1) then
543     if (STREAMICE_umask(i+colx,j+coly,bi,bj).eq.1) then
544     adA_vu(i,j,bi,bj,colx,coly) =
545     & adA_vu(i,j,bi,bj,colx,coly) -
546     & cg_Vin(i,j,bi,bj) *
547     & U_state(i+colx,j+coly,bi,bj)
548     endif
549     if (STREAMICE_vmask(i+colx,j+coly,bi,bj).eq.1) then
550     adA_vv(i,j,bi,bj,colx,coly) =
551     & adA_vv(i,j,bi,bj,colx,coly) -
552     & cg_Vin(i,j,bi,bj) *
553     & V_state(i+colx,j+coly,bi,bj)
554     endif
555     endif
556    
557     enddo
558     enddo
559     ! if (i.eq.1.and.j.eq.1) then
560     ! print *, "adA_uu", adA_uu(i,j,bi,bj,-1,-1),
561     ! & adA_uu(i,j,bi,bj,-1,0),
562     ! & adA_uu(i,j,bi,bj,-1,1),
563     ! & cg_Uin(i,j,bi,bj)
564     ! endif
565     enddo
566     enddo
567     enddo
568     enddo
569    
570     ! print *,"MATRIX 1"
571     ! do i=1,sNx
572     ! do j=1,sNy
573     ! print *, ada_uu(i,j,1,1,-1,-1), ada_uu(i,j,1,1,0,-1),
574     ! & ada_uu(i,j,1,1,1,-1), ada_uu(i,j,1,1,-1,0),
575     ! & ada_uu(i,j,1,1,0,0), ada_uu(i,j,1,1,1,0),
576     ! & ada_uu(i,j,1,1,-1,1), ada_uu(i,j,1,1,0,1),ada_uu(i,j,1,1,1,1)
577     ! enddo
578     ! enddo
579     !
580     ! print *,"MATRIX 2"
581     ! do i=1,sNx
582     ! do j=1,sNy
583     ! print *, ada_uv(i,j,1,1,-1,-1), ada_uv(i,j,1,1,0,-1),
584     ! & ada_uv(i,j,1,1,1,-1), ada_uv(i,j,1,1,-1,0),
585     ! & ada_uv(i,j,1,1,0,0), ada_uv(i,j,1,1,1,0),
586     ! & ada_uv(i,j,1,1,-1,1), ada_uv(i,j,1,1,0,1),ada_uv(i,j,1,1,1,1)
587     ! enddo
588     ! enddo
589     !
590     ! print *,"MATRIX 3"
591     ! do i=1,sNx
592     ! do j=1,sNy
593     ! print *, ada_vu(i,j,1,1,-1,-1), ada_vu(i,j,1,1,0,-1),
594     ! & ada_vu(i,j,1,1,1,-1), ada_vu(i,j,1,1,-1,0),
595     ! & ada_vu(i,j,1,1,0,0), ada_vu(i,j,1,1,1,0),
596     ! & ada_vu(i,j,1,1,-1,1), ada_vu(i,j,1,1,0,1),ada_vu(i,j,1,1,1,1)
597     ! enddo
598     ! enddo
599     !
600     ! print *,"MATRIX 4"
601     ! do i=1,sNx
602     ! do j=1,sNy
603     ! print *, ada_vv(i,j,1,1,-1,-1), ada_vv(i,j,1,1,0,-1),
604     ! & ada_vv(i,j,1,1,1,-1), ada_vv(i,j,1,1,-1,0),
605     ! & ada_vv(i,j,1,1,0,0), ada_vv(i,j,1,1,1,0),
606     ! & ada_vv(i,j,1,1,-1,1), ada_vv(i,j,1,1,0,1),ada_vv(i,j,1,1,1,1)
607     ! enddo
608     ! enddo
609    
610    
611    
612     DO bj = myByLo(myThid), myByHi(myThid)
613     DO bi = myBxLo(myThid), myBxHi(myThid)
614     DO j=1-Oly,sNy+Oly
615     DO i=1-Olx,sNx+Olx
616     if (i.lt.1.or.i.gt.sNx.or.
617     & j.lt.1.or.j.gt.sNy) then
618     cg_Uin (i,j,bi,bj) = 0.0
619     cg_Vin (i,j,bi,bj) = 0.0
620    
621     DO colx=-1,1
622     DO coly=-1,1
623     ada_uu(i,j,bi,bj,colx,coly)=0.0
624     ada_uv(i,j,bi,bj,colx,coly)=0.0
625     ada_vu(i,j,bi,bj,colx,coly)=0.0
626     ada_vv(i,j,bi,bj,colx,coly)=0.0
627     enddo
628     enddo
629    
630    
631     endif
632     cg_Uin (i,j,bi,bj) =
633     & cg_Uin (i,j,bi,bj) +
634     & Utemp (i,j,bi,bj)
635     cg_Vin (i,j,bi,bj) =
636     & cg_Vin (i,j,bi,bj) +
637     & Vtemp (i,j,bi,bj)
638     cg_bu (i,j,bi,bj) = 0.
639     cg_bv (i,j,bi,bj) = 0.
640     U_state (i,j,bi,bj) =
641     & UtempSt (i,j,bi,bj)
642     V_state (i,j,bi,bj) =
643     & VtempSt (i,j,bi,bj)
644     ENDDO
645     ENDDO
646     ENDDO
647     ENDDO
648    
649    
650 dgoldberg 1.3 WRITE(msgBuf,'(A,I5,A)') 'DONE WITH MANUAL CG ADJOINT:',tmpiter,
651     & 'ITERS'
652     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
653     & SQUEEZE_RIGHT , 1)
654    
655 dgoldberg 1.1
656     #endif
657     RETURN
658     END
659    
660    
661    
662    
663    
664    
665    
666    

  ViewVC Help
Powered by ViewVC 1.1.22