SICOPOLIS V3.3
ice_material_properties_m.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Module : i c e _ m a t e r i a l _ p r o p e r t i e s _ m
4 !
5 !> @file
6 !!
7 !! Material properties of ice:
8 !! Rate factor, heat conductivity, specific heat (heat capacity),
9 !! creep function, viscosity.
10 !!
11 !! @section Copyright
12 !!
13 !! Copyright 2009-2017 Ralf Greve
14 !!
15 !! @section License
16 !!
17 !! This file is part of SICOPOLIS.
18 !!
19 !! SICOPOLIS is free software: you can redistribute it and/or modify
20 !! it under the terms of the GNU General Public License as published by
21 !! the Free Software Foundation, either version 3 of the License, or
22 !! (at your option) any later version.
23 !!
24 !! SICOPOLIS is distributed in the hope that it will be useful,
25 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
26 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
27 !! GNU General Public License for more details.
28 !!
29 !! You should have received a copy of the GNU General Public License
30 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
31 !<
32 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
33 
34 !-------------------------------------------------------------------------------
35 !> Material properties of ice:
36 !! Rate factor, heat conductivity, specific heat (heat capacity),
37 !! creep function, viscosity.
38 !<------------------------------------------------------------------------------
40 
41 use sico_types_m
42 
43 implicit none
44 save
45 
46 !> RF(n): Tabulated values for the rate factor of cold ice
47  real(dp), dimension(-256:255), private :: rf
48 
49 !> R_T: Coefficient of the water-content dependence in the rate factor
50 !> for temperate ice
51  real(dp) , private :: r_t
52 
53 !> KAPPA(n): Tabulated values for the heat conductivity of ice
54  real(dp), dimension(-256:255), private :: kappa
55 
56 !> C(n): Tabulated values for the specific heat of ice
57  real(dp), dimension(-256:255), private :: c
58 
59 !> n_temp_min: Lower index limit of properly defined values in RF, KAPPA and C
60 !> (n_temp_min >= -256).
61  integer(i4b), private :: n_temp_min
62 
63 !> n_temp_max: Upper index limit of properly defined values in RF, KAPPA and C
64 !> (n_temp_max <= 255).
65  integer(i4b), private :: n_temp_max
66 
67 !> RHO_I: Density of ice
68  real(dp) , private :: rho_i
69  ! only for the Martian ice caps
70 !> RHO_C: Density of crustal material (dust)
71  real(dp) , private :: rho_c
72  ! only for the Martian ice caps
73 !> KAPPA_C: Heat conductivity of crustal material (dust)
74  real(dp) , private :: kappa_c
75  ! only for the Martian ice caps
76 !> C_C: Specific heat of crustal material (dust)
77  real(dp) , private :: c_c
78  ! only for the Martian ice caps
79 
80 private
81 public :: ice_mat_eqs_pars, &
84 
85 contains
86 
87 !-------------------------------------------------------------------------------
88 !> Setting of required physical parameters.
89 !<------------------------------------------------------------------------------
90 subroutine ice_mat_eqs_pars(RF_table, R_T_val, KAPPA_table, C_table, &
91  n_tmp_min, n_tmp_max, &
92  rho_i_val, rho_c_val, kappa_c_val, c_c_val)
93 
94 implicit none
95 
96 integer(i4b), intent(in) :: n_tmp_min, n_tmp_max
97 real(dp), dimension(n_tmp_min:n_tmp_max), &
98  intent(in) :: RF_table, KAPPA_table, C_table
99 real(dp), intent(in) :: R_T_val
100 real(dp), optional, intent(in) :: RHO_I_val, RHO_C_val, KAPPA_C_val, C_C_val
101 
102 integer(i4b) :: n
103 
104 !-------- Initialisation --------
105 
106 rf = 0.0_dp
107 kappa = 0.0_dp
108 c = 0.0_dp
109 
110 n_temp_min = n_tmp_min
111 n_temp_max = n_tmp_max
112 
113 if ((n_temp_min <= -256).or.(n_temp_max >= 255)) &
114  stop ' >>> ice_mat_eqs_pars: Temperature indices out of allowed range!'
115 
116 !-------- Assignment --------
117 
118 do n=n_temp_min, n_temp_max
119  rf(n) = rf_table(n)
120  kappa(n) = kappa_table(n)
121  c(n) = c_table(n)
122 end do
123 
124 do n=-256, n_temp_min-1
125  rf(n) = rf(n_temp_min) ! dummy values
126  kappa(n) = kappa(n_temp_min) ! dummy values
127  c(n) = c(n_temp_min) ! dummy values
128 end do
129 
130 do n=n_temp_max+1, 255
131  rf(n) = rf(n_temp_max) ! dummy values
132  kappa(n) = kappa(n_temp_max) ! dummy values
133  c(n) = c(n_temp_max) ! dummy values
134 end do
135 
136 r_t = r_t_val
137 
138 !-------- Martian stuff --------
139 
140 if ( present(rho_i_val) ) then
141  rho_i = rho_i_val
142 else
143  rho_i = 0.0_dp ! dummy value
144 end if
145 
146 if ( present(rho_c_val) ) then
147  rho_c = rho_c_val
148 else
149  rho_c = 0.0_dp ! dummy value
150 end if
151 
152 if ( present(kappa_c_val) ) then
153  kappa_c = kappa_c_val
154 else
155  kappa_c = 0.0_dp ! dummy value
156 end if
157 
158 if ( present(c_c_val) ) then
159  c_c = c_c_val
160 else
161  c_c = 0.0_dp ! dummy value
162 end if
163 
164 end subroutine ice_mat_eqs_pars
165 
166 !-------------------------------------------------------------------------------
167 !> Rate factor for cold ice:
168 !! Linear interpolation of tabulated values in RF(.).
169 !<------------------------------------------------------------------------------
170 function ratefac_c(temp_val, temp_m_val)
172 use sico_types_m
174 use sico_vars_m
175 
176 implicit none
177 real(dp) :: ratefac_c
178 real(dp), intent(in) :: temp_val, temp_m_val
179 
180 integer(i4b) :: n_temp_1, n_temp_2
181 real(dp) :: temp_h_val
182 
183 temp_h_val = temp_val-temp_m_val
184 
185 n_temp_1 = floor(temp_h_val)
186 n_temp_1 = max(min(n_temp_1, n_temp_max-1), n_temp_min)
187 n_temp_2 = n_temp_1 + 1
188 
189 ratefac_c = rf(n_temp_1) &
190  + (rf(n_temp_2)-rf(n_temp_1)) &
191  * (temp_h_val-real(n_temp_1,dp)) ! Linear interpolation
192 
193 end function ratefac_c
194 
195 !-------------------------------------------------------------------------------
196 !> Rate factor for temperate ice.
197 !<------------------------------------------------------------------------------
198 function ratefac_t(omega_val)
200 use sico_types_m
202 use sico_vars_m
203 
204 implicit none
205 real(dp) :: ratefac_t
206 real(dp), intent(in) :: omega_val
207 
208 ratefac_t = rf(0)*(1.0_dp+r_t*(omega_val))
209 
210 end function ratefac_t
211 
212 !-------------------------------------------------------------------------------
213 !> Rate factor for cold and temperate ice:
214 !! Combination of ratefac_c and ratefac_t (only for the enthalpy method).
215 !<------------------------------------------------------------------------------
216 function ratefac_c_t(temp_val, omega_val, temp_m_val)
218 use sico_types_m
220 use sico_vars_m
221 
222 implicit none
223 real(dp) :: ratefac_c_t
224 real(dp), intent(in) :: temp_val, temp_m_val, omega_val
225 
226 integer(i4b) :: n_temp_1, n_temp_2
227 real(dp) :: temp_h_val
228 
229 temp_h_val = temp_val-temp_m_val
230 
231 n_temp_1 = floor(temp_h_val)
232 n_temp_1 = max(min(n_temp_1, n_temp_max-1), n_temp_min)
233 n_temp_2 = n_temp_1 + 1
234 
235 ratefac_c_t = ( rf(n_temp_1) &
236  + (rf(n_temp_2)-rf(n_temp_1)) &
237  * (temp_h_val-real(n_temp_1,dp)) ) &
238  * (1.0_dp+r_t*(omega_val))
239 
240 end function ratefac_c_t
241 
242 !-------------------------------------------------------------------------------
243 !> Heat conductivity of ice:
244 !! Linear interpolation of tabulated values in KAPPA(.).
245 !<------------------------------------------------------------------------------
246 function kappa_val(temp_val)
248 use sico_types_m
250 use sico_vars_m
251 
252 implicit none
253 real(dp) :: kappa_val
254 real(dp), intent(in) :: temp_val
255 
256 integer(i4b) :: n_temp_1, n_temp_2
257 real(dp) :: kappa_ice
258 
259 !-------- Heat conductivity of pure ice --------
260 
261 n_temp_1 = floor(temp_val)
262 n_temp_1 = max(min(n_temp_1, n_temp_max-1), n_temp_min)
263 n_temp_2 = n_temp_1 + 1
264 
265 #if defined(FRAC_DUST)
266 kappa_ice = kappa(n_temp_1) &
267  + (kappa(n_temp_2)-kappa(n_temp_1)) &
268  * (temp_val-real(n_temp_1,dp))
269 #else
270 kappa_val = kappa(n_temp_1) &
271  + (kappa(n_temp_2)-kappa(n_temp_1)) &
272  * (temp_val-real(n_temp_1,dp))
273 #endif
274 
275 !-------- If dust is present (polar caps of Mars):
276 ! Heat conductivity of ice-dust mixture --------
277 
278 #if defined(FRAC_DUST)
279 kappa_val = (1.0_dp-frac_dust)*kappa_ice + frac_dust*kappa_c
280 #endif
281 
282 end function kappa_val
283 
284 !-------------------------------------------------------------------------------
285 !> Specific heat of ice:
286 !! Linear interpolation of tabulated values in C(.).
287 !<------------------------------------------------------------------------------
288 function c_val(temp_val)
290 use sico_types_m
292 use sico_vars_m
293 
294 implicit none
295 real(dp) :: c_val
296 real(dp), intent(in) :: temp_val
297 
298 integer(i4b) :: n_temp_1, n_temp_2
299 real(dp) :: c_ice
300 
301 !-------- Specific heat of pure ice --------
302 
303 n_temp_1 = floor(temp_val)
304 n_temp_1 = max(min(n_temp_1, n_temp_max-1), n_temp_min)
305 n_temp_2 = n_temp_1 + 1
306 
307 #if defined(FRAC_DUST)
308 c_ice = c(n_temp_1) &
309  + (c(n_temp_2)-c(n_temp_1)) &
310  * (temp_val-real(n_temp_1,dp))
311 #else
312 c_val = c(n_temp_1) &
313  + (c(n_temp_2)-c(n_temp_1)) &
314  * (temp_val-real(n_temp_1,dp))
315 #endif
316 
317 !-------- If dust is present (polar caps of Mars):
318 ! Specific heat of ice-dust mixture --------
319 
320 #if defined(FRAC_DUST)
321 c_val = rho_inv * ( (1.0_dp-frac_dust)*rho_i*c_ice + frac_dust*rho_c*c_c )
322 #endif
323 
324 end function c_val
325 
326 !-------------------------------------------------------------------------------
327 !> Creep response function for ice.
328 !<------------------------------------------------------------------------------
329 function creep(sigma_val)
331 use sico_types_m
333 use sico_vars_m
334 
335 implicit none
336 
337 real(dp) :: creep
338 real(dp), intent(in) :: sigma_val
339 
340 #if (FLOW_LAW==4)
341 real(dp), parameter :: sm_coeff_1 = 8.5112e-15_dp, & ! s^-1 Pa^-1
342  sm_coeff_2 = 8.1643e-25_dp, & ! s^-1 Pa^-3
343  sm_coeff_3 = 9.2594e-12_dp ! Pa^-2
344 #endif
345 
346 #if (FIN_VISC==1)
347 
348 #if (FLOW_LAW==1)
349 
350 creep = sigma_val*sigma_val
351 ! Glen's flow law (n=3)
352 
353 #elif (FLOW_LAW==2)
354 
355 creep = sigma_val**0.8_dp * gr_size**(-1.4_dp)
356 ! Goldsby-Kohlstedt flow law (n=1.8, d=1.4)
357 
358 #elif (FLOW_LAW==3)
359 
360 creep = sigma_val*sigma_val*sigma_val
361 ! Durham's flow law (n=4)
362 
363 #endif
364 
365 #elif (FIN_VISC==2)
366 
367 #if (FLOW_LAW==1)
368 
369 creep = sigma_val*sigma_val + sigma_res*sigma_res
370 ! Glen's flow (n=3) with additional finite viscosity
371 
372 #elif (FLOW_LAW==2)
373 
374 creep = (sigma_val**0.8_dp + sigma_res**0.8_dp) * gr_size**(-1.4_dp)
375 ! Goldsby-Kohlstedt flow law (n=1.8, d=1.4)
376 ! with additional finite viscosity
377 
378 #elif (FLOW_LAW==3)
379 
380 creep = sigma_val*sigma_val*sigma_val + sigma_res*sigma_res*sigma_res
381 ! Durham's flow law (n=4) with additional finite viscosity
382 
383 #endif
384 
385 #endif
386 
387 #if (FLOW_LAW==4)
388 
389 creep = sm_coeff_1 &
390  + sm_coeff_2 * (sigma_val*sigma_val) &
391  * (1.0_dp + sm_coeff_3 * (sigma_val*sigma_val))
392 ! Smith-Morland (polynomial) flow law, normalised to a dimensionless
393 ! rate factor with A(-10C)=1.
394 
395 #endif
396 
397 end function creep
398 
399 !-------------------------------------------------------------------------------
400 !> Ice viscosity as a function of the effective strain rate and the temperature
401 !! (in cold ice) or the water content (in temperate ice) or both (for the
402 !! enthalpy method).
403 !<------------------------------------------------------------------------------
404 function viscosity(de_val, temp_val, temp_m_val, omega_val, enh_val, &
405  i_flag_cold_temp)
407 use sico_types_m
409 use sico_vars_m
410 
411 implicit none
412 
413 real(dp) :: viscosity
414 real(dp) , intent(in) :: de_val
415 real(dp) , intent(in) :: temp_val, temp_m_val
416 real(dp) , intent(in) :: omega_val
417 real(dp) , intent(in) :: enh_val
418 integer(i2b), intent(in) :: i_flag_cold_temp
419 
420 real(dp) :: ratefac_val
421 real(dp) :: de_val_m
422 
423 real(dp) :: n_power_law, inv_n_power_law, n_grain_size
424 
425 real(dp), parameter :: de_min = 1.0e-30_dp ! minimum value for the
426  ! effective strain rate
427 
428 #if (FLOW_LAW==4)
429 real(dp), parameter :: sm_coeff_1 = 8.5112e-15_dp, & ! s^-1 Pa^-1
430  sm_coeff_2 = 8.1643e-25_dp, & ! s^-1 Pa^-3
431  sm_coeff_3 = 9.2594e-12_dp ! Pa^-2
432 #endif
433 
434 !-------- Rate factor and effective strain rate --------
435 
436 if (i_flag_cold_temp == 0_i2b) then ! cold ice
437  ratefac_val = ratefac_c(temp_val, temp_m_val)
438 else if (i_flag_cold_temp == 1_i2b) then ! temperate ice
439  ratefac_val = ratefac_t(omega_val)
440 else ! enthalpy method
441  ratefac_val = ratefac_c_t(temp_val, omega_val, temp_m_val)
442 end if
443 
444 de_val_m = max(de_val, de_min)
445 
446 #if (FIN_VISC==1)
447 
448 #if (FLOW_LAW==1)
449 
450 !-------- Glen's flow law (n=3) --------
451 
452 inv_n_power_law = 0.333333333333333_dp ! 1/3
453 
454 viscosity = 0.5_dp * de_val_m**(inv_n_power_law-1.0_dp) &
455  * (enh_val*ratefac_val)**(-inv_n_power_law)
456 
457 #elif (FLOW_LAW==2)
458 
459 !-------- Goldsby-Kohlstedt flow law (n=1.8, d=1.4) --------
460 
461 inv_n_power_law = 0.555555555555555_dp ! 1/1.8
462 n_grain_size = 1.4_dp
463 
464 viscosity = 0.5_dp * de_val_m**(inv_n_power_law-1.0_dp) &
465  * gr_size**(n_grain_size*inv_n_power_law) &
466  * (enh_val*ratefac_val)**(-inv_n_power_law)
467 
468 #elif (FLOW_LAW==3)
469 
470 !-------- Durham's flow law (n=4) --------
471 
472 inv_n_power_law = 0.25_dp ! 1/4
473 
474 viscosity = 0.5_dp * de_val_m**(inv_n_power_law-1.0_dp) &
475  * (enh_val*ratefac_val)**(-inv_n_power_law)
476 
477 #endif
478 
479 #elif (FIN_VISC==2)
480 
481 #if (FLOW_LAW==1)
482 
483 !-------- Glen's flow (n=3) with additional finite viscosity --------
484 
485 n_power_law = 3.0_dp
486 
487 viscosity = visc_iter(de_val_m, ratefac_val, enh_val, n_power_law, sigma_res)
488 
489 #elif (FLOW_LAW==2)
490 
491 !-------- Goldsby-Kohlstedt flow law (n=1.8, d=1.4)
492 ! with additional finite viscosity --------
493 
494 n_power_law = 1.8_dp
495 n_grain_size = 1.4_dp
496 
497 write(6,'(a)') ' >>> viscosity: Computation of the viscosity as a function'
498 write(6,'(a)') ' of the effective strain rate not yet implemented'
499 write(6,'(a)') ' for grain-size-dependent finite viscosity flow laws!'
500 stop
501 
502 #elif (FLOW_LAW==3)
503 
504 !-------- Durham's flow law (n=4) with additional finite viscosity --------
505 
506 n_power_law = 4.0_dp
507 
508 viscosity = visc_iter(de_val_m, ratefac_val, enh_val, n_power_law, sigma_res)
509 
510 #endif
511 
512 #endif
513 
514 #if (FLOW_LAW==4)
515 
516 !-------- Smith-Morland (polynomial) flow law --------
517 
518 viscosity = visc_iter_sm(de_val_m, ratefac_val, enh_val, &
519  sm_coeff_1, sm_coeff_2, sm_coeff_3)
520 
521 #endif
522 
523 end function viscosity
524 
525 #if (FIN_VISC==2)
526 
527 !-------------------------------------------------------------------------------
528 !> Iterative computation of the viscosity by solving equation (4.28)
529 !! by Greve and Blatter (Springer, 2009).
530 !<------------------------------------------------------------------------------
531 function visc_iter(de_val_m, ratefac_val, enh_val, n_power_law, sigma_res)
532 
533 use sico_types_m
534 
535 implicit none
536 
537 real(dp) :: visc_iter
538 real(dp), intent(in) :: de_val_m
539 real(dp), intent(in) :: ratefac_val, enh_val
540 real(dp), intent(in) :: n_power_law, sigma_res
541 
542 real(dp) :: visc_val, res
543 
544 real(dp), parameter :: eps = 1.0e-05_dp ! convergence parameter
545 
546 !-------- Determination of the order of magnitude --------
547 
548 visc_val = 1.0e+10_dp ! initial guess (very low value)
549 
550 do
551 
552  res = fct_visc(de_val_m, ratefac_val, enh_val, visc_val, &
553  n_power_law, sigma_res)
554 
555  if (res < 0.0_dp) then
556  visc_val = 10.0_dp*visc_val
557  else
558  exit
559  end if
560 
561 end do
562 
563 !-------- Newton's method --------
564 
565 do
566 
567  visc_val = visc_val &
568  - res &
569  /fct_visc_deriv(de_val_m, ratefac_val, enh_val, visc_val, &
570  n_power_law, sigma_res)
571 
572  res = fct_visc(de_val_m, ratefac_val, enh_val, visc_val, &
573  n_power_law, sigma_res)
574 
575  if (abs(res) < eps) exit
576 
577 end do
578 
579 visc_iter = visc_val
580 
581 end function visc_iter
582 
583 !-------------------------------------------------------------------------------
584 !> Viscosity polynomial
585 !! [equation (4.28) by Greve and Blatter (Springer, 2009)].
586 !<------------------------------------------------------------------------------
587 function fct_visc(de_val_m, ratefac_val, enh_val, visc_val, &
588  n_power_law, sigma_res)
589 
590 use sico_types_m
591 
592 implicit none
593 
594 real(dp) :: fct_visc
595 real(dp), intent(in) :: de_val_m
596 real(dp), intent(in) :: ratefac_val, enh_val
597 real(dp), intent(in) :: visc_val
598 real(dp), intent(in) :: n_power_law, sigma_res
599 
600 fct_visc = 2.0_dp**n_power_law &
601  *enh_val*ratefac_val &
602  *de_val_m**(n_power_law-1.0_dp) &
603  *visc_val**n_power_law &
604  + 2.0_dp*enh_val*ratefac_val &
605  *sigma_res**(n_power_law-1.0_dp) &
606  *visc_val &
607  - 1.0_dp
608 
609 end function fct_visc
610 
611 !-------------------------------------------------------------------------------
612 !> Derivative of the viscosity polynomial
613 !! [equation (4.28) by Greve and Blatter (Springer, 2009)].
614 !<------------------------------------------------------------------------------
615 function fct_visc_deriv(de_val_m, ratefac_val, enh_val, visc_val, &
616  n_power_law, sigma_res)
617 
618 use sico_types_m
619 
620 implicit none
621 
622 real(dp) :: fct_visc_deriv
623 real(dp), intent(in) :: de_val_m
624 real(dp), intent(in) :: ratefac_val, enh_val
625 real(dp), intent(in) :: visc_val
626 real(dp), intent(in) :: n_power_law, sigma_res
627 
628 fct_visc_deriv = 2.0_dp**n_power_law*n_power_law &
629  *enh_val*ratefac_val &
630  *de_val_m**(n_power_law-1.0_dp) &
631  *visc_val**(n_power_law-1.0_dp) &
632  + 2.0_dp*enh_val*ratefac_val &
633  *sigma_res**(n_power_law-1.0_dp)
634 
635 end function fct_visc_deriv
636 
637 #endif
638 
639 #if (FLOW_LAW==4)
640 
641 !-------------------------------------------------------------------------------
642 !> Iterative computation of the viscosity by solving equation (4.33)
643 !! [analogous to (4.28)] by Greve and Blatter (Springer, 2009).
644 !<------------------------------------------------------------------------------
645 function visc_iter_sm(de_val_m, ratefac_val, enh_val, &
646  sm_coeff_1, sm_coeff_2, sm_coeff_3)
647 
648 use sico_types_m
649 
650 implicit none
651 
652 real(dp) :: visc_iter_sm
653 real(dp), intent(in) :: de_val_m
654 real(dp), intent(in) :: ratefac_val, enh_val
655 real(dp), intent(in) :: sm_coeff_1, sm_coeff_2, sm_coeff_3
656 
657 real(dp) :: visc_val, res
658 
659 real(dp), parameter :: eps = 1.0e-05_dp ! convergence parameter
660 
661 !-------- Determination of the order of magnitude --------
662 
663 visc_val = 1.0e+10_dp ! initial guess (very low value)
664 
665 do
666 
667  res = fct_visc_sm(de_val_m, ratefac_val, enh_val, visc_val, &
668  sm_coeff_1, sm_coeff_2, sm_coeff_3)
669 
670  if (res < 0.0_dp) then
671  visc_val = 10.0_dp*visc_val
672  else
673  exit
674  end if
675 
676 end do
677 
678 !-------- Newton's method --------
679 
680 do
681 
682  visc_val = visc_val &
683  - res &
684  /fct_visc_sm_deriv(de_val_m, ratefac_val, enh_val, visc_val, &
685  sm_coeff_1, sm_coeff_2, sm_coeff_3)
686 
687  res = fct_visc_sm(de_val_m, ratefac_val, enh_val, visc_val, &
688  sm_coeff_1, sm_coeff_2, sm_coeff_3)
689 
690  if (abs(res) < eps) exit
691 
692 end do
693 
694 visc_iter_sm = visc_val
695 
696 end function visc_iter_sm
697 
698 !-------------------------------------------------------------------------------
699 !> Viscosity polynomial
700 !! [equation (4.33) by Greve and Blatter (Springer, 2009)].
701 !<------------------------------------------------------------------------------
702 function fct_visc_sm(de_val_m, ratefac_val, enh_val, visc_val, &
703  sm_coeff_1, sm_coeff_2, sm_coeff_3)
704 
705 use sico_types_m
706 
707 implicit none
708 
709 real(dp) :: fct_visc_sm
710 real(dp), intent(in) :: de_val_m
711 real(dp), intent(in) :: ratefac_val, enh_val
712 real(dp), intent(in) :: visc_val
713 real(dp), intent(in) :: sm_coeff_1, sm_coeff_2, sm_coeff_3
714 
715 real(dp) :: de_visc_factor
716 
717 de_visc_factor = de_val_m*de_val_m*visc_val*visc_val
718 
719 fct_visc_sm = 2.0_dp*enh_val*ratefac_val*visc_val &
720  * ( sm_coeff_1 &
721  + 4.0_dp*sm_coeff_2*de_visc_factor &
722  * ( 1.0_dp + 4.0_dp*sm_coeff_3*de_visc_factor ) ) &
723  - 1.0_dp
724 
725 end function fct_visc_sm
726 
727 !-------------------------------------------------------------------------------
728 !> Derivative of the viscosity polynomial
729 !! [equation (4.33) by Greve and Blatter (Springer, 2009)].
730 !<------------------------------------------------------------------------------
731 function fct_visc_sm_deriv(de_val_m, ratefac_val, enh_val, visc_val, &
732  sm_coeff_1, sm_coeff_2, sm_coeff_3)
733 
734 use sico_types_m
735 
736 implicit none
737 
738 real(dp) :: fct_visc_sm_deriv
739 real(dp), intent(in) :: de_val_m
740 real(dp), intent(in) :: ratefac_val, enh_val
741 real(dp), intent(in) :: visc_val
742 real(dp), intent(in) :: sm_coeff_1, sm_coeff_2, sm_coeff_3
743 
744 real(dp) :: de_visc_factor
745 
746 real(dp), parameter :: twenty_over_three = 6.666666666666667_dp
747 
748 de_visc_factor = de_val_m*de_val_m*visc_val*visc_val
749 
750 fct_visc_sm_deriv = 2.0_dp*sm_coeff_1*enh_val*ratefac_val &
751  + 24.0_dp*sm_coeff_2*enh_val*ratefac_val*de_visc_factor &
752  * ( 1.0_dp + twenty_over_three*sm_coeff_3*de_visc_factor )
753 
754 end function fct_visc_sm_deriv
755 
756 #endif
757 
758 !-------------------------------------------------------------------------------
759 
760 end module ice_material_properties_m
761 !
Material properties of ice: Rate factor, heat conductivity, specific heat (heat capacity), creep function, viscosity.
real(dp) function, public kappa_val(temp_val)
Heat conductivity of ice: Linear interpolation of tabulated values in KAPPA(.).
real(dp) function, public creep(sigma_val)
Creep response function for ice.
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars_m.F90:35
real(dp) function, public ratefac_c(temp_val, temp_m_val)
Rate factor for cold ice: Linear interpolation of tabulated values in RF(.).
Declarations of kind types for SICOPOLIS.
real(dp) function, public ratefac_t(omega_val)
Rate factor for temperate ice.
real(dp) function, public c_val(temp_val)
Specific heat of ice: Linear interpolation of tabulated values in C(.).
real(dp) rho_inv
rho_inv: Inverse of the density of ice-dust mixture
Definition: sico_vars_m.F90:85
subroutine, public ice_mat_eqs_pars(RF_table, R_T_val, KAPPA_table, C_table, n_tmp_min, n_tmp_max, RHO_I_val, RHO_C_val, KAPPA_C_val, C_C_val)
Setting of required physical parameters.
real(dp) function, public 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...
integer, parameter dp
Double-precision reals.
Declarations of global variables for SICOPOLIS.