tcreateSetup_flowline.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
---
tcreateSetup_flowline.py (7260B)
---
1 #!/usr/bin/env python3
2 # Copyright (C) 2011, 2013, 2014, 2016, 2018 Torsten Albrecht and Moritz Huetten
3
4 # ./createSetup_flowline.py -a 0.0 -r 10.0
5
6 import sys
7 import getopt
8 import math
9 import numpy as np
10 try:
11 from netCDF4 import Dataset as NC
12 except:
13 print("netCDF4 is not installed!")
14 sys.exit(1)
15
16 # geometry setup flowline
17
18 WRIT_FILE = 'flowline_setup.nc' # default name
19 accumrate = 0.5 # accumulation rate in m/a
20
21 #### command line arguments ####
22 try:
23 opts, args = getopt.getopt(sys.argv[1:], "a:r:", ["accumrate=", "resolution="])
24 for opt, arg in opts:
25 if opt in ("-a", "--accumulation"):
26 accumrate = arg
27 if opt in ("-r", "--resolution"): # resolution in km
28 boxWidth = arg
29 except getopt.GetoptError:
30 print('Incorrect command line arguments')
31 sys.exit(2)
32
33
34 # geometry setup
35 boxWidth = float(boxWidth)
36 accumrate = float(accumrate)
37 WRIT_FILE = 'flowline_setup_' + str(int(boxWidth)) + 'km.nc'
38
39 ### CONSTANTS ###
40 secpera = 31556926.
41 yExtent = 2 * boxWidth # in km
42 xExtent = 500 # in km
43 standard_gravity = 9.81 # g
44 B0 = 1.9e8 # ice hardness
45 rho_ice = 910.0 # in kg/m^3
46 rho_ocean = 1028.0 # in kg/m^3
47
48 # create config overwrite file as used in T. Albrecht, M. A. Martin, R. Winkelmann, M. Haseloff, A. Levermann; Parameterization for subgrid-scale motion of ice-shelf calving fronts; (2011), The Cryosphere 5, 35-44, DOI:10.5194/tc-5-35-2011.
49 '''Generates a config file containing flags and parameters
50 for a particular experiment and step.
51
52 This takes care of flags and parameters that *cannot* be set using
53 command-line options. (We try to use command-line options whenever we can.)
54 '''
55
56 filename = "flowline_config.nc"
57 nc = NC(filename, 'w', format="NETCDF3_CLASSIC")
58 var = nc.createVariable("pism_overrides", 'i')
59 attrs = {"ocean.always_grounded": "no",
60 "geometry.update.use_basal_melt_rate": "no",
61 "stress_balance.ssa.compute_surface_gradient_inward": "no",
62 "flow_law.isothermal_Glen.ice_softness": (B0) ** -3,
63 "constants.ice.density": rho_ice,
64 "constants.sea_water.density": rho_ocean,
65 "bootstrapping.defaults.geothermal_flux": 0.0,
66 "stress_balance.ssa.Glen_exponent": 3.0,
67 "constants.standard_gravity": standard_gravity,
68 "ocean.sub_shelf_heat_flux_into_ice": 0.0,
69 "stress_balance.sia.bed_smoother.range": 0.0,
70 }
71
72 for name, value in attrs.items():
73 var.setncattr(name, value)
74 nc.close()
75
76
77 # grid size: # of boxes
78
79 ny = int(np.floor(yExtent / boxWidth / 2) * 2 + 1) # make it an odd number
80 nx = int(np.floor(xExtent / boxWidth / 2) * 2 + 1) # make it an odd number
81
82 # grid size: extent in km's, origin (0,0) to the left of the domain
83
84 x = np.linspace(0, xExtent, nx) * 1000.0
85 y = np.linspace(-yExtent / 2, yExtent / 2, ny) * 1000.0
86
87 nxcenter = int(np.floor(0.5 * nx))
88 nycenter = int(np.floor(0.5 * ny))
89
90 thk = np.zeros((ny, nx))
91 topg = np.zeros((ny, nx))
92 ice_surface_temp = np.zeros((ny, nx))
93 precip = np.zeros((ny, nx))
94
95
96 print("Informations from createSetup_flowline.py:")
97 print("grid size (pixel):")
98 print(ny)
99 print(nx)
100 print("grid size center:")
101 print(nxcenter)
102 print(nycenter)
103
104 print("domain range in meters:")
105 print("x-dir:")
106 print(x[0])
107 print(x[nx - 1])
108 print("y-dir:")
109 print(y[0])
110 print(y[ny - 1])
111
112
113 thickness = 600.0 # initial ice thickness in meters
114 xfront = 700.0 # x-position of fixed calving front in km
115
116 # bc
117 distbc = 50.0 # km
118 igl = int(np.floor(distbc / boxWidth))
119 vel_bc = 300 # m/yr
120 bc_mask = np.zeros((ny, nx))
121 ubar = np.zeros((ny, nx))
122 vbar = np.zeros((ny, nx))
123
124 # "typical constant ice parameter" as defined in the paper and in Van der
125 # Veen's "Fundamentals of Glacier Dynamics", 1999
126 #C = (rho_ice * standard_gravity * (1.0 - rho_ice/rho_ocean) / (4 * B0))**3
127 # H0=thickness
128 # Q0=vel_bc*thickness/secpera
129
130 nxfront = int(xfront / boxWidth)
131
132 # define bedrock geometry topg:
133
134 for i in range(0, nx):
135 for j in range(0, ny):
136 topg[j, i] = -2000.0
137
138 # define constant initial ice-thickness and extent:
139
140 # for i in range(nxcenter-nxfront,nxcenter+nxfront):
141 for i in range(0, nx):
142 for j in range(0, ny):
143 if i <= igl:
144 #thk[j,i] = (4.0 * C / Q0 * (i-igl) + 1 / H0**4)**(-0.25)
145 thk[j, i] = thickness
146
147 # define precipitation field:
148
149 for i in range(0, nx):
150 for j in range(0, ny):
151 precip[j, i] = accumrate / secpera
152
153 # defining dummy temperature:
154
155 for i in range(0, nx):
156 for j in range(0, ny):
157 ice_surface_temp[j, i] = 247.0
158
159 # defining boundary conditions:
160
161 for i in range(0, nx):
162 for j in range(0, ny):
163 if i <= igl:
164 ubar[j, i] = vel_bc / secpera
165 vbar[j, i] = 0.0
166 bc_mask[j, i] = 1.0
167
168
169 ##### define dimensions in NetCDF file #####
170 ncfile = NC(WRIT_FILE, 'w', format='NETCDF3_CLASSIC')
171 xdim = ncfile.createDimension('x', nx)
172 ydim = ncfile.createDimension('y', ny)
173
174 ##### define variables, set attributes, write data #####
175 # format: ['units', 'long_name', 'standard_name', '_FillValue', array]
176
177 vars = {'y': ['m',
178 'y-coordinate in Cartesian system',
179 'projection_y_coordinate',
180 None,
181 y],
182 'x': ['m',
183 'x-coordinate in Cartesian system',
184 'projection_x_coordinate',
185 None,
186 x],
187 'thk': ['m',
188 'floating ice shelf thickness',
189 'land_ice_thickness',
190 1.0,
191 thk],
192 'topg': ['m',
193 'bedrock surface elevation',
194 'bedrock_altitude',
195 -600.0,
196 topg],
197 'ice_surface_temp': ['K',
198 'annual mean air temperature at ice surface',
199 'surface_temperature',
200 248.0,
201 ice_surface_temp],
202 'climatic_mass_balance': ['kg m-2 year-1',
203 'mean annual net ice equivalent accumulation rate',
204 'land_ice_surface_specific_mass_balance_flux',
205 0.2 * 910.0,
206 precip],
207 'bc_mask': ['',
208 'bc_mask',
209 'bc_mask',
210 0.0,
211 bc_mask],
212 'u_ssa_bc': ['m s-1',
213 'X-component of the SSA velocity boundary conditions',
214 'ubar',
215 0.0,
216 ubar],
217 'v_ssa_bc': ['m s-1',
218 'Y-component of the SSA velocity boundary conditions',
219 'vbar',
220 0.0,
221 vbar],
222 }
223
224 for name in list(vars.keys()):
225 [_, _, _, fill_value, data] = vars[name]
226 if name in ['x', 'y']:
227 var = ncfile.createVariable(name, 'f4', (name,))
228 else:
229 var = ncfile.createVariable(name, 'f4', ('y', 'x'), fill_value=fill_value)
230 for each in zip(['units', 'long_name', 'standard_name'], vars[name]):
231 if each[1]:
232 setattr(var, each[0], each[1])
233 var[:] = data
234
235 # finish up
236 ncfile.close()
237 print("NetCDF file ", WRIT_FILE, " created")
238 print('')