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

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

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


Revision 1.2 - (show annotations) (download)
Wed Jul 25 23:11:27 2012 UTC (13 years ago) by dgoldberg
Branch: MAIN
Changes since 1.1: +11 -2 lines
added overlap update for parallel mode

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

  ViewVC Help
Powered by ViewVC 1.1.22