37 function viscosity(de_val, temp_val, temp_m_val, omega_val, enh_val, &
46 real(dp),
intent(in) :: de_val
47 real(dp),
intent(in) :: temp_val, temp_m_val
48 real(dp),
intent(in) :: omega_val
49 real(dp),
intent(in) :: enh_val
50 logical,
intent(in) :: flag_cold_temp
52 real(dp) :: ratefac_val
55 real(dp) :: n_power_law, inv_n_power_law, n_grain_size
57 real(dp),
parameter :: de_min = 1.0e-30_dp
61 real(dp),
parameter :: sm_coeff_1 = 8.5112e-15_dp, &
62 sm_coeff_2 = 8.1643e-25_dp, &
63 sm_coeff_3 = 9.2594e-12_dp
68 if (flag_cold_temp)
then
69 ratefac_val =
ratefac(temp_val,temp_m_val)
74 de_val_m = max(de_val, de_min)
82 inv_n_power_law = 0.333333333333333_dp
84 viscosity = 0.5_dp * de_val_m**(inv_n_power_law-1.0_dp) &
85 * (enh_val*ratefac_val)**(-inv_n_power_law)
91 inv_n_power_law = 0.555555555555555_dp
94 viscosity = 0.5_dp * de_val_m**(inv_n_power_law-1.0_dp) &
95 * gr_size**(n_grain_size*inv_n_power_law) &
96 * (enh_val*ratefac_val)**(-inv_n_power_law)
102 inv_n_power_law = 0.25_dp
104 viscosity = 0.5_dp * de_val_m**(inv_n_power_law-1.0_dp) &
105 * (enh_val*ratefac_val)**(-inv_n_power_law)
117 viscosity = visc_iter(de_val_m, ratefac_val, enh_val, n_power_law, sigma_res)
125 n_grain_size = 1.4_dp
127 write(6,
'(a)')
' viscosity: Computation of the viscosity as a function'
128 write(6,
'(a)')
' of the effective strain rate not yet implemented'
129 write(6,
'(a)')
' for grain-size-dependent finite viscosity flow laws!'
138 viscosity = visc_iter(de_val_m, ratefac_val, enh_val, n_power_law, sigma_res)
148 viscosity = visc_iter_sm(de_val_m, ratefac_val, enh_val, &
149 sm_coeff_1, sm_coeff_2, sm_coeff_3)
161 function visc_iter(de_val_m, ratefac_val, enh_val, n_power_law, sigma_res)
167 real(dp) :: visc_iter
168 real(dp),
intent(in) :: de_val_m
169 real(dp),
intent(in) :: ratefac_val, enh_val
170 real(dp),
intent(in) :: n_power_law, sigma_res
172 real(dp) :: visc_val, res
174 real(dp),
parameter :: eps = 1.0e-05_dp
178 visc_val = 1.0e+10_dp
182 res = fct_visc(de_val_m, ratefac_val, enh_val, visc_val, &
183 n_power_law, sigma_res)
185 if (res < 0.0_dp)
then
186 visc_val = 10.0_dp*visc_val
197 visc_val = visc_val &
199 /fct_visc_deriv(de_val_m, ratefac_val, enh_val, visc_val, &
200 n_power_law, sigma_res)
202 res = fct_visc(de_val_m, ratefac_val, enh_val, visc_val, &
203 n_power_law, sigma_res)
205 if (abs(res) < eps)
exit
211 end function visc_iter
217 function fct_visc(de_val_m, ratefac_val, enh_val, visc_val, &
218 n_power_law, sigma_res)
225 real(dp),
intent(in) :: de_val_m
226 real(dp),
intent(in) :: ratefac_val, enh_val
227 real(dp),
intent(in) :: visc_val
228 real(dp),
intent(in) :: n_power_law, sigma_res
230 fct_visc = 2.0_dp**n_power_law &
231 *enh_val*ratefac_val &
232 *de_val_m**(n_power_law-1.0_dp) &
233 *visc_val**n_power_law &
234 + 2.0_dp*enh_val*ratefac_val &
235 *sigma_res**(n_power_law-1.0_dp) &
239 end function fct_visc
245 function fct_visc_deriv(de_val_m, ratefac_val, enh_val, visc_val, &
246 n_power_law, sigma_res)
252 real(dp) :: fct_visc_deriv
253 real(dp),
intent(in) :: de_val_m
254 real(dp),
intent(in) :: ratefac_val, enh_val
255 real(dp),
intent(in) :: visc_val
256 real(dp),
intent(in) :: n_power_law, sigma_res
258 fct_visc_deriv = 2.0_dp**n_power_law*n_power_law &
259 *enh_val*ratefac_val &
260 *de_val_m**(n_power_law-1.0_dp) &
261 *visc_val**(n_power_law-1.0_dp) &
262 + 2.0_dp*enh_val*ratefac_val &
263 *sigma_res**(n_power_law-1.0_dp)
265 end function fct_visc_deriv
275 function visc_iter_sm(de_val_m, ratefac_val, enh_val, &
276 sm_coeff_1, sm_coeff_2, sm_coeff_3)
282 real(dp) :: visc_iter_sm
283 real(dp),
intent(in) :: de_val_m
284 real(dp),
intent(in) :: ratefac_val, enh_val
285 real(dp),
intent(in) :: sm_coeff_1, sm_coeff_2, sm_coeff_3
287 real(dp) :: visc_val, res
289 real(dp),
parameter :: eps = 1.0e-05_dp
293 visc_val = 1.0e+10_dp
297 res = fct_visc_sm(de_val_m, ratefac_val, enh_val, visc_val, &
298 sm_coeff_1, sm_coeff_2, sm_coeff_3)
300 if (res < 0.0_dp)
then
301 visc_val = 10.0_dp*visc_val
312 visc_val = visc_val &
314 /fct_visc_sm_deriv(de_val_m, ratefac_val, enh_val, visc_val, &
315 sm_coeff_1, sm_coeff_2, sm_coeff_3)
317 res = fct_visc_sm(de_val_m, ratefac_val, enh_val, visc_val, &
318 sm_coeff_1, sm_coeff_2, sm_coeff_3)
320 if (abs(res) < eps)
exit
324 visc_iter_sm = visc_val
326 end function visc_iter_sm
332 function fct_visc_sm(de_val_m, ratefac_val, enh_val, visc_val, &
333 sm_coeff_1, sm_coeff_2, sm_coeff_3)
339 real(dp) :: fct_visc_sm
340 real(dp),
intent(in) :: de_val_m
341 real(dp),
intent(in) :: ratefac_val, enh_val
342 real(dp),
intent(in) :: visc_val
343 real(dp),
intent(in) :: sm_coeff_1, sm_coeff_2, sm_coeff_3
345 real(dp) :: de_visc_factor
347 de_visc_factor = de_val_m*de_val_m*visc_val*visc_val
349 fct_visc_sm = 2.0_dp*enh_val*ratefac_val*visc_val &
351 + 4.0_dp*sm_coeff_2*de_visc_factor &
352 * ( 1.0_dp + 4.0_dp*sm_coeff_3*de_visc_factor ) ) &
355 end function fct_visc_sm
361 function fct_visc_sm_deriv(de_val_m, ratefac_val, enh_val, visc_val, &
362 sm_coeff_1, sm_coeff_2, sm_coeff_3)
368 real(dp) :: fct_visc_sm_deriv
369 real(dp),
intent(in) :: de_val_m
370 real(dp),
intent(in) :: ratefac_val, enh_val
371 real(dp),
intent(in) :: visc_val
372 real(dp),
intent(in) :: sm_coeff_1, sm_coeff_2, sm_coeff_3
374 real(dp) :: de_visc_factor
376 real(dp),
parameter :: twenty_over_three = 6.666666666666667_dp
378 de_visc_factor = de_val_m*de_val_m*visc_val*visc_val
380 fct_visc_sm_deriv = 2.0_dp*sm_coeff_1*enh_val*ratefac_val &
381 + 24.0_dp*sm_coeff_2*enh_val*ratefac_val*de_visc_factor &
382 * ( 1.0_dp + twenty_over_three*sm_coeff_3*de_visc_factor )
384 end function fct_visc_sm_deriv