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()