texactV.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
---
texactV.py (5278B)
---
1 #!/usr/bin/env python3
2
3 from pylab import figure, plot, xlabel, ylabel, title, show, axis, linspace, hold, subplot, grid, step
4 import numpy as np
5 import subprocess
6 import shlex
7 import sys
8
9 from netCDF4 import Dataset as NC
10
11 # Setup
12
13 secpera = 3.15569259747e7 # seconds per year
14 rho_sw = 1028.0 # sea water density
15 rho_ice = 910.0 # ice density
16 standard_gravity = 9.81 # g
17 B0 = 1.9e8 # ice hardness
18
19 # "typical constant ice parameter" as defined in the paper and in Van der
20 # Veen's "Fundamentals of Glacier Dynamics", 1999
21 C = (rho_ice * standard_gravity * (1.0 - rho_ice / rho_sw) / (4 * B0)) ** 3
22
23 # upstream ice thickness
24 H0 = 600.0 # meters
25 # upstream ice velocity
26 v0 = 300.0 / secpera # 300 meters/year
27 # upstream ice flux
28 Q0 = H0 * v0
29
30
31 def H(x):
32 """Ice thickness."""
33 return (4 * C / Q0 * x + 1 / H0 ** 4) ** (-0.25)
34
35
36 def v(x):
37 """Ice velocity."""
38 return Q0 / H(x)
39
40
41 def x_c(t):
42 """Location of the calving front."""
43 return Q0 / (4 * C) * ((3 * C * t + 1 / H0 ** 3) ** (4.0 / 3.0) - 1 / H0 ** 4)
44
45
46 def plot_xc(t_years):
47 """Plot the location of the calving front."""
48 x = x_c(t_years * secpera) / 1000.0 # convert to km
49 _, _, y_min, y_max = axis()
50
51 hold(True)
52 plot([x, x], [y_min, y_max], '--g')
53
54
55 def run_pismv(Mx, run_length, options, output):
56 command = "pismv -test V -y %f -Mx %d %s -o %s" % (run_length, Mx, options, output)
57 print("Running %s" % command)
58 subprocess.call(shlex.split(command))
59
60
61 def plot_pism_results(filename, figure_title, color, same_figure=False):
62 nc = NC(filename)
63
64 time = nc.variables['time'][0] / secpera # convert to years
65
66 thk = nc.variables['thk'][0, 1, 2:]
67 ubar_ssa = nc.variables['velbar_mag'][0, 1, 2:]
68 x = nc.variables['x'][:]
69 dx = x[1] - x[0]
70 Lx = (x[-1] - x[0]) / 2.0
71 x_nc = (x[2:] + Lx - 2 * dx) / 1000.0
72
73 hold(True)
74
75 if same_figure == False:
76 figure(1)
77
78 subplot(211)
79 title(figure_title)
80 plotter(x_nc, H(x_nc * 1000.0), color='black', linestyle='dashed')
81 plotter(x_nc, thk, color=color, linewidth=2)
82 plot_xc(time)
83 ylabel("Ice thickness, m")
84 axis(xmin=0, xmax=400, ymax=600)
85 grid(True)
86
87 subplot(212)
88 plotter(x_nc, v(x_nc * 1000.0) * secpera, color='black', linestyle='dashed')
89 plotter(x_nc, ubar_ssa, color=color, linewidth=2)
90 plot_xc(time)
91 axis(xmin=0, xmax=400, ymax=1000)
92 xlabel("km")
93 ylabel("ice velocity, m/year")
94 grid(True)
95
96 nc.close()
97
98
99 import argparse
100 # Set up the option parser
101 parser = argparse.ArgumentParser(formatter_class=argparse.RawDescriptionHelpFormatter)
102 parser.description = """Manages PISM runs reproducing Figure 6 in Albrecht et al
103 'Parameterization for subgrid-scale motion of ice-shelf calving fronts', 2011"""
104
105 parser.epilog = """Model "variants":
106
107 0: no subgrid parameterization, no stress boundary condition at the calving front
108 1: -cfbc -part_grid
109 2: -cfbc -part_grid -part_grid_reduce_frontal_thickness
110
111 Here -part_grid_reduce_frontal_thickness adjusts the thickness
112 threshold used to decide when a 'partially filled' cell becomes full.
113 This is done to try and match the van der Veen profile this particular
114 setup is based on. Don't use it."""
115
116 parser.add_argument("-v", dest="variant", type=int,
117 help="choose the 'model variant', choose from 0, 1, 2", default=2)
118 parser.add_argument("-Mx", dest="Mx", type=int,
119 help="number of grid points", default=201)
120 parser.add_argument("-y", dest="y", type=float,
121 help="run length", default=300)
122 parser.add_argument("-s", dest="step_plot", action='store_true',
123 help="use 'plt.step()' to plot")
124 options = parser.parse_args()
125
126 Mx = options.Mx
127 x = linspace(0, 400e3, Mx)
128 run_length = options.y
129
130 opt = "-ssa_method fd -Lx 250 -o_order zyx"
131 extras = " -extra_file ex.nc -extra_vars flux_mag,thk,nuH,flux_divergence,velbar -extra_times 1"
132
133 if options.step_plot:
134 plotter = step
135 else:
136 plotter = plot
137
138 if options.variant == 0:
139 run_pismv(Mx, run_length, opt, "out.nc")
140 plot_pism_results("out.nc", "Figure 6 (a-b) (control)", 'blue')
141
142 opt = opt + extras
143 run_pismv(Mx, run_length, opt + " -max_dt 1", "out.nc")
144 plot_pism_results("out.nc", "Figure 6 (a-b) (control)", 'green', same_figure=True)
145 elif options.variant == 1:
146 opt += " -part_grid -cfbc"
147 run_pismv(Mx, run_length, opt, "out.nc")
148 plot_pism_results("out.nc", "Figure 6 (c-d) (-part_grid)", 'blue')
149
150 opt = opt + extras
151 run_pismv(Mx, run_length, opt + " -max_dt 1", "out.nc")
152 plot_pism_results("out.nc", "Figure 6 (c-d) (-part_grid)", 'green', same_figure=True)
153 elif options.variant == 2:
154 opt += " -cfbc -part_grid -part_grid_reduce_frontal_thickness"
155 run_pismv(Mx, run_length, opt, "out.nc")
156 plot_pism_results("out.nc", "Figure 6 (e-f) (-part_grid, reduce frontal thickness)", 'blue')
157
158 opt = opt + extras
159 run_pismv(Mx, run_length, opt + " -max_dt 1", "out.nc")
160 plot_pism_results("out.nc", "Figure 6 (e-f) (-part_grid)", 'green', same_figure=True)
161 else:
162 print("Wrong variant number. Choose one of 0, 1, 2.")
163 sys.exit(1)
164
165 show()