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

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

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


Revision 1.2 - (hide annotations) (download)
Wed Jun 12 22:19:07 2013 UTC (10 years, 11 months ago) by dgoldberg
Branch: MAIN
Changes since 1.1: +2 -2 lines
more digits in output

1 dgoldberg 1.2 C $Header: /u/gcmpack/MITgcm/pkg/streamice/streamice_vel_solve.F,v 1.1 2013/06/12 21:30:22 dgoldberg Exp $
2 dgoldberg 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     #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 dgoldberg 1.2 WRITE(msgBuf,'(A,1PE23.14)') 'STREAMICE_FP_ERROR',
421 dgoldberg 1.1 & 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