tgenerate_input.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 --- tgenerate_input.py (2180B) --- 1 #!/usr/bin/env python3 2 from PISMNC import PISMDataset as NC 3 import numpy as np 4 import argparse 5 6 parser = argparse.ArgumentParser() 7 parser.description = "Generates an input file for the 'tongues' experiment." 8 9 parser.add_argument("-M", dest="M", type=int, help="grid size", default=101) 10 parser.add_argument("-L", dest="L", type=float, help="domain size, meters", default=1e5) 11 parser.add_argument("-o", dest="output", help="output file name", default="tongues.nc") 12 options = parser.parse_args() 13 14 L = options.L 15 M = options.M 16 17 # grid 18 x = np.linspace(-L, L, M) 19 y = np.linspace(-L, L, M) 20 xx, yy = np.meshgrid(x, y) 21 22 23 def tongue(xx, x0, width): 24 "create one ice tongue" 25 result = np.zeros_like(xx) 26 result[:, x0:x0 + width] = 100.0 27 return result 28 29 30 thk = np.zeros_like(xx) 31 32 x0 = 3 33 width = 1 34 spacing = 5 35 while x0 + width < M - 1: 36 thk += tongue(xx, x0, width) 37 x0 += width + spacing 38 width += 1 39 40 # make tongues shorter 41 thk[5:, :] = 0 42 43 bc_mask = np.zeros_like(thk) 44 bc_mask[thk > 0] = 1 45 46 # make the bed deep everywhere except in icy areas, where it is barely 47 # grounded 48 z = np.zeros_like(xx) - 1000.0 49 z[thk > 0] = -(910.0 / 1028.0) * 100.0 + 1 50 51 # Velocity Dirichlet B.C.: 52 ubar = np.zeros_like(thk) 53 vbar = np.zeros_like(thk) 54 vbar[bc_mask == 1] = 100.0 55 56 try: 57 nc = NC(options.output, 'w') 58 59 nc.create_dimensions(x, y, time_dependent=False) 60 61 nc.define_2d_field("topg", attrs={"units": "m", 62 "long_name": "bedrock topography"}) 63 nc.define_2d_field("thk", attrs={"units": "m", 64 "long_name": "ice thickness"}) 65 66 nc.define_2d_field("climatic_mass_balance", attrs={"units": "kg m-2 year-1"}) 67 nc.define_2d_field("ice_surface_temp", attrs={"units": "Celsius"}) 68 69 nc.define_2d_field("u_ssa_bc", attrs={"units": "m/year"}) 70 nc.define_2d_field("v_ssa_bc", attrs={"units": "m/year"}) 71 except: 72 nc = NC(options.output, 'a') 73 74 nc.write("topg", z) 75 nc.write("thk", thk) 76 nc.write("climatic_mass_balance", np.zeros_like(xx)) 77 nc.write("ice_surface_temp", np.zeros_like(xx) - 30.0) # irrelevant 78 nc.write("u_ssa_bc", ubar) 79 nc.write("v_ssa_bc", vbar) 80 nc.write("bc_mask", bc_mask) 81 82 nc.close()