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 (2035B)
       ---
            1 #!/usr/bin/env python3
            2 
            3 from PISMNC import PISMDataset as NC
            4 import numpy as np
            5 
            6 import argparse
            7 
            8 parser = argparse.ArgumentParser()
            9 parser.description = "Generates an input file for the 'flood' experiment."
           10 
           11 parser.add_argument("-M", dest="M", type=int,
           12                     help="grid size", default=301)
           13 parser.add_argument("-L", dest="L", type=float,
           14                     help="domain size, meters", default=1e6)
           15 parser.add_argument("-o", dest="output", help="output file name",
           16                     default="flood.nc")
           17 options = parser.parse_args()
           18 L = options.L
           19 M = options.M
           20 
           21 x = np.linspace(-L, L, M)
           22 y = np.linspace(-L, L, M)
           23 
           24 xx, yy = np.meshgrid(x, y)
           25 
           26 z = np.zeros_like(xx)
           27 
           28 thk = np.zeros_like(z)
           29 for phi in np.linspace(0, 2 * np.pi, 4):
           30     R = 0.5 * L
           31     r = 0.2 * L
           32     center_x = R * np.cos(phi)
           33     center_y = R * np.sin(phi)
           34     T = phi / (2.0 * np.pi) * 1500.0
           35     thk += (T / r) * np.sqrt(np.maximum(r ** 2 - (xx - center_x) ** 2 - (yy - center_y) ** 2, 0))
           36 
           37 try:
           38     nc = NC(options.output, 'w')
           39 except:
           40     nc = NC(options.output, 'a')
           41 
           42 try:
           43     nc.create_dimensions(x, y, time_dependent=False)
           44 
           45     nc.define_2d_field("topg", attrs={"units": "m",
           46                                       "long_name": "bedrock topography"})
           47     nc.define_2d_field("thk", attrs={"units": "m",
           48                                      "long_name": "ice thickness"})
           49 
           50     nc.define_2d_field("climatic_mass_balance", attrs={"units": "kg m-2 year-1"})
           51     nc.define_2d_field("ice_surface_temp", attrs={"units": "Kelvin"})
           52 except:
           53     pass
           54 
           55 nc.write("topg", z)
           56 nc.write("thk", thk)
           57 nc.write("climatic_mass_balance", np.zeros_like(xx))
           58 nc.write("ice_surface_temp", np.zeros_like(xx) + 273.15 - 30.0)
           59 
           60 # generate sea level data
           61 time = np.linspace(0, 1000, 1001)
           62 sea_level = np.linspace(0, 700, 1001)
           63 
           64 nc.createDimension("time")
           65 nc.define_timeseries("time", attrs={"units": "years since 1-1-1"})
           66 nc.write_timeseries("time", time)
           67 
           68 nc.define_timeseries("delta_SL", attrs={"units": "meters"})
           69 nc.write_timeseries("delta_SL", sea_level)
           70 
           71 nc.close()