tIceModel.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
---
tIceModel.cc (37569B)
---
1 // Copyright (C) 2004-2020 Jed Brown, Ed Bueler 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 <cmath>
20 #include <cstring>
21 #include <algorithm>
22 #include <petscsys.h>
23
24 #include "pism/pism_config.hh"
25
26 #include "IceModel.hh"
27
28 #include "pism/basalstrength/YieldStress.hh"
29 #include "pism/basalstrength/basal_resistance.hh"
30 #include "pism/frontretreat/util/IcebergRemover.hh"
31 #include "pism/frontretreat/calving/CalvingAtThickness.hh"
32 #include "pism/frontretreat/calving/EigenCalving.hh"
33 #include "pism/frontretreat/calving/FloatKill.hh"
34 #include "pism/frontretreat/calving/HayhurstCalving.hh"
35 #include "pism/frontretreat/calving/vonMisesCalving.hh"
36 #include "pism/energy/BedThermalUnit.hh"
37 #include "pism/hydrology/Hydrology.hh"
38 #include "pism/stressbalance/StressBalance.hh"
39 #include "pism/util/IceGrid.hh"
40 #include "pism/util/Mask.hh"
41 #include "pism/util/ConfigInterface.hh"
42 #include "pism/util/Diagnostic.hh"
43 #include "pism/util/error_handling.hh"
44 #include "pism/util/pism_options.hh"
45 #include "pism/coupler/SeaLevel.hh"
46 #include "pism/coupler/OceanModel.hh"
47 #include "pism/coupler/SurfaceModel.hh"
48 #include "pism/earth/BedDef.hh"
49 #include "pism/util/EnthalpyConverter.hh"
50 #include "pism/util/pism_signal.h"
51 #include "pism/util/Vars.hh"
52 #include "pism/util/Profiling.hh"
53 #include "pism/util/pism_utilities.hh"
54 #include "pism/age/AgeModel.hh"
55 #include "pism/energy/EnergyModel.hh"
56 #include "pism/util/io/File.hh"
57 #include "pism/util/iceModelVec2T.hh"
58 #include "pism/fracturedensity/FractureDensity.hh"
59 #include "pism/coupler/util/options.hh" // ForcingOptions
60
61 namespace pism {
62
63 IceModel::IceModel(IceGrid::Ptr g, Context::Ptr context)
64 : m_grid(g),
65 m_config(context->config()),
66 m_ctx(context),
67 m_sys(context->unit_system()),
68 m_log(context->log()),
69 m_time(context->time()),
70 m_output_global_attributes("PISM_GLOBAL", m_sys),
71 m_run_stats("run_stats", m_sys),
72 m_geometry(m_grid),
73 m_new_bed_elevation(true),
74 m_thickness_change(g),
75 m_ts_times(new std::vector<double>()),
76 m_extra_bounds("time_bounds", m_config->get_string("time.dimension_name"), m_sys),
77 m_timestamp("timestamp", m_config->get_string("time.dimension_name"), m_sys) {
78
79 // time-independent info
80 {
81 m_run_stats.set_string("source", std::string("PISM ") + pism::revision);
82 m_run_stats.set_string("long_name", "Run statistics");
83 }
84
85 m_extra_bounds.set_string("units", m_time->units_string());
86
87 m_timestamp.set_string("units", "hours");
88 m_timestamp.set_string("long_name", "wall-clock time since the beginning of the run");
89
90 pism_signal = 0;
91 signal(SIGTERM, pism_signal_handler);
92 signal(SIGUSR1, pism_signal_handler);
93 signal(SIGUSR2, pism_signal_handler);
94
95 m_surface = nullptr;
96 m_ocean = nullptr;
97 m_beddef = nullptr;
98
99 m_btu = nullptr;
100 m_energy_model = nullptr;
101
102 m_output_global_attributes.set_string("Conventions", "CF-1.6");
103 m_output_global_attributes.set_string("source", pism::version());
104
105 // Do not save snapshots by default:
106 m_save_snapshots = false;
107 // Do not save time-series by default:
108 m_save_extra = false;
109
110 m_fracture = nullptr;
111
112 reset_counters();
113
114 // allocate temporary storage
115 {
116 const unsigned int WIDE_STENCIL = m_config->get_number("grid.max_stencil_width");
117
118 // various internal quantities
119 // 2d work vectors
120 for (int j = 0; j < m_n_work2d; j++) {
121 char namestr[30];
122 snprintf(namestr, sizeof(namestr), "work_vector_%d", j);
123 m_work2d[j].create(m_grid, namestr, WITH_GHOSTS, WIDE_STENCIL);
124 }
125 }
126
127 auto surface_input_file = m_config->get_string("hydrology.surface_input.file");
128 if (not surface_input_file.empty()) {
129 ForcingOptions surface_input(*m_ctx, "hydrology.surface_input");
130 int buffer_size = m_config->get_number("input.forcing.buffer_size");
131 int evaluations_per_year = m_config->get_number("input.forcing.evaluations_per_year");
132
133 File file(m_grid->com, surface_input.filename, PISM_NETCDF3, PISM_READONLY);
134
135 m_surface_input_for_hydrology = IceModelVec2T::ForcingField(m_grid,
136 file,
137 "water_input_rate",
138 "", // no standard name
139 buffer_size,
140 evaluations_per_year,
141 surface_input.period);
142 m_surface_input_for_hydrology->set_attrs("diagnostic",
143 "water input rate for the subglacial hydrology model",
144 "kg m-2 s-1", "kg m-2 year-1", "", 0);
145 m_surface_input_for_hydrology->metadata().set_number("valid_min", 0.0);
146 }
147 }
148
149 double IceModel::dt() const {
150 return m_dt;
151 }
152
153 void IceModel::reset_counters() {
154 dt_TempAge = 0.0;
155 m_dt = 0.0;
156 m_skip_countdown = 0;
157
158 m_timestep_hit_multiples_last_time = m_time->current();
159 }
160
161
162 IceModel::~IceModel() {
163
164 delete m_beddef;
165
166 delete m_btu;
167 delete m_energy_model;
168 }
169
170
171 //! Allocate all IceModelVecs defined in IceModel.
172 /*!
173 This procedure allocates the memory used to store model state, diagnostic and
174 work vectors and sets metadata.
175
176 Default values should not be set here; please use set_vars_from_options().
177
178 All the memory allocated here is freed by IceModelVecs' destructors.
179 */
180 void IceModel::allocate_storage() {
181
182 const unsigned int WIDE_STENCIL = m_config->get_number("grid.max_stencil_width");
183
184 // FIXME: this should do for now, but we should pass a const reference to Geometry to sub-models
185 // as a function argument.
186 m_grid->variables().add(m_geometry.ice_surface_elevation);
187 m_grid->variables().add(m_geometry.ice_thickness);
188 m_grid->variables().add(m_geometry.cell_type);
189 m_grid->variables().add(m_geometry.sea_level_elevation);
190 m_grid->variables().add(m_geometry.longitude);
191 m_grid->variables().add(m_geometry.latitude);
192
193 if (m_config->get_flag("geometry.grounded_cell_fraction")) {
194 m_grid->variables().add(m_geometry.cell_grounded_fraction);
195 }
196
197 if (m_config->get_flag("geometry.part_grid.enabled")) {
198 m_grid->variables().add(m_geometry.ice_area_specific_volume);
199 }
200
201 // yield stress for basal till (plastic or pseudo-plastic model)
202 {
203 m_basal_yield_stress.create(m_grid, "tauc", WITH_GHOSTS, WIDE_STENCIL);
204 // PROPOSED standard_name = land_ice_basal_material_yield_stress
205 m_basal_yield_stress.set_attrs("diagnostic",
206 "yield stress for basal till (plastic or pseudo-plastic model)",
207 "Pa", "Pa", "", 0);
208 m_grid->variables().add(m_basal_yield_stress);
209 }
210
211 {
212 m_bedtoptemp.create(m_grid, "bedtoptemp", WITHOUT_GHOSTS);
213 m_bedtoptemp.set_attrs("diagnostic",
214 "temperature at the top surface of the bedrock thermal layer",
215 "Kelvin", "Kelvin", "", 0);
216 }
217
218 // basal melt rate
219 m_basal_melt_rate.create(m_grid, "bmelt", WITHOUT_GHOSTS);
220 m_basal_melt_rate.set_attrs("internal",
221 "ice basal melt rate from energy conservation and subshelf melt, in ice thickness per time",
222 "m s-1", "m year-1", "land_ice_basal_melt_rate", 0);
223 m_basal_melt_rate.metadata().set_string("comment", "positive basal melt rate corresponds to ice loss");
224 m_grid->variables().add(m_basal_melt_rate);
225
226 // SSA Dirichlet B.C. locations and values
227 //
228 // The mask m_ssa_dirichlet_bc_mask is also used to prescribe locations of ice thickness Dirichlet
229 // B.C. (FIXME)
230 {
231 m_ssa_dirichlet_bc_mask.create(m_grid, "bc_mask", WITH_GHOSTS, WIDE_STENCIL);
232 m_ssa_dirichlet_bc_mask.set_attrs("model_state", "Dirichlet boundary mask",
233 "", "", "", 0);
234 m_ssa_dirichlet_bc_mask.metadata().set_numbers("flag_values", {0, 1});
235 m_ssa_dirichlet_bc_mask.metadata().set_string("flag_meanings", "no_data bc_condition");
236 m_ssa_dirichlet_bc_mask.metadata().set_output_type(PISM_INT);
237 m_ssa_dirichlet_bc_mask.set_time_independent(true);
238
239 // FIXME: this is used by the inverse modeling code. Do NOT get
240 // this field from m_grid->variables() elsewhere in the code!
241 m_grid->variables().add(m_ssa_dirichlet_bc_mask);
242
243 m_ssa_dirichlet_bc_mask.set(0.0);
244 }
245 // SSA Dirichlet B.C. values
246 {
247 double fill_value = units::convert(m_sys, m_config->get_number("output.fill_value"),
248 "m year-1", "m second-1");
249 double valid_range = units::convert(m_sys, 1e6, "m year-1", "m second-1");
250 // vel_bc
251 m_ssa_dirichlet_bc_values.create(m_grid, "_ssa_bc", WITH_GHOSTS, WIDE_STENCIL); // u_ssa_bc and v_ssa_bc
252 m_ssa_dirichlet_bc_values.set_attrs("model_state",
253 "X-component of the SSA velocity boundary conditions",
254 "m s-1", "m year-1", "", 0);
255 m_ssa_dirichlet_bc_values.set_attrs("model_state",
256 "Y-component of the SSA velocity boundary conditions",
257 "m s-1", "m year-1", "", 1);
258 for (int j = 0; j < 2; ++j) {
259 m_ssa_dirichlet_bc_values.metadata(j).set_numbers("valid_range", {-valid_range, valid_range});
260 m_ssa_dirichlet_bc_values.metadata(j).set_number("_FillValue", fill_value);
261 }
262
263 // FIXME: this is used by the inverse modeling code. Do NOT get
264 // this field from m_grid->variables() elsewhere in the code!
265 m_grid->variables().add(m_ssa_dirichlet_bc_values);
266 }
267
268 // Add some variables to the list of "model state" fields.
269 m_model_state.insert(&m_ssa_dirichlet_bc_mask);
270 m_model_state.insert(&m_ssa_dirichlet_bc_values);
271
272 m_model_state.insert(&m_geometry.latitude);
273 m_model_state.insert(&m_geometry.longitude);
274 m_model_state.insert(&m_geometry.ice_thickness);
275 m_model_state.insert(&m_geometry.ice_area_specific_volume);
276 }
277
278 //! Update the surface elevation and the flow-type mask when the geometry has changed.
279 /*!
280 This procedure should be called whenever necessary to maintain consistency of geometry.
281
282 For instance, it should be called when either ice thickness or bed elevation change.
283 In particular we always want \f$h = H + b\f$ to apply at grounded points, and, on the
284 other hand, we want the mask to reflect that the ice is floating if the flotation
285 criterion applies at a point.
286
287 If `flag == REMOVE_ICEBERGS`, also calls the code which removes icebergs, to avoid
288 stress balance solver problems caused by ice that is not attached to the grounded ice
289 sheet.
290 */
291 void IceModel::enforce_consistency_of_geometry(ConsistencyFlag flag) {
292
293 m_geometry.bed_elevation.copy_from(m_beddef->bed_elevation());
294 m_geometry.sea_level_elevation.copy_from(m_sea_level->elevation());
295
296 if (m_iceberg_remover and flag == REMOVE_ICEBERGS) {
297 // The iceberg remover has to use the same mask as the stress balance code, hence the
298 // stress-balance-related threshold here.
299 m_geometry.ensure_consistency(m_config->get_number("stress_balance.ice_free_thickness_standard"));
300
301 m_iceberg_remover->update(m_ssa_dirichlet_bc_mask,
302 m_geometry.cell_type,
303 m_geometry.ice_thickness);
304 // The call above modifies ice thickness and updates the mask accordingly, but we re-compute the
305 // mask (we need to use a different threshold).
306 }
307
308 // This will ensure that ice area specific volume is zero if ice thickness is greater
309 // than zero, then compute new surface elevation and mask.
310 m_geometry.ensure_consistency(m_config->get_number("geometry.ice_free_thickness_standard"));
311
312 if (flag == REMOVE_ICEBERGS) {
313 // clean up partially-filled cells that are not next to ice
314 IceModelVec::AccessList list{&m_geometry.ice_area_specific_volume,
315 &m_geometry.cell_type};
316
317 for (Points p(*m_grid); p; p.next()) {
318 const int i = p.i(), j = p.j();
319
320 if (m_geometry.ice_area_specific_volume(i, j) > 0.0 and
321 not m_geometry.cell_type.next_to_ice(i, j)) {
322 m_geometry.ice_area_specific_volume(i, j) = 0.0;
323 }
324 }
325 }
326 }
327
328 stressbalance::Inputs IceModel::stress_balance_inputs() {
329 stressbalance::Inputs result;
330 if (m_config->get_flag("geometry.update.use_basal_melt_rate")) {
331 result.basal_melt_rate = &m_basal_melt_rate;
332 }
333
334 result.basal_yield_stress = &m_basal_yield_stress;
335 result.melange_back_pressure = &m_ocean->melange_back_pressure_fraction();
336 result.geometry = &m_geometry;
337 result.new_bed_elevation = m_new_bed_elevation;
338 result.enthalpy = &m_energy_model->enthalpy();
339 result.age = m_age_model ? &m_age_model->age() : nullptr;
340
341 if (m_config->get_flag("stress_balance.ssa.dirichlet_bc")) {
342 result.bc_mask = &m_ssa_dirichlet_bc_mask;
343 result.bc_values = &m_ssa_dirichlet_bc_values;
344 }
345
346 if (m_config->get_flag("fracture_density.enabled")) {
347 result.fracture_density = &m_fracture->density();
348 }
349
350 return result;
351 }
352
353 energy::Inputs IceModel::energy_model_inputs() {
354 energy::Inputs result;
355
356 result.basal_frictional_heating = &m_stress_balance->basal_frictional_heating();
357 result.basal_heat_flux = &m_btu->flux_through_top_surface(); // bedrock thermal layer
358 result.cell_type = &m_geometry.cell_type; // geometry
359 result.ice_thickness = &m_geometry.ice_thickness; // geometry
360 result.shelf_base_temp = &m_ocean->shelf_base_temperature(); // ocean model
361 result.till_water_thickness = &m_subglacial_hydrology->till_water_thickness();
362 result.surface_liquid_fraction = &m_surface->liquid_water_fraction(); // surface model
363 result.surface_temp = &m_surface->temperature(); // surface model
364
365 result.volumetric_heating_rate = &m_stress_balance->volumetric_strain_heating();
366 result.u3 = &m_stress_balance->velocity_u();
367 result.v3 = &m_stress_balance->velocity_v();
368 result.w3 = &m_stress_balance->velocity_w();
369
370 result.check(); // make sure all data members were set
371
372 return result;
373 }
374
375 YieldStressInputs IceModel::yield_stress_inputs() {
376 YieldStressInputs result;
377
378 result.geometry = &m_geometry;
379 result.till_water_thickness = &m_subglacial_hydrology->till_water_thickness();
380 result.subglacial_water_thickness = &m_subglacial_hydrology->subglacial_water_thickness();
381
382 return result;
383 }
384
385 //! The contents of the main PISM time-step.
386 /*!
387 During the time-step we perform the following actions:
388 */
389 void IceModel::step(bool do_mass_continuity,
390 bool do_skip) {
391
392 const Profiling &profiling = m_ctx->profiling();
393
394 double current_time = m_time->current();
395
396 //! \li call pre_step_hook() to let derived classes do more
397 pre_step_hook();
398
399 //! \li update the velocity field; in some cases the whole three-dimensional
400 //! field is updated and in some cases just the vertically-averaged
401 //! horizontal velocity is updated
402
403 // always "update" ice velocity (possibly trivially); only update
404 // SSA and only update velocities at depth if suggested by temp and age
405 // stability criterion; note *lots* of communication is avoided by skipping
406 // SSA (and temp/age)
407
408 const bool updateAtDepth = (m_skip_countdown == 0);
409
410 // Combine basal melt rate in grounded (computed during the energy
411 // step) and floating (provided by an ocean model) areas.
412 //
413 // Basal melt rate may be used by a stress balance model to compute vertical velocity of
414 // ice.
415 {
416 combine_basal_melt_rate(m_geometry,
417 m_ocean->shelf_base_mass_flux(),
418 m_energy_model->basal_melt_rate(),
419 m_basal_melt_rate);
420 }
421
422 try {
423 profiling.begin("stress_balance");
424 m_stress_balance->update(stress_balance_inputs(), updateAtDepth);
425 profiling.end("stress_balance");
426 } catch (RuntimeError &e) {
427 std::string output_file = m_config->get_string("output.file_name");
428
429 if (output_file.empty()) {
430 m_log->message(2, "WARNING: output.file_name is empty. Using unnamed.nc instead.\n");
431 output_file = "unnamed.nc";
432 }
433
434 std::string o_file = filename_add_suffix(output_file,
435 "_stressbalance_failed", "");
436 File file(m_grid->com, o_file,
437 string_to_backend(m_config->get_string("output.format")),
438 PISM_READWRITE_MOVE,
439 m_ctx->pio_iosys_id());
440
441 update_run_stats();
442 write_metadata(file, WRITE_MAPPING, PREPEND_HISTORY);
443
444 save_variables(file, INCLUDE_MODEL_STATE, output_variables("small"),
445 m_time->current());
446
447 e.add_context("performing a time step. (Note: Model state was saved to '%s'.)",
448 o_file.c_str());
449 throw;
450 }
451
452
453 m_stdout_flags += m_stress_balance->stdout_report();
454
455 m_stdout_flags += (updateAtDepth ? "v" : "V");
456
457 //! \li determine the time step according to a variety of stability criteria
458 max_timestep(m_dt, m_skip_countdown);
459
460 //! \li update the yield stress for the plastic till model (if appropriate)
461 if (m_basal_yield_stress_model) {
462 profiling.begin("basal_yield_stress");
463 m_basal_yield_stress_model->update(yield_stress_inputs(), current_time, m_dt);
464 profiling.end("basal_yield_stress");
465 m_basal_yield_stress.copy_from(m_basal_yield_stress_model->basal_material_yield_stress());
466 m_stdout_flags += "y";
467 } else {
468 m_stdout_flags += "$";
469 }
470
471 dt_TempAge += m_dt;
472
473 //! \li update the age of the ice (if appropriate)
474 if (m_age_model and updateAtDepth) {
475 AgeModelInputs inputs;
476 inputs.ice_thickness = &m_geometry.ice_thickness;
477 inputs.u3 = &m_stress_balance->velocity_u();
478 inputs.v3 = &m_stress_balance->velocity_v();
479 inputs.w3 = &m_stress_balance->velocity_w();
480
481 profiling.begin("age");
482 m_age_model->update(current_time, dt_TempAge, inputs);
483 profiling.end("age");
484 m_stdout_flags += "a";
485 } else {
486 m_stdout_flags += "$";
487 }
488
489 //! \li update the enthalpy (or temperature) field according to the conservation of
490 //! energy model based (especially) on the new velocity field; see
491 //! energy_step()
492 if (updateAtDepth) { // do the energy step
493 profiling.begin("energy");
494 energy_step();
495 profiling.end("energy");
496 m_stdout_flags += "E";
497 } else {
498 m_stdout_flags += "$";
499 }
500
501 //! \li update the fracture density field; see update_fracture_density()
502 if (m_config->get_flag("fracture_density.enabled")) {
503 profiling.begin("fracture_density");
504 update_fracture_density();
505 profiling.end("fracture_density");
506 }
507
508 //! \li update the thickness of the ice according to the mass conservation model and calving
509 //! parameterizations
510
511 // FIXME: thickness B.C. mask should be separate
512 IceModelVec2Int &thickness_bc_mask = m_ssa_dirichlet_bc_mask;
513
514 if (do_mass_continuity) {
515 profiling.begin("mass_transport");
516 {
517 // Note that there are three adaptive time-stepping criteria. Two of them (using max.
518 // diffusion and 2D CFL) are limiting the mass-continuity time-step and the third (3D
519 // CFL) limits the energy and age time-steps.
520
521 // The mass-continuity time-step is usually smaller, and the skipping mechanism lets us
522 // do several mass-continuity steps for each energy step.
523
524 // When -no_mass is set, mass-continuity-related time-step restrictions are disabled,
525 // making "skipping" unnecessary.
526
527 // This is why the following two lines appear here and are executed only if
528 // do_mass_continuity is true.
529 if (do_skip and m_skip_countdown > 0) {
530 m_skip_countdown--;
531 }
532
533 m_geometry_evolution->flow_step(m_geometry,
534 m_dt,
535 m_stress_balance->advective_velocity(),
536 m_stress_balance->diffusive_flux(),
537 m_ssa_dirichlet_bc_mask,
538 thickness_bc_mask);
539
540 m_geometry_evolution->apply_flux_divergence(m_geometry);
541
542 enforce_consistency_of_geometry(DONT_REMOVE_ICEBERGS);
543 }
544 profiling.end("mass_transport");
545
546 // calving, frontal melt, and discharge accounting
547 profiling.begin("front_retreat");
548 front_retreat_step();
549 profiling.end("front_retreat");
550
551 m_stdout_flags += "h";
552 } else {
553 m_stdout_flags += "$";
554 }
555
556 profiling.begin("sea_level");
557 m_sea_level->update(m_geometry, current_time, m_dt);
558 profiling.end("sea_level");
559
560 profiling.begin("ocean");
561 m_ocean->update(m_geometry, current_time, m_dt);
562 profiling.end("ocean");
563
564 // The sea level elevation might have changed, so we need to update the mask, etc. Note
565 // that THIS MAY PRODUCE ICEBERGS, but we assume that the surface model does not care.
566 enforce_consistency_of_geometry(DONT_REMOVE_ICEBERGS);
567
568 //! \li Update surface and ocean models.
569 profiling.begin("surface");
570 m_surface->update(m_geometry, current_time, m_dt);
571 profiling.end("surface");
572
573
574 if (do_mass_continuity) {
575 // compute and apply effective surface and basal mass balance
576
577 m_geometry_evolution->source_term_step(m_geometry, m_dt,
578 thickness_bc_mask,
579 m_surface->mass_flux(),
580 m_basal_melt_rate);
581 m_geometry_evolution->apply_mass_fluxes(m_geometry);
582
583 // add removed icebergs to discharge due to calving
584 {
585 IceModelVec2S
586 &old_H = m_work2d[0],
587 &old_Href = m_work2d[1];
588
589 {
590 old_H.copy_from(m_geometry.ice_thickness);
591 old_Href.copy_from(m_geometry.ice_area_specific_volume);
592 }
593
594 // the last call has to remove icebergs
595 enforce_consistency_of_geometry(REMOVE_ICEBERGS);
596
597 bool add_values = true;
598 compute_geometry_change(m_geometry.ice_thickness,
599 m_geometry.ice_area_specific_volume,
600 old_H, old_Href,
601 add_values,
602 m_thickness_change.calving);
603 }
604 }
605
606 //! \li update the state variables in the subglacial hydrology model (typically
607 //! water thickness and sometimes pressure)
608 profiling.begin("basal_hydrology");
609 hydrology_step();
610 profiling.end("basal_hydrology");
611
612 //! \li compute the bed deformation, which depends on current thickness, bed elevation,
613 //! and sea level
614 if (m_beddef) {
615 int topg_state_counter = m_beddef->bed_elevation().state_counter();
616
617 profiling.begin("bed_deformation");
618 m_beddef->update(m_geometry.ice_thickness,
619 m_geometry.sea_level_elevation,
620 current_time, m_dt);
621 profiling.end("bed_deformation");
622
623 if (m_beddef->bed_elevation().state_counter() != topg_state_counter) {
624 m_new_bed_elevation = true;
625 } else {
626 m_new_bed_elevation = false;
627 }
628 } else {
629 m_new_bed_elevation = false;
630 }
631
632 if (m_new_bed_elevation) {
633 enforce_consistency_of_geometry(DONT_REMOVE_ICEBERGS);
634 m_stdout_flags += "b";
635 } else {
636 m_stdout_flags += " ";
637 }
638
639 //! \li call post_step_hook() to let derived classes do more
640 post_step_hook();
641
642 // Done with the step; now adopt the new time.
643 m_time->step(m_dt);
644
645 if (updateAtDepth) {
646 t_TempAge = m_time->current();
647 dt_TempAge = 0.0;
648 }
649
650 // Check if the ice thickness exceeded the height of the computational box and stop if it did.
651 const bool thickness_too_high = check_maximum_ice_thickness(m_geometry.ice_thickness);
652
653 if (thickness_too_high) {
654 std::string output_file = m_config->get_string("output.file_name");
655
656 if (output_file.empty()) {
657 m_log->message(2, "WARNING: output.file_name is empty. Using unnamed.nc instead.");
658 output_file = "unnamed.nc";
659 }
660
661 std::string o_file = filename_add_suffix(output_file,
662 "_max_thickness", "");
663 File file(m_grid->com,
664 o_file,
665 string_to_backend(m_config->get_string("output.format")),
666 PISM_READWRITE_MOVE,
667 m_ctx->pio_iosys_id());
668
669 update_run_stats();
670 write_metadata(file, WRITE_MAPPING, PREPEND_HISTORY);
671
672 save_variables(file, INCLUDE_MODEL_STATE, output_variables("small"),
673 m_time->current());
674
675 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
676 "Ice thickness exceeds the height of the computational box (%7.4f m).\n"
677 "The model state was saved to '%s'. To continue this simulation,\n"
678 "run with\n"
679 "-i %s -bootstrap -regrid_file %s -allow_extrapolation -Lz N [other options]\n"
680 "where N > %7.4f.",
681 m_grid->Lz(), o_file.c_str(), o_file.c_str(), o_file.c_str(), m_grid->Lz());
682 }
683
684 // end the flag line
685 m_stdout_flags += " " + m_adaptive_timestep_reason;
686 }
687
688 /*!
689 * Note: don't forget to update IceRegionalModel::hydrology_step() if necessary.
690 */
691 void IceModel::hydrology_step() {
692 hydrology::Inputs inputs;
693
694 IceModelVec2S &sliding_speed = m_work2d[0];
695 sliding_speed.set_to_magnitude(m_stress_balance->advective_velocity());
696
697 inputs.no_model_mask = nullptr;
698 inputs.geometry = &m_geometry;
699 inputs.surface_input_rate = nullptr;
700 inputs.basal_melt_rate = &m_basal_melt_rate;
701 inputs.ice_sliding_speed = &sliding_speed;
702
703 if (m_surface_input_for_hydrology) {
704 m_surface_input_for_hydrology->update(m_time->current(), m_dt);
705 m_surface_input_for_hydrology->average(m_time->current(), m_dt);
706 inputs.surface_input_rate = m_surface_input_for_hydrology.get();
707 } else if (m_config->get_flag("hydrology.surface_input_from_runoff")) {
708 // convert [kg m-2] to [kg m-2 s-1]
709 IceModelVec2S &surface_input_rate = m_work2d[1];
710 surface_input_rate.copy_from(m_surface->runoff());
711 surface_input_rate.scale(1.0 / m_dt);
712 inputs.surface_input_rate = &surface_input_rate;
713 }
714
715 m_subglacial_hydrology->update(m_time->current(), m_dt, inputs);
716 }
717
718 //! Virtual. Does nothing in `IceModel`. Derived classes can do more computation in each time step.
719 void IceModel::pre_step_hook() {
720 // empty
721 }
722
723 //! Virtual. Does nothing in `IceModel`. Derived classes can do more computation in each time step.
724 void IceModel::post_step_hook() {
725 // empty
726 }
727
728 /**
729 * Run the time-stepping loop from the current model time to `time`.
730 *
731 * This should be called by the coupler controlling PISM when it is
732 * running alongside a GCM.
733 *
734 * @param run_end model time (in seconds) to run to
735 *
736 * @return 0 on success
737 */
738 void IceModel::run_to(double run_end) {
739
740 m_time->set_end(run_end);
741
742 run();
743 }
744
745
746 /**
747 * Run the time-stepping loop from the current time until the time
748 * specified by the IceModel::grid::time object.
749 *
750 * This is the method used by PISM in the "standalone" mode.
751 */
752 void IceModel::run() {
753 const Profiling &profiling = m_ctx->profiling();
754
755 bool do_mass_conserve = m_config->get_flag("geometry.update.enabled");
756 bool do_energy = m_config->get_flag("energy.enabled");
757 bool do_skip = m_config->get_flag("time_stepping.skip.enabled");
758
759 int stepcount = m_config->get_flag("time_stepping.count_steps") ? 0 : -1;
760
761 // de-allocate diagnostics that are not needed
762 prune_diagnostics();
763
764 // Enforce consistency *and* remove icebergs. During time-stepping we remove icebergs at
765 // the end of the time step, so we need to ensure that ice geometry is "OK" before the
766 // first step.
767 enforce_consistency_of_geometry(REMOVE_ICEBERGS);
768
769 // Update spatially-variable diagnostics at the beginning of the run.
770 write_extras();
771
772 // Update scalar time series to remember the state at the beginning of the run.
773 // This is needed to compute rates of change of the ice mass, volume, etc.
774 {
775 const double time = m_time->current();
776 for (auto d : m_ts_diagnostics) {
777 d.second->update(time, time);
778 }
779 }
780
781 m_log->message(2, "running forward ...\n");
782
783 m_stdout_flags.erase(); // clear it out
784 print_summary_line(true, do_energy, 0.0, 0.0, 0.0, 0.0, 0.0);
785 m_adaptive_timestep_reason = '$'; // no reason for no timestep
786 print_summary(do_energy); // report starting state
787
788 t_TempAge = m_time->current();
789 dt_TempAge = 0.0;
790
791 // main loop for time evolution
792 // IceModel::step calls Time::step(dt), ensuring that this while loop
793 // will terminate
794 profiling.stage_begin("time-stepping loop");
795 while (m_time->current() < m_time->end()) {
796
797 m_stdout_flags.erase(); // clear it out
798
799 step(do_mass_conserve, do_skip);
800
801 update_diagnostics(m_dt);
802
803 // report a summary for major steps or the last one
804 bool updateAtDepth = m_skip_countdown == 0;
805 bool tempAgeStep = updateAtDepth and (m_age_model or do_energy);
806
807 const bool show_step = tempAgeStep or m_adaptive_timestep_reason == "end of the run";
808 print_summary(show_step);
809
810 // update viewers before writing extras because writing extras resets diagnostics
811 update_viewers();
812
813 // writing these fields here ensures that we do it after the last time-step
814 profiling.begin("io");
815 write_snapshot();
816 write_extras();
817 write_backup();
818 profiling.end("io");
819
820 if (stepcount >= 0) {
821 stepcount++;
822 }
823 if (process_signals() != 0) {
824 break;
825 }
826 } // end of the time-stepping loop
827
828 profiling.stage_end("time-stepping loop");
829
830 if (stepcount >= 0) {
831 m_log->message(1,
832 "count_time_steps: run() took %d steps\n"
833 "average dt = %.6f years\n",
834 stepcount,
835 units::convert(m_sys, m_time->end() - m_time->start(), "seconds", "years")/(double)stepcount);
836 }
837 }
838
839 //! Manage the initialization of the IceModel object.
840 /*!
841 Please see the documenting comments of the functions called below to find
842 explanations of their intended uses.
843 */
844 void IceModel::init() {
845 // Get the start time in seconds and ensure that it is consistent
846 // across all processors.
847 m_start_time = GlobalMax(m_grid->com, get_time());
848
849 const Profiling &profiling = m_ctx->profiling();
850
851 profiling.begin("initialization");
852
853 //! The IceModel initialization sequence is this:
854
855 //! 1) Initialize model time:
856 time_setup();
857
858 //! 2) Process the options:
859 process_options();
860
861 //! 3) Memory allocation:
862 allocate_storage();
863
864 //! 4) Allocate PISM components modeling some physical processes.
865 allocate_submodels();
866
867 //! 6) Initialize coupler models and fill the model state variables
868 //! (from a PISM output file, from a bootstrapping file using some
869 //! modeling choices or using formulas). Calls IceModel::regrid()
870 model_state_setup();
871
872 //! 7) Report grid parameters:
873 m_grid->report_parameters();
874
875 //! 8) Miscellaneous stuff: set up the bed deformation model, initialize the
876 //! basal till model, initialize snapshots. This has to happen *after*
877 //! regridding.
878 misc_setup();
879
880 profiling.end("initialization");
881 }
882
883 const Geometry& IceModel::geometry() const {
884 return m_geometry;
885 }
886
887 const GeometryEvolution& IceModel::geometry_evolution() const {
888 return *m_geometry_evolution;
889 }
890
891 const stressbalance::StressBalance* IceModel::stress_balance() const {
892 return this->m_stress_balance.get();
893 }
894
895 const ocean::OceanModel* IceModel::ocean_model() const {
896 return m_ocean.get();
897 }
898
899 const bed::BedDef* IceModel::bed_model() const {
900 return m_beddef;
901 }
902
903 const energy::BedThermalUnit* IceModel::bedrock_thermal_model() const {
904 return m_btu;
905 }
906
907 const energy::EnergyModel* IceModel::energy_balance_model() const {
908 return m_energy_model;
909 }
910
911 /*!
912 * Return thickness change due to calving (over the last time step).
913 */
914 const IceModelVec2S& IceModel::calving() const {
915 return m_thickness_change.calving;
916 }
917
918 /*!
919 * Return thickness change due to frontal melt (over the last time step).
920 */
921 const IceModelVec2S& IceModel::frontal_melt() const {
922 return m_thickness_change.frontal_melt;
923 }
924
925 /*!
926 * Return thickness change due to forced retreat (over the last time step).
927 */
928 const IceModelVec2S& IceModel::forced_retreat() const {
929 return m_thickness_change.forced_retreat;
930 }
931
932 void warn_about_missing(const Logger &log,
933 const std::set<std::string> &vars,
934 const std::string &type,
935 const std::set<std::string> &available,
936 bool stop) {
937 std::vector<std::string> missing;
938 for (auto v : vars) {
939 if (available.find(v) == available.end()) {
940 missing.push_back(v);
941 }
942 }
943
944 if (not missing.empty()) {
945 size_t N = missing.size();
946 const char
947 *ending = N > 1 ? "s" : "",
948 *verb = N > 1 ? "are" : "is";
949 if (stop) {
950 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
951 "%s variable%s %s %s not available!\n"
952 "Available variables:\n- %s",
953 type.c_str(),
954 ending,
955 join(missing, ",").c_str(),
956 verb,
957 set_join(available, ",\n- ").c_str());
958 } else {
959 log.message(2,
960 "\nWARNING: %s variable%s %s %s not available!\n\n",
961 type.c_str(),
962 ending,
963 join(missing, ",").c_str(),
964 verb);
965 }
966 }
967 }
968
969 /*!
970 * De-allocate diagnostics that were not requested.
971 *
972 * Checks viewers, -extra_vars, -backup, -save_vars, and regular output.
973 *
974 * FIXME: I need to make sure that these reporting mechanisms are active. It is possible that
975 * variables are on a list, but that list is not actually used.
976 */
977 void IceModel::prune_diagnostics() {
978
979 // get the list of available diagnostics
980 std::set<std::string> available;
981 for (auto d : m_diagnostics) {
982 available.insert(d.first);
983 }
984
985 auto m_extra_stop = m_config->get_flag("output.extra.stop_missing");
986 warn_about_missing(*m_log, m_output_vars, "output", available, false);
987 warn_about_missing(*m_log, m_snapshot_vars, "snapshot", available, false);
988 warn_about_missing(*m_log, m_backup_vars, "backup", available, false);
989 warn_about_missing(*m_log, m_extra_vars, "diagnostic", available, m_extra_stop);
990
991 // get the list of requested diagnostics
992 auto requested = set_split(m_config->get_string("output.runtime.viewer.variables"), ',');
993 requested = combine(requested, m_output_vars);
994 requested = combine(requested, m_snapshot_vars);
995 requested = combine(requested, m_extra_vars);
996 requested = combine(requested, m_backup_vars);
997
998 // de-allocate diagnostics that were not requested
999 for (auto v : available) {
1000 if (requested.find(v) == requested.end()) {
1001 m_diagnostics.erase(v);
1002 }
1003 }
1004
1005 // scalar time series
1006 std::vector<std::string> missing;
1007 if (not m_ts_filename.empty() and m_ts_vars.empty()) {
1008 // use all diagnostics
1009 } else {
1010 TSDiagnosticList diagnostics;
1011 for (auto v : m_ts_vars) {
1012 if (m_ts_diagnostics.find(v) != m_ts_diagnostics.end()) {
1013 diagnostics[v] = m_ts_diagnostics[v];
1014 } else {
1015 missing.push_back(v);
1016 }
1017 }
1018 // replace m_ts_diagnostics with requested diagnostics, de-allocating the rest
1019 m_ts_diagnostics = diagnostics;
1020 }
1021
1022 if (not missing.empty()) {
1023 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
1024 "requested scalar diagnostics %s are not available",
1025 join(missing, ",").c_str());
1026 }
1027 }
1028
1029 /*!
1030 * Update diagnostics.
1031 *
1032 * This usually involves accumulating data needed to computed time-averaged quantities.
1033 *
1034 * Call this after prune_diagnostics() to avoid unnecessary work.
1035 */
1036 void IceModel::update_diagnostics(double dt) {
1037 for (auto d : m_diagnostics) {
1038 d.second->update(dt);
1039 }
1040
1041 const double time = m_time->current();
1042 for (auto d : m_ts_diagnostics) {
1043 d.second->update(time - dt, time);
1044 }
1045 }
1046
1047 /*!
1048 * Reset accumulators in diagnostics that compute time-averaged quantities.
1049 */
1050 void IceModel::reset_diagnostics() {
1051 for (auto d : m_diagnostics) {
1052 d.second->reset();
1053 }
1054 }
1055
1056 IceModel::ThicknessChanges::ThicknessChanges(IceGrid::ConstPtr grid)
1057 : calving(grid, "thickness_change_due_to_calving", WITHOUT_GHOSTS),
1058 frontal_melt(grid, "thickness_change_due_to_frontal_melt", WITHOUT_GHOSTS),
1059 forced_retreat(grid, "thickness_change_due_to_forced_retreat", WITHOUT_GHOSTS) {
1060 // empty
1061 }
1062
1063 } // end of namespace pism