Introduction
Introduction Statistics Contact Development Disclaimer Help
tMISMIP.py - pism-exp-gsw - ice stream and sediment transport experiments
git clone git://src.adamsgaard.dk/pism-exp-gsw
Log
Files
Refs
README
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 func…
5 computing theoretical steady state profiles corresponding to various MIS…
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, …
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, …
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 …
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…
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 …
182 - theta * s * b_x / h_f
183 - A(experiment, step) * (rho_g * (1.0 - r) / 4.0) ** n() * h…
184
185
186 def x_g(experiment, step, theta=0.0):
187 """Computes the theoretical grounding line location using Newton's m…
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 c…
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 **…
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 …
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…
241
242 # now 'result' contains thickness in reverse order, including the gr…
243 # line point (which is not on the grid); discard it and reverse the …
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) *…
257 * ((q_0 + a() * (x_floating - xg)) ** (n() + 1) - q_0 ** (n() + …
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,…
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='b…
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 corresp…
303 parser.add_option("-e", "--experiment", dest="experiment", type="str…
304 default='1a',
305 help="MISMIP experiments (one of '1a', '1b', '2a',…
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)
You are viewing proxied material from mx1.adamsgaard.dk. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.