/[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.6 - (show annotations) (download)
Wed Mar 26 17:46:48 2014 UTC (10 years, 1 month ago) by dgoldberg
Branch: MAIN
CVS Tags: checkpoint64v
Changes since 1.5: +7 -7 lines
remove unnecessary i/o

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

  ViewVC Help
Powered by ViewVC 1.1.22