SICOPOLIS V3.0
 All Classes Files Functions Variables Macros
calc_vxy_b_sia.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : c a l c _ v x y _ b _s i a
4 !
5 !> @file
6 !!
7 !! Computation of the basal horizontal velocity vx_b, vy_b in the shallow ice
8 !! approximation.
9 !!
10 !! @section Copyright
11 !!
12 !! Copyright 2009-2013 Ralf Greve
13 !!
14 !! @section License
15 !!
16 !! This file is part of SICOPOLIS.
17 !!
18 !! SICOPOLIS is free software: you can redistribute it and/or modify
19 !! it under the terms of the GNU General Public License as published by
20 !! the Free Software Foundation, either version 3 of the License, or
21 !! (at your option) any later version.
22 !!
23 !! SICOPOLIS is distributed in the hope that it will be useful,
24 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
25 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
26 !! GNU General Public License for more details.
27 !!
28 !! You should have received a copy of the GNU General Public License
29 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
30 !<
31 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
32 
33 !-------------------------------------------------------------------------------
34 !> Computation of the basal horizontal velocity vx_b, vy_b in the shallow ice
35 !! approximation.
36 !<------------------------------------------------------------------------------
37 subroutine calc_vxy_b_sia(time, z_sl)
38 
39 use sico_types
40 use sico_variables
41 
42 implicit none
43 
44 real(dp), intent(in) :: time, z_sl
45 
46 integer(i4b) :: i, j
47 integer(i4b) :: r_smw_1, s_smw_1, r_smw_2, s_smw_2
48 integer(i4b), dimension(0:JMAX,0:IMAX) :: p_weert, q_weert
49 real(dp), dimension(0:JMAX,0:IMAX) :: c_slide, gamma_slide_inv
50 real(dp) :: gamma_slide_inv_1, gamma_slide_inv_2
51 real(dp) :: smw_coeff_1, smw_coeff_2
52 real(dp), dimension(0:JMAX,0:IMAX) :: p_b, p_b_red, tau_b
53 real(dp) :: cvxy1, cvxy1a
54 real(dp) :: temp_diff
55 real(dp) :: vh_max
56 real(dp) :: year_sec_inv
57 logical, dimension(0:JMAX,0:IMAX) :: sub_melt_flag
58 
59 !-------- Sliding-law coefficients --------
60 
61 year_sec_inv = 1.0_dp/year_sec
62 
63 #if ( defined(ANT) && ( !defined(SEDI_SLIDE) || SEDI_SLIDE==1 ) )
64  ! Antarctica without soft sediment,
65  ! one sliding law for the entire modelling domain
66 
67 gamma_slide_inv_1 = 1.0_dp/max(gamma_slide, eps)
68 
69 do i=0, imax
70 do j=0, jmax
71  c_slide(j,i) = c_slide*year_sec_inv
72  p_weert(j,i) = p_weert
73  q_weert(j,i) = q_weert
74  gamma_slide_inv(j,i) = gamma_slide_inv_1
75  sub_melt_flag(j,i) = (gamma_slide >= eps)
76 end do
77 end do
78 
79 #elif ( defined(ANT) && SEDI_SLIDE==2 )
80  ! Antarctica with soft sediment,
81  ! different sliding laws for hard rock and soft sediment
82 
83 gamma_slide_inv_1 = 1.0_dp/max(gamma_slide, eps)
84 gamma_slide_inv_2 = 1.0_dp/max(gamma_slide_sedi, eps)
85 
86 do i=0, imax
87 do j=0, jmax
88 
89  if (maske_sedi(j,i) /= 7) then ! no soft sediment
90  c_slide(j,i) = c_slide*year_sec_inv
91  p_weert(j,i) = p_weert
92  q_weert(j,i) = q_weert
93  gamma_slide_inv(j,i) = gamma_slide_inv_1
94  sub_melt_flag(j,i) = (gamma_slide >= eps)
95  else ! maske_sedi(j,i) == 7, soft sediment
96  c_slide(j,i) = c_slide_sedi*year_sec_inv
97  p_weert(j,i) = p_weert_sedi
98  q_weert(j,i) = q_weert_sedi
99  gamma_slide_inv(j,i) = gamma_slide_inv_2
100  sub_melt_flag(j,i) = (gamma_slide_sedi >= eps)
101  end if
102 
103 end do
104 end do
105 
106 #elif ( defined(GRL) && ICE_STREAM==1 \
107  || defined(asf) && ice_stream==1 )
108  ! Greenland or Austfonna without ice streams,
109  ! one sliding law for the entire modelling domain,
110  ! with surface meltwater effect
111 
112 gamma_slide_inv_1 = 1.0_dp/max(gamma_slide, eps)
113 
114 r_smw_1 = r_smw
115 s_smw_1 = s_smw
116 smw_coeff_1 = smw_coeff*year_sec**s_smw_1
117 
118 do i=0, imax
119 do j=0, jmax
120  c_slide(j,i) = c_slide*year_sec_inv &
121  *(1.0_dp+smw_coeff_1*runoff(j,i)**s_smw_1 &
122  /max((h_c(j,i)+h_t(j,i))**r_smw_1, eps))
123  p_weert(j,i) = p_weert
124  q_weert(j,i) = q_weert
125  gamma_slide_inv(j,i) = gamma_slide_inv_1
126  sub_melt_flag(j,i) = (gamma_slide >= eps)
127 end do
128 end do
129 
130 #elif ( defined(GRL) && ICE_STREAM==2 \
131  || defined(asf) && ice_stream==2 )
132  ! Greenland or Austfonna with ice streams,
133  ! different sliding laws for hard rock and ice streams
134 
135 gamma_slide_inv_1 = 1.0_dp/max(gamma_slide, eps)
136 gamma_slide_inv_2 = 1.0_dp/max(gamma_slide_sedi, eps)
137 
138 r_smw_1 = r_smw
139 s_smw_1 = s_smw
140 smw_coeff_1 = smw_coeff *year_sec**s_smw_1
141 r_smw_2 = r_smw_sedi
142 s_smw_2 = s_smw_sedi
143 smw_coeff_2 = smw_coeff_sedi*year_sec**s_smw_2
144 
145 do i=0, imax
146 do j=0, jmax
147 
148  if (maske_sedi(j,i) /= 7) then ! no ice stream
149  c_slide(j,i) = c_slide*year_sec_inv &
150  *(1.0_dp+smw_coeff_1*runoff(j,i)**s_smw_1 &
151  /max((h_c(j,i)+h_t(j,i))**r_smw_1, eps))
152  p_weert(j,i) = p_weert
153  q_weert(j,i) = q_weert
154  gamma_slide_inv(j,i) = gamma_slide_inv_1
155  sub_melt_flag(j,i) = (gamma_slide >= eps)
156  else ! maske_sedi(j,i) == 7, ice stream
157  c_slide(j,i) = c_slide_sedi*year_sec_inv &
158  *(1.0_dp+smw_coeff_2*runoff(j,i)**s_smw_2 &
159  /max((h_c(j,i)+h_t(j,i))**r_smw_2, eps))
160  p_weert(j,i) = p_weert_sedi
161  q_weert(j,i) = q_weert_sedi
162  gamma_slide_inv(j,i) = gamma_slide_inv_2
163  sub_melt_flag(j,i) = (gamma_slide_sedi >= eps)
164  end if
165 
166 end do
167 end do
168 
169 #elif ( defined(HEINO) )
170  ! ISMIP HEINO,
171  ! different sliding laws for hard rock and soft sediment
172 
173 gamma_slide_inv_1 = 1.0_dp/max(gamma_slide, eps)
174 gamma_slide_inv_2 = 1.0_dp/max(gamma_slide_sedi, eps)
175 
176 do i=0, imax
177 do j=0, jmax
178 
179  if (maske_sedi(j,i) /= 7) then ! no soft sediment
180  c_slide(j,i) = c_slide*year_sec_inv
181  p_weert(j,i) = p_weert
182  q_weert(j,i) = q_weert
183  gamma_slide_inv(j,i) = gamma_slide_inv_1
184  sub_melt_flag(j,i) = (gamma_slide >= eps)
185  else ! maske_sedi(j,i) == 7, soft sediment
186  c_slide(j,i) = c_slide_sedi*year_sec_inv
187  p_weert(j,i) = p_weert_sedi
188  q_weert(j,i) = q_weert_sedi
189  gamma_slide_inv(j,i) = gamma_slide_inv_2
190  sub_melt_flag(j,i) = (gamma_slide_sedi >= eps)
191  end if
192 
193 end do
194 end do
195 
196 #elif ( defined(SCAND) \
197  || defined(nhem) \
198  || defined(tibet) \
199  || defined(nmars) \
200  || defined(smars) \
201  || defined(emtp2sge) )
202  ! one sliding law for the entire modelling domain
203 
204 gamma_slide_inv_1 = 1.0_dp/max(gamma_slide, eps)
205 
206 do i=0, imax
207 do j=0, jmax
208  c_slide(j,i) = c_slide*year_sec_inv
209  p_weert(j,i) = p_weert
210  q_weert(j,i) = q_weert
211  gamma_slide_inv(j,i) = gamma_slide_inv_1
212  sub_melt_flag(j,i) = (gamma_slide >= eps)
213 end do
214 end do
215 
216 #else
217 
218 stop ' calc_vxy_b_sia: No valid domain specified!'
219 
220 #endif
221 
222 !-------- Computation of basal stresses --------
223 
224 ! ------ Basal pressure p_b, basal water pressure p_b_w,
225 ! reduced pressure p_b_red
226 
227 do i=0, imax
228 do j=0, jmax
229 
230  if (maske(j,i) == 0) then ! glaciated land
231 
232  p_b(j,i) = rho*g*(h_c(j,i)+h_t(j,i))
233  p_b_w(j,i) = rho_sw*g*max((z_sl-zb(j,i)), 0.0_dp)
234  p_b_red(j,i) = p_b(j,i)-p_b_w(j,i)
235  p_b_red(j,i) = max(p_b_red(j,i), red_pres_limit_fact*p_b(j,i))
236  ! in order to avoid very small values, which may lead to
237  ! huge sliding velocities
238 
239  else ! maske(j,i) == (1 or 2)
240 
241  p_b(j,i) = 0.0_dp
242  p_b_w(j,i) = 0.0_dp
243  p_b_red(j,i) = 0.0_dp
244 
245  end if
246 
247 end do
248 end do
249 
250 ! ------ Absolute value of the basal shear stress, tau_b
251 
252 do i=0, imax
253 do j=0, jmax
254  tau_b(j,i) = p_b(j,i)*(dzs_dxi_g(j,i)**2+dzs_deta_g(j,i)**2)**0.5_dp
255 end do
256 end do
257 
258 !-------- Computation of d_help_b (defined on the grid points (i,j)) --------
259 
260 do i=0, imax
261 do j=0, jmax
262 
263  if (maske(j,i) == 0) then ! glaciated land
264 
265 ! ------ Abbreviations
266 
267 #if SLIDE_LAW==1
268  cvxy1 = c_slide(j,i) &
269  * (tau_b(j,i)**(p_weert(j,i)-1)/p_b(j,i)**q_weert(j,i)) &
270  * p_b(j,i)
271  ! Basal sliding at pressure melting
272 #elif SLIDE_LAW==2
273  cvxy1 = c_slide(j,i) &
274  * (tau_b(j,i)**(p_weert(j,i)-1)/p_b_red(j,i)**q_weert(j,i)) &
275  * p_b(j,i)
276  ! Basal sliding at pressure melting
277 #endif
278 
279 ! ------ d_help_b
280 
281  if (n_cts(j,i) == -1) then ! cold ice base
282 
283  if (sub_melt_flag(j,i)) then
284  temp_diff = max((temp_c_m(0,j,i)-temp_c(0,j,i)), 0.0_dp)
285  cvxy1a = exp(-gamma_slide_inv(j,i)*temp_diff) ! sub-melt sliding
286  else
287  cvxy1a = 0.0_dp ! no sub-melt sliding
288  end if
289 
290  d_help_b(j,i) = cvxy1*cvxy1a
291 
292  else if (n_cts(j,i) == 0) then ! temperate ice base
293 
294  d_help_b(j,i) = cvxy1 ! basal sliding
295  ! (pressure-melting conditions)
296 
297  else ! n_cts(j,i) == 1, temperate ice layer
298 
299  d_help_b(j,i) = cvxy1 ! basal sliding
300  ! (pressure-melting conditions)
301 
302  end if
303 
304  else ! maske(j,i) == (1 oder 2); ice-free land or sea
305 
306  d_help_b(j,i) = 0.0_dp
307 
308  end if
309 
310 end do
311 end do
312 
313 !-------- Computation of vx_b (defined at (i+1/2,j)) --------
314 
315 do i=0, imax-1
316 do j=1, jmax-1
317  vx_b(j,i) = -0.5_dp*(d_help_b(j,i)+d_help_b(j,i+1))*dzs_dxi(j,i)
318 end do
319 end do
320 
321 !-------- Computation of vy_b (defined at (i,j+1/2)) --------
322 
323 do i=1, imax-1
324 do j=0, jmax-1
325  vy_b(j,i) = -0.5_dp*(d_help_b(j,i)+d_help_b(j+1,i))*dzs_deta(j,i)
326 end do
327 end do
328 
329 !-------- Computation of vx_b_g and vy_b_g (defined at (i,j)) --------
330 
331 do i=0, imax
332 do j=0, jmax
333  vx_b_g(j,i) = -d_help_b(j,i)*dzs_dxi_g(j,i)
334  vy_b_g(j,i) = -d_help_b(j,i)*dzs_deta_g(j,i)
335 end do
336 end do
337 
338 !-------- Limitation of computed vx_b, vy_b, vx_b_g, vy_b_g to the interval
339 ! [-VH_MAX, VH_MAX] --------
340 
341 vh_max = vh_max/year_sec
342 
343 do i=0, imax
344 do j=0, jmax
345  vx_b(j,i) = max(vx_b(j,i), -vh_max)
346  vx_b(j,i) = min(vx_b(j,i), vh_max)
347  vy_b(j,i) = max(vy_b(j,i), -vh_max)
348  vy_b(j,i) = min(vy_b(j,i), vh_max)
349  vx_b_g(j,i) = max(vx_b_g(j,i), -vh_max)
350  vx_b_g(j,i) = min(vx_b_g(j,i), vh_max)
351  vy_b_g(j,i) = max(vy_b_g(j,i), -vh_max)
352  vy_b_g(j,i) = min(vy_b_g(j,i), vh_max)
353 end do
354 end do
355 
356 end subroutine calc_vxy_b_sia
357 !
Generated on Mon May 20 2013 10:18:41 for SICOPOLIS V3.0 by   doxygen 1.8.1