SICOPOLIS V5-dev  Revision 1368
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-2018 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 use error_m
43 
44 implicit none
45 save
46 
47 #ifndef ALLOW_OPENAD /* Normal */
48 
49 !> RF(n): Tabulated values for the rate factor of cold ice
50  real(dp), dimension(-256:255), private :: rf
51 
52 !> R_T: Coefficient of the water-content dependence in the rate factor
53 !> for temperate ice
54  real(dp) , private :: r_t
55 
56 !> KAPPA(n): Tabulated values for the heat conductivity of ice
57  real(dp), dimension(-256:255), private :: kappa
58 
59 !> C(n): Tabulated values for the specific heat of ice
60  real(dp), dimension(-256:255), private :: c
61 
62 !> n_temp_min: Lower index limit of properly defined values in RF, KAPPA and C
63 !> (n_temp_min >= -256).
64  integer(i4b), private :: n_temp_min
65 
66 !> n_temp_max: Upper index limit of properly defined values in RF, KAPPA and C
67 !> (n_temp_max <= 255).
68  integer(i4b), private :: n_temp_max
69 
70 !> RHO_I: Density of ice
71  real(dp) , private :: rho_i
72  ! only for the Martian ice caps
73 !> RHO_C: Density of crustal material (dust)
74  real(dp) , private :: rho_c
75  ! only for the Martian ice caps
76 !> KAPPA_C: Heat conductivity of crustal material (dust)
77  real(dp) , private :: kappa_c
78  ! only for the Martian ice caps
79 !> C_C: Specific heat of crustal material (dust)
80  real(dp) , private :: c_c
81  ! only for the Martian ice caps
82 
83 private
84 public :: ice_mat_eqs_pars, &
87 
88 #else /* OpenAD */
89 
90 real(dp), dimension(-256:255), public :: rf_imq
91 real(dp) , public :: r_t_imq
92 real(dp), dimension(-256:255), public :: kappa_imq
93 real(dp), dimension(-256:255), public :: c_imq
94 integer(i4b), public :: n_temp_min_imq
95 integer(i4b), public :: n_temp_max_imq
96 real(dp) , public :: rho_i_imq
97 real(dp) , public :: rho_c_imq
98 real(dp) , public :: kappa_c_imq
99 real(dp) , public :: c_c_imq
100 
101 public :: ice_mat_eqs_pars, &
104 private :: myfloor_local
105 
106 #endif /* Normal vs. OpenAD */
107 
108 contains
109 
110 !-------------------------------------------------------------------------------
111 !> Setting of required physical parameters.
112 !<------------------------------------------------------------------------------
113 subroutine ice_mat_eqs_pars(RF_table, R_T_val, KAPPA_table, C_table, &
114  n_tmp_min, n_tmp_max, &
115  rho_i_val, rho_c_val, kappa_c_val, c_c_val)
117 implicit none
118 
119 integer(i4b), intent(in) :: n_tmp_min, n_tmp_max
120 real(dp), dimension(n_tmp_min:n_tmp_max), &
121  intent(in) :: RF_table, KAPPA_table, C_table
122 real(dp), intent(in) :: R_T_val
123 real(dp), optional, intent(in) :: RHO_I_val, RHO_C_val, KAPPA_C_val, C_C_val
124 
125 integer(i4b) :: n
126 character(len=256) :: errormsgg
127 
128 !-------- Initialisation --------
129 
130 #ifndef ALLOW_OPENAD /* Normal */
131 
132 rf = 0.0_dp
133 kappa = 0.0_dp
134 c = 0.0_dp
135 n_temp_min = n_tmp_min
136 n_temp_max = n_tmp_max
137 
138 #else /* OpenAD */
139 
140 rf_imq = 0.0_dp
141 kappa_imq = 0.0_dp
142 c_imq = 0.0_dp
143 n_temp_min_imq = n_tmp_min
144 n_temp_max_imq = n_tmp_max
145 
146 #endif /* Normal vs. OpenAD */
147 
148 #ifndef ALLOW_OPENAD /* Normal */
149 if ((n_temp_min <= -256).or.(n_temp_max >= 255)) then
150 #else /* OpenAD */
151 if ((n_temp_min_imq <= -256).or.(n_temp_max_imq >= 255)) then
152 #endif /* Normal vs. OpenAD */
153  errormsgg = ' >>> ice_mat_eqs_pars: ' &
154  //'Temperature indices out of allowed range!'
155  call error(errormsgg)
156 end if
157 
158 !-------- Assignment --------
159 
160 #ifndef ALLOW_OPENAD /* Normal */
161 
162 do n=n_temp_min, n_temp_max
163  rf(n) = rf_table(n)
164  kappa(n) = kappa_table(n)
165  c(n) = c_table(n)
166 end do
167 
168 #else /* OpenAD */
169 
170 do n=n_temp_min_imq, n_temp_max_imq
171  rf_imq(n) = rf_table(n)
172  kappa_imq(n) = kappa_table(n)
173  c_imq(n) = c_table(n)
174 end do
175 
176 #endif /* Normal vs. OpenAD */
177 
178 #ifndef ALLOW_OPENAD /* Normal */
179 
180 do n=-256, n_temp_min-1
181  rf(n) = rf(n_temp_min) ! dummy values
182  kappa(n) = kappa(n_temp_min) ! dummy values
183  c(n) = c(n_temp_min) ! dummy values
184 end do
185 
186 #else /* OpenAD */
187 
188 do n=-256, n_temp_min_imq-1
189  rf_imq(n) = rf_imq(n_temp_min_imq) ! dummy values
190  kappa_imq(n) = kappa_imq(n_temp_min_imq) ! dummy values
191  c_imq(n) = c_imq(n_temp_min_imq) ! dummy values
192 end do
193 
194 #endif /* Normal vs. OpenAD */
195 
196 #ifndef ALLOW_OPENAD /* Normal */
197 
198 do n=n_temp_max+1, 255
199  rf(n) = rf(n_temp_max) ! dummy values
200  kappa(n) = kappa(n_temp_max) ! dummy values
201  c(n) = c(n_temp_max) ! dummy values
202 end do
203 
204 #else /* OpenAD */
205 
206 do n=n_temp_max_imq+1, 255
207  rf_imq(n) = rf_imq(n_temp_max_imq) ! dummy values
208  kappa_imq(n) = kappa_imq(n_temp_max_imq) ! dummy values
209  c_imq(n) = c_imq(n_temp_max_imq) ! dummy values
210 end do
211 
212 #endif /* Normal vs. OpenAD */
213 
214 #ifndef ALLOW_OPENAD /* Normal */
215 r_t = r_t_val
216 #else /* OpenAD */
217 r_t_imq = r_t_val
218 #endif /* Normal vs. OpenAD */
219 
220 !-------- Martian stuff --------
221 
222 #ifndef ALLOW_OPENAD /* Normal */
223 
224 if ( present(rho_i_val) ) then
225  rho_i = rho_i_val
226 else
227  rho_i = 0.0_dp ! dummy value
228 end if
229 
230 #else /* OpenAD */
231 
232 if ( present(rho_i_val) ) then
233  rho_i_imq = rho_i_val
234 else
235  rho_i_imq = 0.0_dp ! dummy value
236 end if
237 
238 #endif /* Normal vs. OpenAD */
239 
240 #ifndef ALLOW_OPENAD /* Normal */
241 
242 if ( present(rho_c_val) ) then
243  rho_c = rho_c_val
244 else
245  rho_c = 0.0_dp ! dummy value
246 end if
247 
248 #else /* OpenAD */
249 
250 if ( present(rho_c_val) ) then
251  rho_c_imq = rho_c_val
252 else
253  rho_c_imq = 0.0_dp ! dummy value
254 end if
255 
256 #endif /* Normal vs. OpenAD */
257 
258 #ifndef ALLOW_OPENAD /* Normal */
259 
260 if ( present(kappa_c_val) ) then
261  kappa_c = kappa_c_val
262 else
263  kappa_c = 0.0_dp ! dummy value
264 end if
265 
266 #else /* OpenAD */
267 
268 if ( present(kappa_c_val) ) then
269  kappa_c_imq = kappa_c_val
270 else
271  kappa_c_imq = 0.0_dp ! dummy value
272 end if
273 
274 #endif /* Normal vs. OpenAD */
275 
276 #ifndef ALLOW_OPENAD /* Normal */
277 
278 if ( present(c_c_val) ) then
279  c_c = c_c_val
280 else
281  c_c = 0.0_dp ! dummy value
282 end if
283 
284 #else /* OpenAD */
285 
286 if ( present(c_c_val) ) then
287  c_c_imq = c_c_val
288 else
289  c_c_imq = 0.0_dp ! dummy value
290 end if
291 
292 #endif /* Normal vs. OpenAD */
293 
294 end subroutine ice_mat_eqs_pars
295 
296 !-------------------------------------------------------------------------------
297 !> Rate factor for cold ice:
298 !! Linear interpolation of tabulated values in RF(.).
299 !<------------------------------------------------------------------------------
300 function ratefac_c(temp_val, temp_m_val)
303 use sico_vars_m
304 
305 implicit none
306 real(dp) :: ratefac_c
307 real(dp), intent(in) :: temp_val, temp_m_val
308 
309 integer(i4b) :: n_temp_1, n_temp_2
310 real(dp) :: temp_h_val
311 #ifdef ALLOW_OPENAD /* OpenAD */
312 integer(i4b) :: itemp_h_val
313 #endif /* OpenAD */
314 
315 temp_h_val = temp_val-temp_m_val
316 
317 #ifndef ALLOW_OPENAD /* Normal */
318 
319 n_temp_1 = floor(temp_h_val)
320 n_temp_1 = max(min(n_temp_1, n_temp_max-1), n_temp_min)
321 n_temp_2 = n_temp_1 + 1
322 
323 ratefac_c = rf(n_temp_1) &
324  + (rf(n_temp_2)-rf(n_temp_1)) &
325  * (temp_h_val-real(n_temp_1,dp)) ! Linear interpolation
326 
327 #else /* OpenAD */
328 
329 n_temp_1 = myfloor_local(temp_h_val)
330 n_temp_1 = max(min(n_temp_1, n_temp_max_imq-1), n_temp_min_imq)
331 n_temp_2 = n_temp_1 + 1
332 
333 ratefac_c = rf_imq(n_temp_1) &
334  +(rf_imq(n_temp_2)-rf_imq(n_temp_1)) &
335  *(temp_h_val-real(n_temp_1,dp))
336 
337 #endif /* Normal vs. OpenAD */
338 
339 end function ratefac_c
340 
341 !-------------------------------------------------------------------------------
342 !> Rate factor for temperate ice.
343 !<------------------------------------------------------------------------------
344 function ratefac_t(omega_val)
347 use sico_vars_m
348 
349 implicit none
350 real(dp) :: ratefac_t
351 real(dp), intent(in) :: omega_val
352 
353 #ifndef ALLOW_OPENAD /* Normal */
354 ratefac_t = rf(0)*(1.0_dp+r_t*(omega_val))
355 #else /* OpenAD */
356 ratefac_t = rf_imq(0)*(1.0_dp+r_t_imq*(omega_val))
357 #endif /* Normal vs. OpenAD */
358 
359 end function ratefac_t
360 
361 !-------------------------------------------------------------------------------
362 !> Rate factor for cold and temperate ice:
363 !! Combination of ratefac_c and ratefac_t (only for the enthalpy method).
364 !<------------------------------------------------------------------------------
365 function ratefac_c_t(temp_val, omega_val, temp_m_val)
368 use sico_vars_m
369 
370 implicit none
371 real(dp) :: ratefac_c_t
372 real(dp), intent(in) :: temp_val, temp_m_val, omega_val
373 
374 integer(i4b) :: n_temp_1, n_temp_2
375 real(dp) :: temp_h_val
376 #ifdef ALLOW_OPENAD /* OpenAD */
377 integer(i4b) :: itemp_h_val
378 #endif /* OpenAD */
379 
380 temp_h_val = temp_val-temp_m_val
381 
382 #ifndef ALLOW_OPENAD /* Normal */
383 
384 n_temp_1 = floor(temp_h_val)
385 n_temp_1 = max(min(n_temp_1, n_temp_max-1), n_temp_min)
386 n_temp_2 = n_temp_1 + 1
387 
388 ratefac_c_t = ( rf(n_temp_1) &
389  + (rf(n_temp_2)-rf(n_temp_1)) &
390  * (temp_h_val-real(n_temp_1,dp)) ) &
391  *(1.0_dp+r_t*(omega_val))
392 
393 #else /* OpenAD */
394 
395 n_temp_1 = myfloor_local(temp_h_val)
396 n_temp_1 = max(min(n_temp_1, n_temp_max_imq-1), n_temp_min_imq)
397 n_temp_2 = n_temp_1 + 1
398 
399 ratefac_c_t = ( rf_imq(n_temp_1) &
400  +(rf_imq(n_temp_2)-rf_imq(n_temp_1)) &
401  *(temp_h_val-real(n_temp_1,dp)) ) &
402  *(1.0_dp+r_t_imq*(omega_val))
403 
404 #endif /* Normal vs. OpenAD */
405 
406 end function ratefac_c_t
407 
408 !-------------------------------------------------------------------------------
409 !> Heat conductivity of ice:
410 !! Linear interpolation of tabulated values in KAPPA(.).
411 !<------------------------------------------------------------------------------
412 function kappa_val(temp_val)
415 use sico_vars_m
416 
417 implicit none
418 real(dp) :: kappa_val
419 real(dp), intent(in) :: temp_val
420 
421 integer(i4b) :: n_temp_1, n_temp_2
422 real(dp) :: kappa_ice
423 #ifdef ALLOW_OPENAD /* OpenAD */
424 integer(i4b) :: itemp_val
425 #endif /* OpenAD */
426 
427 !-------- Heat conductivity of pure ice --------
428 
429 #ifndef ALLOW_OPENAD /* Normal */
430 
431 n_temp_1 = floor(temp_val)
432 n_temp_1 = max(min(n_temp_1, n_temp_max-1), n_temp_min)
433 n_temp_2 = n_temp_1 + 1
434 
435 #if defined(FRAC_DUST)
436 kappa_ice = kappa(n_temp_1) &
437  + (kappa(n_temp_2)-kappa(n_temp_1)) &
438  * (temp_val-real(n_temp_1,dp))
439 #else
440 kappa_val = kappa(n_temp_1) &
441  + (kappa(n_temp_2)-kappa(n_temp_1)) &
442  * (temp_val-real(n_temp_1,dp))
443 #endif
444 
445 #else /* OpenAD */
446 
447 n_temp_1 = myfloor_local(temp_val)
448 n_temp_1 = max(min(n_temp_1, n_temp_max_imq-1), n_temp_min_imq)
449 n_temp_2 = n_temp_1 + 1
450 
451 #if defined(FRAC_DUST)
452 kappa_ice = kappa_imq(n_temp_1) &
453  + (kappa_imq(n_temp_2)-kappa_imq(n_temp_1)) &
454  * (temp_val-real(n_temp_1,dp))
455 #else
456 kappa_val = kappa_imq(n_temp_1) &
457  + (kappa_imq(n_temp_2)-kappa_imq(n_temp_1)) &
458  * (temp_val-real(n_temp_1,dp))
459 #endif
460 
461 #endif /* Normal vs. OpenAD */
462 
463 !-------- If dust is present (polar caps of Mars):
464 ! Heat conductivity of ice-dust mixture --------
465 
466 #if defined(FRAC_DUST)
467 #ifndef ALLOW_OPENAD /* Normal */
468 kappa_val = (1.0_dp-frac_dust)*kappa_ice + frac_dust*kappa_c
469 #else /* OpenAD */
470 kappa_val = (1.0_dp-frac_dust)*kappa_ice + frac_dust*kappa_c_imq
471 #endif /* Normal vs. OpenAD */
472 #endif
473 
474 end function kappa_val
475 
476 !-------------------------------------------------------------------------------
477 !> Specific heat of ice:
478 !! Linear interpolation of tabulated values in C(.).
479 !<------------------------------------------------------------------------------
480 function c_val(temp_val)
483 use sico_vars_m
484 
485 implicit none
486 real(dp) :: c_val
487 real(dp), intent(in) :: temp_val
488 
489 integer(i4b) :: n_temp_1, n_temp_2
490 real(dp) :: c_ice
491 #ifdef ALLOW_OPENAD /* OpenAD */
492 integer(i4b) :: itemp_val
493 #endif /* OpenAD */
494 
495 !-------- Specific heat of pure ice --------
496 
497 #ifndef ALLOW_OPENAD /* Normal */
498 
499 n_temp_1 = floor(temp_val)
500 n_temp_1 = max(min(n_temp_1, n_temp_max-1), n_temp_min)
501 n_temp_2 = n_temp_1 + 1
502 
503 #if defined(FRAC_DUST)
504 c_ice = c(n_temp_1) &
505  + (c(n_temp_2)-c(n_temp_1)) &
506  * (temp_val-real(n_temp_1,dp))
507 #else
508 c_val = c(n_temp_1) &
509  + (c(n_temp_2)-c(n_temp_1)) &
510  * (temp_val-real(n_temp_1,dp))
511 #endif
512 
513 #else /* OpenAD */
514 
515 n_temp_1 = myfloor_local(temp_val)
516 n_temp_1 = max(min(n_temp_1, n_temp_max_imq-1), n_temp_min_imq)
517 n_temp_2 = n_temp_1 + 1
518 
519 #if defined(FRAC_DUST)
520 c_ice = c_imq(n_temp_1) &
521  + (c_imq(n_temp_2)-c_imq(n_temp_1)) &
522  * (temp_val-real(n_temp_1,dp))
523 #else
524 c_val = c_imq(n_temp_1) &
525  + (c_imq(n_temp_2)-c_imq(n_temp_1)) &
526  * (temp_val-real(n_temp_1,dp))
527 #endif
528 
529 #endif /* Normal vs. OpenAD */
530 
531 !-------- If dust is present (polar caps of Mars):
532 ! Specific heat of ice-dust mixture --------
533 
534 #if defined(FRAC_DUST)
535 #ifndef ALLOW_OPENAD /* Normal */
536 c_val = rho_inv * ( (1.0_dp-frac_dust)*rho_i*c_ice + frac_dust*rho_c*c_c )
537 #else /* OpenAD */
538 c_val = rho_inv * ( (1.0_dp-frac_dust)*rho_i_imq*c_ice +
539 frac_dust*rho_c_imq*c_c_imq )
540 #endif /* Normal vs. OpenAD */
541 #endif
542 
543 end function c_val
544 
545 !-------------------------------------------------------------------------------
546 !> Creep response function for ice.
547 !<------------------------------------------------------------------------------
548 function creep(sigma_val)
551 use sico_vars_m
552 
553 implicit none
554 
555 real(dp) :: creep
556 real(dp), intent(in) :: sigma_val
557 
558 #if (FLOW_LAW==4)
559 real(dp), parameter :: sm_coeff_1 = 8.5112e-15_dp, & ! s^-1 Pa^-1
560  sm_coeff_2 = 8.1643e-25_dp, & ! s^-1 Pa^-3
561  sm_coeff_3 = 9.2594e-12_dp ! Pa^-2
562 #endif
563 
564 #if (FIN_VISC==1)
565 
566 #if (FLOW_LAW==1)
567 
568 creep = sigma_val*sigma_val
569 ! Glen's flow law (n=3)
570 
571 #elif (FLOW_LAW==2)
572 
573 creep = sigma_val**0.8_dp * gr_size**(-1.4_dp)
574 ! Goldsby-Kohlstedt flow law (n=1.8, d=1.4)
575 
576 #elif (FLOW_LAW==3)
577 
578 creep = sigma_val*sigma_val*sigma_val
579 ! Durham's flow law (n=4)
580 
581 #endif
582 
583 #elif (FIN_VISC==2)
584 
585 #if (FLOW_LAW==1)
586 
587 creep = sigma_val*sigma_val + sigma_res*sigma_res
588 ! Glen's flow (n=3) with additional finite viscosity
589 
590 #elif (FLOW_LAW==2)
591 
592 creep = (sigma_val**0.8_dp + sigma_res**0.8_dp) * gr_size**(-1.4_dp)
593 ! Goldsby-Kohlstedt flow law (n=1.8, d=1.4)
594 ! with additional finite viscosity
595 
596 #elif (FLOW_LAW==3)
597 
598 creep = sigma_val*sigma_val*sigma_val + sigma_res*sigma_res*sigma_res
599 ! Durham's flow law (n=4) with additional finite viscosity
600 
601 #endif
602 
603 #endif
604 
605 #if (FLOW_LAW==4)
606 
607 creep = sm_coeff_1 &
608  + sm_coeff_2 * (sigma_val*sigma_val) &
609  * (1.0_dp + sm_coeff_3 * (sigma_val*sigma_val))
610 ! Smith-Morland (polynomial) flow law, normalised to a dimensionless
611 ! rate factor with A(-10C)=1.
612 
613 #endif
614 
615 end function creep
616 
617 !-------------------------------------------------------------------------------
618 !> Ice viscosity as a function of the effective strain rate and the temperature
619 !! (in cold ice) or the water content (in temperate ice) or both (for the
620 !! enthalpy method).
621 !<------------------------------------------------------------------------------
622 function viscosity(de_val, temp_val, temp_m_val, omega_val, enh_val, &
623  i_flag_cold_temp)
626 use sico_vars_m
627 
628 implicit none
629 
630 real(dp) :: viscosity
631 real(dp) , intent(in) :: de_val
632 real(dp) , intent(in) :: temp_val, temp_m_val
633 real(dp) , intent(in) :: omega_val
634 real(dp) , intent(in) :: enh_val
635 integer(i1b), intent(in) :: i_flag_cold_temp
636 
637 real(dp) :: ratefac_val
638 real(dp) :: de_val_m
639 
640 real(dp) :: n_power_law, inv_n_power_law, n_grain_size
641 
642 real(dp), parameter :: de_min = 1.0e-30_dp ! minimum value for the
643  ! effective strain rate
644 
645 #if (FLOW_LAW==4)
646 real(dp), parameter :: sm_coeff_1 = 8.5112e-15_dp, & ! s^-1 Pa^-1
647  sm_coeff_2 = 8.1643e-25_dp, & ! s^-1 Pa^-3
648  sm_coeff_3 = 9.2594e-12_dp ! Pa^-2
649 #endif
650 
651 !-------- Rate factor and effective strain rate --------
652 
653 if (i_flag_cold_temp == 0_i1b) then ! cold ice
654  ratefac_val = ratefac_c(temp_val, temp_m_val)
655 else if (i_flag_cold_temp == 1_i1b) then ! temperate ice
656  ratefac_val = ratefac_t(omega_val)
657 else ! enthalpy method
658  ratefac_val = ratefac_c_t(temp_val, omega_val, temp_m_val)
659 end if
660 
661 de_val_m = max(de_val, de_min)
662 
663 #if (FIN_VISC==1)
664 
665 #if (FLOW_LAW==1)
666 
667 !-------- Glen's flow law (n=3) --------
668 
669 inv_n_power_law = 0.333333333333333_dp ! 1/3
670 
671 viscosity = 0.5_dp * de_val_m**(inv_n_power_law-1.0_dp) &
672  * (enh_val*ratefac_val)**(-inv_n_power_law)
673 
674 #elif (FLOW_LAW==2)
675 
676 !-------- Goldsby-Kohlstedt flow law (n=1.8, d=1.4) --------
677 
678 inv_n_power_law = 0.555555555555555_dp ! 1/1.8
679 n_grain_size = 1.4_dp
680 
681 viscosity = 0.5_dp * de_val_m**(inv_n_power_law-1.0_dp) &
682  * gr_size**(n_grain_size*inv_n_power_law) &
683  * (enh_val*ratefac_val)**(-inv_n_power_law)
684 
685 #elif (FLOW_LAW==3)
686 
687 !-------- Durham's flow law (n=4) --------
688 
689 inv_n_power_law = 0.25_dp ! 1/4
690 
691 viscosity = 0.5_dp * de_val_m**(inv_n_power_law-1.0_dp) &
692  * (enh_val*ratefac_val)**(-inv_n_power_law)
693 
694 #endif
695 
696 #elif (FIN_VISC==2)
697 
698 #if (FLOW_LAW==1)
699 
700 !-------- Glen's flow (n=3) with additional finite viscosity --------
701 
702 n_power_law = 3.0_dp
703 
704 viscosity = visc_iter(de_val_m, ratefac_val, enh_val, n_power_law, sigma_res)
705 
706 #elif (FLOW_LAW==2)
707 
708 !-------- Goldsby-Kohlstedt flow law (n=1.8, d=1.4)
709 ! with additional finite viscosity --------
710 
711 n_power_law = 1.8_dp
712 n_grain_size = 1.4_dp
713 
714 errormsg = ' >>> viscosity: Computation of the viscosity as a function' &
715  // end_of_line &
716  //' of the effective strain rate not yet implemented' &
717  // end_of_line &
718  //' for grain-size-dependent finite viscosity flow laws!'
719 call error(errormsg)
720 
721 #elif (FLOW_LAW==3)
722 
723 !-------- Durham's flow law (n=4) with additional finite viscosity --------
724 
725 n_power_law = 4.0_dp
726 
727 viscosity = visc_iter(de_val_m, ratefac_val, enh_val, n_power_law, sigma_res)
728 
729 #endif
730 
731 #endif
732 
733 #if (FLOW_LAW==4)
734 
735 !-------- Smith-Morland (polynomial) flow law --------
736 
737 viscosity = visc_iter_sm(de_val_m, ratefac_val, enh_val, &
738  sm_coeff_1, sm_coeff_2, sm_coeff_3)
739 
740 #endif
741 
742 end function viscosity
743 
744 #if (FIN_VISC==2)
745 
746 !-------------------------------------------------------------------------------
747 !> Iterative computation of the viscosity by solving equation (4.28)
748 !! by Greve and Blatter (Springer, 2009).
749 !<------------------------------------------------------------------------------
750 function visc_iter(de_val_m, ratefac_val, enh_val, n_power_law, sigma_res)
751 
752 implicit none
753 
754 real(dp) :: visc_iter
755 real(dp), intent(in) :: de_val_m
756 real(dp), intent(in) :: ratefac_val, enh_val
757 real(dp), intent(in) :: n_power_law, sigma_res
758 
759 integer(i4b) :: n
760 integer(i4b) :: max_iters
761 real(dp) :: visc_val, res
762 logical :: flag_rescheck1, flag_rescheck2
763 
764 real(dp), parameter :: eps = 1.0e-05_dp ! convergence parameter
765 
766 !-------- Determination of the order of magnitude --------
767 
768 visc_val = 1.0e+10_dp ! initial guess (very low value)
769 
770 flag_rescheck1 = .false.
771 n = 0
772 max_iters = 30
773 
774 do while ((.not.flag_rescheck1).and.(n <= max_iters))
775 
776  n = n+1
777 
778  res = fct_visc(de_val_m, ratefac_val, enh_val, visc_val, &
779  n_power_law, sigma_res)
780 
781  if (res < 0.0_dp) then
782  visc_val = 10.0_dp*visc_val
783  else
784  flag_rescheck1 = .true.
785  end if
786 
787 end do
788 
789 !-------- Newton's method --------
790 
791 if (flag_rescheck1) then
792  ! only if order of magnitude could be detected successfully
793 
794  flag_rescheck2 = .false.
795  n = 0
796  max_iters = 1000
797 
798  do while ((.not.flag_rescheck2).and.(n <= max_iters))
799 
800  n = n+1
801 
802  visc_val = visc_val &
803  - res &
804  /fct_visc_deriv(de_val_m, ratefac_val, enh_val, visc_val, &
805  n_power_law, sigma_res)
806 
807  res = fct_visc(de_val_m, ratefac_val, enh_val, visc_val, &
808  n_power_law, sigma_res)
809 
810  if (abs(res) < eps) then
811  flag_rescheck2 = .true.
812  end if
813 
814  end do
815 
816 end if
817 
818 visc_iter = visc_val
819 
820 end function visc_iter
821 
822 !-------------------------------------------------------------------------------
823 !> Viscosity polynomial
824 !! [equation (4.28) by Greve and Blatter (Springer, 2009)].
825 !<------------------------------------------------------------------------------
826 function fct_visc(de_val_m, ratefac_val, enh_val, visc_val, &
827  n_power_law, sigma_res)
828 
829 implicit none
830 
831 real(dp) :: fct_visc
832 real(dp), intent(in) :: de_val_m
833 real(dp), intent(in) :: ratefac_val, enh_val
834 real(dp), intent(in) :: visc_val
835 real(dp), intent(in) :: n_power_law, sigma_res
836 
837 fct_visc = 2.0_dp**n_power_law &
838  *enh_val*ratefac_val &
839  *de_val_m**(n_power_law-1.0_dp) &
840  *visc_val**n_power_law &
841  + 2.0_dp*enh_val*ratefac_val &
842  *sigma_res**(n_power_law-1.0_dp) &
843  *visc_val &
844  - 1.0_dp
845 
846 end function fct_visc
847 
848 !-------------------------------------------------------------------------------
849 !> Derivative of the viscosity polynomial
850 !! [equation (4.28) by Greve and Blatter (Springer, 2009)].
851 !<------------------------------------------------------------------------------
852 function fct_visc_deriv(de_val_m, ratefac_val, enh_val, visc_val, &
853  n_power_law, sigma_res)
854 
855 implicit none
856 
857 real(dp) :: fct_visc_deriv
858 real(dp), intent(in) :: de_val_m
859 real(dp), intent(in) :: ratefac_val, enh_val
860 real(dp), intent(in) :: visc_val
861 real(dp), intent(in) :: n_power_law, sigma_res
862 
863 fct_visc_deriv = 2.0_dp**n_power_law*n_power_law &
864  *enh_val*ratefac_val &
865  *de_val_m**(n_power_law-1.0_dp) &
866  *visc_val**(n_power_law-1.0_dp) &
867  + 2.0_dp*enh_val*ratefac_val &
868  *sigma_res**(n_power_law-1.0_dp)
869 
870 end function fct_visc_deriv
871 
872 #endif /* FIN_VISC==2 */
873 
874 #if (FLOW_LAW==4)
875 
876 !-------------------------------------------------------------------------------
877 !> Iterative computation of the viscosity by solving equation (4.33)
878 !! [analogous to (4.28)] by Greve and Blatter (Springer, 2009).
879 !<------------------------------------------------------------------------------
880 function visc_iter_sm(de_val_m, ratefac_val, enh_val, &
881  sm_coeff_1, sm_coeff_2, sm_coeff_3)
882 
883 implicit none
884 
885 real(dp) :: visc_iter_sm
886 real(dp), intent(in) :: de_val_m
887 real(dp), intent(in) :: ratefac_val, enh_val
888 real(dp), intent(in) :: sm_coeff_1, sm_coeff_2, sm_coeff_3
889 
890 integer(i4b) :: n
891 integer(i4b) :: max_iters
892 real(dp) :: visc_val, res
893 logical :: flag_rescheck1, flag_rescheck2
894 
895 real(dp), parameter :: eps = 1.0e-05_dp ! convergence parameter
896 
897 !-------- Determination of the order of magnitude --------
898 
899 visc_val = 1.0e+10_dp ! initial guess (very low value)
900 
901 flag_rescheck1 = .false.
902 n = 0
903 max_iters = 30
904 
905 do while ((.not.flag_rescheck1).and.(n <= max_iters))
906 
907  n = n+1
908 
909  res = fct_visc_sm(de_val_m, ratefac_val, enh_val, visc_val, &
910  sm_coeff_1, sm_coeff_2, sm_coeff_3)
911 
912  if (res < 0.0_dp) then
913  visc_val = 10.0_dp*visc_val
914  else
915  flag_rescheck1 = .true.
916  end if
917 
918 end do
919 
920 !-------- Newton's method --------
921 
922 if (flag_rescheck1) then
923  ! only if order of magnitude could be detected successfully
924 
925  flag_rescheck2 = .false.
926  n = 0
927  max_iters = 1000
928 
929  do while ((.not.flag_rescheck2).and.(n <= max_iters))
930 
931  n = n+1
932 
933  visc_val = visc_val &
934  - res &
935  /fct_visc_sm_deriv(de_val_m, ratefac_val, &
936  enh_val, visc_val, &
937  sm_coeff_1, sm_coeff_2, sm_coeff_3)
938 
939  res = fct_visc_sm(de_val_m, ratefac_val, enh_val, visc_val, &
940  sm_coeff_1, sm_coeff_2, sm_coeff_3)
941 
942  if (abs(res) < eps) then
943  flag_rescheck2 = .true.
944  end if
945 
946  end do
947 
948 end if
949 
950 visc_iter_sm = visc_val
951 
952 end function visc_iter_sm
953 
954 !-------------------------------------------------------------------------------
955 !> Viscosity polynomial
956 !! [equation (4.33) by Greve and Blatter (Springer, 2009)].
957 !<------------------------------------------------------------------------------
958 function fct_visc_sm(de_val_m, ratefac_val, enh_val, visc_val, &
959  sm_coeff_1, sm_coeff_2, sm_coeff_3)
960 
961 implicit none
962 
963 real(dp) :: fct_visc_sm
964 real(dp), intent(in) :: de_val_m
965 real(dp), intent(in) :: ratefac_val, enh_val
966 real(dp), intent(in) :: visc_val
967 real(dp), intent(in) :: sm_coeff_1, sm_coeff_2, sm_coeff_3
968 
969 real(dp) :: de_visc_factor
970 
971 de_visc_factor = de_val_m*de_val_m*visc_val*visc_val
972 
973 fct_visc_sm = 2.0_dp*enh_val*ratefac_val*visc_val &
974  * ( sm_coeff_1 &
975  + 4.0_dp*sm_coeff_2*de_visc_factor &
976  * ( 1.0_dp + 4.0_dp*sm_coeff_3*de_visc_factor ) ) &
977  - 1.0_dp
978 
979 end function fct_visc_sm
980 
981 !-------------------------------------------------------------------------------
982 !> Derivative of the viscosity polynomial
983 !! [equation (4.33) by Greve and Blatter (Springer, 2009)].
984 !<------------------------------------------------------------------------------
985 function fct_visc_sm_deriv(de_val_m, ratefac_val, enh_val, visc_val, &
986  sm_coeff_1, sm_coeff_2, sm_coeff_3)
987 
988 implicit none
989 
990 real(dp) :: fct_visc_sm_deriv
991 real(dp), intent(in) :: de_val_m
992 real(dp), intent(in) :: ratefac_val, enh_val
993 real(dp), intent(in) :: visc_val
994 real(dp), intent(in) :: sm_coeff_1, sm_coeff_2, sm_coeff_3
995 
996 real(dp) :: de_visc_factor
997 
998 real(dp), parameter :: twenty_over_three = 6.666666666666667_dp
999 
1000 de_visc_factor = de_val_m*de_val_m*visc_val*visc_val
1001 
1002 fct_visc_sm_deriv = 2.0_dp*sm_coeff_1*enh_val*ratefac_val &
1003  + 24.0_dp*sm_coeff_2*enh_val*ratefac_val*de_visc_factor &
1004  * ( 1.0_dp + twenty_over_three*sm_coeff_3*de_visc_factor )
1005 
1006 end function fct_visc_sm_deriv
1007 
1008 #endif /* FLOW_LAW==4 */
1009 
1010 !-------------------------------------------------------------------------------
1011 
1012 #ifdef ALLOW_OPENAD /* OpenAD */
1013 
1014  function myfloor_local(num)
1015 
1016  use sico_variables_m
1017 
1018  implicit none
1019 
1020  real(dp), intent(in) :: num
1021 
1022  integer(i4b) :: inum
1023  integer(i4b) :: myfloor_local
1024 
1025  inum = int(num);
1026 
1027  if (abs(num-real(inum,dp)) <= &
1028  (abs(num+real(inum,dp))*myepsilon_dp) ) then
1029  myfloor_local = inum;
1030  else if (num>0) then
1031  myfloor_local = inum;
1032  else if (num<0) then
1033  myfloor_local = inum - 1;
1034  end if
1035 
1036  end function myfloor_local
1037 
1038 #endif /* OpenAD */
1039 
1040 !-------------------------------------------------------------------------------
1041 
1042 end module ice_material_properties_m
1043 !
subroutine error(error_message)
Main routine of error_m: Writing of error messages and stopping execution.
Definition: error_m.F90:47
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.
Writing of error messages and stopping execution.
Definition: error_m.F90:35
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
character, parameter end_of_line
end_of_line: End-of-line string
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...
character(len=256) errormsg
errormsg: Error-message string
integer, parameter dp
Double-precision reals.
Declarations of global variables for SICOPOLIS.