21#include "pism/coupler/ocean/PicoGeometry.hh"
22#include "pism/util/connected_components/label_components.hh"
23#include "pism/util/array/CellType.hh"
24#include "pism/util/pism_utilities.hh"
26#include "pism/coupler/util/options.hh"
27#include "pism/util/Interpolation1D.hh"
28#include "pism/util/Profiling.hh"
35 m_continental_shelf(grid,
"pico_contshelf_mask"),
36 m_boxes(grid,
"pico_box_mask"),
37 m_ice_shelves(grid,
"pico_shelf_mask"),
38 m_basin_mask(m_grid,
"basins"),
39 m_distance_gl(grid,
"pico_distance_gl"),
40 m_distance_cf(grid,
"pico_distance_cf"),
41 m_ocean_mask(grid,
"pico_ocean_mask"),
42 m_lake_mask(grid,
"pico_lake_mask"),
43 m_ice_rises(grid,
"pico_ice_rise_mask"),
44 m_tmp(grid,
"temporary_storage") {
61 "ocean ice_rise continental_ice_sheet, floating_ice";
114 for (
int basin_id = 1; basin_id <
m_n_basins; ++basin_id) {
116 if (not set.empty()) {
117 std::vector<std::string> neighbors;
118 for (
const auto &
n : set) {
121 std::string neighbor_list =
pism::join(neighbors,
", ");
122 m_log->message(3,
"PICO: basin %d neighbors: %s\n", basin_id, neighbor_list.c_str());
127 bool exclude_ice_rises =
m_config->get_flag(
"ocean.pico.exclude_ice_rises");
158 std::vector<int> cfs_in_basins_per_shelf(n_shelves*
m_n_basins, 0);
159 std::vector<int> most_shelf_cells_in_basin(n_shelves, 0);
161 most_shelf_cells_in_basin, cfs_in_basins_per_shelf);
164 most_shelf_cells_in_basin, cfs_in_basins_per_shelf, n_shelves,
167 double continental_shelf_depth =
m_config->get_number(
"ocean.pico.continental_shelf_depth");
173 int n_boxes =
static_cast<int>(
m_config->get_number(
"ocean.pico.number_of_boxes"));
194 auto grid = mask.
grid();
196 int max_index =
static_cast<int>(
array::max(mask));
204 std::vector<double> area(max_index + 1, 0.0);
205 std::vector<double> area1(max_index + 1, 0.0);
210 for (
auto p = grid->points(); p; p.next()) {
211 const int i = p.i(), j = p.j();
213 int index = mask.
as_int(i, j);
215 if (index > max_index or index < 0) {
228 GlobalSum(grid->com, area.data(), area1.data(),
static_cast<int>(area.size()));
233 for (
unsigned int k = 0;
k < area.size(); ++
k) {
234 area[
k] = grid->cell_area() * area[
k];
240 int biggest_component = 0;
241 for (
unsigned int k = 0;
k < area.size(); ++
k) {
242 if (area[
k] > area[biggest_component]) {
243 biggest_component =
static_cast<int>(
k);
248 for (
auto p = grid->points(); p; p.next()) {
249 const int i = p.i(), j = p.j();
251 int component_index = mask.
as_int(i, j);
253 if (component_index == biggest_component) {
255 }
else if (component_index > 0) {
262 for (
auto p = grid->points(); p; p.next()) {
263 const int i = p.i(), j = p.j();
265 int component_index = mask.
as_int(i, j);
267 if (area[component_index] > threshold) {
269 }
else if (component_index > 0) {
295 int reachable_from_domain_edge = 2;
300 for (
auto p =
grid->points(); p; p.next()) {
301 const int i = p.i(), j = p.j();
303 if (cell_type.
ocean(i, j)) {
304 m_tmp(i, j) = floating;
307 m_tmp(i, j) = reachable_from_domain_edge;
310 m_tmp(i, j) = background;
342 for (
auto p =
m_grid->points(); p; p.next()) {
343 const int i = p.i(), j = p.j();
353 if (exclude_ice_rises) {
362 for (
auto p =
m_grid->points(); p; p.next()) {
363 const int i = p.i(), j = p.j();
365 if (
m_tmp(i, j) == 0.0 and cell_type.
icy(i, j)) {
386 double bed_elevation_threshold,
392 for (
auto p =
m_grid->points(); p; p.next()) {
393 const int i = p.i(), j = p.j();
397 if (bed_elevation(i, j) > bed_elevation_threshold) {
415 for (
auto p =
m_grid->points(); p; p.next()) {
416 const int i = p.i(), j = p.j();
418 if (
m_tmp(i, j) > 0.0) {
448 for (
auto p =
m_grid->points(); p; p.next()) {
449 const int i = p.i(), j = p.j();
465 for (
auto p =
m_grid->points(); p; p.next()) {
466 const int i = p.i(), j = p.j();
494 for (
auto p =
m_grid->points(); p; p.next()) {
495 const int i = p.i(), j = p.j();
532 auto mark_as_neighbors = [&](
int b1,
int b2) {
538 auto adjacent = [&](
int b1,
int b2) {
539 return adjacency_matrix[b1 *
m_n_basins + b2] > 0;
544 for (
auto p =
m_grid->points(); p; p.next()) {
545 const int i = p.i(), j = p.j();
552 if (B.c == 0 or not next_to_icefront) {
560 B.n *=
static_cast<int>(j < (
int)
m_grid->My() - 1);
561 B.e *=
static_cast<int>(i < (
int)
m_grid->Mx() - 1);
562 B.s *=
static_cast<int>(j > 0);
563 B.w *=
static_cast<int>(i > 0);
568 if (ice_free_ocean(M.n)) {
569 mark_as_neighbors(B.c, B.n);
572 if (ice_free_ocean(M.s)) {
573 mark_as_neighbors(B.c, B.s);
576 if (ice_free_ocean(M.e)) {
577 mark_as_neighbors(B.c, B.e);
580 if (ice_free_ocean(M.w)) {
581 mark_as_neighbors(B.c, B.w);
587 std::vector<int> tmp(adjacency_matrix.size(), 0);
589 static_cast<int>(tmp.size()));
591 adjacency_matrix = tmp;
595 std::vector<std::set<int> > result(
m_n_basins);
597 for (
int b2 = b1 + 1; b2 <
m_n_basins; ++b2) {
598 if (adjacent(b1, b2)) {
599 result[b1].insert(b2);
600 result[b2].insert(b1);
616 std::vector<int> &most_shelf_cells_in_basin,
617 std::vector<int> &cfs_in_basins_per_shelf) {
619 std::vector<int> n_shelf_cells_per_basin(n_shelves *
m_n_basins,0);
621 std::vector<int> n_shelf_cells_per_basinr(n_shelves *
m_n_basins,0);
622 std::vector<int> cfs_in_basins_per_shelfr(n_shelves *
m_n_basins,0);
623 std::vector<int> most_shelf_cells_in_basinr(n_shelves, 0);
628 for (
auto p =
m_grid->points(); p; p.next()) {
629 const int i = p.i(), j = p.j();
630 int s = shelf_mask.
as_int(i, j);
633 n_shelf_cells_per_basin[sb]++;
636 auto M = cell_type.
star(i, j);
641 if (cfs_in_basins_per_shelf[sb] != b) {
642 cfs_in_basins_per_shelf[sb] = b;
649 cfs_in_basins_per_shelfr.data(), n_shelves*
m_n_basins);
651 n_shelf_cells_per_basinr.data(), n_shelves*
m_n_basins);
653 cfs_in_basins_per_shelf = cfs_in_basins_per_shelfr;
654 n_shelf_cells_per_basin = n_shelf_cells_per_basinr;
656 for (
int s = 0; s < n_shelves; s++) {
657 int n_shelf_cells_per_basin_max = 0;
660 if (n_shelf_cells_per_basin[sb] > n_shelf_cells_per_basin_max) {
661 most_shelf_cells_in_basin[s] = b;
662 n_shelf_cells_per_basin_max = n_shelf_cells_per_basin[sb];
675 const std::vector<std::set<int> > &basin_neighbors,
676 const std::vector<int> &most_shelf_cells_in_basin,
677 const std::vector<int> &cfs_in_basins_per_shelf,
685 std::vector<int> n_shelf_cells_to_split(n_shelves *
m_n_basins, 0);
689 for (
auto p =
m_grid->points(); p; p.next()) {
690 const int i = p.i(), j = p.j();
693 int shelf = shelf_mask.
as_int(i, j);
694 int b0 = most_shelf_cells_in_basin[shelf];
698 if (basin !=
b0 and (not neighbors) and
699 cfs_in_basins_per_shelf[shelf *
m_n_basins + basin] > 0) {
700 n_shelf_cells_to_split[shelf *
m_n_basins + basin]++;
706 std::vector<int> tmp(n_shelves *
m_n_basins, 0);
710 n_shelf_cells_to_split = tmp;
714 std::vector<int> add_shelf_instance(n_shelves *
m_n_basins, 0);
715 int n_shelf_numbers_to_add = 0;
716 for (
int s = 0; s < n_shelves; s++) {
717 int b0 = most_shelf_cells_in_basin[s];
719 if (n_shelf_cells_to_split[s *
m_n_basins + b] > 0) {
720 n_shelf_numbers_to_add += 1;
721 add_shelf_instance[s *
m_n_basins + b] = n_shelves + n_shelf_numbers_to_add;
722 m_log->message(3,
"\nPICO, split ice shelf s=%d with bmax=%d "
723 "and b=%d and n=%d and si=%d\n", s,
b0, b,
730 for (
auto p =
m_grid->points(); p; p.next()) {
731 const int i = p.i(), j = p.j();
734 int s = shelf_mask.
as_int(i, j);
735 if (add_shelf_instance[s *
m_n_basins + b] > 0) {
749 bool exclude_ice_rises,
762 for (
auto p =
m_grid->points(); p; p.next()) {
763 const int i = p.i(), j = p.j();
766 ocean_mask.
as_int(i, j) == 1 or
767 (exclude_ice_rises and ice_rises.
as_int(i, j) ==
RISE)) {
771 bool neighbor_to_land =
781 if (neighbor_to_land) {
806 bool exclude_ice_rises,
815 for (
auto p =
m_grid->points(); p; p.next()) {
816 const int i = p.i(), j = p.j();
819 ocean_mask.
as_int(i, j) == 1 or
820 (exclude_ice_rises and ice_rises.
as_int(i, j) ==
RISE)) {
824 auto M = ocean_mask.
star(i, j);
826 if (M.n == 2 or M.e == 2 or M.s == 2 or M.w == 2) {
869 auto grid = mask.
grid();
871 double current_label = 1;
872 double continue_loop = 1;
873 while (continue_loop != 0) {
877 for (
auto p = grid->points(); p; p.next()) {
878 const int i = p.i(), j = p.j();
880 if (mask.
as_int(i, j) == 0) {
882 auto R = mask.
star(i, j);
885 (R.n == current_label or R.s == current_label or
886 R.e == current_label or R.w == current_label)) {
889 mask(i, j) = current_label + 1;
898 continue_loop =
GlobalMax(grid->com, continue_loop);
912 int n_shelves =
static_cast<int>(
array::max(shelf_mask)) + 1;
914 std::vector<double> GL_distance_max(n_shelves, 0.0);
915 std::vector<double> GL_distance_max1(n_shelves, 0.0);
916 std::vector<double> CF_distance_max(n_shelves, 0.0);
917 std::vector<double> CF_distance_max1(n_shelves, 0.0);
919 for (
auto p =
m_grid->points(); p; p.next()) {
920 const int i = p.i(), j = p.j();
922 int shelf_id = shelf_mask.
as_int(i, j);
923 assert(shelf_id >= 0);
924 assert(shelf_id < n_shelves + 1);
931 if (D_gl(i, j) > GL_distance_max[shelf_id]) {
932 GL_distance_max[shelf_id] = D_gl(i, j);
935 if (D_cf(i, j) > CF_distance_max[shelf_id]) {
936 CF_distance_max[shelf_id] = D_cf(i, j);
941 GlobalMax(
m_grid->com, GL_distance_max.data(), GL_distance_max1.data(), n_shelves);
942 GlobalMax(
m_grid->com, CF_distance_max.data(), CF_distance_max1.data(), n_shelves);
944 GL_distance_max = GL_distance_max1;
945 CF_distance_max = CF_distance_max1;
947 double GL_distance_ref = *std::max_element(GL_distance_max.begin(), GL_distance_max.end());
951 std::vector<int> n_boxes(n_shelves, 0);
953 const double zeta = 0.5;
955 for (
int k = 0;
k < n_shelves; ++
k) {
956 n_boxes[
k] = n_min + round(pow((GL_distance_max[
k] / GL_distance_ref), zeta) * (max_number_of_boxes - n_min));
958 n_boxes[
k] = std::min(n_boxes[
k], max_number_of_boxes);
963 for (
auto p =
m_grid->points(); p; p.next()) {
964 const int i = p.i(), j = p.j();
966 int d_gl = D_gl.
as_int(i, j);
967 int d_cf = D_cf.
as_int(i, j);
969 if (shelf_mask.
as_int(i, j) > 0 and d_gl > 0.0 and d_cf > 0.0) {
970 int shelf_id = shelf_mask.
as_int(i, j);
971 int n = n_boxes[shelf_id];
975 double r = d_gl / (
double)(d_gl + d_cf);
977 double C = pow(1.0 - r, 2.0);
978 for (
int k = 0;
k <
n; ++
k) {
979 if ((
n -
k - 1) / (
double)
n <= C and C <= (
n -
k) / (
double)
n) {
980 result(i, j) = std::min(d_gl,
k + 1);
std::shared_ptr< const Grid > grid() const
const Config::ConstPtr m_config
configuration database used by this component
const Logger::ConstPtr m_log
logger (for easy access)
const std::shared_ptr< const Grid > m_grid
grid used by this component
const Profiling & profiling() const
A class defining a common interface for most PISM sub-models.
void failed()
Indicates a failure of a parallel section.
void begin(const char *name) const
void end(const char *name) const
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 copy_from(const Array2D< T > &source)
stencils::Star< T > star(int i, int j) const
void set_interpolation_type(InterpolationType type)
std::shared_ptr< const Grid > grid() const
void set(double c)
Result: v[j] <- c for all j.
void regrid(const std::string &filename, io::Default default_value)
void update_ghosts()
Updates ghost points.
unsigned int stencil_width() const
Get the stencil width of the current Array. Returns 0 if ghosts are not available.
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
bool next_to_ice(int i, int j) const
Ice-free margin (at least one of four neighbors has ice).
stencils::Star< int > star_int(int i, int j) const
bool ice_free_ocean(int i, int j) const
bool ocean(int i, int j) const
bool grounded(int i, int j) const
bool icy(int i, int j) const
"Cell type" mask. Adds convenience methods to array::Scalar.
int as_int(int i, int j) const
stencils::Star< int > star_int(int i, int j) const
const array::Scalar & continental_shelf_mask() const
void split_ice_shelves(const array::CellType &cell_type, const array::Scalar &basin_mask, const std::vector< std::set< int > > &basin_neighbors, const std::vector< int > &most_shelf_cells_in_basin, const std::vector< int > &cfs_in_basins_per_shelf, int n_shelves, array::Scalar &shelf_mask)
array::Scalar m_continental_shelf
void compute_box_mask(const array::Scalar &D_gl, const array::Scalar &D_cf, const array::Scalar &shelf_mask, int max_number_of_boxes, array::Scalar &result)
array::Scalar1 m_ice_rises
array::Scalar m_lake_mask
PicoGeometry(std::shared_ptr< const Grid > grid)
void compute_ocean_mask(const array::CellType &cell_type, array::Scalar &result)
const array::Scalar & ice_shelf_mask() const
array::Scalar1 m_distance_gl
void compute_distances_gl(const array::Scalar &ocean_mask, const array::Scalar1 &ice_rises, bool exclude_ice_rises, array::Scalar1 &result)
void compute_distances_cf(const array::Scalar1 &ocean_mask, const array::Scalar &ice_rises, bool exclude_ice_rises, array::Scalar1 &result)
std::vector< std::set< int > > basin_neighbors(const array::CellType1 &cell_type, const array::Scalar1 &basin_mask)
const array::Scalar & ice_rise_mask() const
void compute_ice_shelf_mask(const array::Scalar &ice_rise_mask, const array::Scalar &lake_mask, array::Scalar &result)
array::Scalar1 m_basin_mask
void compute_ice_rises(const array::CellType &cell_type, bool exclude_ice_rises, array::Scalar &result)
array::Scalar1 m_ocean_mask
std::vector< std::set< int > > m_basin_neighbors
array::Scalar m_ice_shelves
void compute_lakes(const array::CellType &cell_type, array::Scalar &result)
void update(const array::Scalar &bed_elevation, const array::CellType1 &cell_type)
void compute_continental_shelf_mask(const array::Scalar &bed_elevation, const array::Scalar &ice_rise_mask, double bed_elevation_threshold, array::Scalar &result)
const array::Scalar & basin_mask() const
const array::Scalar & box_mask() const
array::Scalar1 m_distance_cf
void identify_calving_front_connection(const array::CellType1 &cell_type, const array::Scalar &basin_mask, const array::Scalar &shelf_mask, int n_shelves, std::vector< int > &most_shelf_cells_in_basin, std::vector< int > &cfs_in_basins_per_shelf)
#define PISM_ERROR_LOCATION
double max(const array::Scalar &input)
Finds maximum over all the values in an array::Scalar object. Ignores ghosts.
void label_isolated(array::Scalar1 &mask, int reachable)
void label(array::Scalar1 &mask)
bool domain_edge(const Grid &grid, int i, int j)
bool ice_free_ocean(int M)
bool ocean(int M)
An ocean cell (floating ice or ice-free).
void eikonal_equation(array::Scalar1 &mask)
static void relabel(RelabelingType type, double threshold, array::Scalar &mask)
void GlobalMax(MPI_Comm comm, double *local, double *result, int count)
std::string printf(const char *format,...)
std::string join(const std::vector< std::string > &strings, const std::string &separator)
Concatenate strings, inserting separator between elements.
void GlobalSum(MPI_Comm comm, double *local, double *result, int count)