39 function viscosity(de_val, temp_val, temp_m_val, omega_val, enh_val, &
50 real(dp) ,
intent(in) :: de_val
51 real(dp) ,
intent(in) :: temp_val, temp_m_val
52 real(dp) ,
intent(in) :: omega_val
53 real(dp) ,
intent(in) :: enh_val
54 integer(i2b),
intent(in) :: i_flag_cold_temp
56 real(dp) :: ratefac_val
59 real(dp) :: n_power_law, inv_n_power_law, n_grain_size
61 real(dp),
parameter :: de_min = 1.0e-30_dp
65 real(dp),
parameter :: sm_coeff_1 = 8.5112e-15_dp, &
66 sm_coeff_2 = 8.1643e-25_dp, &
67 sm_coeff_3 = 9.2594e-12_dp
72 if (i_flag_cold_temp == 0_i2b)
then
73 ratefac_val =
ratefac_c(temp_val, temp_m_val)
74 else if (i_flag_cold_temp == 1_i2b)
then
77 ratefac_val =
ratefac_c_t(temp_val, omega_val, temp_m_val)
80 de_val_m = max(de_val, de_min)
88 inv_n_power_law = 0.333333333333333_dp
90 viscosity = 0.5_dp * de_val_m**(inv_n_power_law-1.0_dp) &
91 * (enh_val*ratefac_val)**(-inv_n_power_law)
97 inv_n_power_law = 0.555555555555555_dp
100 viscosity = 0.5_dp * de_val_m**(inv_n_power_law-1.0_dp) &
101 * gr_size**(n_grain_size*inv_n_power_law) &
102 * (enh_val*ratefac_val)**(-inv_n_power_law)
108 inv_n_power_law = 0.25_dp
110 viscosity = 0.5_dp * de_val_m**(inv_n_power_law-1.0_dp) &
111 * (enh_val*ratefac_val)**(-inv_n_power_law)
123 viscosity = visc_iter(de_val_m, ratefac_val, enh_val, n_power_law, sigma_res)
131 n_grain_size = 1.4_dp
133 write(6,
'(a)')
' viscosity: Computation of the viscosity as a function'
134 write(6,
'(a)')
' of the effective strain rate not yet implemented'
135 write(6,
'(a)')
' for grain-size-dependent finite viscosity flow laws!'
144 viscosity = visc_iter(de_val_m, ratefac_val, enh_val, n_power_law, sigma_res)
154 viscosity = visc_iter_sm(de_val_m, ratefac_val, enh_val, &
155 sm_coeff_1, sm_coeff_2, sm_coeff_3)
167 function visc_iter(de_val_m, ratefac_val, enh_val, n_power_law, sigma_res)
173 real(dp) :: visc_iter
174 real(dp),
intent(in) :: de_val_m
175 real(dp),
intent(in) :: ratefac_val, enh_val
176 real(dp),
intent(in) :: n_power_law, sigma_res
178 real(dp) :: visc_val, res
180 real(dp),
parameter :: eps = 1.0e-05_dp
184 visc_val = 1.0e+10_dp
188 res = fct_visc(de_val_m, ratefac_val, enh_val, visc_val, &
189 n_power_law, sigma_res)
191 if (res < 0.0_dp)
then
192 visc_val = 10.0_dp*visc_val
203 visc_val = visc_val &
205 /fct_visc_deriv(de_val_m, ratefac_val, enh_val, visc_val, &
206 n_power_law, sigma_res)
208 res = fct_visc(de_val_m, ratefac_val, enh_val, visc_val, &
209 n_power_law, sigma_res)
211 if (abs(res) < eps)
exit
217 end function visc_iter
223 function fct_visc(de_val_m, ratefac_val, enh_val, visc_val, &
224 n_power_law, sigma_res)
231 real(dp),
intent(in) :: de_val_m
232 real(dp),
intent(in) :: ratefac_val, enh_val
233 real(dp),
intent(in) :: visc_val
234 real(dp),
intent(in) :: n_power_law, sigma_res
236 fct_visc = 2.0_dp**n_power_law &
237 *enh_val*ratefac_val &
238 *de_val_m**(n_power_law-1.0_dp) &
239 *visc_val**n_power_law &
240 + 2.0_dp*enh_val*ratefac_val &
241 *sigma_res**(n_power_law-1.0_dp) &
245 end function fct_visc
251 function fct_visc_deriv(de_val_m, ratefac_val, enh_val, visc_val, &
252 n_power_law, sigma_res)
258 real(dp) :: fct_visc_deriv
259 real(dp),
intent(in) :: de_val_m
260 real(dp),
intent(in) :: ratefac_val, enh_val
261 real(dp),
intent(in) :: visc_val
262 real(dp),
intent(in) :: n_power_law, sigma_res
264 fct_visc_deriv = 2.0_dp**n_power_law*n_power_law &
265 *enh_val*ratefac_val &
266 *de_val_m**(n_power_law-1.0_dp) &
267 *visc_val**(n_power_law-1.0_dp) &
268 + 2.0_dp*enh_val*ratefac_val &
269 *sigma_res**(n_power_law-1.0_dp)
271 end function fct_visc_deriv
281 function visc_iter_sm(de_val_m, ratefac_val, enh_val, &
282 sm_coeff_1, sm_coeff_2, sm_coeff_3)
288 real(dp) :: visc_iter_sm
289 real(dp),
intent(in) :: de_val_m
290 real(dp),
intent(in) :: ratefac_val, enh_val
291 real(dp),
intent(in) :: sm_coeff_1, sm_coeff_2, sm_coeff_3
293 real(dp) :: visc_val, res
295 real(dp),
parameter :: eps = 1.0e-05_dp
299 visc_val = 1.0e+10_dp
303 res = fct_visc_sm(de_val_m, ratefac_val, enh_val, visc_val, &
304 sm_coeff_1, sm_coeff_2, sm_coeff_3)
306 if (res < 0.0_dp)
then
307 visc_val = 10.0_dp*visc_val
318 visc_val = visc_val &
320 /fct_visc_sm_deriv(de_val_m, ratefac_val, enh_val, visc_val, &
321 sm_coeff_1, sm_coeff_2, sm_coeff_3)
323 res = fct_visc_sm(de_val_m, ratefac_val, enh_val, visc_val, &
324 sm_coeff_1, sm_coeff_2, sm_coeff_3)
326 if (abs(res) < eps)
exit
330 visc_iter_sm = visc_val
332 end function visc_iter_sm
338 function fct_visc_sm(de_val_m, ratefac_val, enh_val, visc_val, &
339 sm_coeff_1, sm_coeff_2, sm_coeff_3)
345 real(dp) :: fct_visc_sm
346 real(dp),
intent(in) :: de_val_m
347 real(dp),
intent(in) :: ratefac_val, enh_val
348 real(dp),
intent(in) :: visc_val
349 real(dp),
intent(in) :: sm_coeff_1, sm_coeff_2, sm_coeff_3
351 real(dp) :: de_visc_factor
353 de_visc_factor = de_val_m*de_val_m*visc_val*visc_val
355 fct_visc_sm = 2.0_dp*enh_val*ratefac_val*visc_val &
357 + 4.0_dp*sm_coeff_2*de_visc_factor &
358 * ( 1.0_dp + 4.0_dp*sm_coeff_3*de_visc_factor ) ) &
361 end function fct_visc_sm
367 function fct_visc_sm_deriv(de_val_m, ratefac_val, enh_val, visc_val, &
368 sm_coeff_1, sm_coeff_2, sm_coeff_3)
374 real(dp) :: fct_visc_sm_deriv
375 real(dp),
intent(in) :: de_val_m
376 real(dp),
intent(in) :: ratefac_val, enh_val
377 real(dp),
intent(in) :: visc_val
378 real(dp),
intent(in) :: sm_coeff_1, sm_coeff_2, sm_coeff_3
380 real(dp) :: de_visc_factor
382 real(dp),
parameter :: twenty_over_three = 6.666666666666667_dp
384 de_visc_factor = de_val_m*de_val_m*visc_val*visc_val
386 fct_visc_sm_deriv = 2.0_dp*sm_coeff_1*enh_val*ratefac_val &
387 + 24.0_dp*sm_coeff_2*enh_val*ratefac_val*de_visc_factor &
388 * ( 1.0_dp + twenty_over_three*sm_coeff_3*de_visc_factor )
390 end function fct_visc_sm_deriv
real(dp) function, public ratefac_t(omega_val)
Rate factor for temperate ice.
Declarations of kind types for SICOPOLIS.
real(dp) function viscosity(de_val, temp_val, temp_m_val, omega_val, enh_val, i_flag_cold_temp)
Ice viscosity as a function of the effective strain rate and the temperature (in cold ice) or the wat...
real(dp) function, public ratefac_c_t(temp_val, omega_val, temp_m_val)
Rate factor for cold and temperate ice: Combination of ratefac_c and ratefac_t (only for the enthalpy...
Declarations of global variables for SICOPOLIS (for the ANT domain).
real(dp) function, public ratefac_c(temp_val, temp_m_val)
Rate factor for cold ice: Linear interpolation of tabulated values in RF(.).
Material quantities of ice: Rate factor, heat conductivity, specific heat (heat capacity).
Declarations of global variables for SICOPOLIS.