tinitialization.cc - 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
---
tinitialization.cc (35224B)
---
1 // Copyright (C) 2009--2020 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 //This file contains various initialization routines. See the IceModel::init()
20 //documentation comment in iceModel.cc for the order in which they are called.
21
22 #include "IceModel.hh"
23 #include "pism/basalstrength/ConstantYieldStress.hh"
24 #include "pism/basalstrength/MohrCoulombYieldStress.hh"
25 #include "pism/basalstrength/basal_resistance.hh"
26 #include "pism/frontretreat/util/IcebergRemover.hh"
27 #include "pism/frontretreat/calving/CalvingAtThickness.hh"
28 #include "pism/frontretreat/calving/EigenCalving.hh"
29 #include "pism/frontretreat/calving/FloatKill.hh"
30 #include "pism/frontretreat/calving/HayhurstCalving.hh"
31 #include "pism/frontretreat/calving/vonMisesCalving.hh"
32 #include "pism/energy/BedThermalUnit.hh"
33 #include "pism/hydrology/NullTransport.hh"
34 #include "pism/hydrology/Routing.hh"
35 #include "pism/hydrology/SteadyState.hh"
36 #include "pism/hydrology/Distributed.hh"
37 #include "pism/stressbalance/StressBalance.hh"
38 #include "pism/stressbalance/sia/SIAFD.hh"
39 #include "pism/stressbalance/ssa/SSAFD.hh"
40 #include "pism/stressbalance/ssa/SSAFEM.hh"
41 #include "pism/util/Mask.hh"
42 #include "pism/util/ConfigInterface.hh"
43 #include "pism/util/Time.hh"
44 #include "pism/util/error_handling.hh"
45 #include "pism/util/io/File.hh"
46 #include "pism/util/pism_options.hh"
47 #include "pism/coupler/OceanModel.hh"
48 #include "pism/coupler/SurfaceModel.hh"
49 #include "pism/coupler/atmosphere/Factory.hh"
50 #include "pism/coupler/ocean/Factory.hh"
51 #include "pism/coupler/ocean/Initialization.hh"
52 #include "pism/coupler/ocean/sea_level/Factory.hh"
53 #include "pism/coupler/ocean/sea_level/Initialization.hh"
54 #include "pism/coupler/surface/Factory.hh"
55 #include "pism/coupler/surface/Initialization.hh"
56 #include "pism/earth/LingleClark.hh"
57 #include "pism/earth/BedDef.hh"
58 #include "pism/util/EnthalpyConverter.hh"
59 #include "pism/util/Vars.hh"
60 #include "pism/util/io/io_helpers.hh"
61 #include "pism/util/projection.hh"
62 #include "pism/util/pism_utilities.hh"
63 #include "pism/age/AgeModel.hh"
64 #include "pism/energy/EnthalpyModel.hh"
65 #include "pism/energy/TemperatureModel.hh"
66 #include "pism/fracturedensity/FractureDensity.hh"
67 #include "pism/frontretreat/FrontRetreat.hh"
68 #include "pism/frontretreat/PrescribedRetreat.hh"
69 #include "pism/coupler/frontalmelt/Factory.hh"
70 #include "pism/coupler/util/options.hh" // ForcingOptions
71
72 namespace pism {
73
74 //! Initialize time from an input file or command-line options.
75 void IceModel::time_setup() {
76 initialize_time(m_grid->com,
77 m_config->get_string("time.dimension_name"),
78 *m_log, *m_time);
79
80 bool use_calendar = m_config->get_flag("output.runtime.time_use_calendar");
81
82 if (use_calendar) {
83 m_log->message(2,
84 "* Run time: [%s, %s] (%s years, using the '%s' calendar)\n",
85 m_time->start_date().c_str(),
86 m_time->end_date().c_str(),
87 m_time->run_length().c_str(),
88 m_time->calendar().c_str());
89 } else {
90 std::string time_units = m_config->get_string("output.runtime.time_unit_name");
91
92 double
93 start = m_time->convert_time_interval(m_time->start(), time_units),
94 end = m_time->convert_time_interval(m_time->end(), time_units),
95 length = end - start;
96
97 m_log->message(2,
98 "* Run time: [%f %s, %f %s] (%f %s)\n",
99 start, time_units.c_str(),
100 end, time_units.c_str(),
101 length, time_units.c_str());
102 }
103 }
104
105 //! Sets the starting values of model state variables.
106 /*!
107 There are two cases:
108
109 1) Initializing from a PISM output file.
110
111 2) Setting the values using command-line options only (verification and
112 simplified geometry runs, for example) or from a bootstrapping file, using
113 heuristics to fill in missing and 3D fields.
114
115 Calls IceModel::regrid().
116
117 This function is called after all the memory allocation is done and all the
118 physical parameters are set.
119
120 Calling this method should be all one needs to set model state variables.
121 Please avoid modifying them in other parts of the initialization sequence.
122
123 Also, please avoid operations that would make it unsafe to call this more
124 than once (memory allocation is one example).
125 */
126 void IceModel::model_state_setup() {
127
128 // Check if we are initializing from a PISM output file:
129 InputOptions input = process_input_options(m_ctx->com(), m_config);
130
131 const bool use_input_file = input.type == INIT_BOOTSTRAP or input.type == INIT_RESTART;
132
133 std::unique_ptr<File> input_file;
134
135 if (use_input_file) {
136 input_file.reset(new File(m_grid->com, input.filename, PISM_GUESS, PISM_READONLY));
137 }
138
139 // Initialize 2D fields owned by IceModel (ice geometry, etc)
140 {
141 switch (input.type) {
142 case INIT_RESTART:
143 restart_2d(*input_file, input.record);
144 break;
145 case INIT_BOOTSTRAP:
146 bootstrap_2d(*input_file);
147 break;
148 case INIT_OTHER:
149 default:
150 initialize_2d();
151 }
152
153 regrid();
154 }
155
156 // Get projection information and compute latitudes and longitudes *before* a component
157 // decides to use them...
158 {
159 if (use_input_file) {
160 std::string mapping_name = m_grid->get_mapping_info().mapping.get_name();
161 MappingInfo info = get_projection_info(*input_file, mapping_name,
162 m_ctx->unit_system());
163
164 if (not info.proj.empty()) {
165 m_log->message(2, "* Got projection parameters \"%s\" from \"%s\".\n",
166 info.proj.c_str(), input.filename.c_str());
167 }
168
169 m_output_global_attributes.set_string("proj", info.proj);
170 m_grid->set_mapping_info(info);
171
172 std::string history = input_file->read_text_attribute("PISM_GLOBAL", "history");
173 m_output_global_attributes.set_string("history",
174 history + m_output_global_attributes.get_string("history"));
175
176 }
177
178 compute_lat_lon();
179 }
180
181 m_sea_level->init(m_geometry);
182
183 // Initialize a bed deformation model.
184 if (m_beddef) {
185 m_beddef->init(input, m_geometry.ice_thickness, m_sea_level->elevation());
186 m_grid->variables().add(m_beddef->bed_elevation());
187 m_grid->variables().add(m_beddef->uplift());
188 }
189
190 // Now ice thickness, bed elevation, and sea level are available, so we can compute the ice
191 // surface elevation and the cell type mask. This also ensures consistency of ice geometry.
192 //
193 // The stress balance code is executed near the beginning of a step and ice geometry is
194 // updated near the end, so during the second time step the stress balance code is
195 // guaranteed not to see "icebergs". Here we make sure that the first time step is OK
196 // too.
197 enforce_consistency_of_geometry(REMOVE_ICEBERGS);
198
199 // By now ice geometry is set (including regridding) and so we can initialize the ocean model,
200 // which may need ice thickness, bed topography, and the cell type mask.
201 {
202 m_ocean->init(m_geometry);
203 }
204
205 // Now surface elevation is initialized, so we can initialize surface models (some use
206 // elevation-based parameterizations of surface temperature and/or mass balance).
207 m_surface->init(m_geometry);
208
209 if (m_subglacial_hydrology) {
210
211 switch (input.type) {
212 case INIT_RESTART:
213 m_subglacial_hydrology->restart(*input_file, input.record);
214 break;
215 case INIT_BOOTSTRAP:
216 m_subglacial_hydrology->bootstrap(*input_file,
217 m_geometry.ice_thickness);
218 break;
219 case INIT_OTHER:
220 {
221 IceModelVec2S
222 &W_till = m_work2d[0],
223 &W = m_work2d[1],
224 &P = m_work2d[2];
225
226 W_till.set(m_config->get_number("bootstrapping.defaults.tillwat"));
227 W.set(m_config->get_number("bootstrapping.defaults.bwat"));
228 P.set(m_config->get_number("bootstrapping.defaults.bwp"));
229
230 m_subglacial_hydrology->init(W_till, W, P);
231 break;
232 }
233 }
234 }
235
236 // basal_yield_stress_model->init() needs till water thickness so this must happen
237 // after subglacial_hydrology->init()
238 if (m_basal_yield_stress_model) {
239 auto inputs = yield_stress_inputs();
240
241 switch (input.type) {
242 case INIT_RESTART:
243 m_basal_yield_stress_model->restart(*input_file, input.record);
244 break;
245 case INIT_BOOTSTRAP:
246 m_basal_yield_stress_model->bootstrap(*input_file, inputs);
247 break;
248 default:
249 case INIT_OTHER:
250 m_basal_yield_stress_model->init(inputs);
251 }
252 m_basal_yield_stress.copy_from(m_basal_yield_stress_model->basal_material_yield_stress());
253 }
254
255 // Initialize the bedrock thermal layer model.
256 //
257 // If
258 // - PISM is bootstrapping and
259 // - we are using a non-zero-thickness thermal layer
260 //
261 // initialization of m_btu requires the temperature at the top of the bedrock. This is a problem
262 // because get_bed_top_temp() uses the enthalpy field that is not initialized until later and
263 // bootstrapping enthalpy uses the flux through the bottom surface of the ice (top surface of the
264 // bedrock) provided by m_btu.
265 //
266 // We get out of this by using the fact that the full state of m_btu is not needed and
267 // bootstrapping of the temperature field can be delayed.
268 //
269 // Note that to bootstrap m_btu we use the steady state solution of the heat equation in columns
270 // of the bedrock (a straight line at each column), so the flux through the top surface of the
271 // bedrock after bootstrapping is the same as the time-independent geothermal flux applied at the
272 // BOTTOM surface of the bedrock layer.
273 //
274 // The code then delays bootstrapping of the thickness field until the first time step.
275 if (m_btu) {
276 m_btu->init(input);
277 }
278
279 if (m_age_model) {
280 m_age_model->init(input);
281 m_grid->variables().add(m_age_model->age());
282 }
283
284 // Initialize the energy balance sub-model.
285 {
286 switch (input.type) {
287 case INIT_RESTART:
288 {
289 m_energy_model->restart(*input_file, input.record);
290 break;
291 }
292 case INIT_BOOTSTRAP:
293 {
294
295 m_energy_model->bootstrap(*input_file,
296 m_geometry.ice_thickness,
297 m_surface->temperature(),
298 m_surface->mass_flux(),
299 m_btu->flux_through_top_surface());
300 break;
301 }
302 case INIT_OTHER:
303 default:
304 {
305 m_basal_melt_rate.set(m_config->get_number("bootstrapping.defaults.bmelt"));
306
307 m_energy_model->initialize(m_basal_melt_rate,
308 m_geometry.ice_thickness,
309 m_surface->temperature(),
310 m_surface->mass_flux(),
311 m_btu->flux_through_top_surface());
312
313 }
314 }
315 m_grid->variables().add(m_energy_model->enthalpy());
316 }
317
318 // this has to go after we add enthalpy to m_grid->variables()
319 if (m_stress_balance) {
320 m_stress_balance->init();
321 }
322
323 // miscellaneous steps
324 {
325 reset_counters();
326
327 auto startstr = pism::printf("PISM (%s) started on %d procs.",
328 pism::revision, (int)m_grid->size());
329 prepend_history(startstr + args_string());
330 }
331 }
332
333 //! Initialize 2D model state fields managed by IceModel from a file (for re-starting).
334 /*!
335 * This method should eventually go away as IceModel turns into a "coupler" and all physical
336 * processes are handled by sub-models.
337 */
338 void IceModel::restart_2d(const File &input_file, unsigned int last_record) {
339 std::string filename = input_file.filename();
340
341 m_log->message(2, "initializing 2D fields from NetCDF file '%s'...\n", filename.c_str());
342
343 for (auto variable : m_model_state) {
344 variable->read(input_file, last_record);
345 }
346 }
347
348 void IceModel::bootstrap_2d(const File &input_file) {
349
350 m_log->message(2, "bootstrapping from file '%s'...\n", input_file.filename().c_str());
351
352 auto usurf = input_file.find_variable("usurf", "surface_altitude");
353
354 bool mask_found = input_file.find_variable("mask");
355
356 // now work through all the 2d variables, regridding if present and otherwise
357 // setting to default values appropriately
358
359 if (mask_found) {
360 m_log->message(2, " WARNING: 'mask' found; IGNORING IT!\n");
361 }
362
363 if (usurf.exists) {
364 m_log->message(2, " WARNING: surface elevation 'usurf' found; IGNORING IT!\n");
365 }
366
367 m_log->message(2, " reading 2D model state variables by regridding ...\n");
368
369 // longitude
370 {
371 m_geometry.longitude.regrid(input_file, OPTIONAL);
372
373 auto lon = input_file.find_variable("lon", "longitude");
374
375 if (not lon.exists) {
376 m_geometry.longitude.metadata().set_string("missing_at_bootstrap", "true");
377 }
378 }
379
380 // latitude
381 {
382 m_geometry.latitude.regrid(input_file, OPTIONAL);
383
384 auto lat = input_file.find_variable("lat", "latitude");
385
386 if (not lat.exists) {
387 m_geometry.latitude.metadata().set_string("missing_at_bootstrap", "true");
388 }
389 }
390
391 m_geometry.ice_thickness.regrid(input_file, OPTIONAL,
392 m_config->get_number("bootstrapping.defaults.ice_thickness"));
393 // check the range of the ice thickness
394 {
395 Range thk_range = m_geometry.ice_thickness.range();
396
397 if (thk_range.max >= m_grid->Lz() + 1e-6) {
398 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "Maximum ice thickness (%f meters)\n"
399 "exceeds the height of the computational domain (%f meters).",
400 thk_range.max, m_grid->Lz());
401 }
402 }
403
404 if (m_config->get_flag("geometry.part_grid.enabled")) {
405 // Read the Href field from an input file. This field is
406 // grid-dependent, so interpolating it from one grid to a
407 // different one does not make sense in general.
408 // (IceModel::Href_cleanup() will take care of the side effects of
409 // such interpolation, though.)
410 //
411 // On the other hand, we need to read it in to be able to re-start
412 // from a PISM output file using the -bootstrap option.
413 m_geometry.ice_area_specific_volume.regrid(input_file, OPTIONAL, 0.0);
414 }
415
416 if (m_config->get_flag("stress_balance.ssa.dirichlet_bc")) {
417 // Do not use Dirichlet B.C. anywhere if bc_mask is not present.
418 m_ssa_dirichlet_bc_mask.regrid(input_file, OPTIONAL, 0.0);
419 // In the absence of u_ssa_bc and v_ssa_bc in the file the only B.C. that make sense are the
420 // zero Dirichlet B.C.
421 m_ssa_dirichlet_bc_values.regrid(input_file, OPTIONAL, 0.0);
422 } else {
423 m_ssa_dirichlet_bc_mask.set(0.0);
424 m_ssa_dirichlet_bc_values.set(0.0);
425 }
426
427 // check if Lz is valid
428 Range thk_range = m_geometry.ice_thickness.range();
429
430 if (thk_range.max > m_grid->Lz()) {
431 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "Max. ice thickness (%3.3f m)\n"
432 "exceeds the height of the computational domain (%3.3f m).",
433 thk_range.max, m_grid->Lz());
434 }
435 }
436
437 void IceModel::initialize_2d() {
438 // This method should NOT have the "noreturn" attribute. (This attribute does not mix with virtual
439 // methods).
440 throw RuntimeError(PISM_ERROR_LOCATION, "cannot initialize IceModel without an input file");
441 }
442
443 //! Manage regridding based on user options.
444 void IceModel::regrid() {
445
446 auto filename = m_config->get_string("input.regrid.file");
447 auto regrid_vars = set_split(m_config->get_string("input.regrid.vars"), ',');
448
449 // Return if no regridding is requested:
450 if (filename.empty()) {
451 return;
452 }
453
454 m_log->message(2, "regridding from file %s ...\n", filename.c_str());
455
456 {
457 File regrid_file(m_grid->com, filename, PISM_GUESS, PISM_READONLY);
458 for (auto v : m_model_state) {
459 if (regrid_vars.find(v->get_name()) != regrid_vars.end()) {
460 v->regrid(regrid_file, CRITICAL);
461 }
462 }
463
464 // Check the range of the ice thickness.
465 {
466 double
467 max_thickness = m_geometry.ice_thickness.range().max,
468 Lz = m_grid->Lz();
469
470 if (max_thickness >= Lz + 1e-6) {
471 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "Maximum ice thickness (%f meters)\n"
472 "exceeds the height of the computational domain (%f meters).",
473 max_thickness, Lz);
474 }
475 }
476 }
477 }
478
479 //! \brief Decide which stress balance model to use.
480 void IceModel::allocate_stressbalance() {
481
482 if (m_stress_balance) {
483 return;
484 }
485
486 m_log->message(2, "# Allocating a stress balance model...\n");
487
488 // false means "not regional"
489 m_stress_balance = stressbalance::create(m_config->get_string("stress_balance.model"),
490 m_grid, false);
491
492 m_submodels["stress balance"] = m_stress_balance.get();
493 }
494
495 void IceModel::allocate_geometry_evolution() {
496 if (m_geometry_evolution) {
497 return;
498 }
499
500 m_log->message(2, "# Allocating the geometry evolution model...\n");
501
502 m_geometry_evolution.reset(new GeometryEvolution(m_grid));
503
504 m_submodels["geometry_evolution"] = m_geometry_evolution.get();
505 }
506
507
508 void IceModel::allocate_iceberg_remover() {
509
510 if (m_iceberg_remover != NULL) {
511 return;
512 }
513
514 m_log->message(2,
515 "# Allocating an iceberg remover (part of a calving model)...\n");
516
517 if (m_config->get_flag("geometry.remove_icebergs")) {
518
519 // this will throw an exception on failure
520 m_iceberg_remover.reset(new calving::IcebergRemover(m_grid));
521
522 // Iceberg Remover does not have a state, so it is OK to
523 // initialize here.
524 m_iceberg_remover->init();
525
526 m_submodels["iceberg remover"] = m_iceberg_remover.get();
527 }
528 }
529
530 void IceModel::allocate_age_model() {
531
532 if (m_age_model) {
533 return;
534 }
535
536 if (m_config->get_flag("age.enabled")) {
537 m_log->message(2, "# Allocating an ice age model...\n");
538
539 if (m_stress_balance == NULL) {
540 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
541 "Cannot allocate an age model: m_stress_balance == NULL.");
542 }
543
544 m_age_model.reset(new AgeModel(m_grid, m_stress_balance.get()));
545 m_submodels["age model"] = m_age_model.get();
546 }
547 }
548
549 void IceModel::allocate_energy_model() {
550
551 if (m_energy_model != NULL) {
552 return;
553 }
554
555 m_log->message(2, "# Allocating an energy balance model...\n");
556
557 if (m_config->get_flag("energy.enabled")) {
558 if (m_config->get_flag("energy.temperature_based")) {
559 m_energy_model = new energy::TemperatureModel(m_grid, m_stress_balance.get());
560 } else {
561 m_energy_model = new energy::EnthalpyModel(m_grid, m_stress_balance.get());
562 }
563 } else {
564 m_energy_model = new energy::DummyEnergyModel(m_grid, m_stress_balance.get());
565 }
566
567 m_submodels["energy balance model"] = m_energy_model;
568 }
569
570
571 //! \brief Decide which bedrock thermal unit to use.
572 void IceModel::allocate_bedrock_thermal_unit() {
573
574 if (m_btu != NULL) {
575 return;
576 }
577
578 m_log->message(2, "# Allocating a bedrock thermal layer model...\n");
579
580 m_btu = energy::BedThermalUnit::FromOptions(m_grid, m_ctx);
581
582 m_submodels["bedrock thermal model"] = m_btu;
583 }
584
585 //! \brief Decide which subglacial hydrology model to use.
586 void IceModel::allocate_subglacial_hydrology() {
587
588 using namespace pism::hydrology;
589
590 std::string hydrology_model = m_config->get_string("hydrology.model");
591
592 if (m_subglacial_hydrology) { // indicates it has already been allocated
593 return;
594 }
595
596 m_log->message(2, "# Allocating a subglacial hydrology model...\n");
597
598 if (hydrology_model == "null") {
599 m_subglacial_hydrology.reset(new NullTransport(m_grid));
600 } else if (hydrology_model == "routing") {
601 m_subglacial_hydrology.reset(new Routing(m_grid));
602 } else if (hydrology_model == "steady") {
603 m_subglacial_hydrology.reset(new SteadyState(m_grid));
604 } else if (hydrology_model == "distributed") {
605 m_subglacial_hydrology.reset(new Distributed(m_grid));
606 } else {
607 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
608 "unknown 'hydrology.model': %s", hydrology_model.c_str());
609 }
610
611 m_submodels["subglacial hydrology"] = m_subglacial_hydrology.get();
612 }
613
614 //! \brief Decide which basal yield stress model to use.
615 void IceModel::allocate_basal_yield_stress() {
616
617 if (m_basal_yield_stress_model) {
618 return;
619 }
620
621 m_log->message(2,
622 "# Allocating a basal yield stress model...\n");
623
624 std::string model = m_config->get_string("stress_balance.model");
625
626 // only these two use the yield stress (so far):
627 if (model == "ssa" || model == "ssa+sia") {
628 std::string yield_stress_model = m_config->get_string("basal_yield_stress.model");
629
630 if (yield_stress_model == "constant") {
631 m_basal_yield_stress_model.reset(new ConstantYieldStress(m_grid));
632 } else if (yield_stress_model == "mohr_coulomb") {
633 m_basal_yield_stress_model.reset(new MohrCoulombYieldStress(m_grid));
634 } else {
635 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "yield stress model '%s' is not supported.",
636 yield_stress_model.c_str());
637 }
638
639 m_submodels["basal yield stress"] = m_basal_yield_stress_model.get();
640 }
641 }
642
643 //! Allocate PISM's sub-models implementing some physical processes.
644 /*!
645 This method is called after memory allocation but before filling any of
646 IceModelVecs because all the physical parameters should be initialized before
647 setting up the coupling or filling model-state variables.
648 */
649 void IceModel::allocate_submodels() {
650
651 allocate_geometry_evolution();
652
653 allocate_iceberg_remover();
654
655 allocate_stressbalance();
656
657 // this has to happen *after* allocate_stressbalance()
658 {
659 allocate_age_model();
660 allocate_energy_model();
661 allocate_subglacial_hydrology();
662 }
663
664 // this has to happen *after* allocate_subglacial_hydrology()
665 allocate_basal_yield_stress();
666
667 allocate_bedrock_thermal_unit();
668
669 allocate_bed_deformation();
670
671 allocate_couplers();
672
673 if (m_config->get_flag("fracture_density.enabled")) {
674 m_fracture.reset(new FractureDensity(m_grid, m_stress_balance->shallow()->flow_law()));
675 m_submodels["fracture_density"] = m_fracture.get();
676 }
677 }
678
679
680 void IceModel::allocate_couplers() {
681 // Initialize boundary models:
682
683 if (not m_surface) {
684
685 m_log->message(2, "# Allocating a surface process model or coupler...\n");
686
687 surface::Factory ps(m_grid, atmosphere::Factory(m_grid).create());
688
689 m_surface.reset(new surface::InitializationHelper(m_grid, ps.create()));
690
691 m_submodels["surface process model"] = m_surface.get();
692 }
693
694 if (not m_sea_level) {
695 m_log->message(2, "# Allocating sea level forcing...\n");
696
697 using namespace ocean::sea_level;
698
699 m_sea_level.reset(new InitializationHelper(m_grid, Factory(m_grid).create()));
700
701 m_submodels["sea level forcing"] = m_sea_level.get();
702 }
703
704 if (not m_ocean) {
705 m_log->message(2, "# Allocating an ocean model or coupler...\n");
706
707 using namespace ocean;
708
709 m_ocean.reset(new InitializationHelper(m_grid, Factory(m_grid).create()));
710
711 m_submodels["ocean model"] = m_ocean.get();
712 }
713 }
714
715 //! Miscellaneous initialization tasks plus tasks that need the fields that can come from regridding.
716 void IceModel::misc_setup() {
717
718 m_log->message(3, "Finishing initialization...\n");
719 InputOptions opts = process_input_options(m_ctx->com(), m_config);
720
721 if (not (opts.type == INIT_OTHER)) {
722 // initializing from a file
723 File file(m_grid->com, opts.filename, PISM_GUESS, PISM_READONLY);
724
725 std::string source = file.read_text_attribute("PISM_GLOBAL", "source");
726
727 if (opts.type == INIT_RESTART) {
728 // If it's missing, print a warning
729 if (source.empty()) {
730 m_log->message(1,
731 "PISM WARNING: file '%s' does not have the 'source' global attribute.\n"
732 " If '%s' is a PISM output file, please run the following to get rid of this warning:\n"
733 " ncatted -a source,global,c,c,PISM %s\n",
734 opts.filename.c_str(), opts.filename.c_str(), opts.filename.c_str());
735 } else if (source.find("PISM") == std::string::npos) {
736 // If the 'source' attribute does not contain the string "PISM", then print
737 // a message and stop:
738 m_log->message(1,
739 "PISM WARNING: '%s' does not seem to be a PISM output file.\n"
740 " If it is, please make sure that the 'source' global attribute contains the string \"PISM\".\n",
741 opts.filename.c_str());
742 }
743 }
744 }
745
746 {
747 // A single record of a time-dependent variable cannot exceed 2^32-4
748 // bytes in size. See the NetCDF User's Guide
749 // <http://www.unidata.ucar.edu/software/netcdf/docs/netcdf.html#g_t64-bit-Offset-Limitations>.
750 // Here we use "long int" to avoid integer overflow.
751 const long int two_to_thirty_two = 4294967296L;
752 const long int
753 Mx = m_grid->Mx(),
754 My = m_grid->My(),
755 Mz = m_grid->Mz();
756 std::string output_format = m_config->get_string("output.format");
757 if (Mx * My * Mz * sizeof(double) > two_to_thirty_two - 4 and
758 (output_format == PISM_NETCDF3)) {
759 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
760 "The computational grid is too big to fit in a NetCDF-3 file.\n"
761 "Each 3D variable requires %lu Mb.\n"
762 "Please use '-o_format pnetcdf or -o_format netcdf4_parallel to proceed.",
763 Mx * My * Mz * sizeof(double) / (1024 * 1024));
764 }
765 }
766
767 m_output_vars = output_variables(m_config->get_string("output.size"));
768
769 #if (Pism_USE_PROJ==1)
770 {
771 std::string proj_string = m_grid->get_mapping_info().proj;
772 if (not proj_string.empty()) {
773 m_output_vars.insert("lon_bnds");
774 m_output_vars.insert("lat_bnds");
775 }
776 }
777 #endif
778
779 init_calving();
780 init_frontal_melt();
781 init_front_retreat();
782 init_diagnostics();
783 init_snapshots();
784 init_backups();
785 init_timeseries();
786 init_extras();
787
788 // a report on whether PISM-PIK modifications of IceModel are in use
789 {
790 std::vector<std::string> pik_methods;
791 if (m_config->get_flag("geometry.part_grid.enabled")) {
792 pik_methods.push_back("part_grid");
793 }
794 if (m_config->get_flag("geometry.remove_icebergs")) {
795 pik_methods.push_back("kill_icebergs");
796 }
797
798 if (not pik_methods.empty()) {
799 m_log->message(2,
800 "* PISM-PIK mass/geometry methods are in use: %s\n",
801 join(pik_methods, ", ").c_str());
802 }
803 }
804
805 // initialize diagnostics
806 {
807 // reset: this gives diagnostics a chance to capture the current state of the model at the
808 // beginning of the run
809 for (auto d : m_diagnostics) {
810 d.second->reset();
811 }
812
813 // read in the state (accumulators) if we are re-starting a run
814 if (opts.type == INIT_RESTART) {
815 File file(m_grid->com, opts.filename, PISM_GUESS, PISM_READONLY);
816 for (auto d : m_diagnostics) {
817 d.second->init(file, opts.record);
818 }
819 }
820 }
821
822 if (m_surface_input_for_hydrology) {
823 ForcingOptions surface_input(*m_ctx, "hydrology.surface_input");
824 m_surface_input_for_hydrology->init(surface_input.filename,
825 surface_input.period,
826 surface_input.reference_time);
827 }
828
829 if (m_fracture) {
830 if (opts.type == INIT_OTHER) {
831 m_fracture->initialize();
832 } else {
833 // initializing from a file
834 File file(m_grid->com, opts.filename, PISM_GUESS, PISM_READONLY);
835
836 if (opts.type == INIT_RESTART) {
837 m_fracture->restart(file, opts.record);
838 } else if (opts.type == INIT_BOOTSTRAP) {
839 m_fracture->bootstrap(file);
840 }
841 }
842 }
843 }
844
845 void IceModel::init_frontal_melt() {
846
847 auto frontal_melt = m_config->get_string("frontal_melt.models");
848
849 if (not frontal_melt.empty()) {
850 if (not m_config->get_flag("geometry.part_grid.enabled")) {
851 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
852 "ERROR: frontal melt models require geometry.part_grid.enabled");
853 }
854
855 m_frontal_melt = frontalmelt::Factory(m_grid).create(frontal_melt);
856
857 m_frontal_melt->init(m_geometry);
858
859 m_submodels["frontal melt"] = m_frontal_melt.get();
860
861 if (not m_front_retreat) {
862 m_front_retreat.reset(new FrontRetreat(m_grid));
863 }
864 }
865 }
866
867 void IceModel::init_front_retreat() {
868 auto front_retreat_file = m_config->get_string("geometry.front_retreat.prescribed.file");
869
870 if (not front_retreat_file.empty()) {
871 m_prescribed_retreat.reset(new PrescribedRetreat(m_grid));
872
873 m_prescribed_retreat->init();
874
875 m_submodels["prescribed front retreat"] = m_prescribed_retreat.get();
876 }
877 }
878
879 //! \brief Initialize calving mechanisms.
880 void IceModel::init_calving() {
881
882 std::set<std::string> methods = set_split(m_config->get_string("calving.methods"), ',');
883 bool allocate_front_retreat = false;
884
885 if (member("thickness_calving", methods)) {
886
887 if (not m_thickness_threshold_calving) {
888 m_thickness_threshold_calving.reset(new calving::CalvingAtThickness(m_grid));
889 }
890
891 m_thickness_threshold_calving->init();
892 methods.erase("thickness_calving");
893
894 m_submodels["thickness threshold calving"] = m_thickness_threshold_calving.get();
895 }
896
897
898 if (member("eigen_calving", methods)) {
899 allocate_front_retreat = true;
900
901 if (not m_eigen_calving) {
902 m_eigen_calving.reset(new calving::EigenCalving(m_grid));
903 }
904
905 m_eigen_calving->init();
906 methods.erase("eigen_calving");
907
908 m_submodels["eigen calving"] = m_eigen_calving.get();
909 }
910
911 if (member("vonmises_calving", methods)) {
912 allocate_front_retreat = true;
913
914 if (not m_vonmises_calving) {
915 m_vonmises_calving.reset(new calving::vonMisesCalving(m_grid,
916 m_stress_balance->shallow()->flow_law()));
917 }
918
919 m_vonmises_calving->init();
920 methods.erase("vonmises_calving");
921
922 m_submodels["von Mises calving"] = m_vonmises_calving.get();
923 }
924
925 if (member("hayhurst_calving", methods)) {
926 allocate_front_retreat = true;
927
928 if (not m_hayhurst_calving) {
929 m_hayhurst_calving.reset(new calving::HayhurstCalving(m_grid));
930 }
931
932 m_hayhurst_calving->init();
933 methods.erase("hayhurst_calving");
934
935 m_submodels["Hayhurst calving"] = m_hayhurst_calving.get();
936 }
937
938 if (member("float_kill", methods)) {
939 if (not m_float_kill_calving) {
940 m_float_kill_calving.reset(new calving::FloatKill(m_grid));
941 }
942
943 m_float_kill_calving->init();
944 methods.erase("float_kill");
945
946 m_submodels["float kill calving"] = m_float_kill_calving.get();
947 }
948
949 if (not methods.empty()) {
950 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
951 "PISM ERROR: calving method(s) [%s] are not supported.\n",
952 set_join(methods, ",").c_str());
953 }
954
955 // allocate front retreat code if necessary
956 if (not m_front_retreat and allocate_front_retreat) {
957 m_front_retreat.reset(new FrontRetreat(m_grid));
958 }
959 }
960
961 void IceModel::allocate_bed_deformation() {
962 std::string model = m_config->get_string("bed_deformation.model");
963
964 if (m_beddef != NULL) {
965 return;
966 }
967
968 m_log->message(2,
969 "# Allocating a bed deformation model...\n");
970
971 if (model == "none") {
972 m_beddef = new bed::Null(m_grid);
973 }
974 else if (model == "iso") {
975 m_beddef = new bed::PointwiseIsostasy(m_grid);
976 }
977 else if (model == "lc") {
978 m_beddef = new bed::LingleClark(m_grid);
979 }
980
981 m_submodels["bed deformation"] = m_beddef;
982 }
983
984 //! Read some runtime (command line) options and alter the
985 //! corresponding parameters or flags as appropriate.
986 void IceModel::process_options() {
987
988 m_log->message(3,
989 "Processing physics-related command-line options...\n");
990
991 set_config_from_options(*m_config);
992
993 // Set global attributes using the config database:
994 m_output_global_attributes.set_string("title", m_config->get_string("run_info.title"));
995 m_output_global_attributes.set_string("institution", m_config->get_string("run_info.institution"));
996 m_output_global_attributes.set_string("command", args_string());
997
998 // warn about some option combinations
999
1000 if (m_config->get_number("time_stepping.maximum_time_step") <= 0) {
1001 throw RuntimeError(PISM_ERROR_LOCATION, "time_stepping.maximum_time_step has to be greater than 0.");
1002 }
1003
1004 if (not m_config->get_flag("geometry.update.enabled") &&
1005 m_config->get_flag("time_stepping.skip.enabled")) {
1006 m_log->message(2,
1007 "PISM WARNING: Both -skip and -no_mass are set.\n"
1008 " -skip only makes sense in runs updating ice geometry.\n");
1009 }
1010
1011 if (m_config->get_string("calving.methods").find("thickness_calving") != std::string::npos &&
1012 not m_config->get_flag("geometry.part_grid.enabled")) {
1013 m_log->message(2,
1014 "PISM WARNING: Calving at certain terminal ice thickness (-calving thickness_calving)\n"
1015 " without application of partially filled grid cell scheme (-part_grid)\n"
1016 " may lead to (incorrect) non-moving ice shelf front.\n");
1017 }
1018 }
1019
1020 //! Assembles a list of diagnostics corresponding to an output file size.
1021 std::set<std::string> IceModel::output_variables(const std::string &keyword) {
1022
1023 std::string variables;
1024
1025 if (keyword == "none" or
1026 keyword == "small") {
1027 variables = "";
1028 } else if (keyword == "medium") {
1029 variables = m_config->get_string("output.sizes.medium");
1030 } else if (keyword == "big_2d") {
1031 variables = (m_config->get_string("output.sizes.medium") + "," +
1032 m_config->get_string("output.sizes.big_2d"));
1033 } else if (keyword == "big") {
1034 variables = (m_config->get_string("output.sizes.medium") + "," +
1035 m_config->get_string("output.sizes.big_2d") + "," +
1036 m_config->get_string("output.sizes.big"));
1037 }
1038
1039 return set_split(variables, ',');
1040 }
1041
1042 void IceModel::compute_lat_lon() {
1043
1044 std::string projection = m_grid->get_mapping_info().proj;
1045
1046 if (m_config->get_flag("grid.recompute_longitude_and_latitude") and
1047 not projection.empty()) {
1048 m_log->message(2,
1049 "* Computing longitude and latitude using projection parameters...\n");
1050
1051 compute_longitude(projection, m_geometry.longitude);
1052 m_geometry.longitude.metadata().set_string("missing_at_bootstrap", "");
1053 compute_latitude(projection, m_geometry.latitude);
1054 m_geometry.latitude.metadata().set_string("missing_at_bootstrap", "");
1055 }
1056 }
1057
1058 } // end of namespace pism