tMISMIP.py - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
(HTM) git clone git://src.adamsgaard.dk/pism
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tMISMIP.py (9438B)
---
1 #!/usr/bin/env python3
2 import numpy as np
3
4 """This module contains MISMIP constants and parameters, as well as functions
5 computing theoretical steady state profiles corresponding to various MISMIP
6 experiments.
7
8 It should not be cluttered with plotting or NetCDF output code.
9 """
10
11
12 def secpera():
13 "Number of seconds per year."
14 return 3.15569259747e7
15
16
17 def L():
18 "The length of the MISMIP domain."
19 return 1800e3
20
21
22 def N(mode):
23 "Number of grid spaces corresponding to a MISMIP 'mode.'"
24 if mode == 1:
25 return 150
26
27 if mode == 2:
28 return 1500
29
30 raise ValueError("invalid mode (%s)" % mode)
31
32
33 def A(experiment, step):
34 """Ice softness parameter for given experiment and step."""
35 A1 = np.array([4.6416e-24, 2.1544e-24, 1.0e-24,
36 4.6416e-25, 2.1544e-25, 1.0e-25,
37 4.6416e-26, 2.1544e-26, 1.0e-26])
38 # Values of A to be used in experiments 1 and 2.
39
40 A3a = np.array([3.0e-25, 2.5e-25, 2.0e-25,
41 1.5e-25, 1.0e-25, 5.0e-26,
42 2.5e-26, 5.0e-26, 1.0e-25,
43 1.5e-25, 2.0e-25, 2.5e-25,
44 3.0e-25])
45 # Values of A to be used in experiment 3a.
46
47 A3b = np.array([1.6e-24, 1.4e-24, 1.2e-24,
48 1.0e-24, 8.0e-25, 6.0e-25,
49 4.0e-25, 2.0e-25, 4.0e-25,
50 6.0e-25, 8.0e-25, 1.0e-24,
51 1.2e-24, 1.4e-24, 1.6e-24])
52 # Values of A to be used in experiment 3b.
53
54 try:
55 if experiment in ("1a", "1b", "2a", "2b"):
56 return A1[step - 1]
57
58 if experiment == "3a":
59 return A3a[step - 1]
60
61 if experiment == "3b":
62 return A3b[step - 1]
63 except:
64 raise ValueError("invalid step (%s) for experiment %s" % (step, experiment))
65
66 raise ValueError("invalid experiment (%s)" % experiment)
67
68
69 def run_length(experiment, step):
70 """Returns the time interval for an experiment 3 step."""
71 T3a = np.array([3.0e4, 1.5e4, 1.5e4,
72 1.5e4, 1.5e4, 3.0e4,
73 3.0e4, 1.5e4, 1.5e4,
74 3.0e4, 3.0e4, 3.0e4,
75 1.5e4])
76 # Time intervals to be used in experiment 3a.
77
78 T3b = np.array([3.0e4, 1.5e4, 1.5e4,
79 1.5e4, 1.5e4, 1.5e4,
80 1.5e4, 3.0e4, 1.5e4,
81 1.5e4, 1.5e4, 1.5e4,
82 1.5e4, 3.0e4, 1.5e4])
83 # Time intervals to be used in experiment 3b.
84
85 try:
86 if experiment == "3a":
87 return T3a[step - 1]
88
89 if experiment == "3b":
90 return T3b[step - 1]
91 except:
92 raise ValueError("invalid step (%s) for experiment %s" % (step, experiment))
93
94 return 3e4
95
96
97 def rho_i():
98 "Ice density"
99 return 900.0
100
101
102 def rho_w():
103 "Water density"
104 return 1000.0
105
106
107 def g():
108 """Acceleration due to gravity. (Table 2 on page 19 of mismip_4.pdf
109 uses this value, i.e. g = 9.8 m s-2.)"""
110 return 9.8
111
112
113 def n():
114 "Glen exponent"
115 return 3.0
116
117
118 def a():
119 "Accumulation rate (m/s)"
120 return 0.3 / secpera()
121
122
123 def m(experiment):
124 "Sliding law exponent"
125 if experiment in ("1a", "2a", "3a"):
126 return 1 / 3.0
127
128 if experiment in ("1b", "2b", "3b"):
129 return 1.0
130
131 raise ValueError("invalid experiment (%s)" % experiment)
132
133
134 def C(experiment):
135 "Sliding law coefficient"
136 if experiment in ("1a", "2a", "3a"):
137 return 7.624e6
138
139 if experiment in ("1b", "2b", "3b"):
140 return 7.2082e10
141
142 raise ValueError("invalid experiment (%s)" % experiment)
143
144
145 def b(experiment, x):
146 "Bed depth below sea level. (-b(x) = topg(x))"
147
148 if experiment in ("1a", "1b", "2a", "2b"):
149 return -720. + 778.5 * (x / 7.5e5)
150
151 if experiment in ("3a", "3b"):
152 xx = x / 7.5e5
153 return -(729. - 2184.8 * xx ** 2. + 1031.72 * xx ** 4. - 151.72 * xx ** 6.)
154
155 raise ValueError("invalid experiment (%s)" % experiment)
156
157
158 def b_slope(experiment, x):
159 """The x-derivative of b(experiment, x)."""
160
161 if experiment in ("1a", "1b", "2a", "2b"):
162 return 778.5 / 7.5e5
163
164 if experiment in ("3a", "3b"):
165 xx = x / 7.5e5
166 return -(- 2184.8 * (2. / 7.5e5) * xx
167 + 1031.72 * (4. / 7.5e5) * xx ** 3.
168 - 151.72 * (6. / 7.5e5) * xx ** 5.)
169
170 raise ValueError("invalid experiment (%s)" % experiment)
171
172
173 def cold_function(experiment, step, x, theta=0.0):
174 """Evaluates function whose zeros define x_g in 'cold' steady marine sheet problem."""
175 r = rho_i() / rho_w()
176 h_f = r ** (-1.) * b(experiment, x)
177 b_x = b_slope(experiment, x)
178 s = a() * x
179 rho_g = rho_i() * g()
180 return (theta * a()
181 + C(experiment) * s ** (m(experiment) + 1.0) / (rho_g * h_f ** (m(experiment) + 2.))
182 - theta * s * b_x / h_f
183 - A(experiment, step) * (rho_g * (1.0 - r) / 4.0) ** n() * h_f ** (n() + 1.0))
184
185
186 def x_g(experiment, step, theta=0.0):
187 """Computes the theoretical grounding line location using Newton's method."""
188
189 # set the initial guess
190 if experiment in ("3a", "3b"):
191 x = 800.0e3
192 else:
193 x = 1270.0e3
194
195 delta_x = 10. # Finite difference step size (metres) for gradient calculation
196 tolf = 1.e-4 # Tolerance for finding zeros
197 eps = np.finfo(float).eps
198 normf = tolf + eps
199 toldelta = 1.e1 # Newton step size tolerance
200 dx = toldelta + 1.0
201
202 # this is just a shortcut
203 def F(x):
204 return cold_function(experiment, step, x, theta)
205
206 while (normf > tolf) or (abs(dx) > toldelta):
207 f = F(x)
208 normf = abs(f)
209 grad = (F(x + delta_x) - f) / delta_x
210 dx = -f / grad
211 x = x + dx
212
213 return x
214
215
216 def thickness(experiment, step, x, theta=0.0):
217 """Compute ice thickness for x > 0.
218 """
219 # compute the grounding line position
220 xg = x_g(experiment, step, theta)
221
222 def surface(h, x):
223 b_x = b_slope(experiment, x)
224 rho_g = rho_i() * g()
225 s = a() * np.abs(x)
226 return b_x - (C(experiment) / rho_g) * s ** m(experiment) / h ** (m(experiment) + 1)
227
228 # extract the grounded part of the grid
229 x_grounded = x[x < xg]
230
231 # We will integrate from the grounding line inland. odeint requires that
232 # the first point in x_grid be the one corresponding to the initial
233 # condition; append it and reverse the order.
234 x_grid = np.append(xg, x_grounded[::-1])
235
236 # use thickness at the grounding line as the initial condition
237 h_f = b(experiment, xg) * rho_w() / rho_i()
238
239 import scipy.integrate
240 thk_grounded = scipy.integrate.odeint(surface, [h_f], x_grid, atol=1.e-9, rtol=1.e-9)
241
242 # now 'result' contains thickness in reverse order, including the grounding
243 # line point (which is not on the grid); discard it and reverse the order.
244 thk_grounded = np.squeeze(thk_grounded)[:0:-1]
245
246 # extract the floating part of the grid
247 x_floating = x[x >= xg]
248
249 # compute the flux through the grounding line
250 q_0 = a() * xg
251
252 # Calculate ice thickness for shelf from van der Veen (1986)
253 r = rho_i() / rho_w()
254 rho_g = rho_i() * g()
255 numer = h_f * (q_0 + a() * (x_floating - xg))
256 base = q_0 ** (n() + 1) + h_f ** (n() + 1) * ((1 - r) * rho_g / 4) ** n() * A(experiment, step) \
257 * ((q_0 + a() * (x_floating - xg)) ** (n() + 1) - q_0 ** (n() + 1)) / a()
258 thk_floating = numer / (base ** (1.0 / (n() + 1)))
259
260 return np.r_[thk_grounded, thk_floating]
261
262
263 def plot_profile(experiment, step, out_file):
264 from pylab import figure, subplot, hold, plot, xlabel, ylabel, text, title, axis, vlines, savefig
265
266 if out_file is None:
267 out_file = "MISMIP_%s_A%d.pdf" % (experiment, step)
268
269 xg = x_g(experiment, step)
270
271 x = np.linspace(0, L(), N(2) + 1)
272 thk = thickness(experiment, step, x)
273 x_grounded, thk_grounded = x[x < xg], thk[x < xg]
274 x_floating, thk_floating = x[x >= xg], thk[x >= xg]
275
276 figure(1)
277 ax = subplot(111)
278 hold(True)
279 plot(x / 1e3, np.zeros_like(x), ls='dotted', color='red')
280 plot(x / 1e3, -b(experiment, x), color='black')
281 plot(x / 1e3, np.r_[thk_grounded - b(experiment, x_grounded),
282 thk_floating * (1 - rho_i() / rho_w())],
283 color='blue')
284 plot(x_floating / 1e3, -thk_floating * (rho_i() / rho_w()), color='blue')
285 _, _, ymin, ymax = axis(xmin=0, xmax=x.max() / 1e3)
286 vlines(xg / 1e3, ymin, ymax, linestyles='dashed', color='black')
287
288 xlabel('distance from the summit, km')
289 ylabel('elevation, m')
290 text(0.6, 0.9, "$x_g$ (theory) = %4.0f km" % (xg / 1e3),
291 color='black', transform=ax.transAxes)
292 title("MISMIP experiment %s, step %d" % (experiment, step))
293 savefig(out_file)
294
295
296 if __name__ == "__main__":
297 from optparse import OptionParser
298
299 parser = OptionParser()
300
301 parser.usage = "%prog [options]"
302 parser.description = "Plots the theoretical geometry profile corresponding to MISMIP experiment and step."
303 parser.add_option("-e", "--experiment", dest="experiment", type="string",
304 default='1a',
305 help="MISMIP experiments (one of '1a', '1b', '2a', '2b', '3a', '3b')")
306 parser.add_option("-s", "--step", dest="step", type="int", default=1,
307 help="MISMIP step number")
308 parser.add_option("-o", dest="out_file", help="output file name")
309
310 (opts, args) = parser.parse_args()
311
312 plot_profile(opts.experiment, opts.step, opts.out_file)