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

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

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


Revision 1.11 - (show annotations) (download)
Wed Jan 9 21:56:18 2013 UTC (12 years, 6 months ago) by dgoldberg
Branch: MAIN
Changes since 1.10: +29 -7 lines
changes to accept real datasets

1 C $Header: /u/gcmpack/MITgcm_contrib/dgoldberg/streamice/streamice_vel_solve.F,v 1.10 2012/09/27 20:29:01 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 STREAMICE_VEL_SOLVE( myThid )
10 C /============================================================\
11 C | SUBROUTINE |
12 C | o |
13 C |============================================================|
14 C | |
15 C \============================================================/
16 IMPLICIT NONE
17
18 C === Global variables ===
19 #include "SIZE.h"
20 #include "EEPARAMS.h"
21 #include "PARAMS.h"
22 #include "STREAMICE.h"
23 #include "STREAMICE_CG.h"
24 #ifdef ALLOW_AUTODIFF_TAMC
25 # include "tamc.h"
26 #endif
27
28 C !INPUT/OUTPUT ARGUMENTS
29 INTEGER myThid
30
31 #ifdef ALLOW_STREAMICE
32
33 C LOCAL VARIABLES
34
35 ! real, dimension(:,:), pointer :: TAUDX, TAUDY, u_prev_iterate, v_prev_iterate, &
36 ! u_bdry_cont, v_bdry_cont, Au, Av, err_u, err_v, &
37 ! geolonq, geolatq, u_last, v_last, float_cond, H_node
38 ! type(ocean_grid_type), pointer :: G
39 ! integer :: conv_flag, i, j, k,l, iter, isym, &
40 ! isdq, iedq, jsdq, jedq, isd, ied, jsd, jed, isumstart, jsumstart, nodefloat, nsub
41 ! real :: err_max, err_tempu, err_tempv, err_init, area, max_vel, tempu, tempv, rhoi, rhow
42
43 INTEGER conv_flag, i, j, k, l, iter, cg_iters, bi, bj
44 INTEGER iter_numconv
45 INTEGER ikey_nl
46 _RL err_max, err_tempu, err_tempv, err_init, area, err_max_fp
47 _RL max_vel, tempu, tempv, err_lastchange, cgtol
48 CHARACTER*(MAX_LEN_MBUF) msgBuf
49 ! _RL taudx (1-OLx:sNx+OLx,1-OLy:sNy+Oly,nSx,nSy)
50 ! _RL taudy (1-OLx:sNx+OLx,1-OLy:sNy+Oly,nSx,nSy)
51
52 IF (STREAMICE_ppm_driving_stress) THEN
53 CALL STREAMICE_DRIVING_STRESS_PPM (myThid)
54 ELSE
55 CALL STREAMICE_DRIVING_STRESS (myThid)
56 ENDIF
57
58 cgtol = streamice_cg_tol
59
60 ! CALL WRITE_FULLARRAY_RL ("taudy_SI",taudy_SI,1,0,0,1,0,myThid)
61
62 _EXCH_XY_RL ( taudx_SI , myThid )
63 _EXCH_XY_RL ( taudy_SI , myThid )
64
65 ! CALL WRITE_FULLARRAY_RL ("taudy_SI_2",taudy_SI,1,0,0,1,0,myThid)
66
67 DO bj = myByLo(myThid), myByHi(myThid)
68 DO bi = myBxLo(myThid), myBxHi(myThid)
69 DO j=1-OLy,sNy+OLy
70 DO i=1-OLx,sNx+OLx
71 u_old_SI (i,j,bi,bj) = U_streamice (i,j,bi,bj)
72 v_old_SI (i,j,bi,bj) = V_streamice (i,j,bi,bj)
73 ENDDO
74 ENDDO
75 ENDDO
76 ENDDO
77
78 #ifdef ALLOW_AUTODIFF_TAMC
79 !$TAF STORE U_streamice = comlev1, key=ikey_dynamics
80 !$TAF STORE V_streamice = comlev1, key=ikey_dynamics
81 #ifdef STREAMICE_HYBRID_STRESS
82 !$TAF STORE streamice_taubx = comlev1, key=ikey_dynamics
83 !$TAF STORE streamice_tauby = comlev1, key=ikey_dynamics
84 !$TAF STORE visc_streamice_full = comlev1, key=ikey_dynamics
85 #endif
86 #endif
87
88 #ifdef STREAMICE_HYBRID_STRESS
89 CALL STREAMICE_VISC_BETA_HYBRID ( myThid )
90 #else
91 CALL STREAMICE_VISC_BETA ( myThid )
92 #endif
93
94 #ifdef STREAMICE_HYBRID_STRESS
95 !$TAF STORE visc_streamice_full = comlev1, key=ikey_dynamics
96 #endif
97
98 _EXCH_XY_RL ( tau_beta_eff_streamice , myThid )
99 _EXCH_XY_RL ( visc_streamice , myThid )
100
101 DO bj = myByLo(myThid), myByHi(myThid)
102 DO bi = myBxLo(myThid), myBxHi(myThid)
103 DO j=1,sNy
104 DO i=1,sNx
105 Au_SI (i,j,bi,bj) = 0. _d 0
106 Av_SI (i,j,bi,bj) = 0. _d 0
107 ubd_SI (i,j,bi,bj) = 0. _d 0
108 vbd_SI (i,j,bi,bj) = 0. _d 0
109 ENDDO
110 ENDDO
111 ENDDO
112 ENDDO
113
114
115 CALL STREAMICE_CG_BOUND_VALS( myThid,
116 O ubd_SI,
117 O vbd_SI)
118
119 ! CALL WRITE_FLD_XY_RL("u_bound_cont","",ubd_SI,0,myThid)
120 ! CALL WRITE_FLD_XY_RL("v_bound_cont","",vbd_SI,0,myThid)
121 ! CALL WRITE_FLD_XY_RL("taudx_u","",taudx_SI,0,myThid)
122 ! CALL WRITE_FLD_XY_RL("taudx_v","",taudy_SI,0,myThid)
123
124 #ifdef ALLOW_AUTODIFF_TAMC
125 !$TAF STORE U_streamice = comlev1, key=ikey_dynamics
126 !$TAF STORE V_streamice = comlev1, key=ikey_dynamics
127 #endif
128
129 CALL STREAMICE_CG_ACTION( myThid,
130 O Au_SI,
131 O Av_SI,
132 I U_streamice,
133 I V_streamice,
134 I 0, sNx+1, 0, sNy+1 )
135
136
137
138 err_init = 0. _d 0
139
140 DO bj = myByLo(myThid), myByHi(myThid)
141 DO bi = myBxLo(myThid), myBxHi(myThid)
142 #ifdef ALLOW_AUTODIFF_TAMC
143 !$TAF STORE err_init = comlev1, key=ikey_dynamics
144 #endif
145 DO j=1,sNy
146 DO i=1,sNx
147 err_tempu = 0. _d 0
148 err_tempv = 0. _d 0
149 IF (STREAMICE_umask(i,j,bi,bj).eq.1) THEN
150 err_tempu =
151 & ABS (Au_SI(i,j,bi,bj)+0*ubd_SI(i,j,bi,bj) -
152 & taudx_SI(i,j,bi,bj))
153 ! print *, "err_temp_u", err_tempu
154 ENDIF
155 IF (STREAMICE_vmask(i,j,bi,bj).eq.1) THEN
156 err_tempv = MAX( err_tempu,
157 & ABS (Av_SI(i,j,bi,bj)+0*vbd_SI(i,j,bi,bj) -
158 & taudy_SI(i,j,bi,bj)))
159 ENDIF
160 IF (err_tempv .ge. err_init) THEN
161 err_init = err_tempv
162 ENDIF
163 ENDDO
164 ENDDO
165 ENDDO
166 ENDDO
167 #ifdef ALLOW_AUTODIFF_TAMC
168 !$TAF STORE err_init = comlev1, key=ikey_dynamics
169 #endif
170
171 CALL GLOBAL_MAX_R8 (err_init, myThid)
172
173 WRITE(msgBuf,'(A,E15.7)') 'initial nonlinear resid (error): ',
174 & err_init
175 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
176 & SQUEEZE_RIGHT , 1)
177
178
179 iter_numconv = 0
180 err_max = err_init
181 err_max_fp = streamice_nonlin_tol_fp * 10.
182 err_lastchange = err_init
183
184 C START NL ITER. LOOP
185 C-------------------------------------------------------------------
186
187 DO iter=1,streamice_max_nl_iter
188
189 C To avoid using "exit", loop goes through all iterations
190 C but after convergence loop does nothing
191
192 #ifdef ALLOW_AUTODIFF_TAMC
193 ikey_nl = (ikey_dynamics-1)*streamice_max_nl + iter
194 #endif
195 #ifdef ALLOW_AUTODIFF_TAMC
196 !$TAF STORE cgtol = comlev1_stream_nl, key=ikey_nl
197 !$TAF STORE err_max = comlev1_stream_nl, key=ikey_nl
198 !$TAF STORE err_max_fp = comlev1_stream_nl, key=ikey_nl
199 !$TAF STORE err_tempu = comlev1_stream_nl, key=ikey_nl
200 !$TAF STORE err_tempv = comlev1_stream_nl, key=ikey_nl
201 !$TAF STORE err_lastchange = comlev1_stream_nl, key=ikey_nl
202 !$TAF STORE ru_old_si = comlev1_stream_nl, key=ikey_nl
203 !$TAF STORE rv_old_si = comlev1_stream_nl, key=ikey_nl
204 !$TAF STORE u_old_si = comlev1_stream_nl, key=ikey_nl
205 !$TAF STORE v_old_si = comlev1_stream_nl, key=ikey_nl
206 !$TAF STORE streamice_cg_a1 = comlev1_stream_nl, key=ikey_nl
207 !$TAF STORE streamice_cg_a2 = comlev1_stream_nl, key=ikey_nl
208 !$TAF STORE streamice_cg_a3 = comlev1_stream_nl, key=ikey_nl
209 !$TAF STORE streamice_cg_a4 = comlev1_stream_nl, key=ikey_nl
210 !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
211 !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
212 !$TAF STORE tau_beta_eff_streamice = comlev1_stream_nl, key=ikey_nl
213 !$TAF STORE visc_streamice = comlev1_stream_nl, key=ikey_nl
214 !$TAF STORE zu_old_si = comlev1_stream_nl, key=ikey_nl
215 !$TAF STORE zv_old_si = comlev1_stream_nl, key=ikey_nl
216 #endif
217
218 IF ((err_max .GT. streamice_nonlin_tol * err_init) .and.
219 & (err_max_fp .GT. streamice_nonlin_tol_fp)) THEN
220
221 iter_numconv = iter_numconv + 1
222
223 #ifdef ALLOW_AUTODIFF_TAMC
224 !$TAF STORE cgtol = comlev1_stream_nl, key=ikey_nl
225 !$TAF STORE err_max = comlev1_stream_nl, key=ikey_nl
226 !$TAF STORE err_tempu = comlev1_stream_nl, key=ikey_nl
227 !$TAF STORE err_tempv = comlev1_stream_nl, key=ikey_nl
228 !$TAF STORE err_lastchange = comlev1_stream_nl, key=ikey_nl
229 !$TAF STORE ru_old_si = comlev1_stream_nl, key=ikey_nl
230 !$TAF STORE rv_old_si = comlev1_stream_nl, key=ikey_nl
231 !$TAF STORE streamice_cg_a1 = comlev1_stream_nl, key=ikey_nl
232 !$TAF STORE streamice_cg_a2 = comlev1_stream_nl, key=ikey_nl
233 !$TAF STORE streamice_cg_a3 = comlev1_stream_nl, key=ikey_nl
234 !$TAF STORE streamice_cg_a4 = comlev1_stream_nl, key=ikey_nl
235 !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
236 !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
237 !$TAF STORE tau_beta_eff_streamice = comlev1_stream_nl, key=ikey_nl
238 !$TAF STORE visc_streamice = comlev1_stream_nl, key=ikey_nl
239 !$TAF STORE zu_old_si = comlev1_stream_nl, key=ikey_nl
240 !$TAF STORE zv_old_si = comlev1_stream_nl, key=ikey_nl
241 #endif
242
243 #ifdef ALLOW_AUTODIFF_TAMC
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 ! U_streamice (i,j,bi,bj) = 0. _d 0
249 ! V_streamice (i,j,bi,bj) = 0. _d 0
250 ! ENDDO
251 ! ENDDO
252 ! ENDDO
253 ! ENDDO
254 #endif
255
256 CALL STREAMICE_CG_WRAPPER(
257 & U_streamice,
258 & V_streamice,
259 & taudx_SI,
260 & taudy_SI,
261 & cgtol,
262 & cg_iters,
263 & myThid )
264
265 #ifdef STREAMICE_HYBRID_STRESS
266 #ifdef ALLOW_AUTODIFF_TAMC
267 !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
268 !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
269 #endif
270 #endif
271
272 #ifdef STREAMICE_HYBRID_STRESS
273 CALL STREAMICE_TAUB (myThid)
274 #endif
275
276 WRITE(msgBuf,'(A,I5,A,I4,A)') 'streamice linear solve number',
277 & iter, " ",
278 & cg_iters,
279 & ' iterations '
280 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
281 & SQUEEZE_RIGHT , 1)
282
283 #ifdef ALLOW_AUTODIFF_TAMC
284 !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
285 !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
286 #ifdef STREAMICE_HYBRID_STRESS
287 !$TAF STORE streamice_taubx = comlev1_stream_nl, key=ikey_nl
288 !$TAF STORE streamice_tauby = comlev1_stream_nl, key=ikey_nl
289 !$TAF STORE visc_streamice_full = comlev1_stream_nl, key=ikey_nl
290 #endif
291 #endif
292
293 #ifdef STREAMICE_HYBRID_STRESS
294 CALL STREAMICE_VISC_BETA_HYBRID ( myThid )
295 #else
296 CALL STREAMICE_VISC_BETA ( myThid )
297 #endif
298
299
300 #ifdef ALLOW_AUTODIFF_TAMC
301 !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
302 !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
303 #endif
304
305 _EXCH_XY_RL ( tau_beta_eff_streamice , myThid )
306 _EXCH_XY_RL ( visc_streamice , myThid )
307
308 DO bj = myByLo(myThid), myByHi(myThid)
309 DO bi = myBxLo(myThid), myBxHi(myThid)
310 DO j=1,sNy
311 DO i=1,sNx
312 Au_SI (i,j,bi,bj) = 0. _d 0
313 Av_SI (i,j,bi,bj) = 0. _d 0
314 ubd_SI (i,j,bi,bj) = 0. _d 0
315 vbd_SI (i,j,bi,bj) = 0. _d 0
316 ENDDO
317 ENDDO
318 ENDDO
319 ENDDO
320
321 CALL STREAMICE_CG_BOUND_VALS( myThid,
322 O ubd_SI,
323 O vbd_SI)
324
325 #ifdef ALLOW_AUTODIFF_TAMC
326 !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
327 !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
328 #endif
329
330 CALL STREAMICE_CG_ACTION( myThid,
331 O Au_SI,
332 O Av_SI,
333 I U_streamice,
334 I V_streamice,
335 I 0, sNx+1, 0, sNy+1 )
336
337 if (iter .eq. streamice_max_nl_iter) then
338 CALL WRITE_FLD_XY_RL("u_bound_cont_A","",ubd_SI,0,myThid)
339 CALL WRITE_FLD_XY_RL("v_bound_cont_A","",vbd_SI,0,myThid)
340 CALL WRITE_FLD_XY_RL("u_bound_cont_B","",Au_SI,0,myThid)
341 CALL WRITE_FLD_XY_RL("v_bound_cont_B","",Av_SI,0,myThid)
342 endif
343
344 #ifdef ALLOW_AUTODIFF_TAMC
345 !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
346 !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
347 #endif
348
349 err_max = 0. _d 0
350 err_max_fp = 0. _d 0
351
352 #ifdef ALLOW_AUTODIFF_TAMC
353 !$TAF STORE err_max = comlev1_stream_nl, key=ikey_nl
354 #endif
355 DO bj = myByLo(myThid), myByHi(myThid)
356 DO bi = myBxLo(myThid), myBxHi(myThid)
357 #ifdef ALLOW_AUTODIFF_TAMC
358 !$TAF STORE err_max = comlev1_stream_nl, key=ikey_nl
359 #endif
360 DO j=1,sNy
361 DO i=1,sNx
362 err_tempu = 0. _d 0
363 err_tempv = 0. _d 0
364 IF (STREAMICE_umask(i,j,bi,bj).eq.1) THEN
365 err_tempu =
366 & ABS (Au_SI(i,j,bi,bj)+0*ubd_SI(i,j,bi,bj) -
367 & taudx_SI(i,j,bi,bj))
368 ENDIF
369 IF (STREAMICE_vmask(i,j,bi,bj).eq.1) THEN
370 err_tempv = MAX( err_tempu,
371 & ABS (Av_SI(i,j,bi,bj)+0*vbd_SI(i,j,bi,bj) -
372 & taudy_SI(i,j,bi,bj)))
373 ENDIF
374 ! if (err_tempu.ge.1.e2.or.err_tempv.ge.1.e2) THEN
375 ! print *, "FOUND MAX ", i,j,err_tempu,err_tempv,
376 ! & ubd_SI(i,j,bi,bj),vbd_SI(i,j,bi,bj)
377 ! endif
378 IF (err_tempv .ge. err_max) THEN
379 err_max = err_tempv
380 ENDIF
381 ENDDO
382 ENDDO
383 ENDDO
384 ENDDO
385
386 CALL GLOBAL_MAX_R8 (err_max, myThid)
387
388 DO bj = myByLo(myThid), myByHi(myThid)
389 DO bi = myBxLo(myThid), myBxHi(myThid)
390 DO j=1,sNy
391 DO i=1,sNx
392 err_tempu = 0. _d 0
393 err_tempv = 0. _d 0
394 IF (STREAMICE_umask(i,j,bi,bj).eq.1) THEN
395 err_tempu =
396 & ABS (U_streamice(i,j,bi,bj)-u_old_SI(i,j,bi,bj))
397 ENDIF
398 IF (STREAMICE_vmask(i,j,bi,bj).eq.1) THEN
399 err_tempv = MAX( err_tempu,
400 & ABS (V_streamice(i,j,bi,bj)-v_old_SI(i,j,bi,bj)))
401 ENDIF
402 IF (err_tempv .ge. err_max_fp) err_max_fp = err_tempv
403 ENDDO
404 ENDDO
405 ENDDO
406 ENDDO
407
408 CALL GLOBAL_MAX_R8 (err_max_fp, myThid)
409 WRITE(msgBuf,'(A,E15.7)') '||x_i-x_{i-1}||_inf',
410 & err_max_fp
411 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
412 & SQUEEZE_RIGHT , 1)
413
414 DO bj = myByLo(myThid), myByHi(myThid)
415 DO bi = myBxLo(myThid), myBxHi(myThid)
416 DO j=1-OLy,sNy+OLy
417 DO i=1-OLx,sNx+OLx
418 u_old_SI (i,j,bi,bj) = U_streamice (i,j,bi,bj)
419 v_old_SI (i,j,bi,bj) = V_streamice (i,j,bi,bj)
420 ENDDO
421 ENDDO
422 ENDDO
423 ENDDO
424
425 WRITE(msgBuf,'(A,E15.7)') 'err/err_init',
426 & err_max/err_init
427 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
428 & SQUEEZE_RIGHT , 1)
429
430 IF (err_max<err_lastchange*1.e-2 .and.
431 & STREAMICE_lower_cg_tol) THEN
432 cgtol = cgtol * 5.e-2
433 err_lastchange = err_max
434 WRITE(msgBuf,'(A,E15.7)') 'new cg tol: ',
435 & cgtol
436 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
437 & SQUEEZE_RIGHT , 1)
438 ENDIF
439
440
441 ENDIF ! (err_max .GT. streamice_nonlin_tol * err_init)
442 ENDDO
443
444 C END NL ITER. LOOP
445 C-------------------------------------------------------------------
446
447 if (iter_numconv .lt. streamice_max_nl_iter) then
448 WRITE(msgBuf,'(A,I5,A)') 'VELOCITY SOLVE CONVERGED, ',
449 & iter_numconv, ' iterations'
450 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
451 & SQUEEZE_RIGHT , 1)
452 else
453 WRITE(msgBuf,'(A,I5,A)') 'VELOCITY SOLVE NOT CONVERGED IN ',
454 & iter_numconv, ' iterations'
455 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
456 & SQUEEZE_RIGHT , 1)
457 endif
458
459 _EXCH_XY_RL (U_streamice, myThid)
460 _EXCH_XY_RL (V_streamice, myThid)
461
462 #endif
463 RETURN
464 END
465

  ViewVC Help
Powered by ViewVC 1.1.22