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)