PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
remove_narrow_tongues.cc
Go to the documentation of this file.
1 /* Copyright (C) 2016, 2017, 2018, 2019, 2020, 2022, 2023 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 "pism/frontretreat/util/remove_narrow_tongues.hh"
21 
22 #include "pism/util/Grid.hh"
23 #include "pism/geometry/Geometry.hh"
24 
25 namespace pism {
26 
27 /** Remove tips of one-cell-wide ice tongues ("noses")..
28  *
29  * The center icy cell in ice tongues like this one (and equivalent)
30  *
31  * @code
32  O O ?
33  X X O
34  O O ?
35  @endcode
36  *
37  * where "O" is ice-free and "?" is any mask value, are removed.
38  * Ice tongues like this one
39  *
40  * @code
41  # O ?
42  X X O
43  # O ?
44  @endcode
45  * where one or two of the "#" cells are ice-filled, are not removed.
46  *
47  * See the code for the precise rule, which uses `ice_free_ocean()` for the "O"
48  * cells if the center cell has grounded ice, and uses `ice_free()` if the
49  * center cell has floating ice.
50  *
51  * @note We use `geometry.cell_type` (and not `ice_thickness`) to make decisions. This
52  * means that we can update `ice_thickness` in place without introducing a dependence on
53  * the grid traversal order.
54  *
55  * @param[in,out] mask cell type mask
56  * @param[in,out] ice_thickness modeled ice thickness
57  *
58  * @return 0 on success
59  */
60 void remove_narrow_tongues(const Geometry &geometry,
61  array::Scalar &ice_thickness) {
62 
63  const auto &mask = geometry.cell_type;
64  const auto &bed = geometry.bed_elevation;
65  const auto &sea_level = geometry.sea_level_elevation;
66 
67  auto grid = mask.grid();
68 
69  array::AccessScope list{&mask, &bed, &sea_level, &ice_thickness};
70 
71  for (auto p = grid->points(); p; p.next()) {
72  const int i = p.i(), j = p.j();
73  if (mask.ice_free(i, j) or
74  (mask.grounded_ice(i, j) and bed(i, j) >= sea_level(i, j))) {
75  continue;
76  }
77 
79  auto M = mask.box_int(i, j);
80 
81  // Note: i,j cannot be ice-free (see the if-block above), so it is either grounded ice
82  // or floating ice
83  if (mask::grounded_ice(M.c)) {
85  // if (i,j) is grounded ice then we will remove it if it has
86  // exclusively ice-free ocean neighbors
87  ice_free.n = ice_free_ocean(M.n);
88  ice_free.e = ice_free_ocean(M.e);
89  ice_free.s = ice_free_ocean(M.s);
90  ice_free.w = ice_free_ocean(M.w);
91  ice_free.ne = ice_free_ocean(M.ne);
92  ice_free.nw = ice_free_ocean(M.nw);
93  ice_free.se = ice_free_ocean(M.se);
94  ice_free.sw = ice_free_ocean(M.sw);
95  } else {
96  // if (i,j) is floating then we will remove it if its neighbors are
97  // ice-free, whether ice-free ocean or ice-free ground
98  ice_free.n = mask::ice_free(M.n);
99  ice_free.e = mask::ice_free(M.e);
100  ice_free.s = mask::ice_free(M.s);
101  ice_free.w = mask::ice_free(M.w);
102  ice_free.ne = mask::ice_free(M.ne);
103  ice_free.nw = mask::ice_free(M.nw);
104  ice_free.se = mask::ice_free(M.se);
105  ice_free.sw = mask::ice_free(M.sw);
106  }
107 
108  if ((not ice_free.w and
109  ice_free.nw and
110  ice_free.sw and
111  ice_free.n and
112  ice_free.s and
113  ice_free.e) or
114  (not ice_free.n and
115  ice_free.nw and
116  ice_free.ne and
117  ice_free.w and
118  ice_free.e and
119  ice_free.s) or
120  (not ice_free.e and
121  ice_free.ne and
122  ice_free.se and
123  ice_free.w and
124  ice_free.s and
125  ice_free.n) or
126  (not ice_free.s and
127  ice_free.sw and
128  ice_free.se and
129  ice_free.w and
130  ice_free.e and
131  ice_free.n)) {
132  ice_thickness(i, j) = 0.0;
133  }
134  }
135 }
136 
137 } // end of namespace pism
array::Scalar1 sea_level_elevation
Definition: Geometry.hh:48
array::CellType2 cell_type
Definition: Geometry.hh:55
array::Scalar2 bed_elevation
Definition: Geometry.hh:47
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
Definition: Array.hh:65
std::shared_ptr< const Grid > grid() const
Definition: Array.cc:132
bool grounded_ice(int M)
Definition: Mask.hh:51
bool ice_free_ocean(int M)
Definition: Mask.hh:61
bool ice_free(int M)
Ice-free cell (grounded or ocean).
Definition: Mask.hh:58
void remove_narrow_tongues(const Geometry &geometry, array::Scalar &ice_thickness)