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

Contents of /MITgcm_contrib/dgoldberg/streamice/streamice_vel_solve.F

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


Revision 1.2 - (show annotations) (download)
Wed May 2 02:36:01 2012 UTC (13 years, 2 months ago) by heimbach
Branch: MAIN
Changes since 1.1: +102 -5 lines
Various updates to streamice code

1 C $Header: /u/gcmpack/MITgcm/pkg/streamice/streamice_init_varia.F,v 1.6 2011/06/29 16:24:10 dng Exp $
2 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_AUTODIFF_TAMC
25 # include "tamc.h"
26 #endif
27
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 INTEGER ikey_nl
46 _RL err_max, err_tempu, err_tempv, err_init, area
47 _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 #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 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 CALL STREAMICE_CG_BOUND_VALS( myThid,
98 O ubd_SI,
99 O vbd_SI)
100
101 #ifdef ALLOW_AUTODIFF_TAMC
102 !$TAF STORE U_streamice = comlev1, key=ikey_dynamics
103 !$TAF STORE V_streamice = comlev1, key=ikey_dynamics
104 #endif
105
106 CALL STREAMICE_CG_ACTION( myThid,
107 O Au_SI,
108 O Av_SI,
109 I U_streamice,
110 I V_streamice,
111 I 0, sNx+1, 0, sNy+1 )
112
113 err_init = 0. _d 0
114
115 DO bj = myByLo(myThid), myByHi(myThid)
116 DO bi = myBxLo(myThid), myBxHi(myThid)
117 #ifdef ALLOW_AUTODIFF_TAMC
118 !$TAF STORE err_init = comlev1, key=ikey_dynamics
119 #endif
120 DO j=1,sNy
121 DO i=1,sNx
122 err_tempu = 0. _d 0
123 err_tempv = 0. _d 0
124 IF (STREAMICE_umask(i,j,bi,bj).eq.1) THEN
125 err_tempu =
126 & ABS (Au_SI(i,j,bi,bj)+ubd_SI(i,j,bi,bj) -
127 & taudx_SI(i,j,bi,bj))
128 ENDIF
129 IF (STREAMICE_vmask(i,j,bi,bj).eq.1) THEN
130 err_tempv = MAX( err_tempu,
131 & ABS (Av_SI(i,j,bi,bj)+vbd_SI(i,j,bi,bj) -
132 & taudy_SI(i,j,bi,bj)))
133 ENDIF
134 IF (err_tempv .ge. err_init) err_init = err_tempv
135 ENDDO
136 ENDDO
137 ENDDO
138 ENDDO
139 #ifdef ALLOW_AUTODIFF_TAMC
140 !$TAF STORE err_init = comlev1, key=ikey_dynamics
141 #endif
142
143 CALL GLOBAL_MAX_R8 (err_init, myThid)
144
145 iter_numconv = 0
146 err_max = err_init
147 err_lastchange = err_init
148
149 C START NL ITER. LOOP
150 C-------------------------------------------------------------------
151
152 DO iter=1,streamice_max_nl_iter
153
154 C To avoid using "exit", loop goes through all iterations
155 C but after convergence loop does nothing
156
157 #ifdef ALLOW_AUTODIFF_TAMC
158 ikey_nl = (ikey_dynamics-1)*streamice_max_nl + iter
159 #endif
160 #ifdef ALLOW_AUTODIFF_TAMC
161 !$TAF STORE cgtol = comlev1_stream_nl, key=ikey_nl
162 !$TAF STORE err_max = comlev1_stream_nl, key=ikey_nl
163 !$TAF STORE err_tempu = comlev1_stream_nl, key=ikey_nl
164 !$TAF STORE err_tempv = comlev1_stream_nl, key=ikey_nl
165 !$TAF STORE err_lastchange = comlev1_stream_nl, key=ikey_nl
166 !$TAF STORE ru_old_si = comlev1_stream_nl, key=ikey_nl
167 !$TAF STORE rv_old_si = comlev1_stream_nl, key=ikey_nl
168 !$TAF STORE streamice_cg_a1 = comlev1_stream_nl, key=ikey_nl
169 !$TAF STORE streamice_cg_a2 = comlev1_stream_nl, key=ikey_nl
170 !$TAF STORE streamice_cg_a3 = comlev1_stream_nl, key=ikey_nl
171 !$TAF STORE streamice_cg_a4 = comlev1_stream_nl, key=ikey_nl
172 !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
173 !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
174 !$TAF STORE tau_beta_eff_streamice = comlev1_stream_nl, key=ikey_nl
175 !$TAF STORE visc_streamice = comlev1_stream_nl, key=ikey_nl
176 !$TAF STORE zu_old_si = comlev1_stream_nl, key=ikey_nl
177 !$TAF STORE zv_old_si = comlev1_stream_nl, key=ikey_nl
178 #endif
179
180 IF (err_max .GT. streamice_nonlin_tol * err_init) THEN
181
182 iter_numconv = iter_numconv + 1
183
184 #ifdef ALLOW_AUTODIFF_TAMC
185 !$TAF STORE cgtol = comlev1_stream_nl, key=ikey_nl
186 !$TAF STORE err_max = comlev1_stream_nl, key=ikey_nl
187 !$TAF STORE err_tempu = comlev1_stream_nl, key=ikey_nl
188 !$TAF STORE err_tempv = comlev1_stream_nl, key=ikey_nl
189 !$TAF STORE err_lastchange = comlev1_stream_nl, key=ikey_nl
190 !$TAF STORE ru_old_si = comlev1_stream_nl, key=ikey_nl
191 !$TAF STORE rv_old_si = comlev1_stream_nl, key=ikey_nl
192 !$TAF STORE streamice_cg_a1 = comlev1_stream_nl, key=ikey_nl
193 !$TAF STORE streamice_cg_a2 = comlev1_stream_nl, key=ikey_nl
194 !$TAF STORE streamice_cg_a3 = comlev1_stream_nl, key=ikey_nl
195 !$TAF STORE streamice_cg_a4 = comlev1_stream_nl, key=ikey_nl
196 !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
197 !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
198 !$TAF STORE tau_beta_eff_streamice = comlev1_stream_nl, key=ikey_nl
199 !$TAF STORE visc_streamice = comlev1_stream_nl, key=ikey_nl
200 !$TAF STORE zu_old_si = comlev1_stream_nl, key=ikey_nl
201 !$TAF STORE zv_old_si = comlev1_stream_nl, key=ikey_nl
202 #endif
203
204 CALL STREAMICE_CG_SOLVE(
205 & U_streamice,
206 & V_streamice,
207 & taudx_SI,
208 & taudy_SI,
209 & cgtol,
210 & cg_iters,
211 & myThid )
212
213
214
215 WRITE(msgBuf,'(A,I5,A,I4,A)') 'streamice linear solve number',
216 & iter, " ",
217 & cg_iters,
218 & ' iterations '
219 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
220 & SQUEEZE_RIGHT , 1)
221
222 #ifdef ALLOW_AUTODIFF_TAMC
223 !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
224 !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
225 #endif
226
227 CALL STREAMICE_VISC_BETA ( myThid )
228
229 #ifdef ALLOW_AUTODIFF_TAMC
230 !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
231 !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
232 #endif
233
234 _EXCH_XY_RL ( tau_beta_eff_streamice , myThid )
235 _EXCH_XY_RL ( visc_streamice , myThid )
236
237 DO bj = myByLo(myThid), myByHi(myThid)
238 DO bi = myBxLo(myThid), myBxHi(myThid)
239 DO j=1,sNy
240 DO i=1,sNx
241 Au_SI (i,j,bi,bj) = 0. _d 0
242 Av_SI (i,j,bi,bj) = 0. _d 0
243 ubd_SI (i,j,bi,bj) = 0. _d 0
244 vbd_SI (i,j,bi,bj) = 0. _d 0
245 ENDDO
246 ENDDO
247 ENDDO
248 ENDDO
249
250 CALL STREAMICE_CG_BOUND_VALS( myThid,
251 O ubd_SI,
252 O vbd_SI)
253
254 #ifdef ALLOW_AUTODIFF_TAMC
255 !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
256 !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
257 #endif
258
259 CALL STREAMICE_CG_ACTION( myThid,
260 O Au_SI,
261 O Av_SI,
262 I U_streamice,
263 I V_streamice,
264 I 0, sNx+1, 0, sNy+1 )
265
266 #ifdef ALLOW_AUTODIFF_TAMC
267 !$TAF STORE U_streamice = comlev1_stream_nl, key=ikey_nl
268 !$TAF STORE V_streamice = comlev1_stream_nl, key=ikey_nl
269 #endif
270
271 err_max = 0. _d 0
272
273 #ifdef ALLOW_AUTODIFF_TAMC
274 !$TAF STORE err_max = comlev1_stream_nl, key=ikey_nl
275 #endif
276 DO bj = myByLo(myThid), myByHi(myThid)
277 DO bi = myBxLo(myThid), myBxHi(myThid)
278 #ifdef ALLOW_AUTODIFF_TAMC
279 !$TAF STORE err_max = comlev1_stream_nl, key=ikey_nl
280 #endif
281 DO j=1,sNy
282 DO i=1,sNx
283 err_tempu = 0. _d 0
284 err_tempv = 0. _d 0
285 IF (STREAMICE_umask(i,j,bi,bj).eq.1) THEN
286 err_tempu =
287 & ABS (Au_SI(i,j,bi,bj)+ubd_SI(i,j,bi,bj) -
288 & taudx_SI(i,j,bi,bj))
289 ENDIF
290 IF (STREAMICE_vmask(i,j,bi,bj).eq.1) THEN
291 err_tempv = MAX( err_tempu,
292 & ABS (Av_SI(i,j,bi,bj)+vbd_SI(i,j,bi,bj) -
293 & taudy_SI(i,j,bi,bj)))
294 ENDIF
295 IF (err_tempv .ge. err_max) err_max = err_tempv
296 ENDDO
297 ENDDO
298 ENDDO
299 ENDDO
300
301 CALL GLOBAL_MAX_R8 (err_max, myThid)
302
303 WRITE(msgBuf,'(A,F11.7)') 'err/err_init',
304 & err_max/err_init
305 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
306 & SQUEEZE_RIGHT , 1)
307
308 IF (err_max<err_lastchange*1.e-2 .and.
309 & STREAMICE_lower_cg_tol) THEN
310 cgtol = cgtol * 5.e-2
311 err_lastchange = err_max
312 WRITE(msgBuf,'(A,F11.7)') 'new cg tol: ',
313 & cgtol
314 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
315 & SQUEEZE_RIGHT , 1)
316 ENDIF
317
318 ENDIF ! (err_max .GT. streamice_nonlin_tol * err_init)
319 ENDDO
320
321 C END NL ITER. LOOP
322 C-------------------------------------------------------------------
323
324 if (iter_numconv .lt. streamice_max_nl_iter) then
325 PRINT *, "VELOCITY SOLVE CONVERGED, ", iter_numconv,
326 & " iterations"
327 else
328 PRINT *, "VELOCITY SOLVE DID NOT CONVERGE IN ",
329 & iter, " iterations"
330 endif
331
332 _EXCH_XY_RL (U_streamice, myThid)
333 _EXCH_XY_RL (V_streamice, myThid)
334
335 #endif
336 RETURN
337 END
338

  ViewVC Help
Powered by ViewVC 1.1.22