/[MITgcm]/MITgcm/pkg/streamice/streamice_vel_solve.F
ViewVC logotype

Contents of /MITgcm/pkg/streamice/streamice_vel_solve.F

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


Revision 1.1 - (show annotations) (download)
Wed Jun 12 21:30:22 2013 UTC (10 years, 11 months ago) by dgoldberg
Branch: MAIN
new package from CONTRIB branch

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

  ViewVC Help
Powered by ViewVC 1.1.22