/[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.12 - (show annotations) (download)
Thu Mar 7 15:23:19 2013 UTC (12 years, 4 months ago) by dgoldberg
Branch: MAIN
Changes since 1.11: +4 -1 lines
bug fixes, GL smoothing, changes for controlling bathym with const surf elev

1 C $Header: /u/gcmpack/MITgcm_contrib/dgoldberg/streamice/streamice_vel_solve.F,v 1.11 2013/01/09 21:56: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 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
120 ! CALL WRITE_FLD_XY_RL("u_bound_cont","",ubd_SI,0,myThid)
121 ! CALL WRITE_FLD_XY_RL("v_bound_cont","",vbd_SI,0,myThid)
122 ! CALL WRITE_FLD_XY_RL("taudx_u","",taudx_SI,0,myThid)
123 ! CALL WRITE_FLD_XY_RL("taudx_v","",taudy_SI,0,myThid)
124
125 #ifdef ALLOW_AUTODIFF_TAMC
126 !$TAF STORE U_streamice = comlev1, key=ikey_dynamics
127 !$TAF STORE V_streamice = comlev1, key=ikey_dynamics
128 #endif
129
130 CALL STREAMICE_CG_ACTION( myThid,
131 O Au_SI,
132 O Av_SI,
133 I U_streamice,
134 I V_streamice,
135 I 0, sNx+1, 0, sNy+1 )
136
137
138
139 err_init = 0. _d 0
140
141 DO bj = myByLo(myThid), myByHi(myThid)
142 DO bi = myBxLo(myThid), myBxHi(myThid)
143 #ifdef ALLOW_AUTODIFF_TAMC
144 !$TAF STORE err_init = comlev1, key=ikey_dynamics
145 #endif
146 DO j=1,sNy
147 DO i=1,sNx
148 err_tempu = 0. _d 0
149 err_tempv = 0. _d 0
150 IF (STREAMICE_umask(i,j,bi,bj).eq.1) THEN
151 err_tempu =
152 & ABS (Au_SI(i,j,bi,bj)+0*ubd_SI(i,j,bi,bj) -
153 & taudx_SI(i,j,bi,bj))
154 ! print *, "err_temp_u", err_tempu
155 ENDIF
156 IF (STREAMICE_vmask(i,j,bi,bj).eq.1) THEN
157 err_tempv = MAX( err_tempu,
158 & ABS (Av_SI(i,j,bi,bj)+0*vbd_SI(i,j,bi,bj) -
159 & taudy_SI(i,j,bi,bj)))
160 ENDIF
161 IF (err_tempv .ge. err_init) THEN
162 err_init = err_tempv
163 ENDIF
164 ENDDO
165 ENDDO
166 ENDDO
167 ENDDO
168 #ifdef ALLOW_AUTODIFF_TAMC
169 !$TAF STORE err_init = comlev1, key=ikey_dynamics
170 #endif
171
172 CALL GLOBAL_MAX_R8 (err_init, myThid)
173
174 WRITE(msgBuf,'(A,E15.7)') 'initial nonlinear resid (error): ',
175 & err_init
176 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
177 & SQUEEZE_RIGHT , 1)
178
179
180 iter_numconv = 0
181 err_max = err_init
182 err_max_fp = streamice_nonlin_tol_fp * 10.
183 err_lastchange = err_init
184
185 C START NL ITER. LOOP
186 C-------------------------------------------------------------------
187
188 DO iter=1,streamice_max_nl_iter
189
190 C To avoid using "exit", loop goes through all iterations
191 C but after convergence loop does nothing
192
193 #ifdef ALLOW_AUTODIFF_TAMC
194 ikey_nl = (ikey_dynamics-1)*streamice_max_nl + iter
195 #endif
196 #ifdef ALLOW_AUTODIFF_TAMC
197 !$TAF STORE cgtol = comlev1_stream_nl, key=ikey_nl
198 !$TAF STORE err_max = comlev1_stream_nl, key=ikey_nl
199 !$TAF STORE err_max_fp = comlev1_stream_nl, key=ikey_nl
200 !$TAF STORE err_tempu = comlev1_stream_nl, key=ikey_nl
201 !$TAF STORE err_tempv = comlev1_stream_nl, key=ikey_nl
202 !$TAF STORE err_lastchange = comlev1_stream_nl, key=ikey_nl
203 !$TAF STORE ru_old_si = comlev1_stream_nl, key=ikey_nl
204 !$TAF STORE rv_old_si = comlev1_stream_nl, key=ikey_nl
205 !$TAF STORE u_old_si = comlev1_stream_nl, key=ikey_nl
206 !$TAF STORE v_old_si = comlev1_stream_nl, key=ikey_nl
207 !$TAF STORE streamice_cg_a1 = comlev1_stream_nl, key=ikey_nl
208 !$TAF STORE streamice_cg_a2 = comlev1_stream_nl, key=ikey_nl
209 !$TAF STORE streamice_cg_a3 = comlev1_stream_nl, key=ikey_nl
210 !$TAF STORE streamice_cg_a4 = comlev1_stream_nl, key=ikey_nl
211 !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
212 !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
213 !$TAF STORE tau_beta_eff_streamice = comlev1_stream_nl, key=ikey_nl
214 !$TAF STORE visc_streamice = comlev1_stream_nl, key=ikey_nl
215 !$TAF STORE zu_old_si = comlev1_stream_nl, key=ikey_nl
216 !$TAF STORE zv_old_si = comlev1_stream_nl, key=ikey_nl
217 #endif
218
219 IF ((err_max .GT. streamice_nonlin_tol * err_init) .and.
220 & (err_max_fp .GT. streamice_nonlin_tol_fp)) THEN
221
222 iter_numconv = iter_numconv + 1
223
224 #ifdef ALLOW_AUTODIFF_TAMC
225 !$TAF STORE cgtol = comlev1_stream_nl, key=ikey_nl
226 !$TAF STORE err_max = comlev1_stream_nl, key=ikey_nl
227 !$TAF STORE err_tempu = comlev1_stream_nl, key=ikey_nl
228 !$TAF STORE err_tempv = comlev1_stream_nl, key=ikey_nl
229 !$TAF STORE err_lastchange = comlev1_stream_nl, key=ikey_nl
230 !$TAF STORE ru_old_si = comlev1_stream_nl, key=ikey_nl
231 !$TAF STORE rv_old_si = comlev1_stream_nl, key=ikey_nl
232 !$TAF STORE streamice_cg_a1 = comlev1_stream_nl, key=ikey_nl
233 !$TAF STORE streamice_cg_a2 = comlev1_stream_nl, key=ikey_nl
234 !$TAF STORE streamice_cg_a3 = comlev1_stream_nl, key=ikey_nl
235 !$TAF STORE streamice_cg_a4 = comlev1_stream_nl, key=ikey_nl
236 !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
237 !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
238 !$TAF STORE tau_beta_eff_streamice = comlev1_stream_nl, key=ikey_nl
239 !$TAF STORE visc_streamice = comlev1_stream_nl, key=ikey_nl
240 !$TAF STORE zu_old_si = comlev1_stream_nl, key=ikey_nl
241 !$TAF STORE zv_old_si = comlev1_stream_nl, key=ikey_nl
242 #endif
243
244 #ifdef ALLOW_AUTODIFF_TAMC
245 ! DO bj = myByLo(myThid), myByHi(myThid)
246 ! DO bi = myBxLo(myThid), myBxHi(myThid)
247 ! DO j=1-OLy,sNy+OLy
248 ! DO i=1-OLx,sNx+OLx
249 ! U_streamice (i,j,bi,bj) = 0. _d 0
250 ! V_streamice (i,j,bi,bj) = 0. _d 0
251 ! ENDDO
252 ! ENDDO
253 ! ENDDO
254 ! ENDDO
255 #endif
256
257 CALL STREAMICE_CG_WRAPPER(
258 & U_streamice,
259 & V_streamice,
260 & taudx_SI,
261 & taudy_SI,
262 & cgtol,
263 & cg_iters,
264 & myThid )
265
266 #ifdef STREAMICE_HYBRID_STRESS
267 #ifdef ALLOW_AUTODIFF_TAMC
268 !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
269 !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
270 #endif
271 #endif
272
273 #ifdef STREAMICE_HYBRID_STRESS
274 CALL STREAMICE_TAUB (myThid)
275 #endif
276
277 WRITE(msgBuf,'(A,I5,A,I4,A)') 'streamice linear solve number',
278 & iter, " ",
279 & cg_iters,
280 & ' iterations '
281 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
282 & SQUEEZE_RIGHT , 1)
283
284 #ifdef ALLOW_AUTODIFF_TAMC
285 !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
286 !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
287 #ifdef STREAMICE_HYBRID_STRESS
288 !$TAF STORE streamice_taubx = comlev1_stream_nl, key=ikey_nl
289 !$TAF STORE streamice_tauby = comlev1_stream_nl, key=ikey_nl
290 !$TAF STORE visc_streamice_full = comlev1_stream_nl, key=ikey_nl
291 #endif
292 #endif
293
294 #ifdef STREAMICE_HYBRID_STRESS
295 CALL STREAMICE_VISC_BETA_HYBRID ( myThid )
296 #else
297 CALL STREAMICE_VISC_BETA ( myThid )
298 #endif
299
300
301 #ifdef ALLOW_AUTODIFF_TAMC
302 !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
303 !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
304 #endif
305
306 _EXCH_XY_RL ( tau_beta_eff_streamice , myThid )
307 _EXCH_XY_RL ( visc_streamice , myThid )
308
309 DO bj = myByLo(myThid), myByHi(myThid)
310 DO bi = myBxLo(myThid), myBxHi(myThid)
311 DO j=1,sNy
312 DO i=1,sNx
313 Au_SI (i,j,bi,bj) = 0. _d 0
314 Av_SI (i,j,bi,bj) = 0. _d 0
315 ubd_SI (i,j,bi,bj) = 0. _d 0
316 vbd_SI (i,j,bi,bj) = 0. _d 0
317 ENDDO
318 ENDDO
319 ENDDO
320 ENDDO
321
322 CALL STREAMICE_CG_BOUND_VALS( myThid,
323 O ubd_SI,
324 O vbd_SI)
325
326 #ifdef ALLOW_AUTODIFF_TAMC
327 !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
328 !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
329 #endif
330
331 CALL STREAMICE_CG_ACTION( myThid,
332 O Au_SI,
333 O Av_SI,
334 I U_streamice,
335 I V_streamice,
336 I 0, sNx+1, 0, sNy+1 )
337
338 if (iter .eq. streamice_max_nl_iter) then
339 CALL WRITE_FLD_XY_RL("u_bound_cont_A","",ubd_SI,0,myThid)
340 CALL WRITE_FLD_XY_RL("v_bound_cont_A","",vbd_SI,0,myThid)
341 CALL WRITE_FLD_XY_RL("u_bound_cont_B","",Au_SI,0,myThid)
342 CALL WRITE_FLD_XY_RL("v_bound_cont_B","",Av_SI,0,myThid)
343 endif
344
345 #ifdef ALLOW_AUTODIFF_TAMC
346 !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
347 !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
348 #endif
349
350 err_max = 0. _d 0
351 err_max_fp = 0. _d 0
352
353 #ifdef ALLOW_AUTODIFF_TAMC
354 !$TAF STORE err_max = comlev1_stream_nl, key=ikey_nl
355 #endif
356 DO bj = myByLo(myThid), myByHi(myThid)
357 DO bi = myBxLo(myThid), myBxHi(myThid)
358 #ifdef ALLOW_AUTODIFF_TAMC
359 !$TAF STORE err_max = comlev1_stream_nl, key=ikey_nl
360 #endif
361 DO j=1,sNy
362 DO i=1,sNx
363 err_tempu = 0. _d 0
364 err_tempv = 0. _d 0
365 IF (STREAMICE_umask(i,j,bi,bj).eq.1) THEN
366 err_tempu =
367 & ABS (Au_SI(i,j,bi,bj)+0*ubd_SI(i,j,bi,bj) -
368 & taudx_SI(i,j,bi,bj))
369 ENDIF
370 IF (STREAMICE_vmask(i,j,bi,bj).eq.1) THEN
371 err_tempv = MAX( err_tempu,
372 & ABS (Av_SI(i,j,bi,bj)+0*vbd_SI(i,j,bi,bj) -
373 & taudy_SI(i,j,bi,bj)))
374 ENDIF
375 ! if (err_tempu.ge.1.e2.or.err_tempv.ge.1.e2) THEN
376 ! print *, "FOUND MAX ", i,j,err_tempu,err_tempv,
377 ! & ubd_SI(i,j,bi,bj),vbd_SI(i,j,bi,bj)
378 ! endif
379 IF (err_tempv .ge. err_max) THEN
380 err_max = err_tempv
381 ENDIF
382 ENDDO
383 ENDDO
384 ENDDO
385 ENDDO
386
387 CALL GLOBAL_MAX_R8 (err_max, myThid)
388
389 DO bj = myByLo(myThid), myByHi(myThid)
390 DO bi = myBxLo(myThid), myBxHi(myThid)
391 DO j=1,sNy
392 DO i=1,sNx
393 err_tempu = 0. _d 0
394 err_tempv = 0. _d 0
395 IF (STREAMICE_umask(i,j,bi,bj).eq.1) THEN
396 err_tempu =
397 & ABS (U_streamice(i,j,bi,bj)-u_old_SI(i,j,bi,bj))
398 ENDIF
399 IF (STREAMICE_vmask(i,j,bi,bj).eq.1) THEN
400 err_tempv = MAX( err_tempu,
401 & ABS (V_streamice(i,j,bi,bj)-v_old_SI(i,j,bi,bj)))
402 ENDIF
403 IF (err_tempv .ge. err_max_fp) err_max_fp = err_tempv
404 ENDDO
405 ENDDO
406 ENDDO
407 ENDDO
408
409 CALL GLOBAL_MAX_R8 (err_max_fp, myThid)
410 WRITE(msgBuf,'(A,E15.7)') '||x_i-x_{i-1}||_inf',
411 & err_max_fp
412 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
413 & SQUEEZE_RIGHT , 1)
414
415 DO bj = myByLo(myThid), myByHi(myThid)
416 DO bi = myBxLo(myThid), myBxHi(myThid)
417 DO j=1-OLy,sNy+OLy
418 DO i=1-OLx,sNx+OLx
419 u_old_SI (i,j,bi,bj) = U_streamice (i,j,bi,bj)
420 v_old_SI (i,j,bi,bj) = V_streamice (i,j,bi,bj)
421 ENDDO
422 ENDDO
423 ENDDO
424 ENDDO
425
426 WRITE(msgBuf,'(A,E15.7)') 'err/err_init',
427 & err_max/err_init
428 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
429 & SQUEEZE_RIGHT , 1)
430
431 IF (err_max<err_lastchange*1.e-2 .and.
432 & STREAMICE_lower_cg_tol) THEN
433 cgtol = cgtol * 5.e-2
434 err_lastchange = err_max
435 WRITE(msgBuf,'(A,E15.7)') 'new cg tol: ',
436 & cgtol
437 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
438 & SQUEEZE_RIGHT , 1)
439 ENDIF
440
441
442 ENDIF ! (err_max .GT. streamice_nonlin_tol * err_init)
443 ENDDO
444
445 C END NL ITER. LOOP
446 C-------------------------------------------------------------------
447
448 if (iter_numconv .lt. streamice_max_nl_iter) then
449 WRITE(msgBuf,'(A,I5,A)') 'VELOCITY SOLVE CONVERGED, ',
450 & iter_numconv, ' iterations'
451 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
452 & SQUEEZE_RIGHT , 1)
453 else
454 WRITE(msgBuf,'(A,I5,A)') 'VELOCITY SOLVE NOT CONVERGED IN ',
455 & iter_numconv, ' iterations'
456 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
457 & SQUEEZE_RIGHT , 1)
458 endif
459
460 _EXCH_XY_RL (U_streamice, myThid)
461 _EXCH_XY_RL (V_streamice, myThid)
462 ! CALL WRITE_FLD_XY_RL("taubx","",streamice_taubx,0,myThid)
463 ! CALL WRITE_FLD_XY_RL("tauby","",streamice_tauby,0,myThid)
464
465 #endif
466 RETURN
467 END
468

  ViewVC Help
Powered by ViewVC 1.1.22