tremove_narrow_tongues.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
---
tremove_narrow_tongues.cc (4450B)
---
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 "remove_narrow_tongues.hh"
21
22 #include "pism/util/IceGrid.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 `pism_mask` (and not ice_thickness) to make decisions.
52 * This means that we can update `ice_thickness` in place without
53 * introducing a dependence on 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 IceModelVec2S &ice_thickness) {
62
63 auto &mask = geometry.cell_type;
64 auto &bed = geometry.bed_elevation;
65 auto &sea_level = geometry.sea_level_elevation;
66
67 IceGrid::ConstPtr grid = mask.grid();
68
69 IceModelVec::AccessList list{&mask, &bed, &sea_level, &ice_thickness};
70
71 for (Points p(*grid); 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
78 bool
79 ice_free_N = false, ice_free_E = false,
80 ice_free_S = false, ice_free_W = false,
81 ice_free_NE = false, ice_free_NW = false,
82 ice_free_SE = false, ice_free_SW = false;
83
84 if (mask.grounded_ice(i,j)) {
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 = mask.ice_free_ocean(i, j + 1);
88 ice_free_E = mask.ice_free_ocean(i + 1, j);
89 ice_free_S = mask.ice_free_ocean(i, j - 1);
90 ice_free_W = mask.ice_free_ocean(i - 1, j);
91 ice_free_NE = mask.ice_free_ocean(i + 1, j + 1);
92 ice_free_NW = mask.ice_free_ocean(i - 1, j + 1);
93 ice_free_SE = mask.ice_free_ocean(i + 1, j - 1);
94 ice_free_SW = mask.ice_free_ocean(i - 1, j - 1);
95 } else if (mask.floating_ice(i,j)) {
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(i, j + 1);
99 ice_free_E = mask.ice_free(i + 1, j);
100 ice_free_S = mask.ice_free(i, j - 1);
101 ice_free_W = mask.ice_free(i - 1, j);
102 ice_free_NE = mask.ice_free(i + 1, j + 1);
103 ice_free_NW = mask.ice_free(i - 1, j + 1);
104 ice_free_SE = mask.ice_free(i + 1, j - 1);
105 ice_free_SW = mask.ice_free(i - 1, j - 1);
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