GopherProxy

	tprepare.py - pism-exp-gsw - ice stream and sediment transport experiments
	git clone git://src.adamsgaard.dk/pism-exp-gsw
	Log
	Files
	Refs
	README
	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_…
	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_…
	48
	49 # we expect x to have an odd number of grid points so that one of th…
	50 # at 0
	51 if x.size % 2 != 1:
	52 raise ValueError("x has to have an odd number of points, got %d"…
	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_nonnegati…
	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_pro…
	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, semia…
	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 …
	104 'standard_name': 'land_ice_surface_specifi…
	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 t…
	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 ic…
	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 wit…
	128 parser.add_option("-o", dest="output_filename",
	129 help="output file")
	130 parser.add_option("-e", "--experiment", dest="experiment", type="str…
	131 help="MISMIP experiment (one of '1a', '1b', '2a', …
	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…
	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=…
	143 help="Calving front location, in meters (e.g. 1600…
	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 %…
	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_profi…
	170
	171 print("done.")