tIceRegionalModel.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
---
tIceRegionalModel.cc (14879B)
---
1 /* Copyright (C) 2015, 2016, 2017, 2018, 2019 PISM Authors
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
20 #include "IceRegionalModel.hh"
21 #include "pism/util/Vars.hh"
22 #include "pism/util/EnthalpyConverter.hh"
23 #include "pism/stressbalance/ShallowStressBalance.hh"
24 #include "SSAFD_Regional.hh"
25 #include "pism/stressbalance/SSB_Modifier.hh"
26 #include "SIAFD_Regional.hh"
27 #include "pism/stressbalance/StressBalance.hh"
28 #include "pism/basalstrength/ConstantYieldStress.hh"
29 #include "RegionalYieldStress.hh"
30 #include "pism/util/io/File.hh"
31 #include "pism/coupler/OceanModel.hh"
32 #include "pism/coupler/SurfaceModel.hh"
33 #include "EnthalpyModel_Regional.hh"
34 #include "pism/energy/CHSystem.hh"
35 #include "pism/energy/BedThermalUnit.hh"
36 #include "pism/energy/utilities.hh"
37 #include "pism/util/iceModelVec2T.hh"
38 #include "pism/hydrology/Hydrology.hh"
39
40 namespace pism {
41
42 IceRegionalModel::IceRegionalModel(IceGrid::Ptr g, Context::Ptr c)
43 : IceModel(g, c) {
44 // empty
45
46 if (m_config->get_flag("energy.ch_warming.enabled")) {
47 m_ch_warming_flux.reset(new IceModelVec3(m_grid, "ch_warming_flux", WITHOUT_GHOSTS));
48 }
49 }
50
51
52 void IceRegionalModel::allocate_storage() {
53
54 IceModel::allocate_storage();
55
56 m_log->message(2,
57 " creating IceRegionalModel vecs ...\n");
58
59 // stencil width of 2 needed by SIAFD_Regional::compute_surface_gradient()
60 m_no_model_mask.create(m_grid, "no_model_mask", WITH_GHOSTS, 2);
61 m_no_model_mask.set_attrs("model_state",
62 "mask: zeros (modeling domain) and ones"
63 " (no-model buffer near grid edges)",
64 "", "", "", 0); // no units and no standard name
65 m_no_model_mask.metadata().set_numbers("flag_values", {0, 1});
66 m_no_model_mask.metadata().set_string("flag_meanings", "normal special_treatment");
67 m_no_model_mask.set_time_independent(true);
68 m_no_model_mask.metadata().set_output_type(PISM_INT);
69 m_no_model_mask.set(0);
70
71 // stencil width of 2 needed for differentiation because GHOSTS=1
72 m_usurf_stored.create(m_grid, "usurfstore", WITH_GHOSTS, 2);
73 m_usurf_stored.set_attrs("model_state",
74 "saved surface elevation for use to keep surface gradient constant"
75 " in no_model strip",
76 "m", "m",
77 "", 0); // no standard name
78
79 // stencil width of 1 needed for differentiation
80 m_thk_stored.create(m_grid, "thkstore", WITH_GHOSTS, 1);
81 m_thk_stored.set_attrs("model_state",
82 "saved ice thickness for use to keep driving stress constant"
83 " in no_model strip",
84 "m", "m",
85 "", 0); // no standard name
86
87 m_model_state.insert(&m_thk_stored);
88 m_model_state.insert(&m_usurf_stored);
89 m_model_state.insert(&m_no_model_mask);
90 }
91
92 void IceRegionalModel::model_state_setup() {
93
94 // initialize the model state (including special fields)
95 IceModel::model_state_setup();
96
97 InputOptions input = process_input_options(m_ctx->com(), m_config);
98
99 // Initialize stored ice thickness and surface elevation. This goes here and not in
100 // bootstrap_2d because bed topography is not initialized at the time bootstrap_2d is
101 // called.
102 if (input.type == INIT_BOOTSTRAP) {
103 if (m_config->get_flag("regional.zero_gradient")) {
104 m_usurf_stored.set(0.0);
105 m_thk_stored.set(0.0);
106 } else {
107 m_usurf_stored.copy_from(m_geometry.ice_surface_elevation);
108 m_thk_stored.copy_from(m_geometry.ice_thickness);
109 }
110 }
111
112 m_geometry_evolution->set_no_model_mask(m_no_model_mask);
113
114 if (m_ch_system) {
115 const bool use_input_file = input.type == INIT_BOOTSTRAP or input.type == INIT_RESTART;
116
117 std::unique_ptr<File> input_file;
118
119 if (use_input_file) {
120 input_file.reset(new File(m_grid->com, input.filename, PISM_GUESS, PISM_READONLY));
121 }
122
123 switch (input.type) {
124 case INIT_RESTART:
125 {
126 m_ch_system->restart(*input_file, input.record);
127 break;
128 }
129 case INIT_BOOTSTRAP:
130 {
131
132 m_ch_system->bootstrap(*input_file,
133 m_geometry.ice_thickness,
134 m_surface->temperature(),
135 m_surface->mass_flux(),
136 m_btu->flux_through_top_surface());
137 break;
138 }
139 case INIT_OTHER:
140 default:
141 {
142 m_basal_melt_rate.set(m_config->get_number("bootstrapping.defaults.bmelt"));
143
144 m_ch_system->initialize(m_basal_melt_rate,
145 m_geometry.ice_thickness,
146 m_surface->temperature(),
147 m_surface->mass_flux(),
148 m_btu->flux_through_top_surface());
149
150 }
151 }
152 }
153 }
154
155 void IceRegionalModel::allocate_geometry_evolution() {
156 if (m_geometry_evolution) {
157 return;
158 }
159
160 m_log->message(2, "# Allocating the geometry evolution model (regional version)...\n");
161
162 m_geometry_evolution.reset(new RegionalGeometryEvolution(m_grid));
163
164 m_submodels["geometry_evolution"] = m_geometry_evolution.get();
165 }
166
167 void IceRegionalModel::allocate_energy_model() {
168
169 if (m_energy_model != NULL) {
170 return;
171 }
172
173 m_log->message(2, "# Allocating an energy balance model...\n");
174
175 if (m_config->get_flag("energy.enabled")) {
176 if (m_config->get_flag("energy.temperature_based")) {
177 throw RuntimeError(PISM_ERROR_LOCATION,
178 "pismr -regional does not support the '-energy cold' mode.");
179 } else {
180 m_energy_model = new energy::EnthalpyModel_Regional(m_grid, m_stress_balance.get());
181 }
182 } else {
183 m_energy_model = new energy::DummyEnergyModel(m_grid, m_stress_balance.get());
184 }
185
186 m_submodels["energy balance model"] = m_energy_model;
187
188 if (m_config->get_flag("energy.ch_warming.enabled") and
189 not m_ch_system) {
190
191 m_log->message(2, "# Allocating the cryo-hydrologic warming model...\n");
192
193 m_ch_system.reset(new energy::CHSystem(m_grid, m_stress_balance.get()));
194 m_submodels["cryo-hydrologic warming"] = m_ch_system.get();
195 }
196 }
197
198 void IceRegionalModel::allocate_stressbalance() {
199
200 if (m_stress_balance) {
201 return;
202 }
203
204 bool regional = true;
205 m_stress_balance = stressbalance::create(m_config->get_string("stress_balance.model"),
206 m_grid, regional);
207
208 m_submodels["stress balance"] = m_stress_balance.get();
209 }
210
211
212 void IceRegionalModel::allocate_basal_yield_stress() {
213
214 if (m_basal_yield_stress_model) {
215 return;
216 }
217
218 IceModel::allocate_basal_yield_stress();
219
220 if (m_basal_yield_stress_model) {
221 // IceModel allocated a basal yield stress model. This means that we are using a
222 // stress balance model that uses it and we need to add regional modifications.
223 m_basal_yield_stress_model.reset(new RegionalYieldStress(m_basal_yield_stress_model));
224
225 m_submodels["basal yield stress"] = m_basal_yield_stress_model.get();
226 }
227 }
228
229 //! Bootstrap a "regional" model.
230 /*!
231 * Need to initialize all the variables managed by IceModel, as well as
232 * - no_model_mask
233 * - usurf_stored
234 * - thk_stored
235 */
236 void IceRegionalModel::bootstrap_2d(const File &input_file) {
237
238 IceModel::bootstrap_2d(input_file);
239
240 // no_model_mask
241 {
242 if (input_file.find_variable(m_no_model_mask.metadata().get_name())) {
243 m_no_model_mask.regrid(input_file, CRITICAL);
244 } else {
245 // set using the no_model_strip parameter
246 double strip_width = m_config->get_number("regional.no_model_strip", "meters");
247 set_no_model_strip(*m_grid, strip_width, m_no_model_mask);
248 }
249
250 // m_no_model_mask was added to m_model_state, so
251 // IceModel::regrid() will take care of regridding it, if necessary.
252 }
253
254 if (m_config->get_flag("stress_balance.ssa.dirichlet_bc")) {
255 IceModelVec::AccessList list{&m_no_model_mask, &m_ssa_dirichlet_bc_mask};
256
257 for (Points p(*m_grid); p; p.next()) {
258 const int i = p.i(), j = p.j();
259
260 if (m_no_model_mask(i, j) > 0.5) {
261 m_ssa_dirichlet_bc_mask(i, j) = 1;
262 }
263 }
264 }
265 }
266
267 stressbalance::Inputs IceRegionalModel::stress_balance_inputs() {
268 stressbalance::Inputs result = IceModel::stress_balance_inputs();
269
270 result.no_model_mask = &m_no_model_mask;
271 result.no_model_ice_thickness = &m_thk_stored;
272 result.no_model_surface_elevation = &m_usurf_stored;
273
274 return result;
275 }
276
277 energy::Inputs IceRegionalModel::energy_model_inputs() {
278 energy::Inputs result = IceModel::energy_model_inputs();
279
280 result.no_model_mask = &m_no_model_mask;
281
282 return result;
283 }
284
285 void IceRegionalModel::energy_step() {
286
287 if (m_ch_system) {
288 bedrock_thermal_model_step();
289
290 energy::Inputs inputs = energy_model_inputs();
291 const IceModelVec3 *strain_heating = inputs.volumetric_heating_rate;
292 inputs.volumetric_heating_rate = m_ch_warming_flux.get();
293
294 energy::cryo_hydrologic_warming_flux(m_config->get_number("constants.ice.thermal_conductivity"),
295 m_config->get_number("energy.ch_warming.average_channel_spacing"),
296 m_geometry.ice_thickness,
297 m_energy_model->enthalpy(),
298 m_ch_system->enthalpy(),
299 *m_ch_warming_flux);
300
301 // Convert to the loss of energy by the CH system:
302 m_ch_warming_flux->scale(-1.0);
303
304 m_ch_system->update(t_TempAge, dt_TempAge, inputs);
305
306 // Add CH warming flux to the strain heating term:
307 m_ch_warming_flux->scale(-1.0);
308 m_ch_warming_flux->add(1.0, *strain_heating);
309
310 m_energy_model->update(t_TempAge, dt_TempAge, inputs);
311
312 m_stdout_flags = m_energy_model->stdout_flags() + m_stdout_flags;
313 } else {
314 IceModel::energy_step();
315 }
316 }
317
318 YieldStressInputs IceRegionalModel::yield_stress_inputs() {
319 YieldStressInputs result = IceModel::yield_stress_inputs();
320
321 result.no_model_mask = &m_no_model_mask;
322
323 return result;
324 }
325
326 const energy::CHSystem* IceRegionalModel::cryo_hydrologic_system() const {
327 return m_ch_system.get();
328 }
329
330 /*! @brief Report temperature of the cryo-hydrologic system */
331 class CHTemperature : public Diag<IceRegionalModel>
332 {
333 public:
334 CHTemperature(const IceRegionalModel *m)
335 : Diag<IceRegionalModel>(m) {
336
337 m_vars = {SpatialVariableMetadata(m_sys, "ch_temp", m_grid->z())};
338
339 set_attrs("temperature of the cryo-hydrologic system", "",
340 "Kelvin", "Kelvin", 0);
341 }
342
343 protected:
344 IceModelVec::Ptr compute_impl() const {
345
346 IceModelVec3::Ptr result(new IceModelVec3(m_grid, "ch_temp", WITHOUT_GHOSTS));
347
348 energy::compute_temperature(model->cryo_hydrologic_system()->enthalpy(),
349 model->geometry().ice_thickness,
350 *result);
351 result->metadata(0) = m_vars[0];
352
353 return result;
354 }
355 };
356
357 /*! @brief Report liquid water fraction in the cryo-hydrologic system */
358 class CHLiquidWaterFraction : public Diag<IceRegionalModel>
359 {
360 public:
361 CHLiquidWaterFraction(const IceRegionalModel *m)
362 : Diag<IceRegionalModel>(m) {
363
364 m_vars = {SpatialVariableMetadata(m_sys, "ch_liqfrac", m_grid->z())};
365
366 set_attrs("liquid water fraction in the cryo-hydrologic system", "",
367 "1", "1", 0);
368 }
369
370 protected:
371 IceModelVec::Ptr compute_impl() const {
372
373 IceModelVec3::Ptr result(new IceModelVec3(m_grid, "ch_liqfrac", WITHOUT_GHOSTS));
374
375 energy::compute_liquid_water_fraction(model->cryo_hydrologic_system()->enthalpy(),
376 model->geometry().ice_thickness,
377 *result);
378 result->metadata(0) = m_vars[0];
379 return result;
380 }
381 };
382
383
384 /*! @brief Report rate of cryo-hydrologic warming */
385 class CHHeatFlux : public Diag<IceRegionalModel>
386 {
387 public:
388 CHHeatFlux(const IceRegionalModel *m)
389 : Diag<IceRegionalModel>(m) {
390
391 m_vars = {SpatialVariableMetadata(m_sys, "ch_heat_flux", m_grid->z())};
392
393 set_attrs("rate of cryo-hydrologic warming", "",
394 "W m-3", "W m-3", 0);
395 }
396
397 protected:
398 IceModelVec::Ptr compute_impl() const {
399
400 IceModelVec3::Ptr result(new IceModelVec3(m_grid, "ch_heat_flux", WITHOUT_GHOSTS));
401 result->metadata(0) = m_vars[0];
402
403 energy::cryo_hydrologic_warming_flux(m_config->get_number("constants.ice.thermal_conductivity"),
404 m_config->get_number("energy.ch_warming.average_channel_spacing"),
405 model->geometry().ice_thickness,
406 model->energy_balance_model()->enthalpy(),
407 model->cryo_hydrologic_system()->enthalpy(),
408 *result);
409 return result;
410 }
411 };
412
413 void IceRegionalModel::init_diagnostics() {
414 IceModel::init_diagnostics();
415
416 if (m_ch_system) {
417 m_diagnostics["ch_temp"] = Diagnostic::Ptr(new CHTemperature(this));
418 m_diagnostics["ch_liqfrac"] = Diagnostic::Ptr(new CHLiquidWaterFraction(this));
419 m_diagnostics["ch_heat_flux"] = Diagnostic::Ptr(new CHHeatFlux(this));
420 }
421 }
422
423 void IceRegionalModel::hydrology_step() {
424 hydrology::Inputs inputs;
425
426 IceModelVec2S &sliding_speed = m_work2d[0];
427 sliding_speed.set_to_magnitude(m_stress_balance->advective_velocity());
428
429 inputs.no_model_mask = &m_no_model_mask;
430 inputs.geometry = &m_geometry;
431 inputs.surface_input_rate = nullptr;
432 inputs.basal_melt_rate = &m_basal_melt_rate;
433 inputs.ice_sliding_speed = &sliding_speed;
434
435 if (m_surface_input_for_hydrology) {
436 m_surface_input_for_hydrology->update(m_time->current(), m_dt);
437 m_surface_input_for_hydrology->average(m_time->current(), m_dt);
438 inputs.surface_input_rate = m_surface_input_for_hydrology.get();
439 } else if (m_config->get_flag("hydrology.surface_input_from_runoff")) {
440 // convert [kg m-2] to [kg m-2 s-1]
441 IceModelVec2S &surface_input_rate = m_work2d[1];
442 surface_input_rate.copy_from(m_surface->runoff());
443 surface_input_rate.scale(1.0 / m_dt);
444 inputs.surface_input_rate = &surface_input_rate;
445 }
446
447 m_subglacial_hydrology->update(m_time->current(), m_dt, inputs);
448 }
449
450
451 } // end of namespace pism