tIBIceModel.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
---
tIBIceModel.cc (18949B)
---
1 #include <iostream>
2
3 #include <pism/energy/BedThermalUnit.hh>
4 #include <pism/util/EnthalpyConverter.hh>
5 #include <pism/util/io/File.hh>
6 #include <pism/util/io/io_helpers.hh>
7 #include "pism/energy/EnergyModel.hh"
8
9 #include "pism/icebin/IBIceModel.hh"
10 #include "pism/icebin/IBSurfaceModel.hh"
11
12 #include "pism/coupler/SeaLevel.hh"
13 #include "pism/coupler/ocean/sea_level/Initialization.hh"
14
15 namespace pism {
16 namespace icebin {
17
18 // ================================
19
20 IBIceModel::IBIceModel(IceGrid::Ptr g, Context::Ptr context, IBIceModel::Params const &_params)
21 : pism::IceModel(g, context), params(_params) {
22 // empty
23 }
24
25 IBIceModel::~IBIceModel() {
26 // empty
27 }
28
29
30 void IBIceModel::allocate_subglacial_hydrology() {
31 printf("BEGIN IBIceModel::allocate_subglacial_hydrology()\n");
32
33 if (pism::IceModel::m_subglacial_hydrology) {
34 return; // indicates it has already been allocated
35 }
36
37 m_subglacial_hydrology.reset(new pism::hydrology::NullTransport(m_grid));
38
39 m_submodels["subglacial hydrology"] = m_subglacial_hydrology.get();
40
41 printf("END IBIceModel::allocate_subglacial_hydrology()\n");
42 }
43
44
45 void IBIceModel::allocate_couplers() {
46 // Initialize boundary models:
47
48 if (not m_surface) {
49
50 m_log->message(2, "# Allocating a surface process model or coupler...\n");
51
52 m_surface.reset(new IBSurfaceModel(m_grid));
53
54 m_submodels["surface process model"] = m_surface.get();
55 }
56
57 if (not m_ocean) {
58 m_log->message(2, "# Allocating an ocean model or coupler...\n");
59
60 ocean::Factory po(m_grid);
61 m_ocean = po.create();
62
63 m_submodels["ocean model"] = m_ocean.get();
64 }
65
66 if (not m_sea_level) {
67 using namespace ocean::sea_level;
68 std::shared_ptr<SeaLevel> sea_level(new SeaLevel(m_grid));
69 m_sea_level.reset(new InitializationHelper(m_grid, sea_level));
70 m_submodels["sea level forcing"] = m_sea_level.get();
71 }
72 }
73
74
75 void IBIceModel::allocate_storage() {
76 super::allocate_storage();
77
78 printf("BEGIN IBIceModel::allocate_storage()\n");
79 base.create(m_grid, "", WITHOUT_GHOSTS);
80 cur.create(m_grid, "", WITHOUT_GHOSTS);
81 rate.create(m_grid, "", WITHOUT_GHOSTS);
82 printf("END IBIceModel::allocate_storage()\n");
83
84 M1.create(m_grid, "M1", pism::WITHOUT_GHOSTS);
85 M2.create(m_grid, "M2", pism::WITHOUT_GHOSTS);
86 M1.create(m_grid, "H1", pism::WITHOUT_GHOSTS);
87 M2.create(m_grid, "H2", pism::WITHOUT_GHOSTS);
88 M1.create(m_grid, "V1", pism::WITHOUT_GHOSTS);
89 M2.create(m_grid, "V2", pism::WITHOUT_GHOSTS);
90
91 std::cout << "IBIceModel Conservation Formulas:" << std::endl;
92 cur.print_formulas(std::cout);
93 }
94
95 void IBIceModel::energy_step() {
96
97 printf("BEGIN IBIceModel::energy_step(t=%f, dt=%f)\n", t_TempAge, dt_TempAge);
98
99 // Enthalpy and mass continuity are stepped with different timesteps.
100 // Fish out the timestep relevant to US.
101 // const double my_t0 = t_TempAge; // Time at beginning of timestep
102 const double my_dt = dt_TempAge;
103
104 // =========== BEFORE Energy Step
105
106 // =========== The Energy Step Itself
107 super::energy_step();
108
109 // =========== AFTER Energy Step
110
111 // We need to integrate over strain_heating and geothermal_flux, which
112 // are given in PISM as rates.
113
114
115 // --------- Upward Geothermal Flux
116 // Use actual geothermal flux, not the long-term average..
117 // See: file:///Users/rpfische/git/pism/build/doc/browser/html/classPISMBedThermalUnit.html#details
118 {
119 cur.upward_geothermal_flux.add(my_dt, m_btu->flux_through_top_surface());
120 }
121
122 // ----------- Geothermal Flux
123 cur.geothermal_flux.add(my_dt, m_btu->flux_through_bottom_surface());
124
125 // ---------- Basal Frictional Heating (see iMenthalpy.cc l. 220)
126 IceModelVec2S const &Rb(m_stress_balance->basal_frictional_heating());
127 cur.basal_frictional_heating.add(my_dt, Rb);
128
129 // NOTE: strain_heating is inf at the coastlines.
130 // See: https://github.com/pism/pism/issues/292
131 // ------------ Volumetric Strain Heating
132 // strain_heating_sum += my_dt * sum_columns(strainheating3p)
133 const IceModelVec3 &strain_heating3(m_stress_balance->volumetric_strain_heating());
134 // cur.strain_heating = cur.strain_heating * 1.0 + my_dt * sum_columns(strain_heating3p)
135 strain_heating3.sumColumns(cur.strain_heating, 1.0, my_dt);
136
137 printf("END IBIceModel::energy_step(time=%f)\n", t_TempAge);
138 }
139
140 void IBIceModel::massContExplicitStep(double dt,
141 const IceModelVec2Stag &diffusive_flux,
142 const IceModelVec2V &advective_velocity) {
143
144 printf("BEGIN IBIceModel::MassContExplicitStep()\n");
145
146 _ice_density = m_config->get_number("constants.ice.density");
147 _meter_per_s_to_kg_per_m2 = dt * _ice_density;
148
149
150 // =========== The Mass Continuity Step Itself
151 // This will call through to accumulateFluxes_massContExplicitStep()
152 // in the inner loop
153 {
154 AccessList access{ &cur.pism_smb, &cur.melt_grounded, &cur.melt_floating,
155 &cur.internal_advection, &cur.href_to_h, &cur.nonneg_rule };
156
157 // FIXME: this is obviously broken now that PISM uses GeometryEvolution instead.
158 (void) diffusive_flux;
159 (void) advective_velocity;
160 // super::massContExplicitStep(dt, diffusive_flux, advective_velocity);
161 }
162
163 // =========== AFTER the Mass Continuity Step
164
165 // ----------- SMB: Pass inputs through to outputs.
166 // They are needed to participate in mass/energy budget
167 IBSurfaceModel *ib_surface = ib_surface_model();
168
169
170 {
171 AccessList access{ &ib_surface->icebin_massxfer, &ib_surface->icebin_enthxfer, &ib_surface->icebin_deltah,
172 &cur.icebin_xfer, &cur.icebin_deltah };
173
174 for (int i = m_grid->xs(); i < m_grid->xs() + m_grid->xm(); ++i) {
175 for (int j = m_grid->ys(); j < m_grid->ys() + m_grid->ym(); ++j) {
176 cur.icebin_xfer.mass(i, j) += m_dt * ib_surface->icebin_massxfer(i, j);
177 cur.icebin_xfer.enth(i, j) += m_dt * ib_surface->icebin_enthxfer(i, j);
178 cur.icebin_deltah(i, j) += m_dt * ib_surface->icebin_deltah(i, j);
179 }
180 }
181 }
182
183 printf("END IBIceModel::MassContExplicitStep()\n");
184 }
185
186
187 /** This is called IMMEDIATELY after ice is gained/lost in
188 iMgeometry.cc (massContExplicitStep()). Here we can record the same
189 values that PISM saw when moving ice around. */
190 void IBIceModel::accumulateFluxes_massContExplicitStep(int i, int j,
191 double surface_mass_balance, // [m s-1] ice equivalent
192 double meltrate_grounded, // [m s-1] ice equivalent
193 double meltrate_floating, // [m s-1] ice equivalent
194 double divQ_SIA, // [m s-1] ice equivalent
195 double divQ_SSA, // [m s-1] ice equivalent
196 double Href_to_H_flux, // [m] ice equivalent
197 double nonneg_rule_flux) // [m] ice equivalent
198 {
199 EnthalpyConverter::Ptr EC = ctx()->enthalpy_converter();
200
201 // -------------- Melting
202 double p_basal = EC->pressure(m_geometry.ice_thickness(i, j));
203 double T = EC->melting_temperature(p_basal);
204 double specific_enth_basal = EC->enthalpy_permissive(T, 1.0, p_basal);
205 double mass;
206
207 // ------- Melting at base of ice sheet
208 mass = -meltrate_grounded * _meter_per_s_to_kg_per_m2;
209 cur.melt_grounded.mass(i, j) += mass;
210 cur.melt_grounded.enth(i, j) += mass * specific_enth_basal;
211
212 // ------- Melting under ice shelf
213 mass = -meltrate_floating * _meter_per_s_to_kg_per_m2;
214 cur.melt_floating.mass(i, j) += mass;
215 cur.melt_floating.enth(i, j) += mass * specific_enth_basal;
216
217
218 // -------------- internal_advection
219 const int ks = m_grid->kBelowHeight(m_geometry.ice_thickness(i, j));
220 // Approximate, we will use the enthalpy of the top layer...
221 double specific_enth_top = m_energy_model->enthalpy().get_column(i, j)[ks];
222
223 mass = -(divQ_SIA + divQ_SSA) * _meter_per_s_to_kg_per_m2;
224
225 cur.internal_advection.mass(i, j) += mass;
226 cur.internal_advection.enth(i, j) += mass * specific_enth_top;
227
228
229 // -------------- Get the easy veriables out of the way...
230 mass = surface_mass_balance * _meter_per_s_to_kg_per_m2;
231 cur.pism_smb.mass(i, j) += mass;
232 cur.pism_smb.enth(i, j) += mass * specific_enth_top;
233 cur.nonneg_rule(i, j) -= nonneg_rule_flux * _ice_density;
234 cur.href_to_h(i, j) += Href_to_H_flux * _ice_density;
235
236
237 // printf("END IBIceModel::accumulateFluxes_MassContExplicitStep()\n");
238 }
239
240
241 void IBIceModel::prepare_nc(std::string const &fname, std::unique_ptr<File> &nc) {
242
243 // nc.reset(new File(m_grid->com, m_grid->ctx()->config()->get_string("output.format")));
244
245 nc.reset(new File(m_grid->com,
246 fname,
247 string_to_backend(m_config->get_string("output.format")),
248 PISM_READWRITE_MOVE,
249 m_grid->ctx()->pio_iosys_id()));
250
251 io::define_time(*nc, m_grid->ctx()->config()->get_string("time.dimension_name"), m_grid->ctx()->time()->calendar(),
252 m_grid->ctx()->time()->CF_units_string(), m_grid->ctx()->unit_system());
253
254 // These are in iMtimseries, but not listed as required in iceModelVec.hh
255 // nc->write_attribute(m_config.get_string("time.dimension_name"),
256 // "bounds", "time_bounds");
257 // write_metadata(nc, true, false);
258 // nc->close():
259 }
260
261 /** @param t0 Time of last time we coupled. */
262 void IBIceModel::set_rate(double dt) {
263
264 printf("BEGIN IBIceModel::set_rate(dt=%f)\n", dt);
265
266 double by_dt = 1.0 / dt;
267
268 compute_enth2(cur.total.enth, cur.total.mass);
269 cur.set_epsilon(m_grid);
270
271 // Compute differences, and set base = cur
272 auto base_ii(base.all_vecs.begin());
273 auto cur_ii(cur.all_vecs.begin());
274 auto rate_ii(rate.all_vecs.begin());
275 for (; base_ii != base.all_vecs.end(); ++base_ii, ++cur_ii, ++rate_ii) {
276 IceModelVec2S &vbase(base_ii->vec);
277 IceModelVec2S &vcur(cur_ii->vec);
278 IceModelVec2S &vrate(rate_ii->vec);
279
280 {
281 AccessList access{ &vbase, &vcur, &vrate };
282 for (int i = m_grid->xs(); i < m_grid->xs() + m_grid->xm(); ++i) {
283 for (int j = m_grid->ys(); j < m_grid->ys() + m_grid->ym(); ++j) {
284 // rate = cur - base: Just for DELTA and EPISLON flagged vectors
285 if (base_ii->flags & (MassEnergyBudget::DELTA | MassEnergyBudget::EPSILON)) {
286 vrate(i, j) = (vcur(i, j) - vbase(i, j)) * by_dt;
287 } else {
288 // Or else just copy the to rate
289 vrate(i, j) = vcur(i, j);
290 }
291 }
292 }
293 }
294 }
295
296 printf("END IBIceModel::set_rate()\n");
297 }
298
299 void IBIceModel::reset_rate() {
300 // Compute differences, and set base = cur
301 auto base_ii(base.all_vecs.begin());
302 auto cur_ii(cur.all_vecs.begin());
303 for (; base_ii != base.all_vecs.end(); ++base_ii, ++cur_ii) {
304 IceModelVec2S &vbase(base_ii->vec);
305 IceModelVec2S &vcur(cur_ii->vec);
306
307 // This cannot go in the loop above with PETSc because
308 // vbase is needed on the RHS of the equations above.
309 AccessList access{ &vbase, &vcur };
310 for (int i = m_grid->xs(); i < m_grid->xs() + m_grid->xm(); ++i) {
311 for (int j = m_grid->ys(); j < m_grid->ys() + m_grid->ym(); ++j) {
312 // base = cur: For ALL vectors
313 vbase(i, j) = vcur(i, j);
314 }
315 }
316 }
317
318
319 #if 0
320 rate.geothermal_flux.begin_access();
321 printf("GG rate.geothermal_flux(%d, %d) = %f (%p)\n", m_grid->xs, m_grid->xs, rate.geothermal_flux(m_grid->xs, m_grid->xs), &rate.geothermal_flux(m_grid->xs, m_grid->xs));
322 rate.geothermal_flux.end_access();
323
324 cur.geothermal_flux.begin_access();
325 printf("GG cur.geothermal_flux(%d, %d) = %f (%p)\n", m_grid->xs, m_grid->xs, cur.geothermal_flux(m_grid->xs, m_grid->xs), &cur.geothermal_flux(m_grid->xs, m_grid->xs));
326 cur.geothermal_flux.end_access();
327
328 base.geothermal_flux.begin_access();
329 printf("GG base.geothermal_flux(%d, %d) = %f (%p)\n", m_grid->xs, m_grid->xs, base.geothermal_flux(m_grid->xs, m_grid->xs), &base.geothermal_flux(m_grid->xs, m_grid->xs));
330 base.geothermal_flux.end_access();
331 #endif
332 }
333
334 /** @param t0 Time of last time we coupled. */
335 void IBIceModel::prepare_outputs(double t0) {
336 printf("BEGIN IBIceModel::prepare_outputs()\n");
337
338 // ------ Difference between now and the last time we were called
339 double t1 = enthalpy_t(); // Current time of the enthalpy portion of ice model.
340 set_rate(t1 - t0);
341
342 // ice_surface_enth & ice_surfac_enth_depth
343 prepare_initial_outputs();
344
345 printf("END IBIceModel::prepare_outputs()\n");
346 }
347
348 void IBIceModel::prepare_initial_outputs() {
349 double ice_density = m_config->get_number("constants.ice.density", "kg m-3");
350
351 const IceModelVec3 &ice_enthalpy = m_energy_model->enthalpy();
352
353 AccessList access{ &ice_enthalpy, &M1, &M2, &H1, &H2, &V1, &V2, &m_geometry.ice_thickness };
354 for (int i = m_grid->xs(); i < m_grid->xs() + m_grid->xm(); ++i) {
355 for (int j = m_grid->ys(); j < m_grid->ys() + m_grid->ym(); ++j) {
356 double const *Enth = ice_enthalpy.get_column(i, j);
357
358 // Top Layer
359 int const ks = m_grid->kBelowHeight(m_geometry.ice_thickness(i, j));
360 V1(i, j) = m_geometry.ice_thickness(i, j) - m_grid->z(ks); // [m^3 m-2]
361 M1(i, j) = V1(i, j) * ice_density; // [kg m-2] = [m^3 m-2] [kg m-3]
362 H1(i, j) = Enth[ks] * M1(i, j); // [J m-2] = [J kg-1] [kg m-2]
363
364 // Second layer
365 int const ks2 = ks - 1;
366 if (ks2 >= 0) {
367 V2(i, j) = m_grid->z(ks) - m_grid->z(ks2);
368 M2(i, j) = V2(i, j) * ice_density;
369 H2(i, j) = Enth[ks2] * M2(i, j);
370 } else {
371 // There is no second layer
372 V2(i, j) = 0;
373 M2(i, j) = 0;
374 H2(i, j) = 0;
375 }
376 }
377 }
378
379 // ====================== Write to the post_energy.nc file (OPTIONAL)
380 }
381
382 void IBIceModel::time_setup() {
383 // super::m_grid_setup() trashes m_time->start(). Now set it correctly.
384 m_time->set_start(params.time_start_s);
385 m_time->set(params.time_start_s);
386
387 m_log->message(2, "* Run time: [%s, %s] (%s years, using the '%s' calendar)\n", m_time->start_date().c_str(),
388 m_time->end_date().c_str(), m_time->run_length().c_str(), m_time->calendar().c_str());
389 }
390
391
392 void IBIceModel::misc_setup() {
393 super::misc_setup();
394
395
396 // ------ Initialize MassEnth structures: base, cur, rate
397 for (auto &ii : cur.all_vecs) {
398 ii.vec.set(0);
399 }
400 compute_enth2(cur.total.enth, cur.total.mass);
401 cur.set_epsilon(m_grid);
402
403 // base = cur
404 auto base_ii(base.all_vecs.begin());
405 auto cur_ii(cur.all_vecs.begin());
406 for (; base_ii != base.all_vecs.end(); ++base_ii, ++cur_ii) {
407 base_ii->vec.copy_from(cur_ii->vec);
408 }
409 }
410
411 /** Sums over columns to compute enthalpy on 2D m_grid->
412
413 NOTE: Unfortunately so far PISM does not keep track of enthalpy in
414 "partially-filled" cells, so Enth2(i,j) is not valid at locations like
415 this one. We need to address this, but so far, it seems to me, the
416 best thing we can do is estimate Enth2(i,j) at partially-filled cells
417 by computing the average over icy neighbors. I think you can re-use
418 the idea from IceModel::get_threshold_thickness(...) (iMpartm_grid->cc). */
419
420
421 void IBIceModel::compute_enth2(pism::IceModelVec2S &enth2, pism::IceModelVec2S &mass2) {
422 // getInternalColumn() is allocated already
423 double ice_density = m_config->get_number("constants.ice.density", "kg m-3");
424
425 const IceModelVec3 *ice_enthalpy = &m_energy_model->enthalpy();
426
427 AccessList access{ &m_geometry.ice_thickness, ice_enthalpy, &enth2, &mass2 };
428 for (int i = m_grid->xs(); i < m_grid->xs() + m_grid->xm(); ++i) {
429 for (int j = m_grid->ys(); j < m_grid->ys() + m_grid->ym(); ++j) {
430 enth2(i, j) = 0;
431 mass2(i, j) = 0;
432
433 // count all ice, including cells that have so little they
434 // are considered "ice-free"
435 if (m_geometry.ice_thickness(i, j) > 0) {
436 const int ks = m_grid->kBelowHeight(m_geometry.ice_thickness(i, j));
437 double const *Enth = ice_enthalpy->get_column(i, j);
438 for (int k = 0; k < ks; ++k) {
439 double dz = (m_grid->z(k + 1) - m_grid->z(k));
440 enth2(i, j) += Enth[k] * dz; // m J / kg
441 }
442
443 // Do the last layer a bit differently
444 double dz = (m_geometry.ice_thickness(i, j) - m_grid->z(ks));
445 enth2(i, j) += Enth[ks] * dz;
446 enth2(i, j) *= ice_density; // --> J/m^2
447 mass2(i, j) = m_geometry.ice_thickness(i, j) * ice_density; // --> kg/m^2
448 }
449 }
450 }
451 }
452
453
454 /** Merges surface temperature derived from the energy balance model into any NaN values
455 in the vector provided.
456 @param deltah IN: Input from Icebin (change in enthalpy of each m_grid
457 cell over the timestep) [W m-2].
458 @param default_val: The value that deltah(i,j) will have if no value
459 is listed for that m_grid cell
460 @param timestep_s: Length of the current coupling timestep [s]
461 @param surface_temp OUT: Resulting surface temperature to use as the Dirichlet B.C.
462 */
463 void IBIceModel::construct_surface_temp(
464 pism::IceModelVec2S &deltah, // IN: Input from Icebin
465 double default_val,
466 double timestep_s, // Length of this coupling interval [s]
467 pism::IceModelVec2S &surface_temp) // OUT: Temperature @ top of ice sheet (to use for Dirichlet B.C.)
468
469 {
470 printf("BEGIN IBIceModel::merge_surface_temp default_val=%g\n", default_val);
471 EnthalpyConverter::Ptr EC = ctx()->enthalpy_converter();
472
473 double ice_density = m_config->get_number("constants.ice.density");
474
475 const IceModelVec3 &ice_enthalpy = m_energy_model->enthalpy();
476
477 {
478 AccessList access{ &ice_enthalpy, &deltah, &m_geometry.ice_thickness, &surface_temp };
479
480 // First time around, set effective_surface_temp to top temperature
481 for (int i = m_grid->xs(); i < m_grid->xs() + m_grid->xm(); ++i) {
482 for (int j = m_grid->ys(); j < m_grid->ys() + m_grid->ym(); ++j) {
483 double &surface_temp_ij(surface_temp(i, j));
484 double const &deltah_ij(deltah(i, j));
485
486 double const *Enth = ice_enthalpy.get_column(i, j);
487
488 // Enthalpy at top of ice sheet
489 const int ks = m_grid->kBelowHeight(m_geometry.ice_thickness(i, j));
490 double spec_enth3 = Enth[ks]; // Specific enthalpy [J kg-1]
491
492 if (deltah_ij != default_val) {
493 // Adjust enthalpy @top by deltah
494 double toplayer_dz = m_geometry.ice_thickness(i, j) - m_grid->z(ks); // [m]
495
496 // [J kg-1] = [J kg-1]
497 // + [J m-2 s-1] * [m^2 m-3] * [m^3 kg-1] * [s]
498 spec_enth3 = spec_enth3 + deltah_ij / (toplayer_dz * ice_density) * timestep_s;
499 }
500
501
502 // Convert specific enthalpy value to surface temperature
503 const double p = 0.0; // Pressure at top of ice sheet
504 surface_temp_ij = EC->temperature(spec_enth3, p); // [K]
505 }
506 }
507 }
508
509 printf("END IBIceModel::merge_surface_temp\n");
510 }
511 }
512 } // namespace icebin::gpism