/[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.8 - (show annotations) (download)
Tue Sep 18 17:06:48 2012 UTC (12 years, 10 months ago) by dgoldberg
Branch: MAIN
Changes since 1.7: +40 -5 lines
changes for periodic boundary conds and hybrid stress balance

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

  ViewVC Help
Powered by ViewVC 1.1.22