/[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.6 - (hide annotations) (download)
Mon Jul 30 19:04:55 2012 UTC (12 years, 11 months ago) by dgoldberg
Branch: MAIN
Changes since 1.5: +38 -3 lines
allow for fixed point iteration for nonlinear solve

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

  ViewVC Help
Powered by ViewVC 1.1.22