MISMIP.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)