/[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.4 - (show 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 C $Header: /u/gcmpack/MITgcm_contrib/dgoldberg/streamice/adstreamice_cg_solve.F,v 1.3 2012/07/26 16:13:18 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 CHARACTER*(MAX_LEN_MBUF) msgBuf
81
82 ! iters = streamice_max_cg_iter
83
84 #ifdef ALLOW_STREAMICE
85
86 WRITE(msgBuf,'(A)') 'CALLING MANUAL CG ADJOINT'
87 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
88 & SQUEEZE_RIGHT , 1)
89
90
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
112
113 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 _EXCH_XY_RL( cg_Bu, myThid )
129 _EXCH_XY_RL( cg_Bv, myThid )
130
131 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 _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 ! 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 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
656 #endif
657 RETURN
658 END
659
660
661
662
663
664
665
666

  ViewVC Help
Powered by ViewVC 1.1.22