SICOPOLIS V3.0
 All Classes Files Functions Variables Macros
calc_vxy_sia.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : c a l c _ v x y _ s i a
4 !
5 !> @file
6 !!
7 !! Computation of the shear stresses txz, tyz, the effective shear stress
8 !! sigma, the depth-averaged fluidity flui_ave_sia, the horizontal
9 !! velocity vx, vy and the horizontal volume flux qx, qy in the shallow ice
10 !! approximation.
11 !!
12 !! @section Copyright
13 !!
14 !! Copyright 2009-2013 Ralf Greve
15 !!
16 !! @section License
17 !!
18 !! This file is part of SICOPOLIS.
19 !!
20 !! SICOPOLIS is free software: you can redistribute it and/or modify
21 !! it under the terms of the GNU General Public License as published by
22 !! the Free Software Foundation, either version 3 of the License, or
23 !! (at your option) any later version.
24 !!
25 !! SICOPOLIS is distributed in the hope that it will be useful,
26 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
27 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
28 !! GNU General Public License for more details.
29 !!
30 !! You should have received a copy of the GNU General Public License
31 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
32 !<
33 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
34 
35 !-------------------------------------------------------------------------------
36 !> Computation of the shear stresses txz, tyz, the effective shear stress
37 !! sigma, the depth-averaged fluidity flui_ave_sia, the horizontal
38 !! velocity vx, vy and the horizontal volume flux qx, qy in the shallow ice
39 !! approximation.
40 !<------------------------------------------------------------------------------
41 subroutine calc_vxy_sia(dzeta_c, dzeta_t)
42 
43 use sico_types
44 use sico_variables
45 
46 implicit none
47 
48 real(dp), intent(in) :: dzeta_c, dzeta_t
49 
50 integer(i4b) :: i, j, kc, kt
51 real(dp) :: avxy3(0:kcmax), aqxy1(0:kcmax)
52 real(dp) :: ctxyz1(0:kcmax,0:jmax,0:imax), &
53  ctxyz2(0:ktmax,0:jmax,0:imax)
54 real(dp) :: flui_t(0:ktmax), flui_c(0:kcmax)
55 real(dp) :: cflui0(0:ktmax), cflui1(0:kcmax)
56 real(dp) :: cvxy2(0:ktmax), cvxy3(0:kcmax)
57 real(dp) :: cqxy0(0:ktmax), cqxy1(0:kcmax)
58 real(dp) :: vh_max
59 real(dp) :: ratefac, ratefac_t, creep
60 
61 !-------- Term abbreviations --------
62 
63 do kc=0, kcmax
64  avxy3(kc) = deform*eaz_c(kc)/(ea-1.0_dp)*dzeta_c
65  aqxy1(kc) = deform/(ea-1.0_dp)*eaz_c(kc)*dzeta_c
66 end do
67 
68 !-------- Computation of stresses --------
69 
70 ! ------ Term abbreviations
71 
72 do i=0, imax
73 do j=0, jmax
74 
75  if (maske(j,i) == 0) then ! grounded ice
76 
77  do kc=0, kcmax
78  ctxyz1(kc,j,i) = rho*g*h_c(j,i)*(1.0_dp-eaz_c_quotient(kc))
79  end do
80 
81  if (n_cts(j,i) == 1) then ! temperate layer
82 
83  do kt=0, ktmax
84  ctxyz2(kt,j,i) = rho*g*h_t(j,i)*(1.0_dp-zeta_t(kt))
85  end do
86 
87  else ! cold base (-1), temperate base (0)
88 
89  do kt=0, ktmax
90  ctxyz2(kt,j,i) = 0.0_dp
91  end do
92 
93  end if
94 
95  else ! maske(j,i) == (1, 2 or 3)
96 
97  do kc=0, kcmax
98  ctxyz1(kc,j,i) = 0.0_dp
99  end do
100 
101  do kt=0, ktmax
102  ctxyz2(kt,j,i) = 0.0_dp
103  end do
104 
105  end if
106 
107 end do
108 end do
109 
110 ! ------ Shear stress txz (defined at (i+1/2,j,kc/t))
111 
112 do i=0, imax-1
113 do j=0, jmax
114 
115  do kc=0, kcmax
116  txz_c(kc,j,i) = -0.5_dp*(ctxyz1(kc,j,i)+ctxyz1(kc,j,i+1)) &
117  *dzs_dxi(j,i)
118  end do
119 
120  do kt=0, ktmax
121  txz_t(kt,j,i) = txz_c(0,j,i) &
122  -0.5_dp*(ctxyz2(kt,j,i)+ctxyz2(kt,j,i+1)) &
123  *dzs_dxi(j,i)
124  end do
125 
126 end do
127 end do
128 
129 ! ------ Shear stress tyz (defined at (i,j+1/2,kc/t))
130 
131 do i=0, imax
132 do j=0, jmax-1
133 
134  do kc=0, kcmax
135  tyz_c(kc,j,i) = -0.5_dp*(ctxyz1(kc,j,i)+ctxyz1(kc,j+1,i)) &
136  *dzs_deta(j,i)
137  end do
138 
139  do kt=0, ktmax
140  tyz_t(kt,j,i) = tyz_c(0,j,i) &
141  -0.5_dp*(ctxyz2(kt,j,i)+ctxyz2(kt,j+1,i)) &
142  *dzs_deta(j,i)
143  end do
144 
145 end do
146 end do
147 
148 ! ------ Effective shear stress sigma (defined at (i,j,kc/t))
149 
150 do i=0, imax
151 do j=0, jmax
152 
153  do kc=0, kcmax
154  sigma_c(kc,j,i) = ctxyz1(kc,j,i) &
155  *(dzs_dxi_g(j,i)**2+dzs_deta_g(j,i)**2)**0.5_dp
156  end do
157 
158  do kt=0, ktmax
159  sigma_t(kt,j,i) = sigma_c(0,j,i) &
160  + ctxyz2(kt,j,i) &
161  *(dzs_dxi_g(j,i)**2+dzs_deta_g(j,i)**2)**0.5_dp
162  end do
163 
164 end do
165 end do
166 
167 !-------- Computation of the depth-averaged fluidity
168 ! (defined on the grid points (i,j)) --------
169 
170 do i=0, imax
171 do j=0, jmax
172 
173  if (maske(j,i) == 0) then ! grounded ice
174 
175 ! ------ Fluidity, abbreviations
176 
177  do kt=0, ktmax
178  flui_t(kt) = 2.0_dp &
179  *enh_t(kt,j,i) &
180  *ratefac_t(omega_t(kt,j,i)) &
181  *creep(sigma_t(kt,j,i))
182  cflui0(kt) = h_t(j,i)*flui_t(kt)*dzeta_t
183  end do
184 
185  do kc=0, kcmax
186  flui_c(kc) = 2.0_dp &
187  *enh_c(kc,j,i) &
188  *ratefac(temp_c(kc,j,i), temp_c_m(kc,j,i)) &
189  *creep(sigma_c(kc,j,i))
190  cflui1(kc) = aqxy1(kc)*h_c(j,i)*flui_c(kc)
191  end do
192 
193 ! ------ Depth average
194 
195  flui_ave_sia(j,i) = 0.0_dp
196 
197  if (n_cts(j,i) == 1) then
198 
199  do kt=0, ktmax-1
200  flui_ave_sia(j,i) = flui_ave_sia(j,i)+0.5_dp*(cflui0(kt+1)+cflui0(kt))
201  end do
202 
203  end if
204 
205  do kc=0, kcmax-1
206  flui_ave_sia(j,i) = flui_ave_sia(j,i)+0.5_dp*(cflui1(kc+1)+cflui1(kc))
207  end do
208 
209  flui_ave_sia(j,i) = flui_ave_sia(j,i)/(h_c(j,i)+h_t(j,i))
210 
211  else ! maske(j,i) == (1, 2 or 3)
212 
213  flui_ave_sia(j,i) = 0.0_dp
214 
215  end if
216 
217 end do
218 end do
219 
220 !-------- Computation of d_help_c/t
221 ! (defined on the grid points (i,j,kc/t)) --------
222 
223 do i=0, imax
224 do j=0, jmax
225 
226  if (maske(j,i) == 0) then ! grounded ice
227 
228 ! ------ Abbreviations
229 
230  do kt=0, ktmax
231  cvxy2(kt) = 2.0_dp*h_t(j,i) &
232  *enh_t(kt,j,i) &
233  *ratefac_t(omega_t(kt,j,i)) &
234  *creep(sigma_t(kt,j,i)) &
235  *(ctxyz1(0,j,i)+ctxyz2(kt,j,i)) &
236  *dzeta_t
237  end do
238 
239  do kc=0, kcmax
240  cvxy3(kc) = 2.0_dp*avxy3(kc)*h_c(j,i) &
241  *enh_c(kc,j,i) &
242  *ratefac(temp_c(kc,j,i), temp_c_m(kc,j,i)) &
243  *creep(sigma_c(kc,j,i)) &
244  *ctxyz1(kc,j,i)
245  end do
246 
247 ! ------ d_help_c, d_help_t
248 
249  if (n_cts(j,i) == -1) then ! cold ice base
250 
251  do kt=0, ktmax
252  d_help_t(kt,j,i) = d_help_b(j,i)
253  end do
254 
255  d_help_c(0,j,i) = d_help_t(ktmax,j,i)
256 
257  do kc=0, kcmax-1
258  d_help_c(kc+1,j,i) = d_help_c(kc,j,i) &
259  +0.5_dp*(cvxy3(kc+1)+cvxy3(kc))
260  end do
261 
262  else if (n_cts(j,i) == 0) then ! temperate ice base
263 
264  do kt=0, ktmax
265  d_help_t(kt,j,i) = d_help_b(j,i)
266  end do
267 
268  d_help_c(0,j,i) = d_help_t(ktmax,j,i)
269 
270  do kc=0, kcmax-1
271  d_help_c(kc+1,j,i) = d_help_c(kc,j,i) &
272  +0.5_dp*(cvxy3(kc+1)+cvxy3(kc))
273  end do
274 
275  else ! n_cts(j,i) == 1, temperate ice layer
276 
277  d_help_t(0,j,i) = d_help_b(j,i)
278 
279  do kt=0, ktmax-1
280  d_help_t(kt+1,j,i) = d_help_t(kt,j,i) &
281  +0.5_dp*(cvxy2(kt+1)+cvxy2(kt))
282  end do
283 
284  d_help_c(0,j,i) = d_help_t(ktmax,j,i)
285 
286  do kc=0, kcmax-1
287  d_help_c(kc+1,j,i) = d_help_c(kc,j,i) &
288  +0.5_dp*(cvxy3(kc+1)+cvxy3(kc))
289  end do
290 
291  end if
292 
293  else ! maske(j,i) == (1, 2 or 3)
294 
295  do kt=0, ktmax
296  d_help_t(kt,j,i) = 0.0_dp
297  end do
298 
299  do kc=0, kcmax
300  d_help_c(kc,j,i) = 0.0_dp
301  end do
302 
303  end if
304 
305 end do
306 end do
307 
308 !-------- Computation of vx_c/t (defined at (i+1/2,j,kc/t)) --------
309 
310 do i=0, imax-1
311 do j=1, jmax-1
312 
313  do kt=0, ktmax
314  vx_t(kt,j,i) = -0.5_dp*(d_help_t(kt,j,i)+d_help_t(kt,j,i+1)) &
315  *dzs_dxi(j,i)
316  end do
317 
318  do kc=0, kcmax
319  vx_c(kc,j,i) = -0.5_dp*(d_help_c(kc,j,i)+d_help_c(kc,j,i+1)) &
320  *dzs_dxi(j,i)
321  end do
322 
323 end do
324 end do
325 
326 !-------- Computation of vy_c/t (defined at (i,j+1/2,kc/t)) --------
327 
328 do i=1, imax-1
329 do j=0, jmax-1
330 
331  do kt=0, ktmax
332  vy_t(kt,j,i) = -0.5_dp*(d_help_t(kt,j,i)+d_help_t(kt,j+1,i)) &
333  *dzs_deta(j,i)
334  end do
335 
336  do kc=0, kcmax
337  vy_c(kc,j,i) = -0.5_dp*(d_help_c(kc,j,i)+d_help_c(kc,j+1,i)) &
338  *dzs_deta(j,i)
339  end do
340 
341 end do
342 end do
343 
344 !-------- Computation of the surface velocities vx_s_g and vy_s_g
345 ! (defined at (i,j)) --------
346 
347 do i=0, imax
348 do j=0, jmax
349  vx_s_g(j,i) = -d_help_c(kcmax,j,i)*dzs_dxi_g(j,i)
350  vy_s_g(j,i) = -d_help_c(kcmax,j,i)*dzs_deta_g(j,i)
351 end do
352 end do
353 
354 !-------- Limitation of computed vx_c/t, vy_c/t, vx_s_g, vy_s_g
355 ! to the interval [-VH_MAX, VH_MAX] --------
356 
357 vh_max = vh_max/year_sec
358 
359 do i=0, imax
360 do j=0, jmax
361  vx_s_g(j,i) = max(vx_s_g(j,i), -vh_max)
362  vx_s_g(j,i) = min(vx_s_g(j,i), vh_max)
363  vy_s_g(j,i) = max(vy_s_g(j,i), -vh_max)
364  vy_s_g(j,i) = min(vy_s_g(j,i), vh_max)
365  do kt=0, ktmax
366  vx_t(kt,j,i) = max(vx_t(kt,j,i), -vh_max)
367  vx_t(kt,j,i) = min(vx_t(kt,j,i), vh_max)
368  vy_t(kt,j,i) = max(vy_t(kt,j,i), -vh_max)
369  vy_t(kt,j,i) = min(vy_t(kt,j,i), vh_max)
370  end do
371  do kc=0, kcmax
372  vx_c(kc,j,i) = max(vx_c(kc,j,i), -vh_max)
373  vx_c(kc,j,i) = min(vx_c(kc,j,i), vh_max)
374  vy_c(kc,j,i) = max(vy_c(kc,j,i), -vh_max)
375  vy_c(kc,j,i) = min(vy_c(kc,j,i), vh_max)
376  end do
377 end do
378 end do
379 
380 !-------- Computation of h_diff
381 ! (defined on the grid points (i,j)) --------
382 
383 do i=0, imax
384 do j=0, jmax
385 
386  if (maske(j,i) == 0) then ! grounded ice
387 
388 ! ------ Abbreviations
389 
390  do kt=0, ktmax
391  cqxy0(kt) = h_t(j,i)*d_help_t(kt,j,i)*dzeta_t
392  end do
393 
394  do kc=0, kcmax
395  cqxy1(kc) = aqxy1(kc)*h_c(j,i)*d_help_c(kc,j,i)
396  end do
397 
398 ! ------ h_diff
399 
400  h_diff(j,i) = 0.0_dp
401 
402  if (n_cts(j,i) == 1) then
403 
404  do kt=0, ktmax-1
405  h_diff(j,i) = h_diff(j,i)+0.5_dp*(cqxy0(kt+1)+cqxy0(kt))
406  end do
407 
408  end if
409 
410  do kc=0, kcmax-1
411  h_diff(j,i) = h_diff(j,i)+0.5_dp*(cqxy1(kc+1)+cqxy1(kc))
412  end do
413 
414 ! ------ Limitation of h_diff
415 
416  if (h_diff(j,i) < hd_min) h_diff(j,i) = 0.0_dp
417  if (h_diff(j,i) > hd_max) h_diff(j,i) = hd_max
418 
419  else ! maske(j,i) == (1, 2 or 3)
420 
421  h_diff(j,i) = 0.0_dp
422 
423  end if
424 
425 end do
426 end do
427 
428 !-------- Computation of the horizontal volume flux
429 ! and the depth-averaged velocity --------
430 
431 do i=0, imax-1
432 do j=0, jmax
433 
434  qx(j,i) = -0.5_dp*(h_diff(j,i)+h_diff(j,i+1))*dzs_dxi(j,i)
435 
436  if ( (maske(j,i)==0).or.(maske(j,i+1)==0) ) then ! at least one neighbour
437  ! point is grounded ice
438 
439  vx_m(j,i) = qx(j,i) / ( 0.5_dp*(h_c(j,i)+h_t(j,i)+h_c(j,i+1)+h_t(j,i+1)) )
440 
441  else
442  vx_m(j,i) = 0.0_dp
443  end if
444 
445 end do
446 end do
447 
448 do i=0, imax
449 do j=0, jmax-1
450 
451  qy(j,i) = -0.5_dp*(h_diff(j,i)+h_diff(j+1,i))*dzs_deta(j,i)
452 
453  if ( (maske(j,i)==0).or.(maske(j+1,i)==0) ) then ! at least one neighbour
454  ! point is grounded ice
455 
456  vy_m(j,i) = qy(j,i) / ( 0.5_dp*(h_c(j,i)+h_t(j,i)+h_c(j+1,i)+h_t(j+1,i)) )
457 
458  else
459  vy_m(j,i) = 0.0_dp
460  end if
461 
462 end do
463 end do
464 
465 !-------- Initialisation of the variable q_gl_g
466 ! [volume flux across the grounding line, to be
467 ! computed in the routine calc_vxy_ssa if MARGIN==3
468 ! (coupled ice-sheet/ice-shelf dynamics)]
469 
470 q_gl_g = 0.0_dp
471 
472 end subroutine calc_vxy_sia
473 !
Generated on Mon May 20 2013 10:18:41 for SICOPOLIS V3.0 by   doxygen 1.8.1