tpreprocess.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
---
tpreprocess.py (6634B)
---
1 #!/usr/bin/env python3
2
3 # Copyright (C) 2013, 2014, 2016, 2018 the PISM Authors
4
5 # This script sets up the bootstrap file.
6 # See also preprocess.sh.
7
8 import sys
9 import time
10 import subprocess
11 import argparse
12 import shlex
13 import numpy as np
14
15 try:
16 from netCDF4 import Dataset as CDF
17 except:
18 print("netCDF4 is not installed!")
19 sys.exit(1)
20
21 parser = argparse.ArgumentParser(description='Preprocess for validation using constant flux experiment from Sayag & Worster (2013). Creates PISM-readable bootstrap file and a configuration overrides file.',
22 formatter_class=argparse.ArgumentDefaultsHelpFormatter)
23 parser.add_argument('-Mx', default=52,
24 help='number of points in each direction on a square grid; note MX -> cell width cases: 52 -> 10mm, 104 -> 5mm, 209 -> 2.5mm, 520 -> 1mm')
25 parser.add_argument('-o', metavar='FILENAME', default='initlab52.nc',
26 help='output file name to create (NetCDF)')
27 args = parser.parse_args()
28
29
30 def create_config():
31 print(" creating PISM-readable config override file gumparams.nc ...")
32 nc = CDF("gumparams.nc", 'w')
33 config = nc.createVariable("pism_overrides", 'i4')
34
35 attrs = {
36 "constants.standard_gravity": 9.81,
37 "constants.standard_gravity_doc": "m s-2; = g",
38
39 "constants.ice.density": 1000.0,
40 "constants.ice.density_doc": "kg m-3; 1% Xanthan gum in water has same density as water",
41
42 "stress_balance.sia.bed_smoother.range": -1.0,
43 "stress_balance.sia.bed_smoother.range_doc": "m; negative value de-activates bed smoother",
44
45 "bootstrapping.defaults.geothermal_flux": 0.0,
46 "bootstrapping.defaults.geothermal_flux_doc": "W m-2; no geothermal",
47
48 "output.runtime.time_unit_name": "second",
49 "output.runtime.time_unit_name_doc": "stdout uses seconds (not years) to show model time",
50
51 "output.runtime.time_use_calendar": "no",
52 "output.runtime.time_use_calendar_doc": "stdout does not use a calendar to show model time",
53
54 "output.runtime.volume_scale_factor_log10": -15,
55 "output.runtime.volume_scale_factor_log10_doc": "; an integer; log base 10 of scale factor to use for volume in summary line to stdout; -15 gives volume in cm^3",
56
57 "output.runtime.area_scale_factor_log10": -10,
58 "output.runtime.area_scale_factor_log10_doc": "; an integer; log base 10 of scale factor to use for area in summary line to stdout; -10 gives area in cm^2",
59
60 "geometry.ice_free_thickness_standard": 1e-8,
61 "geometry.ice_free_thickness_standard_doc": "m; only if the fluid is less than this is a cell marked as ice free",
62
63 "output.ice_free_thickness_standard": 1e-8,
64 "output.ice_free_thickness_standard_doc": "fluid layer exceeding this thickness is included in computations of areas and volumes",
65
66 "time_stepping.adaptive_ratio": 0.08,
67 "time_stepping.adaptive_ratio_doc": "; compare default 0.12; needs to be smaller because gum suspension is more shear-thinning than ice?",
68
69 "stress_balance.sia.Glen_exponent": 5.9,
70 "stress_balance.sia.Glen_exponent_doc": "; : n; Sayag & Worster (2013) give n = 5.9 +- 0.2",
71
72 "flow_law.isothermal_Glen.ice_softness": 9.7316e-09, # vs (e.g.) 4e-25 Pa-3 s-1 for ice
73 "ice_softness_doc": "Pa-n s-1; = A_0 = B_0^(-n) = (2 x 11.4 Pa s^(1/n))^(-n); Sayag & Worster (2013) give B_0/2 = tilde mu = 11.4 +- 0.25 Pa s^(1/n)"
74 }
75
76 keys = list(attrs.keys())
77 keys.sort()
78 for k in keys:
79 config.setncattr(k, attrs[k])
80
81 nc.close()
82
83
84 create_config()
85
86 # lab setup is table with hole in the middle into which is piped the
87 # shear-thinning fluid, which is Xanthan gum 1% solution
88 Lx = 260.0e-3 # m; = 260 mm; maximum observed radius is 25.2 cm so we go out just a bit
89 Ly = Lx # square table
90 flux = 3.8173e-3 # kg s-1; = 3.8173 g s-1; Sayag personal communication
91 pipeR = 8.0e-3 # m; = 8 mm; input pipe has this radius; Sayag personal communication
92 temp = 20.0 # C; fluid is at 20 deg (though it should not matter)
93
94 # set up the grid:
95 Mx = int(args.Mx)
96 My = Mx
97 print(" creating grid of Mx = %d by My = %d points ..." % (Mx, My))
98 dx = (2.0 * Lx) / float(Mx)
99 dy = (2.0 * Ly) / float(My)
100 print(" cells have dimensions dx = %.3f mm by dy = %.3f mm ..." % (dx * 1000.0, dy * 1000.0))
101 x = np.linspace(-Lx - dx / 2.0, Lx + dx / 2.0, Mx)
102 y = np.linspace(-Ly - dy / 2.0, Ly + dy / 2.0, My)
103
104 # create dummy fields
105 [xx, yy] = np.meshgrid(x, y) # if there were "ndgrid" in numpy we would use it
106
107 topg = np.zeros((Mx, My))
108 thk = np.zeros((Mx, My)) # no fluid on table at start
109 artm = np.zeros((Mx, My)) + 273.15 + temp # 20 degrees Celsius
110
111 # smb = flux as m s-1, but scaled so that the total is correct even on a coarse grid
112 smb = np.zeros((Mx, My))
113 smb[xx ** 2 + yy ** 2 <= pipeR ** 2 + 1.0e-10] = 1.0
114 smbpos = sum(sum(smb))
115 if smbpos == 0:
116 print("gridding ERROR: no cells have positive input flux ... ending now")
117 sys.exit(1)
118 else:
119 print(" input flux > 0 at %d cells ..." % smbpos)
120 smb = (flux / (smbpos * dx * dy)) * smb # [flux] = kg s-1 so now [smb] = kg m-2 s-1
121
122 # Write the data:
123 nc = CDF(args.o, "w", format='NETCDF3_CLASSIC') # for netCDF4 module
124
125 # Create dimensions x and y
126 nc.createDimension("x", size=Mx)
127 nc.createDimension("y", size=My)
128
129 x_var = nc.createVariable("x", 'f4', dimensions=("x",))
130 x_var.units = "m"
131 x_var.long_name = "easting"
132 x_var.standard_name = "projection_x_coordinate"
133 x_var[:] = x
134
135 y_var = nc.createVariable("y", 'f4', dimensions=("y",))
136 y_var.units = "m"
137 y_var.long_name = "northing"
138 y_var.standard_name = "projection_y_coordinate"
139 y_var[:] = y
140
141 fill_value = np.nan
142
143
144 def def_var(nc, name, units, fillvalue):
145 # dimension transpose is standard: "float thk(y, x)" in NetCDF file
146 var = nc.createVariable(name, 'f', dimensions=("y", "x"), fill_value=fillvalue)
147 var.units = units
148 return var
149
150
151 bed_var = def_var(nc, "topg", "m", fill_value)
152 bed_var.standard_name = "bedrock_altitude"
153 bed_var[:] = topg
154
155 thk_var = def_var(nc, "thk", "m", fill_value)
156 thk_var.standard_name = "land_ice_thickness"
157 thk_var[:] = thk
158
159 smb_var = def_var(nc, "climatic_mass_balance", "kg m-2 s-1", fill_value)
160 smb_var.standard_name = "land_ice_surface_specific_mass_balance_flux"
161 smb_var[:] = smb
162
163 artm_var = def_var(nc, "ice_surface_temp", "K", fill_value)
164 artm_var[:] = artm
165
166 # set global attributes
167 nc.Conventions = 'CF-1.4'
168 historysep = ' '
169 historystr = time.asctime() + ': ' + historysep.join(sys.argv) + '\n'
170 setattr(nc, 'history', historystr)
171
172 nc.close()
173
174 print(' ... PISM-bootable NetCDF file %s written' % args.o)