tsetup_PXXS.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
---
tsetup_PXXS.py (7182B)
---
1 #!/usr/bin/env python3
2
3 # Copyright (C) 2012, 2014, 2020 Moritz Huetten
4
5 # geometry setup MISMIP3D P75S/P10S-experiment
6
7 from pylab import *
8 import os
9 import sys
10 import getopt
11 import math
12 #from numpy import *
13 import numpy as np
14 try:
15 from netCDF4 import Dataset as NC
16 except:
17 print("netCDF4 is not installed!")
18 sys.exit(1)
19
20 subgl = False
21
22 #### command line arguments ####
23 try:
24 opts, args = getopt.getopt(sys.argv[1:], "si:a:", ["subgl", "inpath=", "amplitude="])
25 for opt, arg in opts:
26 if opt in ("-s", "--subgl"):
27 subgl = True
28 if opt in ("-i", "--inpath"):
29 inpath = arg
30 if opt in ("-a", "--amplitude"):
31 a = arg # pertubation amplitude for tauc
32
33
34 except getopt.GetoptError:
35 print('Incorrect command line arguments')
36 sys.exit(2)
37
38
39 a = float(a) * 100.0
40 a = int(a)
41
42 if a == 75:
43 WRIT_FILE = 'MISMIP3D_P75S_initialSetup.nc'
44 a = 0.75
45 elif a == 10:
46 WRIT_FILE = 'MISMIP3D_P10S_initialSetup.nc'
47 a = 0.1
48 else:
49 WRIT_FILE = 'dummy'
50
51
52 ######## geometry setup (moritz.huetten@pik) #########
53
54
55 ### CONSTANTS ###
56
57 secpera = 31556926.
58 ice_density = 910.0
59
60 yExtent = 2 * 50 # in km
61 xExtent = 2 * 800 # in km
62
63
64 # load data from Stnd-experiment
65
66 try:
67 name = inpath # + '.nc'
68 infile = NC(name, 'r')
69 except Exception:
70 print("file '%s' not found" % name)
71 sys.exit(2)
72 # exit(-1)
73
74
75 x = squeeze(infile.variables["x"][:])
76 nx = len(x)
77
78
79 boxWidth = x[nx - 1] / ((nx - 1) / 2) / 1e3
80
81
82 # grid size: # of boxes
83
84 ny = int(np.floor(yExtent / boxWidth / 2) * 2 + 1) # make it an odd number
85
86
87 # grid size: extent in km's, origin (0,0) in the center of the domain
88
89
90 y = np.linspace(-yExtent / 2, yExtent / 2, ny) * 1000.0
91
92 nxcenter = int(np.floor(0.5 * nx))
93 nycenter = int(np.floor(0.5 * ny))
94
95 thk = np.zeros((ny, nx))
96 topg = np.zeros((ny, nx))
97 ice_surface_temp = np.zeros((ny, nx))
98 precip = np.zeros((ny, nx))
99 tauc = np.zeros((ny, nx))
100
101 print("Informations from createSetup_PXXS.py:")
102 print("grid size:")
103 print(nx)
104 print(ny)
105 print("grid size center:")
106 print(nxcenter)
107 print(nycenter)
108
109 print("domain range in meters:")
110 print("x-dir:")
111 print(x[0])
112 print(x[nx - 1])
113 print("y-dir:")
114 print(y[0])
115 print(y[ny - 1])
116 print("y-boxWidth:")
117 print(y[ny - 1] - y[ny - 2])
118
119
120 # load data from Stnd-result:
121
122 Thk_stnd = squeeze(infile.variables["thk"][:])
123 precip_stnd = squeeze(infile.variables["climatic_mass_balance"][:])
124 Topg_stnd = squeeze(infile.variables["topg"][:])
125 if subgl == True:
126 Gl_mask = squeeze(infile.variables["gl_mask"][:])
127
128 print("number snapshots:")
129 print(len(Thk_stnd[:, 0, 0]))
130 lastslice = len(Thk_stnd[:, 0, 0]) - 1
131
132 thk_stnd = Thk_stnd[lastslice, :, :]
133 precip_stnd = precip_stnd[lastslice, :, :]
134 topg_stnd = Topg_stnd[lastslice, :, :]
135 if subgl == True:
136 gl_mask = Gl_mask[lastslice, :, :]
137
138
139 # load bedrock geometry topg from Stnd-experiment:
140
141 for i in range(0, nx):
142 for j in range(0, ny):
143 topg[j, i] = topg_stnd[0, i]
144
145
146 # load initial ice-thickness:
147
148 for i in range(0, nx):
149 for j in range(0, ny):
150 thk[j, i] = thk_stnd[0, i]
151
152
153 # load precipitation field:
154
155 for i in range(0, nx):
156 for j in range(0, ny):
157 precip[j, i] = precip_stnd[0, 0] / secpera / ice_density
158 print("snow per year in meters")
159 print(precip_stnd[0, 0])
160
161
162 # defining dummy temperature:
163
164 for i in range(0, nx):
165 for j in range(0, ny):
166 ice_surface_temp[j, i] = 268.15
167
168 # create the maximum ice extent mask
169 land_ice_area_fraction_retreat = np.zeros_like(thk)
170 land_ice_area_fraction_retreat[thk > 0] = 1
171 land_ice_area_fraction_retreat[topg > 0] = 1
172
173 # number of grid cells
174 Mx = x.shape[0]
175 middle = (Mx - 1) / 2
176 x1 = x[middle:Mx] / 1000.0 # km
177 dx = x1[1] - x1[0]
178 thk_stnd1 = thk_stnd[1, middle:Mx] # 1D
179 Mask = squeeze(infile.variables["mask"][:])
180 mask = Mask[lastslice, :, :]
181 mask = mask[1, middle:Mx] # 1D
182
183 # find grounding line
184 for i in range(mask.shape[0]):
185 if (thk_stnd1[i] > 0 and mask[i] == 2 and mask[i + 1] == 3):
186 xg = x1[i]
187 if subgl == True:
188 xg_new = xg + dx / 2.0 - (1 - gl_mask[0, i]) * dx + gl_mask[0, i + 1] * dx
189 else:
190 xg_new = xg + dx / 2.0
191
192 print("old grounding line at position:")
193 print(xg, "km")
194
195 print("new grounding line at position:")
196 print(xg_new, "km")
197
198
199 xg = xg_new * 1.0e3
200
201
202 # defining tauc:
203
204 xb = xg
205 yb = 0
206 xc = 150e3
207 yc = 10e3
208 C = 1.0e7
209
210 a = float(a)
211
212 for i in range(nxcenter, nx):
213 for j in range(0, ny):
214 tauc[j, i] = C * (1 - a * math.exp(-(x[i] - xb) ** 2 / (2 * xc ** 2) - (y[j] - yb) ** 2 / (2 * yc ** 2)))
215
216 for i in range(0, nxcenter):
217 for j in range(0, ny):
218 tauc[j, i] = C * (1 - a * math.exp(-(x[i] + xb) ** 2 / (2 * xc ** 2) - (y[j] - yb) ** 2 / (2 * yc ** 2)))
219
220
221 ##### define dimensions in NetCDF file #####
222 ncfile = NC(WRIT_FILE, 'w', format='NETCDF3_CLASSIC')
223 xdim = ncfile.createDimension('x', nx)
224 ydim = ncfile.createDimension('y', ny)
225
226 ##### define variables, set attributes, write data #####
227 # format: ['units', 'long_name', 'standard_name', '_FillValue', array]
228
229 vars = {'y': ['m',
230 'y-coordinate in Cartesian system',
231 'projection_y_coordinate',
232 None,
233 y],
234 'x': ['m',
235 'x-coordinate in Cartesian system',
236 'projection_x_coordinate',
237 None,
238 x],
239 'thk': ['m',
240 'floating ice shelf thickness',
241 'land_ice_thickness',
242 1.0,
243 thk],
244 'topg': ['m',
245 'bedrock surface elevation',
246 'bedrock_altitude',
247 -600.0,
248 topg],
249 'ice_surface_temp': ['K',
250 'annual mean air temperature at ice surface',
251 'surface_temperature',
252 248.0,
253 ice_surface_temp],
254 'climatic_mass_balance': ['kg m-2 year-1',
255 'mean annual net ice equivalent accumulation rate',
256 'land_ice_surface_specific_mass_balance_flux',
257 0.2 * ice_density,
258 precip],
259 'tauc': ['Pa',
260 'yield stress for basal till (plastic or pseudo-plastic model)',
261 'yield_stress_for_basal_till',
262 1e6,
263 tauc],
264 'land_ice_area_fraction_retreat' : ["1",
265 "maximum ice extent mask",
266 "",
267 -1,
268 land_ice_area_fraction_retreat]
269 }
270
271 for name in list(vars.keys()):
272 [_, _, _, fill_value, data] = vars[name]
273 if name in ['x', 'y']:
274 var = ncfile.createVariable(name, 'f4', (name,))
275 else:
276 var = ncfile.createVariable(name, 'f4', ('y', 'x'), fill_value=fill_value)
277 for each in zip(['units', 'long_name', 'standard_name'], vars[name]):
278 if each[1]:
279 setattr(var, each[0], each[1])
280 var[:] = data
281
282 # finish up
283 ncfile.close()
284 print("NetCDF file ", WRIT_FILE, " created")
285 print('')