/[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.10 - (show annotations) (download)
Thu Sep 27 20:29:01 2012 UTC (12 years, 10 months ago) by dgoldberg
Branch: MAIN
Changes since 1.9: +6 -2 lines
various changes

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

  ViewVC Help
Powered by ViewVC 1.1.22