tsg_create_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
---
tsg_create_flowline.py (2842B)
---
1 #!/usr/bin/env python3
2 #
3 # Copyright (C) 2011, 2014, 2018 Andy Aschwanden
4 #
5 # This file is part of PISM.
6 #
7 # PISM is free software; you can redistribute it and/or modify it under the
8 # terms of the GNU General Public License as published by the Free Software
9 # Foundation; either version 3 of the License, or (at your option) any later
10 # version.
11 #
12 # PISM is distributed in the hope that it will be useful, but WITHOUT ANY
13 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14 # FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
15 # details.
16 #
17 # You should have received a copy of the GNU General Public License
18 # along with PISM; if not, write to the Free Software
19 # Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
20
21 import numpy as np
22
23 try:
24 from netCDF4 import Dataset as CDF
25 except:
26 print("netCDF4 is not installed!")
27 sys.exit(1)
28
29 x, topg, thk = np.loadtxt('sg_35m_flowline.txt', unpack=True)
30
31 output = 'storglaciaren_flowline.nc'
32
33 # Write the data:
34 print("Writing the data to '%s'... " % output)
35 nc = CDF(output, "w")
36 nc.createDimension("x", size=len(x))
37
38 x_var = nc.createVariable("x", 'f', dimensions=("x",))
39 x_var.units = "m"
40 x_var[:] = x
41
42 topg_var = nc.createVariable("topg", 'f', dimensions=("x",))
43 topg_var.units = "m"
44 topg_var.standard_name = "bedrock_altitude"
45 topg_var[:] = topg
46
47 thk_var = nc.createVariable("thk", 'f', dimensions=("x",))
48 thk_var.units = "m"
49 thk_var.standard_name = "land_ice_thickness"
50 thk_var[:] = thk
51
52 usurf_var = nc.createVariable("usurf", 'f', dimensions=("x",))
53 usurf_var.units = "m"
54 usurf_var.standard_name = "surface_altitude"
55 usurf_var[:] = topg + thk
56
57 qgeo = 0.042
58 bheatflx_var = nc.createVariable("bheatflx", 'f', dimensions=("x",))
59 bheatflx_var.units = "W m-2"
60 bheatflx_var[:] = qgeo * np.ones_like(x)
61
62 # generate (somewhat) reasonable acab
63 acab_max = 2.5 # m/a
64 acab_min = -3.0 # m/a
65 np.linspace(acab_max, acab_min, 100)
66 acab = np.ones_like(x)
67 acab[:5] = 0
68 acab[5:105] = np.linspace(acab_max, acab_min, 100)
69 acab[105:] = acab_min
70
71 ice_density = 910.0 # kg m-3
72 acab_var = nc.createVariable("climatic_mass_balance", 'f', dimensions=("x",))
73 acab_var.units = "kg m-2 year-1"
74 acab_var.standard_name = "land_ice_surface_specific_mass_balance_flux"
75 # convert from m/year ice equivalent into [kg m-2 year-1]:
76 acab_var[:] = acab * ice_density
77
78 # Set boundary conditions
79 # ------------------------------------------------------------------------------
80 #
81 # (A) Surface temperature for temperature equation bc
82 Tma = -6.0 # degC, mean annual air temperature at Tarfala
83 zcts = 1400 # m a.s.l.; altitude where CTS is at the surface
84
85
86 artm = np.zeros_like(x)
87 artm[(topg + thk) < zcts] = Tma
88 artm_var = nc.createVariable("ice_surface_temp", 'f', dimensions=("x",))
89 artm_var.units = "deg_C"
90 artm_var[:] = artm
91
92 nc.close()