tElevation.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tElevation.cc (9685B)
---
1 // Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018 Andy Aschwanden and Constantine Khroulev
2 //
3 // This file is part of PISM.
4 //
5 // PISM is free software; you can redistribute it and/or modify it under the
6 // terms of the GNU General Public License as published by the Free Software
7 // Foundation; either version 3 of the License, or (at your option) any later
8 // version.
9 //
10 // PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 // details.
14 //
15 // You should have received a copy of the GNU General Public License
16 // along with PISM; if not, write to the Free Software
17 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18
19 #include "Elevation.hh"
20
21 #include "pism/util/iceModelVec.hh"
22 #include "pism/util/io/File.hh"
23 #include "pism/util/Vars.hh"
24 #include "pism/util/IceGrid.hh"
25 #include "pism/util/ConfigInterface.hh"
26 #include "pism/util/error_handling.hh"
27 #include "pism/util/pism_options.hh"
28 #include "pism/util/io/io_helpers.hh"
29 #include "pism/util/MaxTimestep.hh"
30 #include "pism/util/pism_utilities.hh"
31 #include "pism/geometry/Geometry.hh"
32
33 namespace pism {
34 namespace surface {
35
36 ///// Elevation-dependent temperature and surface mass balance.
37 Elevation::Elevation(IceGrid::ConstPtr grid, std::shared_ptr<atmosphere::AtmosphereModel> input)
38 : SurfaceModel(grid) {
39 (void) input;
40
41 m_temperature = allocate_temperature(grid);
42 m_mass_flux = allocate_mass_flux(grid);
43 }
44
45 void Elevation::init_impl(const Geometry &geometry) {
46 (void) geometry;
47
48 bool limits_set = false;
49
50 m_log->message(2,
51 "* Initializing the constant-in-time surface processes model Elevation. Setting...\n");
52
53 units::Converter K(m_sys, "Celsius", "Kelvin");
54 // options
55 {
56 // ice surface temperature
57 {
58 options::RealList IST("-ice_surface_temp", "ice surface temperature parameterization",
59 {-5.0, 0.0, 1325.0, 1350.0});
60
61 if (IST->size() != 4) {
62 throw RuntimeError(PISM_ERROR_LOCATION, "option -ice_surface_temp requires an argument"
63 " (comma-separated list of 4 numbers)");
64 }
65 m_T_min = K(IST[0]);
66 m_T_max = K(IST[1]);
67 m_z_T_min = IST[2];
68 m_z_T_max = IST[3];
69 }
70
71 // climatic mass balance
72 units::Converter meter_per_second(m_sys, "m year-1", "m second-1");
73 {
74 options::RealList CMB("-climatic_mass_balance",
75 "climatic mass balance parameterization",
76 {-3.0, 4.0, 1100.0, 1450.0, 1700.0});
77 if (CMB->size() != 5) {
78 throw RuntimeError(PISM_ERROR_LOCATION, "-climatic_mass_balance requires an argument"
79 " (comma-separated list of 5 numbers)");
80 }
81 m_M_min = meter_per_second(CMB[0]);
82 m_M_max = meter_per_second(CMB[1]);
83 m_z_M_min = CMB[2];
84 m_z_ELA = CMB[3];
85 m_z_M_max = CMB[4];
86 }
87
88 // limits of the climatic mass balance
89 {
90 options::RealList limits("-climatic_mass_balance_limits",
91 "lower and upper limits of the climatic mass balance", {});
92 limits_set = limits.is_set();
93 if (limits.is_set()) {
94 if (limits->size() != 2) {
95 throw RuntimeError(PISM_ERROR_LOCATION, "-climatic_mass_balance_limits requires an argument"
96 " (a comma-separated list of 2 numbers)");
97 }
98 m_M_limit_min = meter_per_second(limits[0]);
99 m_M_limit_max = meter_per_second(limits[1]);
100 } else {
101 m_M_limit_min = m_M_min;
102 m_M_limit_max = m_M_max;
103 }
104 }
105 }
106
107 units::Converter meter_per_year(m_sys, "m second-1", "m year-1");
108 m_log->message(3,
109 " temperature at %.0f m a.s.l. = %.2f deg C\n"
110 " temperature at %.0f m a.s.l. = %.2f deg C\n"
111 " mass balance below %.0f m a.s.l. = %.2f m year-1\n"
112 " mass balance at %.0f m a.s.l. = %.2f m year-1\n"
113 " mass balance at %.0f m a.s.l. = %.2f m year-1\n"
114 " mass balance above %.0f m a.s.l. = %.2f m year-1\n"
115 " equilibrium line altitude z_ELA = %.2f m a.s.l.\n",
116 m_z_T_min, m_T_min,
117 m_z_T_max, m_T_max,
118 m_z_M_min, meter_per_year(m_M_limit_min),
119 m_z_M_min, m_M_min,
120 m_z_M_max, meter_per_year(m_M_max),
121 m_z_M_max, meter_per_year(m_M_limit_max),
122 m_z_ELA);
123
124 // parameterizing the ice surface temperature 'ice_surface_temp'
125 m_log->message(2, " - parameterizing the ice surface temperature 'ice_surface_temp' ... \n");
126 m_log->message(2,
127 " ice temperature at the ice surface (T = ice_surface_temp) is piecewise-linear function\n"
128 " of surface altitude (usurf):\n"
129 " / %2.2f K for usurf < %.f m\n"
130 " T = | %5.2f K + %5.3f * (usurf - %.f m) for %.f m < usurf < %.f m\n"
131 " \\ %5.2f K for %.f m < usurf\n",
132 m_T_min, m_z_T_min,
133 m_T_min, (m_T_max-m_T_min)/(m_z_T_max-m_z_T_min), m_z_T_min, m_z_T_min, m_z_T_max,
134 m_T_max, m_z_T_max);
135
136 // parameterizing the ice surface mass balance 'climatic_mass_balance'
137 m_log->message(2,
138
139 " - parameterizing the ice surface mass balance 'climatic_mass_balance' ... \n");
140
141 if (limits_set) {
142 m_log->message(2,
143 " - option '-climatic_mass_balance_limits' seen, limiting upper and lower bounds ... \n");
144 }
145
146 m_log->message(2,
147 " surface mass balance (M = climatic_mass_balance) is piecewise-linear function\n"
148 " of surface altitue (usurf):\n"
149 " / %5.2f m year-1 for usurf < %3.f m\n"
150 " M = | %5.3f 1/a * (usurf-%.0f m) for %3.f m < usurf < %3.f m\n"
151 " \\ %5.3f 1/a * (usurf-%.0f m) for %3.f m < usurf < %3.f m\n"
152 " \\ %5.2f m year-1 for %3.f m < usurf\n",
153 meter_per_year(m_M_limit_min), m_z_M_min,
154 meter_per_year(-m_M_min)/(m_z_ELA - m_z_M_min), m_z_ELA, m_z_M_min, m_z_ELA,
155 meter_per_year(m_M_max)/(m_z_M_max - m_z_ELA), m_z_ELA, m_z_ELA, m_z_M_max,
156 meter_per_year(m_M_limit_max), m_z_M_max);
157 }
158
159 MaxTimestep Elevation::max_timestep_impl(double t) const {
160 (void) t;
161 return MaxTimestep("surface 'elevation'");
162 }
163
164 void Elevation::update_impl(const Geometry &geometry, double t, double dt) {
165 (void) t;
166 (void) dt;
167
168 compute_mass_flux(geometry.ice_surface_elevation, *m_mass_flux);
169 compute_temperature(geometry.ice_surface_elevation, *m_temperature);
170
171 dummy_accumulation(*m_mass_flux, *m_accumulation);
172 dummy_melt(*m_mass_flux, *m_melt);
173 dummy_runoff(*m_mass_flux, *m_runoff);
174
175 }
176
177 const IceModelVec2S &Elevation::mass_flux_impl() const {
178 return *m_mass_flux;
179 }
180
181 const IceModelVec2S &Elevation::temperature_impl() const {
182 return *m_temperature;
183 }
184
185 const IceModelVec2S &Elevation::accumulation_impl() const {
186 return *m_accumulation;
187 }
188
189 const IceModelVec2S &Elevation::melt_impl() const {
190 return *m_melt;
191 }
192
193 const IceModelVec2S &Elevation::runoff_impl() const {
194 return *m_runoff;
195 }
196
197 void Elevation::compute_mass_flux(const IceModelVec2S &surface, IceModelVec2S &result) const {
198 double dabdz = -m_M_min/(m_z_ELA - m_z_M_min);
199 double dacdz = m_M_max/(m_z_M_max - m_z_ELA);
200
201 IceModelVec::AccessList list{&result, &surface};
202
203 ParallelSection loop(m_grid->com);
204 try {
205 for (Points p(*m_grid); p; p.next()) {
206 const int i = p.i(), j = p.j();
207
208 double z = surface(i, j);
209
210 if (z < m_z_M_min) {
211 result(i, j) = m_M_limit_min;
212 }
213 else if ((z >= m_z_M_min) and (z < m_z_ELA)) {
214 result(i, j) = dabdz * (z - m_z_ELA);
215 }
216 else if ((z >= m_z_ELA) and (z <= m_z_M_max)) {
217 result(i, j) = dacdz * (z - m_z_ELA);
218 }
219 else if (z > m_z_M_max) {
220 result(i, j) = m_M_limit_max;
221 }
222 else {
223 throw RuntimeError(PISM_ERROR_LOCATION,
224 "Elevation::compute_mass_flux: HOW DID I GET HERE?");
225 }
226 }
227 } catch (...) {
228 loop.failed();
229 }
230 loop.check();
231
232 // convert from m second-1 ice equivalent to kg m-2 s-1:
233 result.scale(m_config->get_number("constants.ice.density"));
234 }
235
236 void Elevation::compute_temperature(const IceModelVec2S &surface, IceModelVec2S &result) const {
237
238 IceModelVec::AccessList list{&result, &surface};
239
240 double dTdz = (m_T_max - m_T_min)/(m_z_T_max - m_z_T_min);
241 ParallelSection loop(m_grid->com);
242 try {
243 for (Points p(*m_grid); p; p.next()) {
244 const int i = p.i(), j = p.j();
245
246 double z = surface(i, j);
247
248 if (z <= m_z_T_min) {
249 result(i, j) = m_T_min;
250 }
251 else if ((z > m_z_T_min) and (z < m_z_T_max)) {
252 result(i, j) = m_T_min + dTdz * (z - m_z_T_min);
253 }
254 else if (z >= m_z_T_max) {
255 result(i, j) = m_T_max;
256 }
257 else {
258 throw RuntimeError(PISM_ERROR_LOCATION,
259 "Elevation::compute_temperature: HOW DID I GET HERE?");
260 }
261 }
262 } catch (...) {
263 loop.failed();
264 }
265 loop.check();
266 }
267
268 } // end of namespace surface
269 } // end of namespace pism