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