SICOPOLIS V3.1
 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 
60 !-------- Term abbreviations --------
61 
62 do kc=0, kcmax
63  avxy3(kc) = deform*eaz_c(kc)/(ea-1.0_dp)*dzeta_c
64  aqxy1(kc) = deform/(ea-1.0_dp)*eaz_c(kc)*dzeta_c
65 end do
66 
67 !-------- Computation of stresses --------
68 
69 ! ------ Term abbreviations
70 
71 do i=0, imax
72 do j=0, jmax
73 
74  if (maske(j,i) == 0) then ! grounded ice
75 
76  do kc=0, kcmax
77  ctxyz1(kc,j,i) = rho*g*h_c(j,i)*(1.0_dp-eaz_c_quotient(kc))
78  end do
79 
80  if (n_cts(j,i) == 1) then ! temperate layer
81 
82  do kt=0, ktmax
83  ctxyz2(kt,j,i) = rho*g*h_t(j,i)*(1.0_dp-zeta_t(kt))
84  end do
85 
86  else ! cold base (-1), temperate base (0)
87 
88  do kt=0, ktmax
89  ctxyz2(kt,j,i) = 0.0_dp
90  end do
91 
92  end if
93 
94  else ! maske(j,i) == (1, 2 or 3)
95 
96  do kc=0, kcmax
97  ctxyz1(kc,j,i) = 0.0_dp
98  end do
99 
100  do kt=0, ktmax
101  ctxyz2(kt,j,i) = 0.0_dp
102  end do
103 
104  end if
105 
106 end do
107 end do
108 
109 ! ------ Shear stress txz (defined at (i+1/2,j,kc/t))
110 
111 do i=0, imax-1
112 do j=0, jmax
113 
114  do kc=0, kcmax
115  txz_c(kc,j,i) = -0.5_dp*(ctxyz1(kc,j,i)+ctxyz1(kc,j,i+1)) &
116  *dzs_dxi(j,i)
117  end do
118 
119  do kt=0, ktmax
120  txz_t(kt,j,i) = txz_c(0,j,i) &
121  -0.5_dp*(ctxyz2(kt,j,i)+ctxyz2(kt,j,i+1)) &
122  *dzs_dxi(j,i)
123  end do
124 
125 end do
126 end do
127 
128 ! ------ Shear stress tyz (defined at (i,j+1/2,kc/t))
129 
130 do i=0, imax
131 do j=0, jmax-1
132 
133  do kc=0, kcmax
134  tyz_c(kc,j,i) = -0.5_dp*(ctxyz1(kc,j,i)+ctxyz1(kc,j+1,i)) &
135  *dzs_deta(j,i)
136  end do
137 
138  do kt=0, ktmax
139  tyz_t(kt,j,i) = tyz_c(0,j,i) &
140  -0.5_dp*(ctxyz2(kt,j,i)+ctxyz2(kt,j+1,i)) &
141  *dzs_deta(j,i)
142  end do
143 
144 end do
145 end do
146 
147 ! ------ Effective shear stress sigma (defined at (i,j,kc/t))
148 
149 do i=0, imax
150 do j=0, jmax
151 
152  do kc=0, kcmax
153  sigma_c(kc,j,i) = ctxyz1(kc,j,i) &
154  *(dzs_dxi_g(j,i)**2+dzs_deta_g(j,i)**2)**0.5_dp
155  end do
156 
157  do kt=0, ktmax
158  sigma_t(kt,j,i) = sigma_c(0,j,i) &
159  + ctxyz2(kt,j,i) &
160  *(dzs_dxi_g(j,i)**2+dzs_deta_g(j,i)**2)**0.5_dp
161  end do
162 
163 end do
164 end do
165 
166 !-------- Computation of the depth-averaged fluidity
167 ! (defined on the grid points (i,j)) --------
168 
169 do i=0, imax
170 do j=0, jmax
171 
172  if (maske(j,i) == 0) then ! grounded ice
173 
174 ! ------ Fluidity, abbreviations
175 
176  do kt=0, ktmax
177  flui_t(kt) = 2.0_dp &
178  *enh_t(kt,j,i) &
179  *ratefac_t(omega_t(kt,j,i)) &
180  *creep(sigma_t(kt,j,i))
181  cflui0(kt) = h_t(j,i)*flui_t(kt)*dzeta_t
182  end do
183 
184  do kc=0, kcmax
185  flui_c(kc) = 2.0_dp &
186  *enh_c(kc,j,i) &
187  *ratefac(temp_c(kc,j,i), temp_c_m(kc,j,i)) &
188  *creep(sigma_c(kc,j,i))
189  cflui1(kc) = aqxy1(kc)*h_c(j,i)*flui_c(kc)
190  end do
191 
192 ! ------ Depth average
193 
194  flui_ave_sia(j,i) = 0.0_dp
195 
196  if (n_cts(j,i) == 1) then
197 
198  do kt=0, ktmax-1
199  flui_ave_sia(j,i) = flui_ave_sia(j,i)+0.5_dp*(cflui0(kt+1)+cflui0(kt))
200  end do
201 
202  end if
203 
204  do kc=0, kcmax-1
205  flui_ave_sia(j,i) = flui_ave_sia(j,i)+0.5_dp*(cflui1(kc+1)+cflui1(kc))
206  end do
207 
208  flui_ave_sia(j,i) = flui_ave_sia(j,i)/(h_c(j,i)+h_t(j,i))
209 
210  else ! maske(j,i) == (1, 2 or 3)
211 
212  flui_ave_sia(j,i) = 0.0_dp
213 
214  end if
215 
216 end do
217 end do
218 
219 !-------- Computation of d_help_c/t
220 ! (defined on the grid points (i,j,kc/t)) --------
221 
222 do i=0, imax
223 do j=0, jmax
224 
225  if (maske(j,i) == 0) then ! grounded ice
226 
227 ! ------ Abbreviations
228 
229  do kt=0, ktmax
230  cvxy2(kt) = 2.0_dp*h_t(j,i) &
231  *enh_t(kt,j,i) &
232  *ratefac_t(omega_t(kt,j,i)) &
233  *creep(sigma_t(kt,j,i)) &
234  *(ctxyz1(0,j,i)+ctxyz2(kt,j,i)) &
235  *dzeta_t
236  end do
237 
238  do kc=0, kcmax
239  cvxy3(kc) = 2.0_dp*avxy3(kc)*h_c(j,i) &
240  *enh_c(kc,j,i) &
241  *ratefac(temp_c(kc,j,i), temp_c_m(kc,j,i)) &
242  *creep(sigma_c(kc,j,i)) &
243  *ctxyz1(kc,j,i)
244  end do
245 
246 ! ------ d_help_c, d_help_t
247 
248  if (n_cts(j,i) == -1) then ! cold ice base
249 
250  do kt=0, ktmax
251  d_help_t(kt,j,i) = d_help_b(j,i)
252  end do
253 
254  d_help_c(0,j,i) = d_help_t(ktmax,j,i)
255 
256  do kc=0, kcmax-1
257  d_help_c(kc+1,j,i) = d_help_c(kc,j,i) &
258  +0.5_dp*(cvxy3(kc+1)+cvxy3(kc))
259  end do
260 
261  else if (n_cts(j,i) == 0) then ! temperate ice base
262 
263  do kt=0, ktmax
264  d_help_t(kt,j,i) = d_help_b(j,i)
265  end do
266 
267  d_help_c(0,j,i) = d_help_t(ktmax,j,i)
268 
269  do kc=0, kcmax-1
270  d_help_c(kc+1,j,i) = d_help_c(kc,j,i) &
271  +0.5_dp*(cvxy3(kc+1)+cvxy3(kc))
272  end do
273 
274  else ! n_cts(j,i) == 1, temperate ice layer
275 
276  d_help_t(0,j,i) = d_help_b(j,i)
277 
278  do kt=0, ktmax-1
279  d_help_t(kt+1,j,i) = d_help_t(kt,j,i) &
280  +0.5_dp*(cvxy2(kt+1)+cvxy2(kt))
281  end do
282 
283  d_help_c(0,j,i) = d_help_t(ktmax,j,i)
284 
285  do kc=0, kcmax-1
286  d_help_c(kc+1,j,i) = d_help_c(kc,j,i) &
287  +0.5_dp*(cvxy3(kc+1)+cvxy3(kc))
288  end do
289 
290  end if
291 
292  else ! maske(j,i) == (1, 2 or 3)
293 
294  do kt=0, ktmax
295  d_help_t(kt,j,i) = 0.0_dp
296  end do
297 
298  do kc=0, kcmax
299  d_help_c(kc,j,i) = 0.0_dp
300  end do
301 
302  end if
303 
304 end do
305 end do
306 
307 !-------- Computation of vx_c/t (defined at (i+1/2,j,kc/t)) --------
308 
309 do i=0, imax-1
310 do j=1, jmax-1
311 
312  do kt=0, ktmax
313  vx_t(kt,j,i) = -0.5_dp*(d_help_t(kt,j,i)+d_help_t(kt,j,i+1)) &
314  *dzs_dxi(j,i)
315  end do
316 
317  do kc=0, kcmax
318  vx_c(kc,j,i) = -0.5_dp*(d_help_c(kc,j,i)+d_help_c(kc,j,i+1)) &
319  *dzs_dxi(j,i)
320  end do
321 
322 end do
323 end do
324 
325 !-------- Computation of vy_c/t (defined at (i,j+1/2,kc/t)) --------
326 
327 do i=1, imax-1
328 do j=0, jmax-1
329 
330  do kt=0, ktmax
331  vy_t(kt,j,i) = -0.5_dp*(d_help_t(kt,j,i)+d_help_t(kt,j+1,i)) &
332  *dzs_deta(j,i)
333  end do
334 
335  do kc=0, kcmax
336  vy_c(kc,j,i) = -0.5_dp*(d_help_c(kc,j,i)+d_help_c(kc,j+1,i)) &
337  *dzs_deta(j,i)
338  end do
339 
340 end do
341 end do
342 
343 !-------- Computation of the surface velocities vx_s_g and vy_s_g
344 ! (defined at (i,j)) --------
345 
346 do i=0, imax
347 do j=0, jmax
348  vx_s_g(j,i) = -d_help_c(kcmax,j,i)*dzs_dxi_g(j,i)
349  vy_s_g(j,i) = -d_help_c(kcmax,j,i)*dzs_deta_g(j,i)
350 end do
351 end do
352 
353 !-------- Limitation of computed vx_c/t, vy_c/t, vx_s_g, vy_s_g
354 ! to the interval [-VH_MAX, VH_MAX] --------
355 
356 vh_max = vh_max/year_sec
357 
358 do i=0, imax
359 do j=0, jmax
360  vx_s_g(j,i) = max(vx_s_g(j,i), -vh_max)
361  vx_s_g(j,i) = min(vx_s_g(j,i), vh_max)
362  vy_s_g(j,i) = max(vy_s_g(j,i), -vh_max)
363  vy_s_g(j,i) = min(vy_s_g(j,i), vh_max)
364  do kt=0, ktmax
365  vx_t(kt,j,i) = max(vx_t(kt,j,i), -vh_max)
366  vx_t(kt,j,i) = min(vx_t(kt,j,i), vh_max)
367  vy_t(kt,j,i) = max(vy_t(kt,j,i), -vh_max)
368  vy_t(kt,j,i) = min(vy_t(kt,j,i), vh_max)
369  end do
370  do kc=0, kcmax
371  vx_c(kc,j,i) = max(vx_c(kc,j,i), -vh_max)
372  vx_c(kc,j,i) = min(vx_c(kc,j,i), vh_max)
373  vy_c(kc,j,i) = max(vy_c(kc,j,i), -vh_max)
374  vy_c(kc,j,i) = min(vy_c(kc,j,i), vh_max)
375  end do
376 end do
377 end do
378 
379 !-------- Computation of h_diff
380 ! (defined on the grid points (i,j)) --------
381 
382 do i=0, imax
383 do j=0, jmax
384 
385  if (maske(j,i) == 0) then ! grounded ice
386 
387 ! ------ Abbreviations
388 
389  do kt=0, ktmax
390  cqxy0(kt) = h_t(j,i)*d_help_t(kt,j,i)*dzeta_t
391  end do
392 
393  do kc=0, kcmax
394  cqxy1(kc) = aqxy1(kc)*h_c(j,i)*d_help_c(kc,j,i)
395  end do
396 
397 ! ------ h_diff
398 
399  h_diff(j,i) = 0.0_dp
400 
401  if (n_cts(j,i) == 1) then
402 
403  do kt=0, ktmax-1
404  h_diff(j,i) = h_diff(j,i)+0.5_dp*(cqxy0(kt+1)+cqxy0(kt))
405  end do
406 
407  end if
408 
409  do kc=0, kcmax-1
410  h_diff(j,i) = h_diff(j,i)+0.5_dp*(cqxy1(kc+1)+cqxy1(kc))
411  end do
412 
413 ! ------ Limitation of h_diff
414 
415  if (h_diff(j,i) < hd_min) h_diff(j,i) = 0.0_dp
416  if (h_diff(j,i) > hd_max) h_diff(j,i) = hd_max
417 
418  else ! maske(j,i) == (1, 2 or 3)
419 
420  h_diff(j,i) = 0.0_dp
421 
422  end if
423 
424 end do
425 end do
426 
427 !-------- Computation of the horizontal volume flux
428 ! and the depth-averaged velocity --------
429 
430 do i=0, imax-1
431 do j=0, jmax
432 
433  qx(j,i) = -0.5_dp*(h_diff(j,i)+h_diff(j,i+1))*dzs_dxi(j,i)
434 
435  if ( (maske(j,i)==0).or.(maske(j,i+1)==0) ) then ! at least one neighbour
436  ! point is grounded ice
437 
438  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)) )
439 
440  else
441  vx_m(j,i) = 0.0_dp
442  end if
443 
444 end do
445 end do
446 
447 do i=0, imax
448 do j=0, jmax-1
449 
450  qy(j,i) = -0.5_dp*(h_diff(j,i)+h_diff(j+1,i))*dzs_deta(j,i)
451 
452  if ( (maske(j,i)==0).or.(maske(j+1,i)==0) ) then ! at least one neighbour
453  ! point is grounded ice
454 
455  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)) )
456 
457  else
458  vy_m(j,i) = 0.0_dp
459  end if
460 
461 end do
462 end do
463 
464 !-------- Initialisation of the variable q_gl_g
465 ! [volume flux across the grounding line, to be
466 ! computed in the routine calc_vxy_ssa if MARGIN==3
467 ! (coupled ice-sheet/ice-shelf dynamics)]
468 
469 q_gl_g = 0.0_dp
470 
471 end subroutine calc_vxy_sia
472 !
Generated on Thu Sep 12 2013 09:45:12 for SICOPOLIS V3.1 by   doxygen 1.8.2