/[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.9 - (hide annotations) (download)
Thu Sep 20 02:04:45 2012 UTC (12 years, 10 months ago) by heimbach
Branch: MAIN
Changes since 1.8: +5 -5 lines
Enable working version for #define STREAMICE_HYBRID_STRESS

1 heimbach 1.9 C $Header: /u/gcmpack/MITgcm_contrib/dgoldberg/streamice/streamice_vel_solve.F,v 1.8 2012/09/18 17:06:48 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 dgoldberg 1.8 #ifdef STREAMICE_HYBRID_STRESS
78 heimbach 1.9 !$TAF STORE streamice_taubx = comlev1, key=ikey_dynamics
79     !$TAF STORE streamice_tauby = comlev1, key=ikey_dynamics
80     !$TAF STORE visc_streamice_full = comlev1, key=ikey_dynamics
81 dgoldberg 1.8 #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 heimbach 1.2 #endif
89    
90 dgoldberg 1.8 #ifdef STREAMICE_HYBRID_STRESS
91 heimbach 1.9 !$TAF STORE visc_streamice_full = comlev1, key=ikey_dynamics
92 dgoldberg 1.8 #endif
93 heimbach 1.1
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 dgoldberg 1.3
111 heimbach 1.1 CALL STREAMICE_CG_BOUND_VALS( myThid,
112     O ubd_SI,
113     O vbd_SI)
114    
115 heimbach 1.2 #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 heimbach 1.1 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 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
132     !$TAF STORE err_init = comlev1, key=ikey_dynamics
133     #endif
134 heimbach 1.1 DO j=1,sNy
135     DO i=1,sNx
136 heimbach 1.2 err_tempu = 0. _d 0
137     err_tempv = 0. _d 0
138 heimbach 1.1 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 dgoldberg 1.8 ! print *, "err_temp_u", err_tempu
143 heimbach 1.1 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 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
155     !$TAF STORE err_init = comlev1, key=ikey_dynamics
156     #endif
157 heimbach 1.1
158     CALL GLOBAL_MAX_R8 (err_init, myThid)
159    
160 dgoldberg 1.4 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 heimbach 1.1 iter_numconv = 0
167     err_max = err_init
168 dgoldberg 1.6 err_max_fp = streamice_nonlin_tol_fp * 10.
169 heimbach 1.1 err_lastchange = err_init
170    
171 heimbach 1.2 C START NL ITER. LOOP
172     C-------------------------------------------------------------------
173    
174 heimbach 1.1 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 heimbach 1.2 #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 dgoldberg 1.7 !$TAF STORE err_max_fp = comlev1_stream_nl, key=ikey_nl
186 heimbach 1.2 !$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 dgoldberg 1.7 !$TAF STORE u_old_si = comlev1_stream_nl, key=ikey_nl
192     !$TAF STORE v_old_si = comlev1_stream_nl, key=ikey_nl
193 heimbach 1.2 !$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 dgoldberg 1.6 IF ((err_max .GT. streamice_nonlin_tol * err_init) .and.
206     & (err_max_fp .GT. streamice_nonlin_tol_fp)) THEN
207 heimbach 1.1
208     iter_numconv = iter_numconv + 1
209    
210 heimbach 1.2 #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 dgoldberg 1.4 #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 dgoldberg 1.3
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 heimbach 1.1
252 dgoldberg 1.8 #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 heimbach 1.1
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 heimbach 1.2 #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 dgoldberg 1.8 #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 heimbach 1.2 #endif
279    
280 dgoldberg 1.8 #ifdef STREAMICE_HYBRID_STRESS
281     CALL STREAMICE_VISC_BETA_HYBRID ( myThid )
282     #else
283 heimbach 1.1 CALL STREAMICE_VISC_BETA ( myThid )
284 dgoldberg 1.8 #endif
285 dgoldberg 1.4
286 heimbach 1.1
287 heimbach 1.2 #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 heimbach 1.1 _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 heimbach 1.2 #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 heimbach 1.1 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 heimbach 1.2 #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 heimbach 1.1 err_max = 0. _d 0
330 dgoldberg 1.7 err_max_fp = 0. _d 0
331 heimbach 1.1
332 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
333     !$TAF STORE err_max = comlev1_stream_nl, key=ikey_nl
334     #endif
335 heimbach 1.1 DO bj = myByLo(myThid), myByHi(myThid)
336     DO bi = myBxLo(myThid), myBxHi(myThid)
337 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
338     !$TAF STORE err_max = comlev1_stream_nl, key=ikey_nl
339     #endif
340 heimbach 1.1 DO j=1,sNy
341     DO i=1,sNx
342 heimbach 1.2 err_tempu = 0. _d 0
343     err_tempv = 0. _d 0
344 heimbach 1.1 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 dgoldberg 1.6 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 dgoldberg 1.8 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 dgoldberg 1.7
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 dgoldberg 1.6
399 dgoldberg 1.8 WRITE(msgBuf,'(A,E15.7)') 'err/err_init',
400 heimbach 1.1 & 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 dgoldberg 1.4 WRITE(msgBuf,'(A,E15.7)') 'new cg tol: ',
409 heimbach 1.1 & cgtol
410     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
411     & SQUEEZE_RIGHT , 1)
412     ENDIF
413    
414 dgoldberg 1.4
415 heimbach 1.1 ENDIF ! (err_max .GT. streamice_nonlin_tol * err_init)
416     ENDDO
417    
418 heimbach 1.2 C END NL ITER. LOOP
419     C-------------------------------------------------------------------
420    
421 heimbach 1.1 if (iter_numconv .lt. streamice_max_nl_iter) then
422 dgoldberg 1.5 WRITE(msgBuf,'(A,I5,A)') 'VELOCITY SOLVE CONVERGED, ',
423     & iter_numconv, ' iterations'
424     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
425     & SQUEEZE_RIGHT , 1)
426 heimbach 1.1 else
427 dgoldberg 1.5 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 heimbach 1.1 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