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()