tprepare.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
---
tprepare.py (6087B)
---
1 #!/usr/bin/env python3
2
3 try:
4 from netCDF4 import Dataset as NC
5 except:
6 print("netCDF4 is not installed!")
7 sys.exit(1)
8
9 import MISMIP
10
11 import numpy as np
12
13
14 def bed_topography(experiment, x):
15 """Computes bed elevation as a function of x.
16 experiment can be '1a', '1b', '2a', '2b', '3a', '3b'.
17 """
18
19 return np.tile(-MISMIP.b(experiment, np.abs(x)), (3, 1))
20
21
22 def surface_mass_balance(x):
23 """Computes the surface mass balance."""
24 return np.tile(np.zeros_like(x) + MISMIP.a(), (3, 1)) * MISMIP.rho_i()
25
26
27 def ice_surface_temp(x):
28 """Computes the ice surface temperature (irrelevant)."""
29 return np.tile(np.zeros_like(x) + 273.15, (3, 1))
30
31
32 def x(mismip_mode, N=None):
33 if mismip_mode in (1, 2):
34 return np.linspace(-MISMIP.L(), MISMIP.L(), 2 * MISMIP.N(mismip_mode) + 1)
35
36 return np.linspace(-MISMIP.L(), MISMIP.L(), N)
37
38
39 def y(x):
40 """Computes y coordinates giving the 1:1 aspect ratio.
41 Takes cross-flow grid periodicity into account."""
42 dx = x[1] - x[0]
43 dy = dx
44 return np.array([-dy, 0, dy])
45
46
47 def thickness(experiment, step, x, calving_front=1750e3, semianalytical_profile=True):
48
49 # we expect x to have an odd number of grid points so that one of them is
50 # at 0
51 if x.size % 2 != 1:
52 raise ValueError("x has to have an odd number of points, got %d", x.size)
53
54 x_nonnegative = x[x >= 0]
55 if not semianalytical_profile:
56 thk_nonnegative = np.zeros_like(x_nonnegative) + 10
57 else:
58 thk_nonnegative = MISMIP.thickness(experiment, step, x_nonnegative)
59
60 thk_nonnegative[x_nonnegative > calving_front] = 0
61
62 thk = np.zeros_like(x)
63 thk[x >= 0] = thk_nonnegative
64 thk[x < 0] = thk_nonnegative[:0:-1]
65
66 return np.tile(thk, (3, 1))
67
68
69 def pism_bootstrap_file(filename, experiment, step, mode,
70 calving_front=1750e3, N=None, semianalytical_profile=True):
71 import PISMNC
72
73 xx = x(mode, N)
74 yy = y(xx)
75
76 bed_elevation = bed_topography(experiment, xx)
77 ice_thickness = thickness(experiment, step, xx, calving_front, semianalytical_profile)
78 ice_surface_mass_balance = surface_mass_balance(xx)
79 ice_surface_temperature = ice_surface_temp(xx)
80
81 ice_extent = np.zeros_like(ice_thickness)
82 ice_extent[ice_thickness > 0] = 1
83 ice_extent[bed_elevation > 0] = 1
84
85 nc = PISMNC.PISMDataset(filename, 'w', format="NETCDF3_CLASSIC")
86
87 nc.create_dimensions(xx, yy)
88
89 nc.define_2d_field('topg',
90 attrs={'units': 'm',
91 'long_name': 'bedrock surface elevation',
92 'standard_name': 'bedrock_altitude'})
93 nc.write('topg', bed_elevation)
94
95 nc.define_2d_field('thk',
96 attrs={'units': 'm',
97 'long_name': 'ice thickness',
98 'standard_name': 'land_ice_thickness'})
99 nc.write('thk', ice_thickness)
100
101 nc.define_2d_field('climatic_mass_balance',
102 attrs={'units': 'kg m-2 / s',
103 'long_name': 'ice-equivalent surface mass balance (accumulation/ablation) rate',
104 'standard_name': 'land_ice_surface_specific_mass_balance_flux'})
105 nc.write('climatic_mass_balance', ice_surface_mass_balance)
106
107 nc.define_2d_field('ice_surface_temp',
108 attrs={'units': 'Kelvin',
109 'long_name': 'annual average ice surface temperature, below firn processes'})
110 nc.write('ice_surface_temp', ice_surface_temperature)
111
112 nc.define_2d_field('land_ice_area_fraction_retreat',
113 attrs={'units': '1',
114 'long_name': 'mask defining the maximum ice extent'})
115 nc.write('land_ice_area_fraction_retreat', ice_extent)
116
117 nc.close()
118
119
120 if __name__ == "__main__":
121
122 from optparse import OptionParser
123
124 parser = OptionParser()
125
126 parser.usage = "%prog [options]"
127 parser.description = "Creates a MISMIP boostrapping file for use with PISM."
128 parser.add_option("-o", dest="output_filename",
129 help="output file")
130 parser.add_option("-e", "--experiment", dest="experiment", type="string",
131 help="MISMIP experiment (one of '1a', '1b', '2a', '2b', '3a', '3b')",
132 default="1a")
133 parser.add_option("-s", "--step", dest="step", type="int", default=1,
134 help="MISMIP step")
135 parser.add_option("-u", "--uniform_thickness",
136 action="store_false", dest="semianalytical_profile", default=True,
137 help="Use uniform 10 m ice thickness")
138 parser.add_option("-m", "--mode", dest="mode", type="int", default=2,
139 help="MISMIP grid mode")
140 parser.add_option("-N", dest="N", type="int", default=3601,
141 help="Custom grid size; use with --mode=3")
142 parser.add_option("-c", dest="calving_front", type="float", default=1600e3,
143 help="Calving front location, in meters (e.g. 1600e3)")
144
145 (opts, args) = parser.parse_args()
146
147 experiments = ('1a', '1b', '2a', '2b', '3a', '3b')
148 if opts.experiment not in experiments:
149 print("Invalid experiment %s. Has to be one of %s." % (
150 opts.experiment, experiments))
151 exit(1)
152
153 if not opts.output_filename:
154 output_filename = "MISMIP_%s_%d_%d.nc" % (opts.experiment,
155 opts.step,
156 opts.mode)
157 else:
158 output_filename = opts.output_filename
159
160 print("Creating MISMIP setup for experiment %s, step %s, grid mode %d in %s..." % (
161 opts.experiment, opts.step, opts.mode, output_filename))
162
163 pism_bootstrap_file(output_filename,
164 opts.experiment,
165 opts.step,
166 opts.mode,
167 calving_front=opts.calving_front,
168 N=opts.N,
169 semianalytical_profile=opts.semianalytical_profile)
170
171 print("done.")