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