/[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.7 - (hide annotations) (download)
Tue Sep 4 21:11:44 2012 UTC (12 years, 10 months ago) by dgoldberg
Branch: MAIN
Changes since 1.6: +17 -33 lines
recent changes.. DNG

1 dgoldberg 1.7 C $Header: /u/gcmpack/MITgcm_contrib/dgoldberg/streamice/streamice_vel_solve.F,v 1.6 2012/07/30 19:04:55 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 dgoldberg 1.7 !$TAF STORE err_max_fp = comlev1_stream_nl, key=ikey_nl
172 heimbach 1.2 !$TAF STORE err_tempu = comlev1_stream_nl, key=ikey_nl
173     !$TAF STORE err_tempv = comlev1_stream_nl, key=ikey_nl
174     !$TAF STORE err_lastchange = comlev1_stream_nl, key=ikey_nl
175     !$TAF STORE ru_old_si = comlev1_stream_nl, key=ikey_nl
176     !$TAF STORE rv_old_si = comlev1_stream_nl, key=ikey_nl
177 dgoldberg 1.7 !$TAF STORE u_old_si = comlev1_stream_nl, key=ikey_nl
178     !$TAF STORE v_old_si = comlev1_stream_nl, key=ikey_nl
179 heimbach 1.2 !$TAF STORE streamice_cg_a1 = comlev1_stream_nl, key=ikey_nl
180     !$TAF STORE streamice_cg_a2 = comlev1_stream_nl, key=ikey_nl
181     !$TAF STORE streamice_cg_a3 = comlev1_stream_nl, key=ikey_nl
182     !$TAF STORE streamice_cg_a4 = comlev1_stream_nl, key=ikey_nl
183     !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
184     !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
185     !$TAF STORE tau_beta_eff_streamice = comlev1_stream_nl, key=ikey_nl
186     !$TAF STORE visc_streamice = comlev1_stream_nl, key=ikey_nl
187     !$TAF STORE zu_old_si = comlev1_stream_nl, key=ikey_nl
188     !$TAF STORE zv_old_si = comlev1_stream_nl, key=ikey_nl
189     #endif
190    
191 dgoldberg 1.6 IF ((err_max .GT. streamice_nonlin_tol * err_init) .and.
192     & (err_max_fp .GT. streamice_nonlin_tol_fp)) THEN
193 heimbach 1.1
194     iter_numconv = iter_numconv + 1
195    
196 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
197     !$TAF STORE cgtol = comlev1_stream_nl, key=ikey_nl
198     !$TAF STORE err_max = comlev1_stream_nl, key=ikey_nl
199     !$TAF STORE err_tempu = comlev1_stream_nl, key=ikey_nl
200     !$TAF STORE err_tempv = comlev1_stream_nl, key=ikey_nl
201     !$TAF STORE err_lastchange = comlev1_stream_nl, key=ikey_nl
202     !$TAF STORE ru_old_si = comlev1_stream_nl, key=ikey_nl
203     !$TAF STORE rv_old_si = comlev1_stream_nl, key=ikey_nl
204     !$TAF STORE streamice_cg_a1 = comlev1_stream_nl, key=ikey_nl
205     !$TAF STORE streamice_cg_a2 = comlev1_stream_nl, key=ikey_nl
206     !$TAF STORE streamice_cg_a3 = comlev1_stream_nl, key=ikey_nl
207     !$TAF STORE streamice_cg_a4 = comlev1_stream_nl, key=ikey_nl
208     !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
209     !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
210     !$TAF STORE tau_beta_eff_streamice = comlev1_stream_nl, key=ikey_nl
211     !$TAF STORE visc_streamice = comlev1_stream_nl, key=ikey_nl
212     !$TAF STORE zu_old_si = comlev1_stream_nl, key=ikey_nl
213     !$TAF STORE zv_old_si = comlev1_stream_nl, key=ikey_nl
214     #endif
215    
216 dgoldberg 1.4 #ifdef ALLOW_AUTODIFF_TAMC
217     ! DO bj = myByLo(myThid), myByHi(myThid)
218     ! DO bi = myBxLo(myThid), myBxHi(myThid)
219     ! DO j=1-OLy,sNy+OLy
220     ! DO i=1-OLx,sNx+OLx
221     ! U_streamice (i,j,bi,bj) = 0. _d 0
222     ! V_streamice (i,j,bi,bj) = 0. _d 0
223     ! ENDDO
224     ! ENDDO
225     ! ENDDO
226     ! ENDDO
227     #endif
228 dgoldberg 1.3
229     CALL STREAMICE_CG_WRAPPER(
230     & U_streamice,
231     & V_streamice,
232     & taudx_SI,
233     & taudy_SI,
234     & cgtol,
235     & cg_iters,
236     & myThid )
237 heimbach 1.1
238    
239    
240     WRITE(msgBuf,'(A,I5,A,I4,A)') 'streamice linear solve number',
241     & iter, " ",
242     & cg_iters,
243     & ' iterations '
244     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
245     & SQUEEZE_RIGHT , 1)
246    
247 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
248     !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
249     !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
250     #endif
251    
252 heimbach 1.1 CALL STREAMICE_VISC_BETA ( myThid )
253 dgoldberg 1.4
254 heimbach 1.1
255 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
256     !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
257     !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
258     #endif
259    
260 heimbach 1.1 _EXCH_XY_RL ( tau_beta_eff_streamice , myThid )
261     _EXCH_XY_RL ( visc_streamice , myThid )
262    
263     DO bj = myByLo(myThid), myByHi(myThid)
264     DO bi = myBxLo(myThid), myBxHi(myThid)
265     DO j=1,sNy
266     DO i=1,sNx
267     Au_SI (i,j,bi,bj) = 0. _d 0
268     Av_SI (i,j,bi,bj) = 0. _d 0
269     ubd_SI (i,j,bi,bj) = 0. _d 0
270     vbd_SI (i,j,bi,bj) = 0. _d 0
271     ENDDO
272     ENDDO
273     ENDDO
274     ENDDO
275    
276     CALL STREAMICE_CG_BOUND_VALS( myThid,
277     O ubd_SI,
278     O vbd_SI)
279    
280 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
281     !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
282     !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
283     #endif
284    
285 heimbach 1.1 CALL STREAMICE_CG_ACTION( myThid,
286     O Au_SI,
287     O Av_SI,
288     I U_streamice,
289     I V_streamice,
290     I 0, sNx+1, 0, sNy+1 )
291    
292 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
293     !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
294     !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
295     #endif
296    
297 heimbach 1.1 err_max = 0. _d 0
298 dgoldberg 1.7 err_max_fp = 0. _d 0
299 heimbach 1.1
300 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
301     !$TAF STORE err_max = comlev1_stream_nl, key=ikey_nl
302     #endif
303 heimbach 1.1 DO bj = myByLo(myThid), myByHi(myThid)
304     DO bi = myBxLo(myThid), myBxHi(myThid)
305 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
306     !$TAF STORE err_max = comlev1_stream_nl, key=ikey_nl
307     #endif
308 heimbach 1.1 DO j=1,sNy
309     DO i=1,sNx
310 heimbach 1.2 err_tempu = 0. _d 0
311     err_tempv = 0. _d 0
312 heimbach 1.1 IF (STREAMICE_umask(i,j,bi,bj).eq.1) THEN
313     err_tempu =
314     & ABS (Au_SI(i,j,bi,bj)+ubd_SI(i,j,bi,bj) -
315     & taudx_SI(i,j,bi,bj))
316     ENDIF
317     IF (STREAMICE_vmask(i,j,bi,bj).eq.1) THEN
318     err_tempv = MAX( err_tempu,
319     & ABS (Av_SI(i,j,bi,bj)+vbd_SI(i,j,bi,bj) -
320     & taudy_SI(i,j,bi,bj)))
321     ENDIF
322     IF (err_tempv .ge. err_max) err_max = err_tempv
323     ENDDO
324     ENDDO
325     ENDDO
326     ENDDO
327    
328     CALL GLOBAL_MAX_R8 (err_max, myThid)
329    
330 dgoldberg 1.6 DO bj = myByLo(myThid), myByHi(myThid)
331     DO bi = myBxLo(myThid), myBxHi(myThid)
332     DO j=1,sNy
333     DO i=1,sNx
334     err_tempu = 0. _d 0
335     err_tempv = 0. _d 0
336     IF (STREAMICE_umask(i,j,bi,bj).eq.1) THEN
337     err_tempu =
338     & ABS (U_streamice(i,j,bi,bj)-u_old_SI(i,j,bi,bj))
339     ENDIF
340     IF (STREAMICE_vmask(i,j,bi,bj).eq.1) THEN
341     err_tempv = MAX( err_tempu,
342     & ABS (V_streamice(i,j,bi,bj)-v_old_SI(i,j,bi,bj)))
343     ENDIF
344     IF (err_tempv .ge. err_max_fp) err_max_fp = err_tempv
345     ENDDO
346     ENDDO
347     ENDDO
348     ENDDO
349    
350     CALL GLOBAL_MAX_R8 (err_max_fp, myThid)
351 dgoldberg 1.7 print *, "err_max_fp", err_max_fp
352    
353     DO bj = myByLo(myThid), myByHi(myThid)
354     DO bi = myBxLo(myThid), myBxHi(myThid)
355     DO j=1-OLy,sNy+OLy
356     DO i=1-OLx,sNx+OLx
357     u_old_SI (i,j,bi,bj) = U_streamice (i,j,bi,bj)
358     v_old_SI (i,j,bi,bj) = V_streamice (i,j,bi,bj)
359     ENDDO
360     ENDDO
361     ENDDO
362     ENDDO
363 dgoldberg 1.6
364 heimbach 1.1 WRITE(msgBuf,'(A,F11.7)') 'err/err_init',
365     & err_max/err_init
366     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
367     & SQUEEZE_RIGHT , 1)
368    
369     IF (err_max<err_lastchange*1.e-2 .and.
370     & STREAMICE_lower_cg_tol) THEN
371     cgtol = cgtol * 5.e-2
372     err_lastchange = err_max
373 dgoldberg 1.4 WRITE(msgBuf,'(A,E15.7)') 'new cg tol: ',
374 heimbach 1.1 & cgtol
375     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
376     & SQUEEZE_RIGHT , 1)
377     ENDIF
378    
379 dgoldberg 1.4
380 heimbach 1.1 ENDIF ! (err_max .GT. streamice_nonlin_tol * err_init)
381     ENDDO
382    
383 heimbach 1.2 C END NL ITER. LOOP
384     C-------------------------------------------------------------------
385    
386 heimbach 1.1 if (iter_numconv .lt. streamice_max_nl_iter) then
387 dgoldberg 1.5 WRITE(msgBuf,'(A,I5,A)') 'VELOCITY SOLVE CONVERGED, ',
388     & iter_numconv, ' iterations'
389     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
390     & SQUEEZE_RIGHT , 1)
391 heimbach 1.1 else
392 dgoldberg 1.5 WRITE(msgBuf,'(A,I5,A)') 'VELOCITY SOLVE NOT CONVERGED IN ',
393     & iter_numconv, ' iterations'
394     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
395     & SQUEEZE_RIGHT , 1)
396 heimbach 1.1 endif
397    
398     _EXCH_XY_RL (U_streamice, myThid)
399     _EXCH_XY_RL (V_streamice, myThid)
400    
401     #endif
402     RETURN
403     END
404    

  ViewVC Help
Powered by ViewVC 1.1.22