tsetup_Stnd.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_Stnd.py (5201B)
---
1 #!/usr/bin/env python3
2
3 # Copyright (C) 2012, 2014, 2020 Moritz Huetten
4
5 import sys
6 import getopt
7 import math
8 #from numpy import *
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
17 # geometry setup MISMIP3D Stnd-experiment
18
19
20 WRIT_FILE = 'MISMIP3D_Stnd_initialSetup.nc'
21
22 accumrate = 0.5 # accumulation rate in m/a
23
24 #### command line arguments ####
25 try:
26 opts, args = getopt.getopt(sys.argv[1:], "a:r:", ["accumrate=", "resolution="])
27 for opt, arg in opts:
28 if opt in ("-a", "--accumulation"):
29 accumrate = arg
30 if opt in ("-r", "--resolution"): # resolution in km
31 boxWidth = arg
32 except getopt.GetoptError:
33 print('Incorrect command line arguments')
34 sys.exit(2)
35
36
37 ######## geometry setup (moritz.huetten@pik) #########
38
39
40 boxWidth = float(boxWidth)
41 accumrate = float(accumrate)
42
43 ### CONSTANTS ###
44
45
46 secpera = 31556926.
47 ice_density = 910.0 # [kg m-3]
48
49 yExtent = 2 * boxWidth # in km
50 xExtent = 2 * 800 # in km
51
52
53 # grid size: # of boxes
54
55 ny = int(np.floor(yExtent / boxWidth / 2) * 2 + 1) # make it an odd number
56 nx = int(np.floor(xExtent / boxWidth / 2) * 2 + 1) # make it an odd number
57
58 # grid size: extent in km's, origin (0,0) in the center of the domain
59
60 x = np.linspace(-xExtent / 2, xExtent / 2, nx) * 1000.0
61 y = np.linspace(-yExtent / 2, yExtent / 2, ny) * 1000.0
62
63 nxcenter = int(np.floor(0.5 * nx))
64 nycenter = int(np.floor(0.5 * ny))
65
66 thk = np.zeros((ny, nx))
67 topg = np.zeros((ny, nx))
68 ice_surface_temp = np.zeros((ny, nx))
69 precip = np.zeros((ny, nx))
70
71 # define bed elevation:
72 # linear sloping bed in x-direction with top in the middle, symmetric in y-direction
73
74 # print range(0,int(np.floor(0.5*nx)))
75
76 print("Informations from createSetup_Stnd.py:")
77 print("grid size (pixel):")
78 print(ny)
79 print(nx)
80 print("grid size center:")
81 print(nxcenter)
82 print(nycenter)
83
84 print("domain range in meters:")
85 print("x-dir:")
86 print(x[0])
87 print(x[nx - 1])
88 print("y-dir:")
89 print(y[0])
90 print(y[ny - 1])
91
92
93 # define bedrock geometry topg:
94
95 for i in range(0, nx):
96 for j in range(0, ny):
97 topg[j, i] = -100.0 - abs(x[i]) / 1.0e3
98
99 # define constant initial ice-thickness and extent:
100
101 thickness = 500.0 # initial ice thickness in meters
102 xfront = 700.0 # x-position of fixed calving front in km
103
104 nxfront = int(xfront / boxWidth)
105
106 for i in range(nxcenter - nxfront, nxcenter + nxfront):
107 for j in range(0, ny):
108 thk[j, i] = thickness
109
110
111 # define precipitation field:
112
113 for i in range(0, nx):
114 for j in range(0, ny):
115 precip[j, i] = accumrate / secpera
116
117
118 # defining dummy temperature:
119
120 for i in range(0, nx):
121 for j in range(0, ny):
122 ice_surface_temp[j, i] = 268.15
123
124 # create the maximum ice extent mask
125 land_ice_area_fraction_retreat = np.zeros_like(thk)
126 land_ice_area_fraction_retreat[thk > 0] = 1
127 land_ice_area_fraction_retreat[topg > 0] = 1
128
129 ##### define dimensions in NetCDF file #####
130 ncfile = NC(WRIT_FILE, 'w', format='NETCDF3_CLASSIC')
131 xdim = ncfile.createDimension('x', nx)
132 ydim = ncfile.createDimension('y', ny)
133
134 ##### define variables, set attributes, write data #####
135 # format: ['units', 'long_name', 'standard_name', '_FillValue', array]
136
137 vars = {'y': ['m',
138 'y-coordinate in Cartesian system',
139 'projection_y_coordinate',
140 None,
141 y],
142 'x': ['m',
143 'x-coordinate in Cartesian system',
144 'projection_x_coordinate',
145 None,
146 x],
147 'thk': ['m',
148 'floating ice shelf thickness',
149 'land_ice_thickness',
150 1.0,
151 thk],
152 'topg': ['m',
153 'bedrock surface elevation',
154 'bedrock_altitude',
155 -600.0,
156 topg],
157 'ice_surface_temp': ['K',
158 'annual mean air temperature at ice surface',
159 'surface_temperature',
160 248.0,
161 ice_surface_temp],
162 'climatic_mass_balance': ['kg m-2 year-1',
163 'mean annual net ice equivalent accumulation rate',
164 'land_ice_surface_specific_mass_balance_flux',
165 0.2 * ice_density,
166 precip],
167 'land_ice_area_fraction_retreat' : ["1",
168 "maximum ice extent mask",
169 "",
170 -1,
171 land_ice_area_fraction_retreat]
172 }
173
174 for name in list(vars.keys()):
175 [_, _, _, fill_value, data] = vars[name]
176 if name in ['x', 'y']:
177 var = ncfile.createVariable(name, 'f4', (name,))
178 else:
179 var = ncfile.createVariable(name, 'f4', ('y', 'x'), fill_value=fill_value)
180 for each in zip(['units', 'long_name', 'standard_name'], vars[name]):
181 if each[1]:
182 setattr(var, each[0], each[1])
183 var[:] = data
184
185 # finish up
186 ncfile.close()
187 print("NetCDF file ", WRIT_FILE, " created")
188 print('')