tGeometryEvolution.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
---
tGeometryEvolution.cc (49797B)
---
1 /* Copyright (C) 2016, 2017, 2018, 2019 PISM Authors
2 *
3 * This file is part of PISM.
4 *
5 * PISM is free software; you can redistribute it and/or modify it under the
6 * terms of the GNU General Public License as published by the Free Software
7 * Foundation; either version 3 of the License, or (at your option) any later
8 * version.
9 *
10 * PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 * details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with PISM; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18 */
19
20 #include "GeometryEvolution.hh"
21
22 #include "pism/util/iceModelVec.hh"
23 #include "pism/util/IceGrid.hh"
24 #include "pism/util/Mask.hh"
25
26 #include "pism/geometry/part_grid_threshold_thickness.hh"
27 #include "pism/util/pism_utilities.hh"
28 #include "pism/util/Logger.hh"
29 #include "pism/util/Profiling.hh"
30
31 namespace pism {
32
33 using mask::floating_ice;
34 using mask::grounded_ice;
35 using mask::ice_free;
36 using mask::ice_free_land;
37 using mask::ice_free_ocean;
38 using mask::icy;
39
40 struct GeometryEvolution::Impl {
41 Impl(IceGrid::ConstPtr g);
42
43 const Profiling &profile;
44
45 GeometryCalculator gc;
46
47 double ice_density;
48
49 //! True if the basal melt rate contributes to geometry evolution.
50 bool use_bmr;
51
52 //! True if the part-grid scheme is enabled.
53 bool use_part_grid;
54
55 //! Flux divergence (used to track thickness changes due to flow).
56 IceModelVec2S flux_divergence;
57
58 //! Conservation error due to enforcing non-negativity of ice thickness.
59 IceModelVec2S conservation_error;
60
61 //! Effective surface mass balance.
62 IceModelVec2S effective_SMB;
63
64 //! Effective basal mass balance.
65 IceModelVec2S effective_BMB;
66
67 //! Change in ice thickness due to flow during the last time step.
68 IceModelVec2S thickness_change;
69
70 //! Change in the ice area-specific volume due to flow during the last time step.
71 IceModelVec2S ice_area_specific_volume_change;
72
73 //! Flux through cell interfaces. Ghosted.
74 IceModelVec2Stag flux_staggered;
75
76 // Work space
77 IceModelVec2V input_velocity; // ghosted copy; not modified
78 IceModelVec2S bed_elevation; // ghosted copy; not modified
79 IceModelVec2S sea_level; // ghosted copy; not modified
80 IceModelVec2S ice_thickness; // ghosted; updated in place
81 IceModelVec2S area_specific_volume; // ghosted; updated in place
82 IceModelVec2S surface_elevation; // ghosted; updated to maintain consistency
83 IceModelVec2CellType cell_type; // ghosted; updated to maintain consistency
84 IceModelVec2S residual; // ghosted; temporary storage
85 IceModelVec2S thickness; // ghosted; temporary storage
86 IceModelVec2Int velocity_bc_mask;
87 };
88
89 GeometryEvolution::Impl::Impl(IceGrid::ConstPtr grid)
90 : profile(grid->ctx()->profiling()),
91 gc(*grid->ctx()->config()),
92 flux_divergence(grid, "flux_divergence", WITHOUT_GHOSTS),
93 conservation_error(grid, "conservation_error", WITHOUT_GHOSTS),
94 effective_SMB(grid, "effective_SMB", WITHOUT_GHOSTS),
95 effective_BMB(grid, "effective_BMB", WITHOUT_GHOSTS),
96 thickness_change(grid, "thickness_change", WITHOUT_GHOSTS),
97 ice_area_specific_volume_change(grid, "ice_area_specific_volume_change", WITHOUT_GHOSTS),
98 flux_staggered(grid, "flux_staggered", WITH_GHOSTS),
99 input_velocity(grid, "input_velocity", WITH_GHOSTS),
100 bed_elevation(grid, "bed_elevation", WITH_GHOSTS),
101 sea_level(grid, "sea_level", WITH_GHOSTS),
102 ice_thickness(grid, "ice_thickness", WITH_GHOSTS),
103 area_specific_volume(grid, "area_specific_volume", WITH_GHOSTS),
104 surface_elevation(grid, "surface_elevation", WITH_GHOSTS),
105 cell_type(grid, "cell_type", WITH_GHOSTS),
106 residual(grid, "residual", WITH_GHOSTS),
107 thickness(grid, "thickness", WITH_GHOSTS),
108 velocity_bc_mask(grid, "velocity_bc_mask", WITH_GHOSTS) {
109
110 Config::ConstPtr config = grid->ctx()->config();
111
112 gc.set_icefree_thickness(config->get_number("geometry.ice_free_thickness_standard"));
113
114 // constants
115 {
116 ice_density = config->get_number("constants.ice.density");
117 use_bmr = config->get_flag("geometry.update.use_basal_melt_rate");
118 use_part_grid = config->get_flag("geometry.part_grid.enabled");
119 }
120
121 // reported quantities
122 {
123 // This is the only reported field that is ghosted (we need ghosts to compute flux divergence).
124 flux_staggered.set_attrs("diagnostic", "fluxes through cell interfaces (sides)"
125 " on the staggered grid",
126 "m2 s-1", "m2 year-1", "", 0);
127
128 flux_divergence.set_attrs("diagnostic", "flux divergence", "m s-1", "m year-1", "", 0);
129
130 conservation_error.set_attrs("diagnostic",
131 "conservation error due to enforcing non-negativity of"
132 " ice thickness (over the last time step)",
133 "meters", "meters", "", 0);
134
135 effective_SMB.set_attrs("internal", "effective surface mass balance over the last time step",
136 "meters", "meters", "", 0);
137
138 effective_BMB.set_attrs("internal", "effective basal mass balance over the last time step",
139 "meters", "meters", "", 0);
140
141 thickness_change.set_attrs("internal", "change in thickness due to flow",
142 "meters", "meters", "", 0);
143
144 ice_area_specific_volume_change.set_attrs("interval",
145 "change in area-specific volume due to flow",
146 "meters3 / meters2", "meters3 / meters2", "", 0);
147 }
148
149 // internal storage
150 {
151 input_velocity.set_attrs("internal", "ghosted copy of the input velocity",
152 "meters / second", "meters / second", "", 0);
153
154 bed_elevation.set_attrs("internal", "ghosted copy of the bed elevation",
155 "meters", "meters", "", 0);
156
157 sea_level.set_attrs("internal", "ghosted copy of the sea level elevation",
158 "meters", "meters", "", 0);
159
160 ice_thickness.set_attrs("internal", "working (ghosted) copy of the ice thickness",
161 "meters", "meters", "", 0);
162
163 area_specific_volume.set_attrs("internal", "working (ghosted) copy of the area specific volume",
164 "meters3 / meters2", "meters3 / meters2", "", 0);
165
166 surface_elevation.set_attrs("internal", "working (ghosted) copy of the surface elevation",
167 "meters", "meters", "", 0);
168
169 cell_type.set_attrs("internal", "working (ghosted) copy of the cell type mask",
170 "", "", "", 0);
171
172 residual.set_attrs("internal", "residual area specific volume",
173 "meters3 / meters2", "meters3 / meters2", "", 0);
174
175 thickness.set_attrs("internal", "thickness (temporary storage)",
176 "meters", "meters", "", 0);
177
178 velocity_bc_mask.set_attrs("internal", "ghosted copy of the velocity B.C. mask"
179 " (1 at velocity B.C. location, 0 elsewhere)",
180 "", "", "", 0);
181 }
182 }
183
184 GeometryEvolution::GeometryEvolution(IceGrid::ConstPtr grid)
185 : Component(grid) {
186 m_impl = new Impl(grid);
187 }
188
189 GeometryEvolution::~GeometryEvolution() {
190 delete m_impl;
191 }
192
193 void GeometryEvolution::init(const InputOptions &opts) {
194 this->init_impl(opts);
195 }
196
197 void GeometryEvolution::init_impl(const InputOptions &opts) {
198 (void) opts;
199 // empty: the default implementation has no state
200 }
201
202 const IceModelVec2S& GeometryEvolution::flux_divergence() const {
203 return m_impl->flux_divergence;
204 }
205
206 const IceModelVec2Stag& GeometryEvolution::flux_staggered() const {
207 return m_impl->flux_staggered;
208 }
209
210 const IceModelVec2S& GeometryEvolution::top_surface_mass_balance() const {
211 return m_impl->effective_SMB;
212 }
213
214 const IceModelVec2S& GeometryEvolution::bottom_surface_mass_balance() const {
215 return m_impl->effective_BMB;
216 }
217
218 const IceModelVec2S& GeometryEvolution::thickness_change_due_to_flow() const {
219 return m_impl->thickness_change;
220 }
221
222 const IceModelVec2S& GeometryEvolution::area_specific_volume_change_due_to_flow() const {
223 return m_impl->ice_area_specific_volume_change;
224 }
225
226 const IceModelVec2S& GeometryEvolution::conservation_error() const {
227 return m_impl->conservation_error;
228 }
229
230 /*!
231 * @param[in] geometry ice geometry
232 * @param[in] dt time step, seconds
233 * @param[in] advective_velocity advective (SSA) velocity
234 * @param[in] diffusive_flux diffusive (SIA) flux
235 * @param[in] velocity_bc_mask advective velocity Dirichlet B.C. mask
236 * @param[in] velocity_bc_values advective velocity Dirichlet B.C. values
237 * @param[in] thickness_bc_mask ice thickness Dirichlet B.C. mask
238 *
239 * Results are stored in internal fields accessible using getters.
240 */
241 void GeometryEvolution::flow_step(const Geometry &geometry, double dt,
242 const IceModelVec2V &advective_velocity,
243 const IceModelVec2Stag &diffusive_flux,
244 const IceModelVec2Int &velocity_bc_mask,
245 const IceModelVec2Int &thickness_bc_mask) {
246
247 m_impl->profile.begin("ge.update_ghosted_copies");
248 {
249 // make ghosted copies of input fields
250 m_impl->ice_thickness.copy_from(geometry.ice_thickness);
251 m_impl->area_specific_volume.copy_from(geometry.ice_area_specific_volume);
252 m_impl->sea_level.copy_from(geometry.sea_level_elevation);
253 m_impl->bed_elevation.copy_from(geometry.bed_elevation);
254 m_impl->input_velocity.copy_from(advective_velocity);
255 m_impl->velocity_bc_mask.copy_from(velocity_bc_mask);
256
257 // Compute cell_type and surface_elevation. Ghosts of results are updated.
258 m_impl->gc.compute(m_impl->sea_level, // in (uses ghosts)
259 m_impl->bed_elevation, // in (uses ghosts)
260 m_impl->ice_thickness, // in (uses ghosts)
261 m_impl->cell_type, // out (ghosts are updated)
262 m_impl->surface_elevation); // out (ghosts are updated)
263 }
264 m_impl->profile.end("ge.update_ghosted_copies");
265
266 // Derived classes can include modifications for regional runs.
267 m_impl->profile.begin("ge.interface_fluxes");
268 compute_interface_fluxes(m_impl->cell_type, // in (uses ghosts)
269 m_impl->ice_thickness, // in (uses ghosts)
270 m_impl->input_velocity, // in (uses ghosts)
271 m_impl->velocity_bc_mask, // in (uses ghosts)
272 diffusive_flux, // in
273 m_impl->flux_staggered); // out
274 m_impl->profile.end("ge.interface_fluxes");
275
276 m_impl->flux_staggered.update_ghosts();
277
278 m_impl->profile.begin("ge.flux_divergence");
279 compute_flux_divergence(m_impl->flux_staggered, // in (uses ghosts)
280 thickness_bc_mask, // in
281 m_impl->flux_divergence); // out
282 m_impl->profile.end("ge.flux_divergence");
283
284 // This is where part_grid is implemented.
285 m_impl->profile.begin("ge.update_in_place");
286 update_in_place(dt, // in
287 m_impl->bed_elevation, // in
288 m_impl->sea_level, // in
289 m_impl->flux_divergence, // in
290 m_impl->ice_thickness, // in/out
291 m_impl->area_specific_volume); // in/out
292 m_impl->profile.end("ge.update_in_place");
293
294 // Compute ice thickness and area specific volume changes.
295 m_impl->profile.begin("ge.compute_changes");
296 {
297 m_impl->ice_thickness.add(-1.0, geometry.ice_thickness,
298 m_impl->thickness_change);
299 m_impl->area_specific_volume.add(-1.0, geometry.ice_area_specific_volume,
300 m_impl->ice_area_specific_volume_change);
301 }
302 m_impl->profile.end("ge.compute_changes");
303
304 // Computes the numerical conservation error and corrects ice_thickness_change and
305 // ice_area_specific_volume_change. We can do this here because
306 // compute_surface_and_basal_mass_balance() preserves non-negativity.
307 //
308 // Note that here we use the "old" ice geometry.
309 //
310 // This computation is purely local.
311 m_impl->profile.begin("ge.ensure_nonnegativity");
312 ensure_nonnegativity(geometry.ice_thickness, // in
313 geometry.ice_area_specific_volume, // in
314 m_impl->thickness_change, // in/out
315 m_impl->ice_area_specific_volume_change, // in/out
316 m_impl->conservation_error); // out
317 m_impl->profile.end("ge.ensure_nonnegativity");
318
319 // Now the caller can compute
320 //
321 // H_new = H_old + thickness_change
322 // Href_new = Href_old + ice_area_specific_volume_change.
323
324 // calving is a separate issue
325 }
326
327 void GeometryEvolution::source_term_step(const Geometry &geometry, double dt,
328 const IceModelVec2Int &thickness_bc_mask,
329 const IceModelVec2S &surface_mass_balance_rate,
330 const IceModelVec2S &basal_melt_rate) {
331
332 m_impl->profile.begin("ge.source_terms");
333 compute_surface_and_basal_mass_balance(dt, // in
334 thickness_bc_mask, // in
335 geometry.ice_thickness, // in
336 geometry.cell_type, // in
337 surface_mass_balance_rate, // in
338 basal_melt_rate, // in
339 m_impl->effective_SMB, // out
340 m_impl->effective_BMB); // out
341 m_impl->profile.end("ge.source_terms");
342
343 }
344
345 /*!
346 * Apply changes due to flow to ice geometry and ice area specific volume.
347 */
348 void GeometryEvolution::apply_flux_divergence(Geometry &geometry) const {
349 geometry.ice_thickness.add(1.0, m_impl->thickness_change);
350 geometry.ice_area_specific_volume.add(1.0, m_impl->ice_area_specific_volume_change);
351 }
352
353 /*!
354 * Update geometry by applying changes due to surface and basal mass fluxes.
355 *
356 * Note: This method performs these changes in the same order as the code ensuring
357 * non-negativity. This is important.
358 */
359 void GeometryEvolution::apply_mass_fluxes(Geometry &geometry) const {
360
361 const IceModelVec2S
362 &dH_SMB = top_surface_mass_balance(),
363 &dH_BMB = bottom_surface_mass_balance();
364 IceModelVec2S &H = geometry.ice_thickness;
365
366 IceModelVec::AccessList list{&H, &dH_SMB, &dH_BMB};
367 ParallelSection loop(m_grid->com);
368 try {
369 for (Points p(*m_grid); p; p.next()) {
370 const int i = p.i(), j = p.j();
371
372 // To preserve non-negativity of thickness we need to apply changes in this exact order.
373 // (Recall that floating-point arithmetic is not associative.)
374 const double H_new = (H(i, j) + dH_SMB(i, j)) + dH_BMB(i, j);
375
376 #if (Pism_DEBUG==1)
377 if (H_new < 0.0) {
378 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "H = %f (negative) at i=%d, j=%d",
379 H_new, i, j);
380 }
381 #endif
382
383 H(i, j) = H_new;
384 }
385 } catch (...) {
386 loop.failed();
387 }
388 loop.check();
389 }
390
391
392 /*!
393 * Prevent advective ice flow from floating ice to ice-free land, as well as in the ice-free areas.
394 */
395 static double limit_advective_velocity(int current, int neighbor, double velocity) {
396
397 // Case 1: Flow between grounded_ice and grounded_ice.
398 if (grounded_ice(current) and grounded_ice(neighbor)) {
399 return velocity;
400 }
401
402 // Cases 2 and 3: Flow between grounded_ice and floating_ice.
403 if ((grounded_ice(current) and floating_ice(neighbor)) or
404 (floating_ice(current) and grounded_ice(neighbor))) {
405 return velocity;
406 }
407
408 // Cases 4 and 5: Flow between grounded_ice and ice_free_land.
409 if ((grounded_ice(current) and ice_free_land(neighbor)) or
410 (ice_free_land(current) and grounded_ice(neighbor))) {
411 return velocity;
412 }
413
414 // Cases 6 and 7: Flow between grounded_ice and ice_free_ocean.
415 if ((grounded_ice(current) and ice_free_ocean(neighbor)) or
416 (ice_free_ocean(current) and grounded_ice(neighbor))) {
417 return velocity;
418 }
419
420 // Case 8: Flow between floating_ice and floating_ice.
421 if (floating_ice(current) and floating_ice(neighbor)) {
422 return velocity;
423 }
424
425 // Cases 9 and 10: Flow between floating_ice and ice_free_land.
426 if ((floating_ice(current) and ice_free_land(neighbor)) or
427 (ice_free_land(current) and floating_ice(neighbor))) {
428 // Disable all flow. This ensures that an ice shelf does not climb up a cliff.
429 return 0.0;
430 }
431
432 // Cases 11 and 12: Flow between floating_ice and ice_free_ocean.
433 if ((floating_ice(current) and ice_free_ocean(neighbor)) or
434 (ice_free_ocean(current) and floating_ice(neighbor))) {
435 return velocity;
436 }
437
438 // Case 13: Flow between ice_free_land and ice_free_land.
439 if (ice_free_land(current) and ice_free_land(neighbor)) {
440 return 0.0;
441 }
442
443 // Cases 14 and 15: Flow between ice_free_land and ice_free_ocean.
444 if ((ice_free_land(current) and ice_free_ocean(neighbor)) or
445 (ice_free_ocean(current) and ice_free_land(neighbor))) {
446 return 0.0;
447 }
448
449 // Case 16: Flow between ice_free_ocean and ice_free_ocean.
450 if (ice_free_ocean(current) and ice_free_ocean(neighbor)) {
451 return 0.0;
452 }
453
454 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
455 "cannot handle the case current=%d, neighbor=%d",
456 current, neighbor);
457 }
458
459 /*!
460 * Prevent SIA-driven flow in ice shelves and ice-free areas.
461 */
462 static double limit_diffusive_flux(int current, int neighbor, double flux) {
463
464 // Case 1: Flow between grounded_ice and grounded_ice.
465 if (grounded_ice(current) and grounded_ice(neighbor)) {
466 return flux;
467 }
468
469 // Cases 2 and 3: Flow between grounded_ice and floating_ice.
470 if ((grounded_ice(current) and floating_ice(neighbor)) or
471 (floating_ice(current) and grounded_ice(neighbor))) {
472 return flux;
473 }
474
475 // Cases 4 and 5: Flow between grounded_ice and ice_free_land.
476 if ((grounded_ice(current) and ice_free_land(neighbor)) or
477 (ice_free_land(current) and grounded_ice(neighbor))) {
478 return flux;
479 }
480
481 // Cases 6 and 7: Flow between grounded_ice and ice_free_ocean.
482 if ((grounded_ice(current) and ice_free_ocean(neighbor)) or
483 (ice_free_ocean(current) and grounded_ice(neighbor))) {
484 return flux;
485 }
486
487 // Case 8: Flow between floating_ice and floating_ice.
488 if (floating_ice(current) and floating_ice(neighbor)) {
489 // no diffusive flux in ice shelves
490 return 0.0;
491 }
492
493 // Cases 9 and 10: Flow between floating_ice and ice_free_land.
494 if ((floating_ice(current) and ice_free_land(neighbor)) or
495 (ice_free_land(current) and floating_ice(neighbor))) {
496 // Disable all flow. This ensures that an ice shelf does not climb up a cliff.
497 return 0.0;
498 }
499
500 // Cases 11 and 12: Flow between floating_ice and ice_free_ocean.
501 if ((floating_ice(current) and ice_free_ocean(neighbor)) or
502 (ice_free_ocean(current) and floating_ice(neighbor))) {
503 return 0.0;
504 }
505
506 // Case 13: Flow between ice_free_land and ice_free_land.
507 if (ice_free_land(current) and ice_free_land(neighbor)) {
508 return 0.0;
509 }
510
511 // Cases 14 and 15: Flow between ice_free_land and ice_free_ocean.
512 if ((ice_free_land(current) and ice_free_ocean(neighbor)) or
513 (ice_free_ocean(current) and ice_free_land(neighbor))) {
514 return 0.0;
515 }
516
517 // Case 16: Flow between ice_free_ocean and ice_free_ocean.
518 if (ice_free_ocean(current) and ice_free_ocean(neighbor)) {
519 return 0.0;
520 }
521
522 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
523 "cannot handle the case current=%d, neighbor=%d",
524 current, neighbor);
525 }
526
527 /*!
528 * Combine advective velocity and the diffusive flux on the staggered grid with the ice thickness to
529 * compute the total flux through cell interfaces.
530 *
531 * Uses first-order upwinding to compute the advective flux.
532 *
533 * Limits the diffusive flux to prevent SIA-driven flow in the ocean and ice-free areas.
534 */
535 void GeometryEvolution::compute_interface_fluxes(const IceModelVec2CellType &cell_type,
536 const IceModelVec2S &ice_thickness,
537 const IceModelVec2V &velocity,
538 const IceModelVec2Int &velocity_bc_mask,
539 const IceModelVec2Stag &diffusive_flux,
540 IceModelVec2Stag &output) {
541
542 IceModelVec::AccessList list{&cell_type, &velocity, &velocity_bc_mask, &ice_thickness,
543 &diffusive_flux, &output};
544
545 ParallelSection loop(m_grid->com);
546 try {
547 for (Points p(*m_grid); p; p.next()) {
548 const int
549 i = p.i(),
550 j = p.j(),
551 M = cell_type(i, j),
552 BC = velocity_bc_mask.as_int(i, j);
553
554 const double H = ice_thickness(i, j);
555 const Vector2 V = velocity(i, j);
556
557 for (int n = 0; n < 2; ++n) {
558 const int
559 oi = 1 - n, // offset in the i direction
560 oj = n, // offset in the j direction
561 i_n = i + oi, // i index of a neighbor
562 j_n = j + oj; // j index of a neighbor
563
564 const int M_n = cell_type(i_n, j_n);
565
566 // advective velocity at the current interface
567 double v = 0.0;
568 {
569 const Vector2 V_n = velocity(i_n, j_n);
570
571 // Regular case
572 {
573 if (icy(M) and icy(M_n)) {
574 // Case 1: both sides of the interface are icy
575 v = (n == 0 ? 0.5 * (V.u + V_n.u) : 0.5 * (V.v + V_n.v));
576
577 } else if (icy(M) and ice_free(M_n)) {
578 // Case 2: icy cell next to an ice-free cell
579 v = (n == 0 ? V.u : V.v);
580
581 } else if (ice_free(M) and icy(M_n)) {
582 // Case 3: ice-free cell next to icy cell
583 v = (n == 0 ? V_n.u : V_n.v);
584
585 } else if (ice_free(M) and ice_free(M_n)) {
586 // Case 4: both sides of the interface are ice-free
587 v = 0.0;
588
589 }
590 }
591
592 // The Dirichlet B.C. case:
593 {
594 const int BC_n = velocity_bc_mask.as_int(i_n, j_n);
595
596 if (BC == 1 and BC_n == 1) {
597 // Case 1: both sides of the interface are B.C. locations: average from
598 // the regular grid onto the staggered grid.
599 v = (n == 0 ? 0.5 * (V.u + V_n.u) : 0.5 * (V.v + V_n.v));
600
601 } else if (BC == 1 and BC_n == 0) {
602 // Case 2: at a Dirichlet B.C. location next to a regular location
603 v = (n == 0 ? V.u : V.v);
604
605 } else if (BC == 0 and BC_n == 1) {
606
607 // Case 3: at a regular location next to a Dirichlet B.C. location
608 v = (n == 0 ? V_n.u : V_n.v);
609
610 } else {
611 // Case 4: elsewhere.
612 // No Dirichlet B.C. adjustment here.
613 }
614
615 } // end of the Dirichlet B.C. case
616
617 // finally, limit advective velocities
618 v = limit_advective_velocity(M, M_n, v);
619 }
620
621 // advective flux
622 const double
623 H_n = ice_thickness(i_n, j_n),
624 Q_advective = v * (v > 0.0 ? H : H_n); // first order upwinding
625
626 // diffusive flux
627 const double
628 Q_diffusive = limit_diffusive_flux(M, M_n, diffusive_flux(i, j, n));
629
630 output(i, j, n) = Q_diffusive + Q_advective;
631 } // end of the loop over neighbors (n)
632 }
633 } catch (...) {
634 loop.failed();
635 }
636 loop.check();
637 }
638
639 /*!
640 * Compute flux divergence using cell interface fluxes on the staggered grid.
641 *
642 * The flux divergence at *ice thickness* Dirichlet B.C. locations is set to zero.
643 */
644 void GeometryEvolution::compute_flux_divergence(const IceModelVec2Stag &flux,
645 const IceModelVec2Int &thickness_bc_mask,
646 IceModelVec2S &output) {
647 const double
648 dx = m_grid->dx(),
649 dy = m_grid->dy();
650
651 IceModelVec::AccessList list{&flux, &thickness_bc_mask, &output};
652
653 ParallelSection loop(m_grid->com);
654 try {
655 for (Points p(*m_grid); p; p.next()) {
656 const int i = p.i(), j = p.j();
657
658 if (thickness_bc_mask(i, j) > 0.5) {
659 output(i, j) = 0.0;
660 } else {
661 StarStencil<double> Q = flux.star(i, j);
662
663 output(i, j) = (Q.e - Q.w) / dx + (Q.n - Q.s) / dy;
664 }
665 }
666 } catch (...) {
667 loop.failed();
668 }
669 loop.check();
670 }
671
672 /*!
673 * Update ice thickness and area_specific_volume *in place*.
674 *
675 * It would be better to compute the change in ice thickness and area_specific_volume and then apply
676 * them, but it would require re-writing all the part-grid code from scratch. So, I make copies of
677 * ice thickness and area_specific_volume, use this old code, then compute differences to get changes.
678 * Compute ice thickness changes due to the flow of the ice.
679 *
680 * @param[in] dt time step, seconds
681 * @param[in] bed_elevation bed elevation, meters
682 * @param[in] sea_level sea level elevation
683 * @param[in] ice_thickness ice thickness
684 * @param[in] area_specific_volume area-specific volume (m3/m2)
685 * @param[in] flux_divergence flux divergence
686 * @param[out] thickness_change ice thickness change due to flow
687 * @param[out] area_specific_volume_change area specific volume change due to flow
688 */
689 void GeometryEvolution::update_in_place(double dt,
690 const IceModelVec2S &bed_topography,
691 const IceModelVec2S &sea_level,
692 const IceModelVec2S &flux_divergence,
693 IceModelVec2S &ice_thickness,
694 IceModelVec2S &area_specific_volume) {
695
696 m_impl->gc.compute(sea_level, bed_topography, ice_thickness,
697 m_impl->cell_type, m_impl->surface_elevation);
698
699 IceModelVec::AccessList list{&ice_thickness, &flux_divergence};
700
701 if (m_impl->use_part_grid) {
702 m_impl->residual.set(0.0);
703
704 // Store ice thickness. We need this copy to make sure that modifying ice_thickness in the loop
705 // below does not affect the computation of the threshold thickness. (Note that
706 // part_grid_threshold_thickness uses neighboring values of the mask, ice thickness, and surface
707 // elevation.)
708 m_impl->thickness.copy_from(ice_thickness);
709
710 list.add({&area_specific_volume, &m_impl->residual, &m_impl->thickness,
711 &m_impl->surface_elevation, &bed_topography, &m_impl->cell_type});
712 }
713
714 #if (Pism_DEBUG==1)
715 const double Lz = m_grid->Lz();
716 #endif
717
718 ParallelSection loop(m_grid->com);
719 try {
720 for (Points p(*m_grid); p; p.next()) {
721 const int i = p.i(), j = p.j();
722
723 double divQ = flux_divergence(i, j);
724
725 if (m_impl->use_part_grid) {
726 if (m_impl->cell_type.ice_free_ocean(i, j) and m_impl->cell_type.next_to_ice(i, j)) {
727 // Add the flow contribution to this partially filled cell.
728 area_specific_volume(i, j) += -divQ * dt;
729
730 double threshold = part_grid_threshold_thickness(m_impl->cell_type.int_star(i, j),
731 m_impl->thickness.star(i, j),
732 m_impl->surface_elevation.star(i, j),
733 bed_topography(i, j));
734
735 // if threshold is zero, turn all the area specific volume into ice thickness, with zero
736 // residual
737 if (threshold == 0.0) {
738 threshold = area_specific_volume(i, j);
739 }
740
741 if (area_specific_volume(i, j) >= threshold) {
742 ice_thickness(i, j) += threshold;
743 m_impl->residual(i, j) = area_specific_volume(i, j) - threshold;
744 area_specific_volume(i, j) = 0.0;
745 }
746
747 // In this case the flux goes into the area_specific_volume variable and does not directly
748 // contribute to ice thickness at this location.
749 divQ = 0.0;
750 }
751 } // end of if (use_part_grid)
752
753 ice_thickness(i, j) += - dt * divQ;
754
755 #if (Pism_DEBUG==1)
756 if (ice_thickness(i, j) > Lz) {
757 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "ice thickness exceeds Lz at i=%d, j=%d (H=%f, Lz=%f)",
758 i, j, ice_thickness(i, j), Lz);
759 }
760 #endif
761 }
762 } catch (...) {
763 loop.failed();
764 }
765 loop.check();
766
767 ice_thickness.update_ghosts();
768
769 // Compute the mask corresponding to the new thickness.
770 m_impl->gc.compute_mask(sea_level, bed_topography, ice_thickness, m_impl->cell_type);
771
772 /*
773 Redistribute residual ice mass from subgrid-scale parameterization.
774
775 See [@ref Albrechtetal2011].
776 */
777 if (m_impl->use_part_grid) {
778 const int max_n_iterations = m_config->get_number("geometry.part_grid.max_iterations");
779
780 bool done = false;
781 for (int i = 0; i < max_n_iterations and not done; ++i) {
782 m_log->message(4, "redistribution iteration %d\n", i);
783
784 // this call may set done to true
785 residual_redistribution_iteration(bed_topography,
786 sea_level,
787 m_impl->surface_elevation,
788 ice_thickness,
789 m_impl->cell_type,
790 area_specific_volume,
791 m_impl->residual,
792 done);
793 }
794
795 if (not done) {
796 m_log->message(2,
797 "WARNING: not done redistributing mass after %d iterations, remaining residual: %f m^3.\n",
798 max_n_iterations, m_impl->residual.sum() * m_grid->cell_area());
799
800 // Add residual to ice thickness, preserving total ice mass. (This is not great, but
801 // better than losing mass.)
802 ice_thickness.add(1.0, m_impl->residual);
803 m_impl->residual.set(0.0);
804 }
805 }
806 }
807
808 //! @brief Perform one iteration of the residual mass redistribution.
809 /*!
810 @param[in] bed_topography bed elevation
811 @param[in] sea_level sea level elevation
812 @param[in,out] ice_surface_elevation surface elevation; used as temp. storage
813 @param[in,out] ice_thickness ice thickness; updated
814 @param[in,out] cell_type cell type mask; used as temp. storage
815 @param[in,out] area_specific_volume area specific volume; updated
816 @param[in,out] residual ice volume that still needs to be distributed; updated
817 @param[in,out] done result flag: true if this iteration should be the last one
818 */
819 void GeometryEvolution::residual_redistribution_iteration(const IceModelVec2S &bed_topography,
820 const IceModelVec2S &sea_level,
821 IceModelVec2S &ice_surface_elevation,
822 IceModelVec2S &ice_thickness,
823 IceModelVec2CellType &cell_type,
824 IceModelVec2S &area_specific_volume,
825 IceModelVec2S &residual,
826 bool &done) {
827
828 m_impl->gc.compute_mask(sea_level, bed_topography, ice_thickness, cell_type);
829
830 const Direction directions[4] = {North, East, South, West};
831
832 // First step: distribute residual mass
833 {
834 // will be destroyed at the end of the block
835 IceModelVec::AccessList list{&cell_type, &ice_thickness, &area_specific_volume, &residual};
836
837 for (Points p(*m_grid); p; p.next()) {
838 const int i = p.i(), j = p.j();
839
840 if (residual(i, j) <= 0.0) {
841 continue;
842 }
843
844 StarStencil<int> m = cell_type.int_star(i, j);
845
846 int N = 0; // number of empty or partially filled neighbors
847 for (unsigned int n = 0; n < 4; ++n) {
848 const Direction direction = directions[n];
849 if (ice_free_ocean(m[direction])) {
850 N++;
851 }
852 }
853
854 if (N > 0) {
855 // Remaining ice mass will be redistributed equally among all adjacent
856 // ice-free-ocean cells (is there a more physical way?)
857 residual(i, j) /= N;
858 } else {
859 // Conserve mass, but (possibly) create a "ridge" at the shelf
860 // front
861 ice_thickness(i, j) += residual(i, j);
862 residual(i, j) = 0.0;
863 }
864 }
865
866 residual.update_ghosts();
867
868 // update area_specific_volume using adjusted residuals
869 for (Points p(*m_grid); p; p.next()) {
870 const int i = p.i(), j = p.j();
871
872 if (cell_type.ice_free_ocean(i, j)) {
873 area_specific_volume(i, j) += (residual(i + 1, j) +
874 residual(i - 1, j) +
875 residual(i, j + 1) +
876 residual(i, j - 1));
877 }
878
879 }
880
881 residual.set(0.0);
882 }
883
884 ice_thickness.update_ghosts();
885
886 // Store ice thickness. We need this copy to make sure that modifying ice_thickness in the loop
887 // below does not affect the computation of the threshold thickness. (Note that
888 // part_grid_threshold_thickness uses neighboring values of the mask, ice thickness, and surface
889 // elevation.)
890 m_impl->thickness.copy_from(ice_thickness);
891
892 // The loop above updated ice_thickness, so we need to re-calculate the mask and the
893 // surface elevation:
894 m_impl->gc.compute(sea_level, bed_topography, ice_thickness, cell_type, ice_surface_elevation);
895
896 double remaining_residual = 0.0;
897
898 // Second step: we need to redistribute residual ice volume if
899 // neighbors which gained redistributed ice also become full.
900 {
901 // will be destroyed at the end of the block
902 IceModelVec::AccessList list{&m_impl->thickness, &ice_thickness,
903 &ice_surface_elevation, &bed_topography, &cell_type};
904
905 for (Points p(*m_grid); p; p.next()) {
906 const int i = p.i(), j = p.j();
907
908 if (area_specific_volume(i, j) <= 0.0) {
909 continue;
910 }
911
912 double threshold = part_grid_threshold_thickness(cell_type.int_star(i, j),
913 m_impl->thickness.star(i, j),
914 ice_surface_elevation.star(i, j),
915 bed_topography(i, j));
916
917 // if threshold is zero, turn all the area specific volume into ice thickness, with zero
918 // residual
919 if (threshold == 0.0) {
920 threshold = area_specific_volume(i, j);
921 }
922
923 if (area_specific_volume(i, j) >= threshold) {
924 ice_thickness(i, j) += threshold;
925 residual(i, j) = area_specific_volume(i, j) - threshold;
926 area_specific_volume(i, j) = 0.0;
927
928 remaining_residual += residual(i, j);
929 }
930 }
931 }
932
933 // check if redistribution should be run once more
934 remaining_residual = GlobalSum(m_grid->com, remaining_residual);
935
936 if (remaining_residual > 0.0) {
937 done = false;
938 } else {
939 done = true;
940 }
941
942 ice_thickness.update_ghosts();
943 }
944
945 /*!
946 * Correct `thickness_change` and `area_specific_volume_change` so that applying them will not
947 * result in negative `ice_thickness` and `area_specific_volume`.
948 *
949 * Compute the `conservation_error`, i.e. the amount of ice that is added to preserve
950 * non-negativity.
951 *
952 * @param[in] ice_thickness ice thickness (m)
953 * @param[in] area_specific_volume area-specific volume (m3/m2)
954 * @param[in,out] thickness_change "proposed" thickness change (m)
955 * @param[in,out] area_specific_volume_change "proposed" area-specific volume change (m3/m2)
956 * @param[out] conservation_error computed conservation error (m)
957 *
958 * This computation is purely local.
959 */
960 void GeometryEvolution::ensure_nonnegativity(const IceModelVec2S &ice_thickness,
961 const IceModelVec2S &area_specific_volume,
962 IceModelVec2S &thickness_change,
963 IceModelVec2S &area_specific_volume_change,
964 IceModelVec2S &conservation_error) {
965
966 IceModelVec::AccessList list{&ice_thickness, &area_specific_volume, &thickness_change,
967 &area_specific_volume_change, &conservation_error};
968
969 ParallelSection loop(m_grid->com);
970 try {
971 for (Points p(*m_grid); p; p.next()) {
972 const int i = p.i(), j = p.j();
973
974 conservation_error(i, j) = 0.0;
975
976 const double
977 H = ice_thickness(i, j),
978 dH = thickness_change(i, j);
979
980 // applying thickness_change will lead to negative thickness
981 if (H + dH < 0.0) {
982 thickness_change(i, j) = H;
983 conservation_error(i, j) += - (H + dH);
984 }
985
986 const double
987 V = area_specific_volume(i, j),
988 dV = area_specific_volume_change(i, j);
989
990 if (V + dV < 0.0) {
991 area_specific_volume_change(i, j) = V;
992 conservation_error(i, j) += - (V + dV);
993 }
994 }
995 } catch (...) {
996 loop.failed();
997 }
998 loop.check();
999 }
1000
1001 /*!
1002 * Given ice thickness `H` and the "proposed" change `dH`, compute the corrected change preserving
1003 * non-negativity.
1004 */
1005 static inline double effective_change(double H, double dH) {
1006 if (H + dH <= 0) {
1007 return -H;
1008 } else {
1009 return dH;
1010 }
1011 }
1012
1013 /*!
1014 * Compute effective surface and basal mass balance.
1015 *
1016 * @param[in] dt time step, seconds
1017 * @param[in] thickness_bc_mask mask specifying ice thickness Dirichlet B.C. locations
1018 * @param[in] ice_thickness ice thickness, m
1019 * @param[in] thickness_change thickness change due to flow, m
1020 * @param[in] cell_type cell type mask
1021 * @param[in] smb_rate top surface mass balance rate, kg m-2 s-1
1022 * @param[in] basal_melt_rate basal melt rate, m s-1
1023 * @param[out] effective_smb effective top surface mass balance, m
1024 * @param[out] effective_bmb effective basal mass balance, m
1025 *
1026 * This computation is purely local.
1027 */
1028 void GeometryEvolution::compute_surface_and_basal_mass_balance(double dt,
1029 const IceModelVec2Int &thickness_bc_mask,
1030 const IceModelVec2S &ice_thickness,
1031 const IceModelVec2CellType &cell_type,
1032 const IceModelVec2S &smb_flux,
1033 const IceModelVec2S &basal_melt_rate,
1034 IceModelVec2S &effective_SMB,
1035 IceModelVec2S &effective_BMB) {
1036
1037 IceModelVec::AccessList list{&ice_thickness,
1038 &smb_flux, &basal_melt_rate, &cell_type, &thickness_bc_mask,
1039 &effective_SMB, &effective_BMB};
1040
1041 ParallelSection loop(m_grid->com);
1042 try {
1043 for (Points p(*m_grid); p; p.next()) {
1044 const int i = p.i(), j = p.j();
1045
1046 // Don't modify ice thickness at Dirichlet B.C. locations and in the ice-free ocean.
1047 if (thickness_bc_mask.as_int(i, j) == 1 or cell_type.ice_free_ocean(i, j)) {
1048 effective_SMB(i, j) = 0.0;
1049 effective_BMB(i, j) = 0.0;
1050 continue;
1051 }
1052
1053 const double H = ice_thickness(i, j);
1054
1055 // Thickness change due to the surface mass balance
1056 //
1057 // Note that here we convert surface mass balance from [kg m-2 s-1] to [m s-1].
1058 double dH_SMB = effective_change(H, dt * smb_flux(i, j) / m_impl->ice_density);
1059
1060 // Thickness change due to the basal mass balance
1061 //
1062 // Note that basal_melt_rate is in [m s-1]. Here the negative sign converts the melt rate into
1063 // mass balance.
1064 double dH_BMB = effective_change(H + dH_SMB,
1065 dt * (m_impl->use_bmr ? -basal_melt_rate(i, j) : 0.0));
1066
1067 effective_SMB(i, j) = dH_SMB;
1068 effective_BMB(i, j) = dH_BMB;
1069 }
1070 } catch (...) {
1071 loop.failed();
1072 }
1073 loop.check();
1074 }
1075
1076 namespace diagnostics {
1077
1078 /*! @brief Report the divergence of the ice flux. */
1079 class FluxDivergence : public Diag<GeometryEvolution>
1080 {
1081 public:
1082 FluxDivergence(const GeometryEvolution *m)
1083 : Diag<GeometryEvolution>(m) {
1084 m_vars = {model->flux_divergence().metadata()};
1085 }
1086 protected:
1087 IceModelVec::Ptr compute_impl() const {
1088 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "flux_divergence", WITHOUT_GHOSTS));
1089 result->metadata(0) = m_vars[0];
1090
1091 result->copy_from(model->flux_divergence());
1092
1093 return result;
1094 }
1095 };
1096
1097 /*! @brief Report mass flux on the staggered grid. */
1098 class FluxStaggered : public Diag<GeometryEvolution>
1099 {
1100 public:
1101 FluxStaggered(const GeometryEvolution *m)
1102 : Diag<GeometryEvolution>(m) {
1103 m_vars = {model->flux_staggered().metadata()};
1104 }
1105 protected:
1106 IceModelVec::Ptr compute_impl() const {
1107 IceModelVec2Stag::Ptr result(new IceModelVec2Stag(m_grid, "flux_staggered", WITHOUT_GHOSTS));
1108 result->metadata(0) = m_vars[0];
1109
1110 const IceModelVec2Stag &input = model->flux_staggered();
1111 IceModelVec2Stag &output = *result.get();
1112
1113 // FIXME: implement IceModelVec2Stag::copy_from()
1114
1115 IceModelVec::AccessList list{&input, &output};
1116
1117 ParallelSection loop(m_grid->com);
1118 try {
1119 for (Points p(*m_grid); p; p.next()) {
1120 const int i = p.i(), j = p.j();
1121
1122 output(i, j, 0) = input(i, j, 0);
1123 output(i, j, 1) = input(i, j, 1);
1124 }
1125 } catch (...) {
1126 loop.failed();
1127 }
1128 loop.check();
1129
1130 return result;
1131 }
1132 };
1133
1134 } // end of namespace diagnostics
1135
1136 DiagnosticList GeometryEvolution::diagnostics_impl() const {
1137 using namespace diagnostics;
1138 typedef Diagnostic::Ptr Ptr;
1139
1140 std::map<std::string, Ptr> result;
1141 result = {
1142 {"flux_staggered", Ptr(new FluxStaggered(this))},
1143 {"flux_divergence", Ptr(new FluxDivergence(this))},
1144 };
1145 return result;
1146 }
1147
1148 RegionalGeometryEvolution::RegionalGeometryEvolution(IceGrid::ConstPtr grid)
1149 : GeometryEvolution(grid) {
1150
1151 m_no_model_mask.create(m_grid, "no_model_mask", WITH_GHOSTS);
1152 m_no_model_mask.set_attrs("model_mask", "'no model' mask", "", "", "", 0);
1153 }
1154
1155 void RegionalGeometryEvolution::set_no_model_mask_impl(const IceModelVec2Int &mask) {
1156 m_no_model_mask.copy_from(mask);
1157 }
1158
1159 /*!
1160 * Disable ice flow in "no model" areas.
1161 */
1162 void RegionalGeometryEvolution::compute_interface_fluxes(const IceModelVec2CellType &cell_type,
1163 const IceModelVec2S &ice_thickness,
1164 const IceModelVec2V &velocity,
1165 const IceModelVec2Int &velocity_bc_mask,
1166 const IceModelVec2Stag &diffusive_flux,
1167 IceModelVec2Stag &output) {
1168
1169 GeometryEvolution::compute_interface_fluxes(cell_type, ice_thickness,
1170 velocity, velocity_bc_mask, diffusive_flux,
1171 output);
1172
1173 IceModelVec::AccessList list{&m_no_model_mask, &output};
1174
1175 ParallelSection loop(m_grid->com);
1176 try {
1177 for (Points p(*m_grid); p; p.next()) {
1178 const int i = p.i(), j = p.j();
1179
1180 const int M = m_no_model_mask.as_int(i, j);
1181
1182 for (unsigned int n = 0; n < 2; ++n) {
1183 const int
1184 oi = 1 - n, // offset in the i direction
1185 oj = n, // offset in the j direction
1186 i_n = i + oi, // i index of a neighbor
1187 j_n = j + oj; // j index of a neighbor
1188
1189 const int M_n = m_no_model_mask.as_int(i_n, j_n);
1190
1191 if (not (M == 0 and M_n == 0)) {
1192 output(i, j, n) = 0.0;
1193 }
1194 }
1195 }
1196 } catch (...) {
1197 loop.failed();
1198 }
1199 loop.check();
1200 }
1201
1202 /*!
1203 * Set surface and basal mass balance to zero in "no model" areas.
1204 */
1205 void RegionalGeometryEvolution::compute_surface_and_basal_mass_balance(double dt,
1206 const IceModelVec2Int &thickness_bc_mask,
1207 const IceModelVec2S &ice_thickness,
1208 const IceModelVec2CellType &cell_type,
1209 const IceModelVec2S &surface_mass_flux,
1210 const IceModelVec2S &basal_melt_rate,
1211 IceModelVec2S &effective_SMB,
1212 IceModelVec2S &effective_BMB) {
1213 GeometryEvolution::compute_surface_and_basal_mass_balance(dt,
1214 thickness_bc_mask,
1215 ice_thickness,
1216 cell_type,
1217 surface_mass_flux,
1218 basal_melt_rate,
1219 effective_SMB,
1220 effective_BMB);
1221
1222 IceModelVec::AccessList list{&m_no_model_mask, &effective_SMB, &effective_BMB};
1223
1224 ParallelSection loop(m_grid->com);
1225 try {
1226 for (Points p(*m_grid); p; p.next()) {
1227 const int i = p.i(), j = p.j();
1228
1229 if (m_no_model_mask(i, j) > 0.5) {
1230 effective_SMB(i, j) = 0.0;
1231 effective_BMB(i, j) = 0.0;
1232 }
1233 }
1234 } catch (...) {
1235 loop.failed();
1236 }
1237 loop.check();
1238 }
1239
1240 void GeometryEvolution::set_no_model_mask(const IceModelVec2Int &mask) {
1241 this->set_no_model_mask_impl(mask);
1242 }
1243
1244 void GeometryEvolution::set_no_model_mask_impl(const IceModelVec2Int &mask) {
1245 (void) mask;
1246 // the default implementation is a no-op
1247 }
1248
1249 void grounding_line_flux(const IceModelVec2CellType &cell_type,
1250 const IceModelVec2Stag &flux,
1251 double dt,
1252 InsertMode flag,
1253 IceModelVec2S &output) {
1254
1255 using mask::grounded;
1256
1257 auto grid = output.grid();
1258
1259 const double
1260 dx = grid->dx(),
1261 dy = grid->dy();
1262
1263 auto cell_area = grid->cell_area();
1264
1265 auto ice_density = grid->ctx()->config()->get_number("constants.ice.density");
1266
1267 IceModelVec::AccessList list{&cell_type, &flux, &output};
1268
1269 ParallelSection loop(grid->com);
1270 try {
1271 for (Points p(*grid); p; p.next()) {
1272 const int i = p.i(), j = p.j();
1273
1274 double result = 0.0;
1275
1276 if (cell_type.ocean(i ,j)) {
1277 auto M = cell_type.int_star(i, j);
1278 auto Q = flux.star(i, j);
1279
1280 if (grounded(M.n) and Q.n <= 0.0) {
1281 result += Q.n * dx;
1282 }
1283
1284 if (grounded(M.e) and Q.e <= 0.0) {
1285 result += Q.e * dy;
1286 }
1287
1288 if (grounded(M.s) and Q.s >= 0.0) {
1289 result -= Q.s * dx;
1290 }
1291
1292 if (grounded(M.w) and Q.w >= 0.0) {
1293 result -= Q.w * dy;
1294 }
1295
1296 // convert from "m^3 / s" to "kg / m^2"
1297 result *= dt * (ice_density / cell_area);
1298 }
1299
1300 if (flag == ADD_VALUES) {
1301 output(i, j) += result;
1302 } else {
1303 output(i, j) = result;
1304 }
1305 }
1306 } catch (...) {
1307 loop.failed();
1308 }
1309 loop.check();
1310 }
1311
1312 /*!
1313 * Compute the total grounding line flux over a time step, in kg.
1314 */
1315 double total_grounding_line_flux(const IceModelVec2CellType &cell_type,
1316 const IceModelVec2Stag &flux,
1317 double dt) {
1318 using mask::grounded;
1319
1320 auto grid = cell_type.grid();
1321
1322 const double
1323 dx = grid->dx(),
1324 dy = grid->dy();
1325
1326 auto ice_density = grid->ctx()->config()->get_number("constants.ice.density");
1327
1328 double total_flux = 0.0;
1329
1330 IceModelVec::AccessList list{&cell_type, &flux};
1331
1332 ParallelSection loop(grid->com);
1333 try {
1334 for (Points p(*grid); p; p.next()) {
1335 const int i = p.i(), j = p.j();
1336
1337 double volume_flux = 0.0;
1338
1339 if (cell_type.ocean(i ,j)) {
1340 auto M = cell_type.int_star(i, j);
1341 auto Q = flux.star(i, j); // m^2 / s
1342
1343 if (grounded(M.n) and Q.n <= 0.0) {
1344 volume_flux += Q.n * dx;
1345 }
1346
1347 if (grounded(M.e) and Q.e <= 0.0) {
1348 volume_flux += Q.e * dy;
1349 }
1350
1351 if (grounded(M.s) and Q.s >= 0.0) {
1352 volume_flux -= Q.s * dx;
1353 }
1354
1355 if (grounded(M.w) and Q.w >= 0.0) {
1356 volume_flux -= Q.w * dy;
1357 }
1358 }
1359
1360 // convert from "m^3 / s" to "kg" and sum up
1361 total_flux += volume_flux * dt * ice_density;
1362 }
1363 } catch (...) {
1364 loop.failed();
1365 }
1366 loop.check();
1367
1368 return GlobalSum(grid->com, total_flux);
1369 }
1370
1371 } // end of namespace pism