SICOPOLIS V5-dev  Revision 1279
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 real(dp) :: visc_val, res
769 
770 real(dp), parameter :: eps = 1.0e-05_dp ! convergence parameter
771 
772 #ifdef ALLOW_OPENAD /* OpenAD */
773 logical :: rescheck1, rescheck2
774 #endif /* OpenAD */
775 
776 !-------- Determination of the order of magnitude --------
777 
778 visc_val = 1.0e+10_dp ! initial guess (very low value)
779 
780 #ifndef ALLOW_OPENAD /* Normal */
781 
782 do
783 
784  res = fct_visc(de_val_m, ratefac_val, enh_val, visc_val, &
785  n_power_law, sigma_res)
786 
787  if (res < 0.0_dp) then
788  visc_val = 10.0_dp*visc_val
789  else
790  exit
791  end if
792 
793 end do
794 
795 #else /* OpenAD */
796 
797 rescheck1 = .false.
798 do while (rescheck1.eqv..false.)
799 
800  res = fct_visc(de_val_m, ratefac_val, enh_val, visc_val, &
801  n_power_law, sigma_res)
802 
803  if (res < 0.0_dp) then
804  visc_val = 10.0_dp*visc_val
805  else
806  rescheck1 = .true.
807  end if
808 
809 end do
810 
811 #endif /* Normal vs. OpenAD */
812 
813 !-------- Newton's method --------
814 
815 #ifndef ALLOW_OPENAD /* Normal */
816 
817 do
818 
819  visc_val = visc_val &
820  - res &
821  /fct_visc_deriv(de_val_m, ratefac_val, enh_val, visc_val, &
822  n_power_law, sigma_res)
823 
824  res = fct_visc(de_val_m, ratefac_val, enh_val, visc_val, &
825  n_power_law, sigma_res)
826 
827  if (abs(res) < eps) exit
828 
829 end do
830 
831 #else /* OpenAD */
832 
833 rescheck2 = .false.
834 do while (rescheck2.eqv..false.)
835 
836  visc_val = visc_val &
837  - res &
838  /fct_visc_deriv(de_val_m, ratefac_val, enh_val, visc_val, &
839  n_power_law, sigma_res)
840 
841  res = fct_visc(de_val_m, ratefac_val, enh_val, visc_val, &
842  n_power_law, sigma_res)
843 
844  if (abs(res) < eps) then
845  rescheck2 = .true.
846  end if
847 
848 end do
849 
850 #endif /* Normal vs. OpenAD */
851 
852 visc_iter = visc_val
853 
854 end function visc_iter
855 
856 !-------------------------------------------------------------------------------
857 !> Viscosity polynomial
858 !! [equation (4.28) by Greve and Blatter (Springer, 2009)].
859 !<------------------------------------------------------------------------------
860 function fct_visc(de_val_m, ratefac_val, enh_val, visc_val, &
861  n_power_law, sigma_res)
862 
863 use sico_types_m
864 
865 implicit none
866 
867 real(dp) :: fct_visc
868 real(dp), intent(in) :: de_val_m
869 real(dp), intent(in) :: ratefac_val, enh_val
870 real(dp), intent(in) :: visc_val
871 real(dp), intent(in) :: n_power_law, sigma_res
872 
873 fct_visc = 2.0_dp**n_power_law &
874  *enh_val*ratefac_val &
875  *de_val_m**(n_power_law-1.0_dp) &
876  *visc_val**n_power_law &
877  + 2.0_dp*enh_val*ratefac_val &
878  *sigma_res**(n_power_law-1.0_dp) &
879  *visc_val &
880  - 1.0_dp
881 
882 end function fct_visc
883 
884 !-------------------------------------------------------------------------------
885 !> Derivative of the viscosity polynomial
886 !! [equation (4.28) by Greve and Blatter (Springer, 2009)].
887 !<------------------------------------------------------------------------------
888 function fct_visc_deriv(de_val_m, ratefac_val, enh_val, visc_val, &
889  n_power_law, sigma_res)
890 
891 use sico_types_m
892 
893 implicit none
894 
895 real(dp) :: fct_visc_deriv
896 real(dp), intent(in) :: de_val_m
897 real(dp), intent(in) :: ratefac_val, enh_val
898 real(dp), intent(in) :: visc_val
899 real(dp), intent(in) :: n_power_law, sigma_res
900 
901 fct_visc_deriv = 2.0_dp**n_power_law*n_power_law &
902  *enh_val*ratefac_val &
903  *de_val_m**(n_power_law-1.0_dp) &
904  *visc_val**(n_power_law-1.0_dp) &
905  + 2.0_dp*enh_val*ratefac_val &
906  *sigma_res**(n_power_law-1.0_dp)
907 
908 end function fct_visc_deriv
909 
910 #endif
911 
912 #if (FLOW_LAW==4)
913 
914 !-------------------------------------------------------------------------------
915 !> Iterative computation of the viscosity by solving equation (4.33)
916 !! [analogous to (4.28)] by Greve and Blatter (Springer, 2009).
917 !<------------------------------------------------------------------------------
918 function visc_iter_sm(de_val_m, ratefac_val, enh_val, &
919  sm_coeff_1, sm_coeff_2, sm_coeff_3)
920 
921 use sico_types_m
922 
923 implicit none
924 
925 real(dp) :: visc_iter_sm
926 real(dp), intent(in) :: de_val_m
927 real(dp), intent(in) :: ratefac_val, enh_val
928 real(dp), intent(in) :: sm_coeff_1, sm_coeff_2, sm_coeff_3
929 
930 real(dp) :: visc_val, res
931 
932 real(dp), parameter :: eps = 1.0e-05_dp ! convergence parameter
933 #ifdef ALLOW_OPENAD /* OpenAD */
934 logical:: rescheck1, rescheck2
935 #endif /* OpenAD */
936 
937 !-------- Determination of the order of magnitude --------
938 
939 visc_val = 1.0e+10_dp ! initial guess (very low value)
940 
941 #ifndef ALLOW_OPENAD /* Normal */
942 
943 do
944 
945  res = fct_visc_sm(de_val_m, ratefac_val, enh_val, visc_val, &
946  sm_coeff_1, sm_coeff_2, sm_coeff_3)
947 
948  if (res < 0.0_dp) then
949  visc_val = 10.0_dp*visc_val
950  else
951  exit
952  end if
953 
954 end do
955 
956 #else /* OpenAD */
957 
958 rescheck1 = .false.
959 do while (rescheck1.eqv..false.)
960 
961  res = fct_visc_sm(de_val_m, ratefac_val, enh_val, visc_val, &
962  sm_coeff_1, sm_coeff_2, sm_coeff_3)
963 
964  if (res < 0.0_dp) then
965  visc_val = 10.0_dp*visc_val
966  else
967  rescheck1 = .true.
968  end if
969 
970 end do
971 
972 #endif /* Normal vs. OpenAD */
973 
974 !-------- Newton's method --------
975 
976 #ifndef ALLOW_OPENAD /* Normal */
977 
978 do
979 
980  visc_val = visc_val &
981  - res &
982  /fct_visc_sm_deriv(de_val_m, ratefac_val, enh_val, visc_val, &
983  sm_coeff_1, sm_coeff_2, sm_coeff_3)
984 
985  res = fct_visc_sm(de_val_m, ratefac_val, enh_val, visc_val, &
986  sm_coeff_1, sm_coeff_2, sm_coeff_3)
987 
988  if (abs(res) < eps) exit
989 
990 end do
991 
992 #else /* OpenAD */
993 
994 rescheck2 = .false.
995 
996 do while (rescheck2 .eqv. .false.) !--LL
997 
998  visc_val = visc_val &
999  - res &
1000  /fct_visc_sm_deriv(de_val_m, ratefac_val, enh_val, visc_val, &
1001  sm_coeff_1, sm_coeff_2, sm_coeff_3)
1002 
1003  res = fct_visc_sm(de_val_m, ratefac_val, enh_val, visc_val, &
1004  sm_coeff_1, sm_coeff_2, sm_coeff_3)
1005 
1006  if (abs(res) < eps)
1007  rescheck2 = .true.
1008 
1009 end do
1010 
1011 #endif /* Normal vs. OpenAD */
1012 
1013 visc_iter_sm = visc_val
1014 
1015 end function visc_iter_sm
1016 
1017 !-------------------------------------------------------------------------------
1018 !> Viscosity polynomial
1019 !! [equation (4.33) by Greve and Blatter (Springer, 2009)].
1020 !<------------------------------------------------------------------------------
1021 function fct_visc_sm(de_val_m, ratefac_val, enh_val, visc_val, &
1022  sm_coeff_1, sm_coeff_2, sm_coeff_3)
1023 
1024 use sico_types_m
1025 
1026 implicit none
1027 
1028 real(dp) :: fct_visc_sm
1029 real(dp), intent(in) :: de_val_m
1030 real(dp), intent(in) :: ratefac_val, enh_val
1031 real(dp), intent(in) :: visc_val
1032 real(dp), intent(in) :: sm_coeff_1, sm_coeff_2, sm_coeff_3
1033 
1034 real(dp) :: de_visc_factor
1035 
1036 de_visc_factor = de_val_m*de_val_m*visc_val*visc_val
1037 
1038 fct_visc_sm = 2.0_dp*enh_val*ratefac_val*visc_val &
1039  * ( sm_coeff_1 &
1040  + 4.0_dp*sm_coeff_2*de_visc_factor &
1041  * ( 1.0_dp + 4.0_dp*sm_coeff_3*de_visc_factor ) ) &
1042  - 1.0_dp
1043 
1044 end function fct_visc_sm
1045 
1046 !-------------------------------------------------------------------------------
1047 !> Derivative of the viscosity polynomial
1048 !! [equation (4.33) by Greve and Blatter (Springer, 2009)].
1049 !<------------------------------------------------------------------------------
1050 function fct_visc_sm_deriv(de_val_m, ratefac_val, enh_val, visc_val, &
1051  sm_coeff_1, sm_coeff_2, sm_coeff_3)
1052 
1053 use sico_types_m
1054 
1055 implicit none
1056 
1057 real(dp) :: fct_visc_sm_deriv
1058 real(dp), intent(in) :: de_val_m
1059 real(dp), intent(in) :: ratefac_val, enh_val
1060 real(dp), intent(in) :: visc_val
1061 real(dp), intent(in) :: sm_coeff_1, sm_coeff_2, sm_coeff_3
1062 
1063 real(dp) :: de_visc_factor
1064 
1065 real(dp), parameter :: twenty_over_three = 6.666666666666667_dp
1066 
1067 de_visc_factor = de_val_m*de_val_m*visc_val*visc_val
1068 
1069 fct_visc_sm_deriv = 2.0_dp*sm_coeff_1*enh_val*ratefac_val &
1070  + 24.0_dp*sm_coeff_2*enh_val*ratefac_val*de_visc_factor &
1071  * ( 1.0_dp + twenty_over_three*sm_coeff_3*de_visc_factor )
1072 
1073 end function fct_visc_sm_deriv
1074 
1075 #endif /* FLOW_LAW==4 */
1076 
1077 !-------------------------------------------------------------------------------
1078 
1079 #ifdef ALLOW_OPENAD /* OpenAD */
1080 
1081  function myfloor_local(num)
1082 
1083  use sico_types_m
1084  use sico_variables_m
1085 
1086  implicit none
1087 
1088  real(dp), intent(in) :: num
1089 
1090  integer(i4b) :: inum
1091  integer(i4b) :: myfloor_local
1092 
1093  inum = int(num);
1094 
1095  if (abs(num-real(inum,dp)) <= &
1096  (abs(num+real(inum,dp))*myepsilon_dp) ) then
1097  myfloor_local = inum;
1098  else if (num>0) then
1099  myfloor_local = inum;
1100  else if (num<0) then
1101  myfloor_local = inum - 1;
1102  end if
1103 
1104  end function myfloor_local
1105 
1106 #endif /* OpenAD */
1107 
1108 !-------------------------------------------------------------------------------
1109 
1110 end module ice_material_properties_m
1111 !
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.