SICOPOLIS V3.0
 All Classes Files Functions Variables Macros
viscosity.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Function : v i s c o s i t y
4 !
5 !> @file
6 !!
7 !! Ice viscosity as a function of the effective strain rate and the temperature
8 !! (in cold ice) or the water content (in temperate ice).
9 !!
10 !! @section Copyright
11 !!
12 !! Copyright 2009-2013 Ralf Greve
13 !!
14 !! @section License
15 !!
16 !! This file is part of SICOPOLIS.
17 !!
18 !! SICOPOLIS is free software: you can redistribute it and/or modify
19 !! it under the terms of the GNU General Public License as published by
20 !! the Free Software Foundation, either version 3 of the License, or
21 !! (at your option) any later version.
22 !!
23 !! SICOPOLIS is distributed in the hope that it will be useful,
24 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
25 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
26 !! GNU General Public License for more details.
27 !!
28 !! You should have received a copy of the GNU General Public License
29 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
30 !<
31 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
32 
33 !-------------------------------------------------------------------------------
34 !> Ice viscosity as a function of the effective strain rate and the temperature
35 !! (in cold ice) or the water content (in temperate ice).
36 !<------------------------------------------------------------------------------
37 function viscosity(de_val, temp_val, temp_m_val, omega_val, enh_val, &
38  flag_cold_temp)
39 
40 use sico_types
41 use sico_variables
42 
43 implicit none
44 
45 real(dp) :: viscosity
46 real(dp), intent(in) :: de_val
47 real(dp), intent(in) :: temp_val, temp_m_val
48 real(dp), intent(in) :: omega_val
49 real(dp), intent(in) :: enh_val
50 logical, intent(in) :: flag_cold_temp
51 
52 real(dp) :: ratefac, ratefac_t
53 real(dp) :: visc_iter, visc_iter_sm
54 real(dp) :: ratefac_val
55 real(dp) :: de_val_m
56 
57 real(dp) :: n_power_law, inv_n_power_law, n_grain_size
58 
59 real(dp), parameter :: de_min = 1.0e-30_dp ! minimum value for the
60  ! effective strain rate
61 
62 #if (FLOW_LAW==4)
63 real(dp), parameter :: sm_coeff_1 = 8.5112e-15_dp, & ! s^-1 Pa^-1
64  sm_coeff_2 = 8.1643e-25_dp, & ! s^-1 Pa^-3
65  sm_coeff_3 = 9.2594e-12_dp ! Pa^-2
66 #endif
67 
68 !-------- Rate factor and effective strain rate --------
69 
70 if (flag_cold_temp) then ! cold ice
71  ratefac_val = ratefac(temp_val,temp_m_val)
72 else ! temperate ice
73  ratefac_val = ratefac_t(omega_val)
74 end if
75 
76 de_val_m = max(de_val, de_min)
77 
78 #if (FIN_VISC==1)
79 
80 #if (FLOW_LAW==1)
81 
82 !-------- Glen's flow law (n=3) --------
83 
84 inv_n_power_law = 0.333333333333333_dp ! 1/3
85 
86 viscosity = 0.5_dp * de_val_m**(inv_n_power_law-1.0_dp) &
87  * (enh_val*ratefac_val)**(-inv_n_power_law)
88 
89 #elif (FLOW_LAW==2)
90 
91 !-------- Goldsby-Kohlstedt flow law (n=1.8, d=1.4) --------
92 
93 inv_n_power_law = 0.555555555555555_dp ! 1/1.8
94 n_grain_size = 1.4_dp
95 
96 viscosity = 0.5_dp * de_val_m**(inv_n_power_law-1.0_dp) &
97  * gr_size**(n_grain_size*inv_n_power_law) &
98  * (enh_val*ratefac_val)**(-inv_n_power_law)
99 
100 #elif (FLOW_LAW==3)
101 
102 !-------- Durham's flow law (n=4) --------
103 
104 inv_n_power_law = 0.25_dp ! 1/4
105 
106 viscosity = 0.5_dp * de_val_m**(inv_n_power_law-1.0_dp) &
107  * (enh_val*ratefac_val)**(-inv_n_power_law)
108 
109 #endif
110 
111 #elif (FIN_VISC==2)
112 
113 #if (FLOW_LAW==1)
114 
115 !-------- Glen's flow (n=3) with additional finite viscosity --------
116 
117 n_power_law = 3.0_dp
118 
119 viscosity = visc_iter(de_val_m, ratefac_val, enh_val, n_power_law, sigma_res)
120 
121 #elif (FLOW_LAW==2)
122 
123 !-------- Goldsby-Kohlstedt flow law (n=1.8, d=1.4)
124 ! with additional finite viscosity --------
125 
126 n_power_law = 1.8_dp
127 n_grain_size = 1.4_dp
128 
129 write(6,'(a)') ' viscosity: Computation of the viscosity as a function'
130 write(6,'(a)') ' of the effective strain rate not yet implemented'
131 write(6,'(a)') ' for grain-size-dependent finite viscosity flow laws!'
132 stop
133 
134 #elif (FLOW_LAW==3)
135 
136 !-------- Durham's flow law (n=4) with additional finite viscosity --------
137 
138 n_power_law = 4.0_dp
139 
140 viscosity = visc_iter(de_val_m, ratefac_val, enh_val, n_power_law, sigma_res)
141 
142 #endif
143 
144 #endif
145 
146 #if (FLOW_LAW==4)
147 
148 !-------- Smith-Morland (polynomial) flow law --------
149 
150 viscosity = visc_iter_sm(de_val_m, ratefac_val, enh_val, &
151  sm_coeff_1, sm_coeff_2, sm_coeff_3)
152 
153 #endif
154 
155 end function viscosity
156 
157 #if (FIN_VISC==2)
158 
159 !-------------------------------------------------------------------------------
160 !> Iterative computation of the viscosity by solving equation (4.28)
161 !! by Greve and Blatter (Springer, 2009).
162 !<------------------------------------------------------------------------------
163 function visc_iter(de_val_m, ratefac_val, enh_val, n_power_law, sigma_res)
164 
165 use sico_types
166 
167 implicit none
168 
169 real(dp) :: visc_iter
170 real(dp), intent(in) :: de_val_m
171 real(dp), intent(in) :: ratefac_val, enh_val
172 real(dp), intent(in) :: n_power_law, sigma_res
173 
174 real(dp) :: visc_val, res
175 real(dp) :: fct_visc, fct_visc_deriv
176 
177 real(dp), parameter :: eps = 1.0e-05_dp ! convergence parameter
178 
179 !-------- Determination of the order of magnitude --------
180 
181 visc_val = 1.0e+10_dp ! initial guess (very low value)
182 
183 do
184 
185  res = fct_visc(de_val_m, ratefac_val, enh_val, visc_val, &
186  n_power_law, sigma_res)
187 
188  if (res < 0.0_dp) then
189  visc_val = 10.0_dp*visc_val
190  else
191  exit
192  end if
193 
194 end do
195 
196 !-------- Newton's method --------
197 
198 do
199 
200  visc_val = visc_val &
201  - res &
202  /fct_visc_deriv(de_val_m, ratefac_val, enh_val, visc_val, &
203  n_power_law, sigma_res)
204 
205  res = fct_visc(de_val_m, ratefac_val, enh_val, visc_val, &
206  n_power_law, sigma_res)
207 
208  if (abs(res) < eps) exit
209 
210 end do
211 
212 visc_iter = visc_val
213 
214 end function visc_iter
215 
216 !-------------------------------------------------------------------------------
217 !> Viscosity polynomial
218 !! [equation (4.28) by Greve and Blatter (Springer, 2009)].
219 !<------------------------------------------------------------------------------
220 function fct_visc(de_val_m, ratefac_val, enh_val, visc_val, &
221  n_power_law, sigma_res)
222 
223 use sico_types
224 
225 implicit none
226 
227 real(dp) :: fct_visc
228 real(dp), intent(in) :: de_val_m
229 real(dp), intent(in) :: ratefac_val, enh_val
230 real(dp), intent(in) :: visc_val
231 real(dp), intent(in) :: n_power_law, sigma_res
232 
233 fct_visc = 2.0_dp**n_power_law &
234  *enh_val*ratefac_val &
235  *de_val_m**(n_power_law-1.0_dp) &
236  *visc_val**n_power_law &
237  + 2.0_dp*enh_val*ratefac_val &
238  *sigma_res**(n_power_law-1.0_dp) &
239  *visc_val &
240  - 1.0_dp
241 
242 end function fct_visc
243 
244 !-------------------------------------------------------------------------------
245 !> Derivative of the viscosity polynomial
246 !! [equation (4.28) by Greve and Blatter (Springer, 2009)].
247 !<------------------------------------------------------------------------------
248 function fct_visc_deriv(de_val_m, ratefac_val, enh_val, visc_val, &
249  n_power_law, sigma_res)
250 
251 use sico_types
252 
253 implicit none
254 
255 real(dp) :: fct_visc_deriv
256 real(dp), intent(in) :: de_val_m
257 real(dp), intent(in) :: ratefac_val, enh_val
258 real(dp), intent(in) :: visc_val
259 real(dp), intent(in) :: n_power_law, sigma_res
260 
261 fct_visc_deriv = 2.0_dp**n_power_law*n_power_law &
262  *enh_val*ratefac_val &
263  *de_val_m**(n_power_law-1.0_dp) &
264  *visc_val**(n_power_law-1.0_dp) &
265  + 2.0_dp*enh_val*ratefac_val &
266  *sigma_res**(n_power_law-1.0_dp)
267 
268 end function fct_visc_deriv
269 
270 #endif
271 
272 #if (FLOW_LAW==4)
273 
274 !-------------------------------------------------------------------------------
275 !> Iterative computation of the viscosity by solving equation (4.33)
276 !! [analogous to (4.28)] by Greve and Blatter (Springer, 2009).
277 !<------------------------------------------------------------------------------
278 function visc_iter_sm(de_val_m, ratefac_val, enh_val, &
279  sm_coeff_1, sm_coeff_2, sm_coeff_3)
280 
281 use sico_types
282 
283 implicit none
284 
285 real(dp) :: visc_iter_sm
286 real(dp), intent(in) :: de_val_m
287 real(dp), intent(in) :: ratefac_val, enh_val
288 real(dp), intent(in) :: sm_coeff_1, sm_coeff_2, sm_coeff_3
289 
290 real(dp) :: visc_val, res
291 real(dp) :: fct_visc_sm, fct_visc_sm_deriv
292 
293 real(dp), parameter :: eps = 1.0e-05_dp ! convergence parameter
294 
295 !-------- Determination of the order of magnitude --------
296 
297 visc_val = 1.0e+10_dp ! initial guess (very low value)
298 
299 do
300 
301  res = fct_visc_sm(de_val_m, ratefac_val, enh_val, visc_val, &
302  sm_coeff_1, sm_coeff_2, sm_coeff_3)
303 
304  if (res < 0.0_dp) then
305  visc_val = 10.0_dp*visc_val
306  else
307  exit
308  end if
309 
310 end do
311 
312 !-------- Newton's method --------
313 
314 do
315 
316  visc_val = visc_val &
317  - res &
318  /fct_visc_sm_deriv(de_val_m, ratefac_val, enh_val, visc_val, &
319  sm_coeff_1, sm_coeff_2, sm_coeff_3)
320 
321  res = fct_visc_sm(de_val_m, ratefac_val, enh_val, visc_val, &
322  sm_coeff_1, sm_coeff_2, sm_coeff_3)
323 
324  if (abs(res) < eps) exit
325 
326 end do
327 
328 visc_iter_sm = visc_val
329 
330 end function visc_iter_sm
331 
332 !-------------------------------------------------------------------------------
333 !> Viscosity polynomial
334 !! [equation (4.33) by Greve and Blatter (Springer, 2009)].
335 !<------------------------------------------------------------------------------
336 function fct_visc_sm(de_val_m, ratefac_val, enh_val, visc_val, &
337  sm_coeff_1, sm_coeff_2, sm_coeff_3)
338 
339 use sico_types
340 
341 implicit none
342 
343 real(dp) :: fct_visc_sm
344 real(dp), intent(in) :: de_val_m
345 real(dp), intent(in) :: ratefac_val, enh_val
346 real(dp), intent(in) :: visc_val
347 real(dp), intent(in) :: sm_coeff_1, sm_coeff_2, sm_coeff_3
348 
349 real(dp) :: de_visc_factor
350 
351 de_visc_factor = de_val_m*de_val_m*visc_val*visc_val
352 
353 fct_visc_sm = 2.0_dp*enh_val*ratefac_val*visc_val &
354  * ( sm_coeff_1 &
355  + 4.0_dp*sm_coeff_2*de_visc_factor &
356  * ( 1.0_dp + 4.0_dp*sm_coeff_3*de_visc_factor ) ) &
357  - 1.0_dp
358 
359 end function fct_visc_sm
360 
361 !-------------------------------------------------------------------------------
362 !> Derivative of the viscosity polynomial
363 !! [equation (4.33) by Greve and Blatter (Springer, 2009)].
364 !<------------------------------------------------------------------------------
365 function fct_visc_sm_deriv(de_val_m, ratefac_val, enh_val, visc_val, &
366  sm_coeff_1, sm_coeff_2, sm_coeff_3)
367 
368 use sico_types
369 
370 implicit none
371 
372 real(dp) :: fct_visc_sm_deriv
373 real(dp), intent(in) :: de_val_m
374 real(dp), intent(in) :: ratefac_val, enh_val
375 real(dp), intent(in) :: visc_val
376 real(dp), intent(in) :: sm_coeff_1, sm_coeff_2, sm_coeff_3
377 
378 real(dp) :: de_visc_factor
379 
380 real(dp), parameter :: twenty_over_three = 6.666666666666667_dp
381 
382 de_visc_factor = de_val_m*de_val_m*visc_val*visc_val
383 
384 fct_visc_sm_deriv = 2.0_dp*sm_coeff_1*enh_val*ratefac_val &
385  + 24.0_dp*sm_coeff_2*enh_val*ratefac_val*de_visc_factor &
386  * ( 1.0_dp + twenty_over_three*sm_coeff_3*de_visc_factor )
387 
388 end function fct_visc_sm_deriv
389 
390 #endif
391 !
Generated on Mon May 20 2013 10:18:43 for SICOPOLIS V3.0 by   doxygen 1.8.1