/[MITgcm]/MITgcm_contrib/dgoldberg/streamice/streamice_vel_solve.F
ViewVC logotype

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

1 dgoldberg 1.10 C $Header: /u/gcmpack/MITgcm_contrib/dgoldberg/streamice/streamice_vel_solve.F,v 1.9 2012/09/20 02:04:45 heimbach Exp $
2 heimbach 1.1 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 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
25     # include "tamc.h"
26     #endif
27 heimbach 1.1
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 heimbach 1.2 INTEGER ikey_nl
46 dgoldberg 1.6 _RL err_max, err_tempu, err_tempv, err_init, area, err_max_fp
47 heimbach 1.1 _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 dgoldberg 1.10 IF (STREAMICE_ppm_driving_stress) THEN
53     CALL STREAMICE_DRIVING_STRESS_PPM (myThid)
54     ELSE
55     CALL STREAMICE_DRIVING_STRESS (myThid)
56     ENDIF
57 heimbach 1.1
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 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
79     !$TAF STORE U_streamice = comlev1, key=ikey_dynamics
80     !$TAF STORE V_streamice = comlev1, key=ikey_dynamics
81 dgoldberg 1.8 #ifdef STREAMICE_HYBRID_STRESS
82 heimbach 1.9 !$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 dgoldberg 1.8 #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 heimbach 1.2 #endif
93    
94 dgoldberg 1.8 #ifdef STREAMICE_HYBRID_STRESS
95 heimbach 1.9 !$TAF STORE visc_streamice_full = comlev1, key=ikey_dynamics
96 dgoldberg 1.8 #endif
97 heimbach 1.1
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 dgoldberg 1.3
115 heimbach 1.1 CALL STREAMICE_CG_BOUND_VALS( myThid,
116     O ubd_SI,
117     O vbd_SI)
118    
119 heimbach 1.2 #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 heimbach 1.1 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 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
136     !$TAF STORE err_init = comlev1, key=ikey_dynamics
137     #endif
138 heimbach 1.1 DO j=1,sNy
139     DO i=1,sNx
140 heimbach 1.2 err_tempu = 0. _d 0
141     err_tempv = 0. _d 0
142 heimbach 1.1 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 dgoldberg 1.8 ! print *, "err_temp_u", err_tempu
147 heimbach 1.1 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 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
159     !$TAF STORE err_init = comlev1, key=ikey_dynamics
160     #endif
161 heimbach 1.1
162     CALL GLOBAL_MAX_R8 (err_init, myThid)
163    
164 dgoldberg 1.4 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 heimbach 1.1 iter_numconv = 0
171     err_max = err_init
172 dgoldberg 1.6 err_max_fp = streamice_nonlin_tol_fp * 10.
173 heimbach 1.1 err_lastchange = err_init
174    
175 heimbach 1.2 C START NL ITER. LOOP
176     C-------------------------------------------------------------------
177    
178 heimbach 1.1 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 heimbach 1.2 #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 dgoldberg 1.7 !$TAF STORE err_max_fp = comlev1_stream_nl, key=ikey_nl
190 heimbach 1.2 !$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 dgoldberg 1.7 !$TAF STORE u_old_si = comlev1_stream_nl, key=ikey_nl
196     !$TAF STORE v_old_si = comlev1_stream_nl, key=ikey_nl
197 heimbach 1.2 !$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 dgoldberg 1.6 IF ((err_max .GT. streamice_nonlin_tol * err_init) .and.
210     & (err_max_fp .GT. streamice_nonlin_tol_fp)) THEN
211 heimbach 1.1
212     iter_numconv = iter_numconv + 1
213    
214 heimbach 1.2 #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 dgoldberg 1.4 #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 dgoldberg 1.3
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 heimbach 1.1
256 dgoldberg 1.8 #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 heimbach 1.1
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 heimbach 1.2 #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 dgoldberg 1.8 #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 heimbach 1.2 #endif
283    
284 dgoldberg 1.8 #ifdef STREAMICE_HYBRID_STRESS
285     CALL STREAMICE_VISC_BETA_HYBRID ( myThid )
286     #else
287 heimbach 1.1 CALL STREAMICE_VISC_BETA ( myThid )
288 dgoldberg 1.8 #endif
289 dgoldberg 1.4
290 heimbach 1.1
291 heimbach 1.2 #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 heimbach 1.1 _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 heimbach 1.2 #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 heimbach 1.1 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 heimbach 1.2 #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 heimbach 1.1 err_max = 0. _d 0
334 dgoldberg 1.7 err_max_fp = 0. _d 0
335 heimbach 1.1
336 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
337     !$TAF STORE err_max = comlev1_stream_nl, key=ikey_nl
338     #endif
339 heimbach 1.1 DO bj = myByLo(myThid), myByHi(myThid)
340     DO bi = myBxLo(myThid), myBxHi(myThid)
341 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
342     !$TAF STORE err_max = comlev1_stream_nl, key=ikey_nl
343     #endif
344 heimbach 1.1 DO j=1,sNy
345     DO i=1,sNx
346 heimbach 1.2 err_tempu = 0. _d 0
347     err_tempv = 0. _d 0
348 heimbach 1.1 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 dgoldberg 1.6 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 dgoldberg 1.8 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 dgoldberg 1.7
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 dgoldberg 1.6
403 dgoldberg 1.8 WRITE(msgBuf,'(A,E15.7)') 'err/err_init',
404 heimbach 1.1 & 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 dgoldberg 1.4 WRITE(msgBuf,'(A,E15.7)') 'new cg tol: ',
413 heimbach 1.1 & cgtol
414     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
415     & SQUEEZE_RIGHT , 1)
416     ENDIF
417    
418 dgoldberg 1.4
419 heimbach 1.1 ENDIF ! (err_max .GT. streamice_nonlin_tol * err_init)
420     ENDDO
421    
422 heimbach 1.2 C END NL ITER. LOOP
423     C-------------------------------------------------------------------
424    
425 heimbach 1.1 if (iter_numconv .lt. streamice_max_nl_iter) then
426 dgoldberg 1.5 WRITE(msgBuf,'(A,I5,A)') 'VELOCITY SOLVE CONVERGED, ',
427     & iter_numconv, ' iterations'
428     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
429     & SQUEEZE_RIGHT , 1)
430 heimbach 1.1 else
431 dgoldberg 1.5 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 heimbach 1.1 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