22#include "pism/geometry/Geometry.hh"
24#include "pism/util/array/CellType.hh"
25#include "pism/util/Mask.hh"
26#include "pism/util/pism_utilities.hh"
27#include "pism/geometry/grounded_cell_fraction.hh"
28#include "pism/util/Context.hh"
29#include "pism/util/VariableMetadata.hh"
30#include "pism/util/io/File.hh"
31#include "pism/util/io/io_helpers.hh"
32#include "pism/util/io/IO_Flags.hh"
39 : latitude(grid,
"lat"),
40 longitude(grid,
"lon"),
41 bed_elevation(grid,
"topg"),
42 sea_level_elevation(grid,
"sea_level"),
43 ice_thickness(grid,
"thk"),
44 ice_area_specific_volume(grid,
"ice_area_specific_volume"),
45 cell_type(grid,
"mask"),
46 cell_grounded_fraction(grid,
"cell_grounded_fraction"),
47 ice_surface_elevation(grid,
"usurf") {
51 .
units(
"degree_north")
71 .
long_name(
"sea level elevation above reference ellipsoid")
73 .
standard_name(
"sea_surface_height_above_reference_ellipsoid");
82 .
long_name(
"ice-volume-per-area in partially-filled grid cells")
85 "this variable represents the amount of ice in a partially-filled cell and not "
86 "the corresponding geometry, so thinking about it as 'thickness' is not helpful";
89 .
long_name(
"ice-type (ice-free/grounded/floating/ocean) integer mask")
94 "ice_free_bedrock grounded_ice floating_ice ice_free_ocean";
97 "fractional grounded/floating mask (floating=0, grounded=1)");
100 .
long_name(
"ice upper surface elevation")
126 for (
auto p = grid->points(); p; p.next()) {
127 const int i = p.i(), j = p.j();
131 "H = %e (negative) at point i=%d, j=%d",
153 for (
auto p = grid->points(); p; p.next()) {
154 const int i = p.i(), j = p.j();
173 ice_density = config->get_number(
"constants.ice.density"),
174 ocean_density = config->get_number(
"constants.sea_water.density");
184 e.
add_context(
"computing the grounded cell fraction");
186 std::string output_file = config->get_string(
"output.file");
188 "_grounded_cell_fraction_failed",
"");
190 dump(o_file.c_str());
198 File file(grid->com, filename,
220 auto grid = result.
grid();
221 auto config = grid->ctx()->config();
224 ice_density = config->get_number(
"constants.ice.density"),
225 water_density = config->get_number(
"constants.sea_water.density"),
226 alpha = ice_density / water_density;
236 for (
auto p = grid->points(); p; p.next()) {
237 const int i = p.i(), j = p.j();
240 b_grounded = bed_elevation(i, j),
241 b_floating = sea_level(i, j) - alpha * ice_thickness(i, j);
243 result(i, j) = std::max(b_grounded, b_floating);
256 auto config = grid->ctx()->config();
262 auto cell_area = grid->cell_area();
265 for (
auto p = grid->points(); p; p.next()) {
266 const int i = p.i(), j = p.j();
275 if (config->get_flag(
"geometry.part_grid.enabled")) {
277 for (
auto p = grid->points(); p; p.next()) {
278 const int i = p.i(), j = p.j();
288 double thickness_threshold) {
290 auto config = grid->ctx()->config();
293 sea_water_density = config->get_number(
"constants.sea_water.density"),
294 ice_density = config->get_number(
"constants.ice.density"),
295 cell_area = grid->cell_area();
302 for (
auto p = grid->points(); p; p.next()) {
303 const int i = p.i(), j = p.j();
311 double max_floating_thickness =
312 std::max(sea_level - bed, 0.0) * (sea_water_density / ice_density);
313 volume += cell_area * (thickness - max_floating_thickness);
321 double cell_area = grid.cell_area();
324 for (
auto p = grid.points(); p; p.next()) {
325 const int i = p.i(), j = p.j();
327 if (condition(i, j)) {
339 return geometry.ice_thickness(i, j) >= thickness_threshold;
347 return (geometry.cell_type.grounded(i, j) and
348 geometry.ice_thickness(i, j) >= thickness_threshold);
356 return (geometry.cell_type.ocean(i, j) and geometry.ice_thickness(i, j) >= thickness_threshold);
366 water_density = config->get_number(
"constants.fresh_water.density"),
367 ice_density = config->get_number(
"constants.ice.density"),
368 ocean_area = config->get_number(
"constants.global_ocean_area");
372 thickness_threshold),
373 additional_water_volume = (ice_density / water_density) * volume,
374 sea_level_change = additional_water_volume / ocean_area;
376 return sea_level_change;
392 for (
auto p = grid.points(); p; p.next()) {
393 const int i = p.i(), j = p.j();
std::shared_ptr< const Config > ConstPtr
High-level PISM I/O class.
void set_icefree_thickness(double threshold)
void compute(const array::Scalar &sea_level, const array::Scalar &bed, const array::Scalar &thickness, array::Scalar &out_mask, array::Scalar &out_surface) const
array::Scalar1 sea_level_elevation
array::Scalar cell_grounded_fraction
void ensure_consistency(double ice_free_thickness_threshold)
array::Scalar2 ice_surface_elevation
array::Scalar1 ice_area_specific_volume
array::CellType2 cell_type
void dump(const char *filename) const
array::Scalar2 ice_thickness
Geometry(const std::shared_ptr< const Grid > &grid)
array::Scalar2 bed_elevation
Describes the PISM grid and the distribution of data across processors.
void failed()
Indicates a failure of a parallel section.
void add_context(const std::string &message)
Add a message providing some context. This way we can (sort of) get a stack trace even though C++ exc...
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
void add(double alpha, const Array2D< T > &x)
void write(const std::string &filename) const
std::shared_ptr< const Grid > grid() const
void set(double c)
Result: v[j] <- c for all j.
void update_ghosts()
Updates ghost points.
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
bool grounded(int i, int j) const
#define PISM_ERROR_LOCATION
bool in_null_strip(const Grid &grid, int i, int j, double strip_width)
Check if a point (i,j) is in the strip of stripwidth meters around the edge of the computational doma...
void append_time(const File &file, const Config &config, double time_seconds)
Prepare a file for output.
@ PISM_READWRITE_CLOBBER
create a file for writing, overwrite if present
void define_time(const File &file, const Context &ctx)
Prepare a file for output.
double ice_volume_not_displacing_seawater(const Geometry &geometry, double thickness_threshold)
io::Backend string_to_backend(const std::string &backend)
double ice_area(const Geometry &geometry, double thickness_threshold)
Computes ice area, in m^2.
void set_no_model_strip(const Grid &grid, double width, array::Scalar &result)
Set no_model_mask variable to have value 1 in strip of width 'strip' m around edge of computational d...
void compute_grounded_cell_fraction(double ice_density, double ocean_density, const array::Scalar1 &sea_level, const array::Scalar1 &ice_thickness, const array::Scalar1 &bed_topography, array::Scalar &result)
double ice_area_floating(const Geometry &geometry, double thickness_threshold)
Computes floating ice area, in m^2.
double sea_level_rise_potential(const Geometry &geometry, double thickness_threshold)
Computes the sea level rise that would result if all the ice were melted.
std::string filename_add_suffix(const std::string &filename, const std::string &separator, const std::string &suffix)
Adds a suffix to a filename.
double ice_volume(const Geometry &geometry, double thickness_threshold)
Computes the ice volume, in m^3.
void ice_bottom_surface(const Geometry &geometry, array::Scalar &result)
static double compute_area(const Grid &grid, std::function< bool(int, int)> condition)
double ice_area_grounded(const Geometry &geometry, double thickness_threshold)
Computes grounded ice area, in m^2.
void GlobalSum(MPI_Comm comm, double *local, double *result, int count)