tpiktests_utils.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
---
tpiktests_utils.py (3649B)
---
1 #!/usr/bin/env python3
2
3 import numpy as np
4 import argparse
5
6 # constants
7
8
9 class Parameters:
10 secpera = 3.15569259747e7 # seconds per year
11 standard_gravity = 9.81 # m/s^2
12 B0 = 1.9e8 # ice hardness
13
14 rho_ice = 910.0 # in kg/m^3
15 rho_ocean = 1028.0 # in kg/m^3
16
17 vel_bc = 300 # m/year
18 accumulation_rate = 0.3 # m/year
19 air_temperature = 247.0 # Kelvin
20 domain_size = 1000.0 # km
21 topg_min = -3000.0 # m
22
23 # "typical constant ice parameter" as defined in the paper and in Van der
24 # Veen's "Fundamentals of Glacier Dynamics", 1999
25 C = (rho_ice * standard_gravity * (1.0 - rho_ice / rho_ocean) / (4 * B0)) ** 3
26 H0 = 600 # ice thickness at the grounding line
27 Q0 = (vel_bc / secpera) * H0 # flux at the grounding line
28
29 r_gl = 0.25e6 # grounding line radius, meters
30 r_cf = 0.45e6 # calving front radius, meters
31
32
33 def process_options(default_filename, domain_size):
34 parser = argparse.ArgumentParser()
35 parser.add_argument("-o", dest="output_filename",
36 default=default_filename)
37 parser.add_argument("-Mx", dest="Mx", type=int, default=301)
38 parser.add_argument("-My", dest="My", type=int, default=301)
39 parser.add_argument("-domain_size", dest="domain_size", type=float,
40 default=domain_size)
41 parser.add_argument("-shelf", dest="shelf", action="store_true")
42 parser.add_argument("-square", dest="square", action="store_true")
43 return parser.parse_args()
44
45
46 def create_grid(options):
47 L = options.domain_size
48 dx = L / (options.Mx - 1) * 1000.0 # meters
49 dy = L / (options.My - 1) * 1000.0 # meters
50 print("dx = %.2f km, dy = %.2f km" % (dx / 1000.0, dy / 1000.0))
51 x = np.linspace(-L / 2 * 1000.0, L / 2 * 1000.0, options.Mx)
52 y = np.linspace(-L / 2 * 1000.0, L / 2 * 1000.0, options.My)
53
54 return (dx, dy, x, y)
55
56
57 def prepare_output_file(nc, x, y, include_vel_bc=True):
58 nc.create_dimensions(x, y)
59
60 attrs = {'long_name': "ice thickness",
61 "units": "m",
62 "standard_name": "land_ice_thickness",
63 "_FillValue": 1.0}
64 nc.define_2d_field("thk", attrs=attrs)
65
66 attrs = {'long_name': "bedrock surface elevation",
67 "units": "m",
68 "standard_name": "bedrock_altitude",
69 "_FillValue": -600.0}
70 nc.define_2d_field("topg", attrs=attrs)
71
72 attrs = {'long_name': "mean annual temperature at ice surface",
73 "units": "Kelvin",
74 "_FillValue": 248.0}
75 nc.define_2d_field("ice_surface_temp", attrs=attrs)
76
77 ice_density = 910.0
78 attrs = {'long_name': "mean annual net ice equivalent accumulation rate",
79 "units": "kg m-2 year-1",
80 "standard_name": "land_ice_surface_specific_mass_balance_flux",
81 "_FillValue": 0.2 * ice_density} # 0.2 m/year
82 nc.define_2d_field("climatic_mass_balance", attrs=attrs)
83
84 if include_vel_bc == False:
85 return
86
87 attrs = {'long_name': "Dirichlet boundary condition locations",
88 "units": "1",
89 "_FillValue": 0}
90 nc.define_2d_field("bc_mask", attrs=attrs)
91
92 attrs = {'long_name': "X-component of the SSA velocity boundary conditions",
93 "units": "m/year",
94 "_FillValue": 0.0}
95 nc.define_2d_field("u_ssa_bc", attrs=attrs)
96
97 attrs['long_name'] = "Y-component of the SSA velocity boundary conditions"
98 nc.define_2d_field("v_ssa_bc", attrs=attrs)
99
100
101 def write_data(nc, variables):
102 for name in list(variables.keys()):
103 nc.write(name, variables[name])