SICOPOLIS V5-dev  Revision 1288
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)
302 use sico_types_m
304 use sico_vars_m
305 
306 implicit none
307 real(dp) :: ratefac_c
308 real(dp), intent(in) :: temp_val, temp_m_val
309 
310 integer(i4b) :: n_temp_1, n_temp_2
311 real(dp) :: temp_h_val
312 #ifdef ALLOW_OPENAD /* OpenAD */
313 integer(i4b) :: itemp_h_val
314 #endif /* OpenAD */
315 
316 temp_h_val = temp_val-temp_m_val
317 
318 #ifndef ALLOW_OPENAD /* Normal */
319 
320 n_temp_1 = floor(temp_h_val)
321 n_temp_1 = max(min(n_temp_1, n_temp_max-1), n_temp_min)
322 n_temp_2 = n_temp_1 + 1
323 
324 ratefac_c = rf(n_temp_1) &
325  + (rf(n_temp_2)-rf(n_temp_1)) &
326  * (temp_h_val-real(n_temp_1,dp)) ! Linear interpolation
327 
328 #else /* OpenAD */
329 
330 n_temp_1 = myfloor_local(temp_h_val)
331 n_temp_1 = max(min(n_temp_1, n_temp_max_imq-1), n_temp_min_imq)
332 n_temp_2 = n_temp_1 + 1
333 
334 ratefac_c = rf_imq(n_temp_1) &
335  +(rf_imq(n_temp_2)-rf_imq(n_temp_1)) &
336  *(temp_h_val-real(n_temp_1,dp))
337 
338 #endif /* Normal vs. OpenAD */
339 
340 end function ratefac_c
341 
342 !-------------------------------------------------------------------------------
343 !> Rate factor for temperate ice.
344 !<------------------------------------------------------------------------------
345 function ratefac_t(omega_val)
347 use sico_types_m
349 use sico_vars_m
350 
351 implicit none
352 real(dp) :: ratefac_t
353 real(dp), intent(in) :: omega_val
354 
355 #ifndef ALLOW_OPENAD /* Normal */
356 ratefac_t = rf(0)*(1.0_dp+r_t*(omega_val))
357 #else /* OpenAD */
358 ratefac_t = rf_imq(0)*(1.0_dp+r_t_imq*(omega_val))
359 #endif /* Normal vs. OpenAD */
360 
361 end function ratefac_t
362 
363 !-------------------------------------------------------------------------------
364 !> Rate factor for cold and temperate ice:
365 !! Combination of ratefac_c and ratefac_t (only for the enthalpy method).
366 !<------------------------------------------------------------------------------
367 function ratefac_c_t(temp_val, omega_val, temp_m_val)
369 use sico_types_m
371 use sico_vars_m
372 
373 implicit none
374 real(dp) :: ratefac_c_t
375 real(dp), intent(in) :: temp_val, temp_m_val, omega_val
376 
377 integer(i4b) :: n_temp_1, n_temp_2
378 real(dp) :: temp_h_val
379 #ifdef ALLOW_OPENAD /* OpenAD */
380 integer(i4b) :: itemp_h_val
381 #endif /* OpenAD */
382 
383 temp_h_val = temp_val-temp_m_val
384 
385 #ifndef ALLOW_OPENAD /* Normal */
386 
387 n_temp_1 = floor(temp_h_val)
388 n_temp_1 = max(min(n_temp_1, n_temp_max-1), n_temp_min)
389 n_temp_2 = n_temp_1 + 1
390 
391 ratefac_c_t = ( rf(n_temp_1) &
392  + (rf(n_temp_2)-rf(n_temp_1)) &
393  * (temp_h_val-real(n_temp_1,dp)) ) &
394  *(1.0_dp+r_t*(omega_val))
395 
396 #else /* OpenAD */
397 
398 n_temp_1 = myfloor_local(temp_h_val)
399 n_temp_1 = max(min(n_temp_1, n_temp_max_imq-1), n_temp_min_imq)
400 n_temp_2 = n_temp_1 + 1
401 
402 ratefac_c_t = ( rf_imq(n_temp_1) &
403  +(rf_imq(n_temp_2)-rf_imq(n_temp_1)) &
404  *(temp_h_val-real(n_temp_1,dp)) ) &
405  *(1.0_dp+r_t_imq*(omega_val))
406 
407 #endif /* Normal vs. OpenAD */
408 
409 end function ratefac_c_t
410 
411 !-------------------------------------------------------------------------------
412 !> Heat conductivity of ice:
413 !! Linear interpolation of tabulated values in KAPPA(.).
414 !<------------------------------------------------------------------------------
415 function kappa_val(temp_val)
417 use sico_types_m
419 use sico_vars_m
420 
421 implicit none
422 real(dp) :: kappa_val
423 real(dp), intent(in) :: temp_val
424 
425 integer(i4b) :: n_temp_1, n_temp_2
426 real(dp) :: kappa_ice
427 #ifdef ALLOW_OPENAD /* OpenAD */
428 integer(i4b) :: itemp_val
429 #endif /* OpenAD */
430 
431 !-------- Heat conductivity of pure ice --------
432 
433 #ifndef ALLOW_OPENAD /* Normal */
434 
435 n_temp_1 = floor(temp_val)
436 n_temp_1 = max(min(n_temp_1, n_temp_max-1), n_temp_min)
437 n_temp_2 = n_temp_1 + 1
438 
439 #if defined(FRAC_DUST)
440 kappa_ice = kappa(n_temp_1) &
441  + (kappa(n_temp_2)-kappa(n_temp_1)) &
442  * (temp_val-real(n_temp_1,dp))
443 #else
444 kappa_val = kappa(n_temp_1) &
445  + (kappa(n_temp_2)-kappa(n_temp_1)) &
446  * (temp_val-real(n_temp_1,dp))
447 #endif
448 
449 #else /* OpenAD */
450 
451 n_temp_1 = myfloor_local(temp_val)
452 n_temp_1 = max(min(n_temp_1, n_temp_max_imq-1), n_temp_min_imq)
453 n_temp_2 = n_temp_1 + 1
454 
455 #if defined(FRAC_DUST)
456 kappa_ice = 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 #else
460 kappa_val = kappa_imq(n_temp_1) &
461  + (kappa_imq(n_temp_2)-kappa_imq(n_temp_1)) &
462  * (temp_val-real(n_temp_1,dp))
463 #endif
464 
465 #endif /* Normal vs. OpenAD */
466 
467 !-------- If dust is present (polar caps of Mars):
468 ! Heat conductivity of ice-dust mixture --------
469 
470 #if defined(FRAC_DUST)
471 #ifndef ALLOW_OPENAD /* Normal */
472 kappa_val = (1.0_dp-frac_dust)*kappa_ice + frac_dust*kappa_c
473 #else /* OpenAD */
474 kappa_val = (1.0_dp-frac_dust)*kappa_ice + frac_dust*kappa_c_imq
475 #endif /* Normal vs. OpenAD */
476 #endif
477 
478 end function kappa_val
479 
480 !-------------------------------------------------------------------------------
481 !> Specific heat of ice:
482 !! Linear interpolation of tabulated values in C(.).
483 !<------------------------------------------------------------------------------
484 function c_val(temp_val)
486 use sico_types_m
488 use sico_vars_m
489 
490 implicit none
491 real(dp) :: c_val
492 real(dp), intent(in) :: temp_val
493 
494 integer(i4b) :: n_temp_1, n_temp_2
495 real(dp) :: c_ice
496 #ifdef ALLOW_OPENAD /* OpenAD */
497 integer(i4b) :: itemp_val
498 #endif /* OpenAD */
499 
500 !-------- Specific heat of pure ice --------
501 
502 #ifndef ALLOW_OPENAD /* Normal */
503 
504 n_temp_1 = floor(temp_val)
505 n_temp_1 = max(min(n_temp_1, n_temp_max-1), n_temp_min)
506 n_temp_2 = n_temp_1 + 1
507 
508 #if defined(FRAC_DUST)
509 c_ice = c(n_temp_1) &
510  + (c(n_temp_2)-c(n_temp_1)) &
511  * (temp_val-real(n_temp_1,dp))
512 #else
513 c_val = c(n_temp_1) &
514  + (c(n_temp_2)-c(n_temp_1)) &
515  * (temp_val-real(n_temp_1,dp))
516 #endif
517 
518 #else /* OpenAD */
519 
520 n_temp_1 = myfloor_local(temp_val)
521 n_temp_1 = max(min(n_temp_1, n_temp_max_imq-1), n_temp_min_imq)
522 n_temp_2 = n_temp_1 + 1
523 
524 #if defined(FRAC_DUST)
525 c_ice = c_imq(n_temp_1) &
526  + (c_imq(n_temp_2)-c_imq(n_temp_1)) &
527  * (temp_val-real(n_temp_1,dp))
528 #else
529 c_val = c_imq(n_temp_1) &
530  + (c_imq(n_temp_2)-c_imq(n_temp_1)) &
531  * (temp_val-real(n_temp_1,dp))
532 #endif
533 
534 #endif /* Normal vs. OpenAD */
535 
536 !-------- If dust is present (polar caps of Mars):
537 ! Specific heat of ice-dust mixture --------
538 
539 #if defined(FRAC_DUST)
540 #ifndef ALLOW_OPENAD /* Normal */
541 c_val = rho_inv * ( (1.0_dp-frac_dust)*rho_i*c_ice + frac_dust*rho_c*c_c )
542 #else /* OpenAD */
543 c_val = rho_inv * ( (1.0_dp-frac_dust)*rho_i_imq*c_ice +
544 frac_dust*rho_c_imq*c_c_imq )
545 #endif /* Normal vs. OpenAD */
546 #endif
547 
548 end function c_val
549 
550 !-------------------------------------------------------------------------------
551 !> Creep response function for ice.
552 !<------------------------------------------------------------------------------
553 function creep(sigma_val)
555 use sico_types_m
557 use sico_vars_m
558 
559 implicit none
560 
561 real(dp) :: creep
562 real(dp), intent(in) :: sigma_val
563 
564 #if (FLOW_LAW==4)
565 real(dp), parameter :: sm_coeff_1 = 8.5112e-15_dp, & ! s^-1 Pa^-1
566  sm_coeff_2 = 8.1643e-25_dp, & ! s^-1 Pa^-3
567  sm_coeff_3 = 9.2594e-12_dp ! Pa^-2
568 #endif
569 
570 #if (FIN_VISC==1)
571 
572 #if (FLOW_LAW==1)
573 
574 creep = sigma_val*sigma_val
575 ! Glen's flow law (n=3)
576 
577 #elif (FLOW_LAW==2)
578 
579 creep = sigma_val**0.8_dp * gr_size**(-1.4_dp)
580 ! Goldsby-Kohlstedt flow law (n=1.8, d=1.4)
581 
582 #elif (FLOW_LAW==3)
583 
584 creep = sigma_val*sigma_val*sigma_val
585 ! Durham's flow law (n=4)
586 
587 #endif
588 
589 #elif (FIN_VISC==2)
590 
591 #if (FLOW_LAW==1)
592 
593 creep = sigma_val*sigma_val + sigma_res*sigma_res
594 ! Glen's flow (n=3) with additional finite viscosity
595 
596 #elif (FLOW_LAW==2)
597 
598 creep = (sigma_val**0.8_dp + sigma_res**0.8_dp) * gr_size**(-1.4_dp)
599 ! Goldsby-Kohlstedt flow law (n=1.8, d=1.4)
600 ! with additional finite viscosity
601 
602 #elif (FLOW_LAW==3)
603 
604 creep = sigma_val*sigma_val*sigma_val + sigma_res*sigma_res*sigma_res
605 ! Durham's flow law (n=4) with additional finite viscosity
606 
607 #endif
608 
609 #endif
610 
611 #if (FLOW_LAW==4)
612 
613 creep = sm_coeff_1 &
614  + sm_coeff_2 * (sigma_val*sigma_val) &
615  * (1.0_dp + sm_coeff_3 * (sigma_val*sigma_val))
616 ! Smith-Morland (polynomial) flow law, normalised to a dimensionless
617 ! rate factor with A(-10C)=1.
618 
619 #endif
620 
621 end function creep
622 
623 !-------------------------------------------------------------------------------
624 !> Ice viscosity as a function of the effective strain rate and the temperature
625 !! (in cold ice) or the water content (in temperate ice) or both (for the
626 !! enthalpy method).
627 !<------------------------------------------------------------------------------
628 function viscosity(de_val, temp_val, temp_m_val, omega_val, enh_val, &
629  i_flag_cold_temp)
631 use sico_types_m
633 use sico_vars_m
634 
635 implicit none
636 
637 real(dp) :: viscosity
638 real(dp) , intent(in) :: de_val
639 real(dp) , intent(in) :: temp_val, temp_m_val
640 real(dp) , intent(in) :: omega_val
641 real(dp) , intent(in) :: enh_val
642 integer(i2b), intent(in) :: i_flag_cold_temp
643 
644 real(dp) :: ratefac_val
645 real(dp) :: de_val_m
646 
647 real(dp) :: n_power_law, inv_n_power_law, n_grain_size
648 
649 real(dp), parameter :: de_min = 1.0e-30_dp ! minimum value for the
650  ! effective strain rate
651 
652 #if (FLOW_LAW==4)
653 real(dp), parameter :: sm_coeff_1 = 8.5112e-15_dp, & ! s^-1 Pa^-1
654  sm_coeff_2 = 8.1643e-25_dp, & ! s^-1 Pa^-3
655  sm_coeff_3 = 9.2594e-12_dp ! Pa^-2
656 #endif
657 
658 !-------- Rate factor and effective strain rate --------
659 
660 if (i_flag_cold_temp == 0_i2b) then ! cold ice
661  ratefac_val = ratefac_c(temp_val, temp_m_val)
662 else if (i_flag_cold_temp == 1_i2b) then ! temperate ice
663  ratefac_val = ratefac_t(omega_val)
664 else ! enthalpy method
665  ratefac_val = ratefac_c_t(temp_val, omega_val, temp_m_val)
666 end if
667 
668 de_val_m = max(de_val, de_min)
669 
670 #if (FIN_VISC==1)
671 
672 #if (FLOW_LAW==1)
673 
674 !-------- Glen's flow law (n=3) --------
675 
676 inv_n_power_law = 0.333333333333333_dp ! 1/3
677 
678 viscosity = 0.5_dp * de_val_m**(inv_n_power_law-1.0_dp) &
679  * (enh_val*ratefac_val)**(-inv_n_power_law)
680 
681 #elif (FLOW_LAW==2)
682 
683 !-------- Goldsby-Kohlstedt flow law (n=1.8, d=1.4) --------
684 
685 inv_n_power_law = 0.555555555555555_dp ! 1/1.8
686 n_grain_size = 1.4_dp
687 
688 viscosity = 0.5_dp * de_val_m**(inv_n_power_law-1.0_dp) &
689  * gr_size**(n_grain_size*inv_n_power_law) &
690  * (enh_val*ratefac_val)**(-inv_n_power_law)
691 
692 #elif (FLOW_LAW==3)
693 
694 !-------- Durham's flow law (n=4) --------
695 
696 inv_n_power_law = 0.25_dp ! 1/4
697 
698 viscosity = 0.5_dp * de_val_m**(inv_n_power_law-1.0_dp) &
699  * (enh_val*ratefac_val)**(-inv_n_power_law)
700 
701 #endif
702 
703 #elif (FIN_VISC==2)
704 
705 #if (FLOW_LAW==1)
706 
707 !-------- Glen's flow (n=3) with additional finite viscosity --------
708 
709 n_power_law = 3.0_dp
710 
711 viscosity = visc_iter(de_val_m, ratefac_val, enh_val, n_power_law, sigma_res)
712 
713 #elif (FLOW_LAW==2)
714 
715 !-------- Goldsby-Kohlstedt flow law (n=1.8, d=1.4)
716 ! with additional finite viscosity --------
717 
718 n_power_law = 1.8_dp
719 n_grain_size = 1.4_dp
720 
721 errormsg = ' >>> viscosity: Computation of the viscosity as a function' &
722  // end_of_line &
723  //' of the effective strain rate not yet implemented' &
724  // end_of_line &
725  //' for grain-size-dependent finite viscosity flow laws!'
726 call error(errormsg)
727 
728 #elif (FLOW_LAW==3)
729 
730 !-------- Durham's flow law (n=4) with additional finite viscosity --------
731 
732 n_power_law = 4.0_dp
733 
734 viscosity = visc_iter(de_val_m, ratefac_val, enh_val, n_power_law, sigma_res)
735 
736 #endif
737 
738 #endif
739 
740 #if (FLOW_LAW==4)
741 
742 !-------- Smith-Morland (polynomial) flow law --------
743 
744 viscosity = visc_iter_sm(de_val_m, ratefac_val, enh_val, &
745  sm_coeff_1, sm_coeff_2, sm_coeff_3)
746 
747 #endif
748 
749 end function viscosity
750 
751 #if (FIN_VISC==2)
752 
753 !-------------------------------------------------------------------------------
754 !> Iterative computation of the viscosity by solving equation (4.28)
755 !! by Greve and Blatter (Springer, 2009).
756 !<------------------------------------------------------------------------------
757 function visc_iter(de_val_m, ratefac_val, enh_val, n_power_law, sigma_res)
758 
759 use sico_types_m
760 
761 implicit none
762 
763 real(dp) :: visc_iter
764 real(dp), intent(in) :: de_val_m
765 real(dp), intent(in) :: ratefac_val, enh_val
766 real(dp), intent(in) :: n_power_law, sigma_res
767 
768 integer(i4b) :: n
769 integer(i4b) :: max_iters
770 real(dp) :: visc_val, res
771 logical :: flag_rescheck1, flag_rescheck2
772 
773 real(dp), parameter :: eps = 1.0e-05_dp ! convergence parameter
774 
775 !-------- Determination of the order of magnitude --------
776 
777 visc_val = 1.0e+10_dp ! initial guess (very low value)
778 
779 flag_rescheck1 = .false.
780 n = 0
781 max_iters = 30
782 
783 do while ((.not.flag_rescheck1).and.(n <= max_iters))
784 
785  n = n+1
786 
787  res = fct_visc(de_val_m, ratefac_val, enh_val, visc_val, &
788  n_power_law, sigma_res)
789 
790  if (res < 0.0_dp) then
791  visc_val = 10.0_dp*visc_val
792  else
793  flag_rescheck1 = .true.
794  end if
795 
796 end do
797 
798 !-------- Newton's method --------
799 
800 if (flag_rescheck1) then
801  ! only if order of magnitude could be detected successfully
802 
803  flag_rescheck2 = .false.
804  n = 0
805  max_iters = 1000
806 
807  do while ((.not.flag_rescheck2).and.(n <= max_iters))
808 
809  n = n+1
810 
811  visc_val = visc_val &
812  - res &
813  /fct_visc_deriv(de_val_m, ratefac_val, enh_val, visc_val, &
814  n_power_law, sigma_res)
815 
816  res = fct_visc(de_val_m, ratefac_val, enh_val, visc_val, &
817  n_power_law, sigma_res)
818 
819  if (abs(res) < eps) then
820  flag_rescheck2 = .true.
821  end if
822 
823  end do
824 
825 end if
826 
827 visc_iter = visc_val
828 
829 end function visc_iter
830 
831 !-------------------------------------------------------------------------------
832 !> Viscosity polynomial
833 !! [equation (4.28) by Greve and Blatter (Springer, 2009)].
834 !<------------------------------------------------------------------------------
835 function fct_visc(de_val_m, ratefac_val, enh_val, visc_val, &
836  n_power_law, sigma_res)
837 
838 use sico_types_m
839 
840 implicit none
841 
842 real(dp) :: fct_visc
843 real(dp), intent(in) :: de_val_m
844 real(dp), intent(in) :: ratefac_val, enh_val
845 real(dp), intent(in) :: visc_val
846 real(dp), intent(in) :: n_power_law, sigma_res
847 
848 fct_visc = 2.0_dp**n_power_law &
849  *enh_val*ratefac_val &
850  *de_val_m**(n_power_law-1.0_dp) &
851  *visc_val**n_power_law &
852  + 2.0_dp*enh_val*ratefac_val &
853  *sigma_res**(n_power_law-1.0_dp) &
854  *visc_val &
855  - 1.0_dp
856 
857 end function fct_visc
858 
859 !-------------------------------------------------------------------------------
860 !> Derivative of the viscosity polynomial
861 !! [equation (4.28) by Greve and Blatter (Springer, 2009)].
862 !<------------------------------------------------------------------------------
863 function fct_visc_deriv(de_val_m, ratefac_val, enh_val, visc_val, &
864  n_power_law, sigma_res)
865 
866 use sico_types_m
867 
868 implicit none
869 
870 real(dp) :: fct_visc_deriv
871 real(dp), intent(in) :: de_val_m
872 real(dp), intent(in) :: ratefac_val, enh_val
873 real(dp), intent(in) :: visc_val
874 real(dp), intent(in) :: n_power_law, sigma_res
875 
876 fct_visc_deriv = 2.0_dp**n_power_law*n_power_law &
877  *enh_val*ratefac_val &
878  *de_val_m**(n_power_law-1.0_dp) &
879  *visc_val**(n_power_law-1.0_dp) &
880  + 2.0_dp*enh_val*ratefac_val &
881  *sigma_res**(n_power_law-1.0_dp)
882 
883 end function fct_visc_deriv
884 
885 #endif /* FIN_VISC==2 */
886 
887 #if (FLOW_LAW==4)
888 
889 !-------------------------------------------------------------------------------
890 !> Iterative computation of the viscosity by solving equation (4.33)
891 !! [analogous to (4.28)] by Greve and Blatter (Springer, 2009).
892 !<------------------------------------------------------------------------------
893 function visc_iter_sm(de_val_m, ratefac_val, enh_val, &
894  sm_coeff_1, sm_coeff_2, sm_coeff_3)
895 
896 use sico_types_m
897 
898 implicit none
899 
900 real(dp) :: visc_iter_sm
901 real(dp), intent(in) :: de_val_m
902 real(dp), intent(in) :: ratefac_val, enh_val
903 real(dp), intent(in) :: sm_coeff_1, sm_coeff_2, sm_coeff_3
904 
905 integer(i4b) :: n
906 integer(i4b) :: max_iters
907 real(dp) :: visc_val, res
908 logical :: flag_rescheck1, flag_rescheck2
909 
910 real(dp), parameter :: eps = 1.0e-05_dp ! convergence parameter
911 
912 !-------- Determination of the order of magnitude --------
913 
914 visc_val = 1.0e+10_dp ! initial guess (very low value)
915 
916 flag_rescheck1 = .false.
917 n = 0
918 max_iters = 30
919 
920 do while ((.not.flag_rescheck1).and.(n <= max_iters))
921 
922  n = n+1
923 
924  res = fct_visc_sm(de_val_m, ratefac_val, enh_val, visc_val, &
925  sm_coeff_1, sm_coeff_2, sm_coeff_3)
926 
927  if (res < 0.0_dp) then
928  visc_val = 10.0_dp*visc_val
929  else
930  flag_rescheck1 = .true.
931  end if
932 
933 end do
934 
935 !-------- Newton's method --------
936 
937 if (flag_rescheck1) then
938  ! only if order of magnitude could be detected successfully
939 
940  flag_rescheck2 = .false.
941  n = 0
942  max_iters = 1000
943 
944  do while ((.not.flag_rescheck2).and.(n <= max_iters))
945 
946  n = n+1
947 
948  visc_val = visc_val &
949  - res &
950  /fct_visc_sm_deriv(de_val_m, ratefac_val, &
951  enh_val, visc_val, &
952  sm_coeff_1, sm_coeff_2, sm_coeff_3)
953 
954  res = fct_visc_sm(de_val_m, ratefac_val, enh_val, visc_val, &
955  sm_coeff_1, sm_coeff_2, sm_coeff_3)
956 
957  if (abs(res) < eps) then
958  flag_rescheck2 = .true.
959  end if
960 
961  end do
962 
963 end if
964 
965 visc_iter_sm = visc_val
966 
967 end function visc_iter_sm
968 
969 !-------------------------------------------------------------------------------
970 !> Viscosity polynomial
971 !! [equation (4.33) by Greve and Blatter (Springer, 2009)].
972 !<------------------------------------------------------------------------------
973 function fct_visc_sm(de_val_m, ratefac_val, enh_val, visc_val, &
974  sm_coeff_1, sm_coeff_2, sm_coeff_3)
975 
976 use sico_types_m
977 
978 implicit none
979 
980 real(dp) :: fct_visc_sm
981 real(dp), intent(in) :: de_val_m
982 real(dp), intent(in) :: ratefac_val, enh_val
983 real(dp), intent(in) :: visc_val
984 real(dp), intent(in) :: sm_coeff_1, sm_coeff_2, sm_coeff_3
985 
986 real(dp) :: de_visc_factor
987 
988 de_visc_factor = de_val_m*de_val_m*visc_val*visc_val
989 
990 fct_visc_sm = 2.0_dp*enh_val*ratefac_val*visc_val &
991  * ( sm_coeff_1 &
992  + 4.0_dp*sm_coeff_2*de_visc_factor &
993  * ( 1.0_dp + 4.0_dp*sm_coeff_3*de_visc_factor ) ) &
994  - 1.0_dp
995 
996 end function fct_visc_sm
997 
998 !-------------------------------------------------------------------------------
999 !> Derivative of the viscosity polynomial
1000 !! [equation (4.33) by Greve and Blatter (Springer, 2009)].
1001 !<------------------------------------------------------------------------------
1002 function fct_visc_sm_deriv(de_val_m, ratefac_val, enh_val, visc_val, &
1003  sm_coeff_1, sm_coeff_2, sm_coeff_3)
1004 
1005 use sico_types_m
1006 
1007 implicit none
1008 
1009 real(dp) :: fct_visc_sm_deriv
1010 real(dp), intent(in) :: de_val_m
1011 real(dp), intent(in) :: ratefac_val, enh_val
1012 real(dp), intent(in) :: visc_val
1013 real(dp), intent(in) :: sm_coeff_1, sm_coeff_2, sm_coeff_3
1014 
1015 real(dp) :: de_visc_factor
1016 
1017 real(dp), parameter :: twenty_over_three = 6.666666666666667_dp
1018 
1019 de_visc_factor = de_val_m*de_val_m*visc_val*visc_val
1020 
1021 fct_visc_sm_deriv = 2.0_dp*sm_coeff_1*enh_val*ratefac_val &
1022  + 24.0_dp*sm_coeff_2*enh_val*ratefac_val*de_visc_factor &
1023  * ( 1.0_dp + twenty_over_three*sm_coeff_3*de_visc_factor )
1024 
1025 end function fct_visc_sm_deriv
1026 
1027 #endif /* FLOW_LAW==4 */
1028 
1029 !-------------------------------------------------------------------------------
1030 
1031 #ifdef ALLOW_OPENAD /* OpenAD */
1032 
1033  function myfloor_local(num)
1034 
1035  use sico_types_m
1036  use sico_variables_m
1037 
1038  implicit none
1039 
1040  real(dp), intent(in) :: num
1041 
1042  integer(i4b) :: inum
1043  integer(i4b) :: myfloor_local
1044 
1045  inum = int(num);
1046 
1047  if (abs(num-real(inum,dp)) <= &
1048  (abs(num+real(inum,dp))*myepsilon_dp) ) then
1049  myfloor_local = inum;
1050  else if (num>0) then
1051  myfloor_local = inum;
1052  else if (num<0) then
1053  myfloor_local = inum - 1;
1054  end if
1055 
1056  end function myfloor_local
1057 
1058 #endif /* OpenAD */
1059 
1060 !-------------------------------------------------------------------------------
1061 
1062 end module ice_material_properties_m
1063 !
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.