tGeometry.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
(HTM) git clone git://src.adamsgaard.dk/pism
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tGeometry.cc (13528B)
---
1 /* Copyright (C) 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 "Geometry.hh"
21
22 #include "pism/util/iceModelVec.hh"
23 #include "pism/util/IceModelVec2CellType.hh"
24 #include "pism/util/Mask.hh"
25 #include "pism/util/pism_utilities.hh"
26 #include "pism/geometry/grounded_cell_fraction.hh"
27
28 namespace pism {
29
30 Geometry::Geometry(IceGrid::ConstPtr grid)
31 // FIXME: ideally these fields should be "global", i.e. without ghosts.
32 // (However this may increase communication costs...)
33 : m_stencil_width(grid->ctx()->config()->get_number("grid.max_stencil_width")),
34 latitude(grid, "lat", WITHOUT_GHOSTS),
35 longitude(grid, "lon", WITHOUT_GHOSTS),
36 bed_elevation(grid, "topg", WITH_GHOSTS, m_stencil_width),
37 sea_level_elevation(grid, "sea_level", WITH_GHOSTS),
38 ice_thickness(grid, "thk", WITH_GHOSTS, m_stencil_width),
39 ice_area_specific_volume(grid, "ice_area_specific_volume", WITH_GHOSTS),
40 cell_type(grid, "mask", WITH_GHOSTS, m_stencil_width),
41 cell_grounded_fraction(grid, "cell_grounded_fraction", WITHOUT_GHOSTS),
42 ice_surface_elevation(grid, "usurf", WITH_GHOSTS, m_stencil_width) {
43
44 latitude.set_attrs("mapping", "latitude", "degree_north", "degree_north", "latitude", 0);
45 latitude.set_time_independent(true);
46 latitude.metadata().set_string("grid_mapping", "");
47 latitude.metadata().set_numbers("valid_range", {-90.0, 90.0});
48
49 longitude.set_attrs("mapping", "longitude", "degree_east", "degree_east", "longitude", 0);
50 longitude.set_time_independent(true);
51 longitude.metadata().set_string("grid_mapping", "");
52 longitude.metadata().set_numbers("valid_range", {-180.0, 180.0});
53
54 bed_elevation.set_attrs("model_state", "bedrock surface elevation",
55 "m", "m", "bedrock_altitude", 0);
56
57 sea_level_elevation.set_attrs("model_state",
58 "sea level elevation above reference ellipsoid", "meters", "meters",
59 "sea_surface_height_above_reference_ellipsoid", 0);
60
61 ice_thickness.set_attrs("model_state", "land ice thickness",
62 "m", "m", "land_ice_thickness", 0);
63 ice_thickness.metadata().set_number("valid_min", 0.0);
64
65 ice_area_specific_volume.set_attrs("model_state",
66 "ice-volume-per-area in partially-filled grid cells",
67 "m3/m2", "m3/m2", "", 0);
68 ice_area_specific_volume.metadata().set_string("comment",
69 "this variable represents the amount of ice "
70 "in a partially-filled cell and not "
71 "the corresponding geometry, so thinking "
72 "about it as 'thickness' is not helpful");
73
74 cell_type.set_attrs("diagnostic", "ice-type (ice-free/grounded/floating/ocean) integer mask",
75 "", "", "", 0);
76 cell_type.metadata().set_numbers("flag_values", {MASK_ICE_FREE_BEDROCK, MASK_GROUNDED,
77 MASK_FLOATING, MASK_ICE_FREE_OCEAN});
78 cell_type.metadata().set_string("flag_meanings",
79 "ice_free_bedrock grounded_ice floating_ice ice_free_ocean");
80 cell_type.metadata().set_output_type(PISM_INT);
81
82 cell_grounded_fraction.set_attrs("internal",
83 "fractional grounded/floating mask (floating=0, grounded=1)",
84 "", "", "", 0);
85
86 ice_surface_elevation.set_attrs("diagnostic", "ice upper surface elevation",
87 "m", "m", "surface_altitude", 0);
88
89 // make sure all the fields are initialized
90 latitude.set(0.0);
91 longitude.set(0.0);
92 bed_elevation.set(0.0);
93 sea_level_elevation.set(0.0);
94 ice_thickness.set(0.0);
95 ice_area_specific_volume.set(0.0);
96 ensure_consistency(0.0);
97 }
98
99 void check_minimum_ice_thickness(const IceModelVec2S &ice_thickness) {
100 IceGrid::ConstPtr grid = ice_thickness.grid();
101
102 IceModelVec::AccessList list(ice_thickness);
103
104 ParallelSection loop(grid->com);
105 try {
106 for (Points p(*grid); p; p.next()) {
107 const int i = p.i(), j = p.j();
108
109 if (ice_thickness(i, j) < 0.0) {
110 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
111 "H = %e (negative) at point i=%d, j=%d",
112 ice_thickness(i, j), i, j);
113 }
114 }
115 } catch (...) {
116 loop.failed();
117 }
118 loop.check();
119 }
120
121 void Geometry::ensure_consistency(double ice_free_thickness_threshold) {
122 IceGrid::ConstPtr grid = ice_thickness.grid();
123 Config::ConstPtr config = grid->ctx()->config();
124
125 check_minimum_ice_thickness(ice_thickness);
126
127 IceModelVec::AccessList list{&sea_level_elevation, &bed_elevation,
128 &ice_thickness, &ice_area_specific_volume,
129 &cell_type, &ice_surface_elevation};
130
131 // first ensure that ice_area_specific_volume is 0 if ice_thickness > 0.
132 {
133 ParallelSection loop(grid->com);
134 try {
135 for (Points p(*grid); p; p.next()) {
136 const int i = p.i(), j = p.j();
137
138 if (ice_thickness(i, j) > 0.0 and ice_area_specific_volume(i, j) > 0.0) {
139 ice_thickness(i, j) += ice_area_specific_volume(i, j);
140 ice_area_specific_volume(i, j) = 0.0;
141 }
142 }
143 } catch (...) {
144 loop.failed();
145 }
146 loop.check();
147 }
148
149 // compute cell type and surface elevation
150 {
151 GeometryCalculator gc(*config);
152 gc.set_icefree_thickness(ice_free_thickness_threshold);
153
154 ParallelSection loop(grid->com);
155 try {
156 for (Points p(*grid); p; p.next()) {
157 const int i = p.i(), j = p.j();
158
159 int mask = 0;
160 gc.compute(sea_level_elevation(i, j), bed_elevation(i, j), ice_thickness(i, j),
161 &mask, &ice_surface_elevation(i, j));
162 cell_type(i, j) = mask;
163 }
164 } catch (...) {
165 loop.failed();
166 }
167 loop.check();
168 }
169
170 ice_thickness.update_ghosts();
171 ice_area_specific_volume.update_ghosts();
172 cell_type.update_ghosts();
173 ice_surface_elevation.update_ghosts();
174
175 const double
176 ice_density = config->get_number("constants.ice.density"),
177 ocean_density = config->get_number("constants.sea_water.density");
178
179 compute_grounded_cell_fraction(ice_density,
180 ocean_density,
181 sea_level_elevation,
182 ice_thickness,
183 bed_elevation,
184 cell_grounded_fraction);
185 }
186
187 /*! Compute the elevation of the bottom surface of the ice.
188 */
189 void ice_bottom_surface(const Geometry &geometry, IceModelVec2S &result) {
190
191 auto grid = result.grid();
192 auto config = grid->ctx()->config();
193
194 double
195 ice_density = config->get_number("constants.ice.density"),
196 water_density = config->get_number("constants.sea_water.density"),
197 alpha = ice_density / water_density;
198
199 const IceModelVec2S &ice_thickness = geometry.ice_thickness;
200 const IceModelVec2S &bed_elevation = geometry.bed_elevation;
201 const IceModelVec2S &sea_level = geometry.sea_level_elevation;
202
203 IceModelVec::AccessList list{&ice_thickness, &bed_elevation, &sea_level, &result};
204
205 ParallelSection loop(grid->com);
206 try {
207 for (Points p(*grid); p; p.next()) {
208 const int i = p.i(), j = p.j();
209
210 double
211 b_grounded = bed_elevation(i, j),
212 b_floating = sea_level(i, j) - alpha * ice_thickness(i, j);
213
214 result(i, j) = std::max(b_grounded, b_floating);
215 }
216 } catch (...) {
217 loop.failed();
218 }
219 loop.check();
220
221 result.update_ghosts();
222 }
223
224 //! Computes the ice volume, in m^3.
225 double ice_volume(const Geometry &geometry, double thickness_threshold) {
226 auto grid = geometry.ice_thickness.grid();
227 auto config = grid->ctx()->config();
228
229 IceModelVec::AccessList list{&geometry.ice_thickness};
230
231 double volume = 0.0;
232
233 auto cell_area = grid->cell_area();
234
235 {
236 for (Points p(*grid); p; p.next()) {
237 const int i = p.i(), j = p.j();
238
239 if (geometry.ice_thickness(i,j) >= thickness_threshold) {
240 volume += geometry.ice_thickness(i,j) * cell_area;
241 }
242 }
243 }
244
245 // Add the volume of the ice in Href:
246 if (config->get_flag("geometry.part_grid.enabled")) {
247 list.add(geometry.ice_area_specific_volume);
248 for (Points p(*grid); p; p.next()) {
249 const int i = p.i(), j = p.j();
250
251 volume += geometry.ice_area_specific_volume(i,j) * cell_area;
252 }
253 }
254
255 return GlobalSum(grid->com, volume);
256 }
257
258 double ice_volume_not_displacing_seawater(const Geometry &geometry,
259 double thickness_threshold) {
260 auto grid = geometry.ice_thickness.grid();
261 auto config = grid->ctx()->config();
262
263 const double
264 sea_water_density = config->get_number("constants.sea_water.density"),
265 ice_density = config->get_number("constants.ice.density"),
266 cell_area = grid->cell_area();
267
268 IceModelVec::AccessList list{&geometry.cell_type, &geometry.ice_thickness,
269 &geometry.bed_elevation, &geometry.sea_level_elevation};
270
271 double volume = 0.0;
272
273 for (Points p(*grid); p; p.next()) {
274 const int i = p.i(), j = p.j();
275
276 const double
277 bed = geometry.bed_elevation(i, j),
278 thickness = geometry.ice_thickness(i, j),
279 sea_level = geometry.sea_level_elevation(i, j);
280
281 if (geometry.cell_type.grounded(i, j) and thickness > thickness_threshold) {
282 const double cell_ice_volume = thickness * cell_area;
283 if (bed > sea_level) {
284 volume += cell_ice_volume;
285 } else {
286 const double max_floating_volume = (sea_level - bed) * cell_area * (sea_water_density / ice_density);
287 volume += cell_ice_volume - max_floating_volume;
288 }
289 }
290 } // end of the loop over grid points
291
292 return GlobalSum(grid->com, volume);
293 }
294
295 //! Computes ice area, in m^2.
296 double ice_area(const Geometry &geometry, double thickness_threshold) {
297 auto grid = geometry.ice_thickness.grid();
298
299 double area = 0.0;
300
301 auto cell_area = grid->cell_area();
302
303 IceModelVec::AccessList list{&geometry.ice_thickness};
304 for (Points p(*grid); p; p.next()) {
305 const int i = p.i(), j = p.j();
306
307 if (geometry.ice_thickness(i, j) >= thickness_threshold) {
308 area += cell_area;
309 }
310 }
311
312 return GlobalSum(grid->com, area);
313 }
314
315 //! Computes grounded ice area, in m^2.
316 double ice_area_grounded(const Geometry &geometry, double thickness_threshold) {
317 auto grid = geometry.ice_thickness.grid();
318
319 double area = 0.0;
320
321 auto cell_area = grid->cell_area();
322
323 IceModelVec::AccessList list{&geometry.cell_type, &geometry.ice_thickness};
324 for (Points p(*grid); p; p.next()) {
325 const int i = p.i(), j = p.j();
326
327 if (geometry.cell_type.grounded(i, j) and
328 geometry.ice_thickness(i, j) >= thickness_threshold) {
329 area += cell_area;
330 }
331 }
332
333 return GlobalSum(grid->com, area);
334 }
335
336 //! Computes floating ice area, in m^2.
337 double ice_area_floating(const Geometry &geometry, double thickness_threshold) {
338 auto grid = geometry.ice_thickness.grid();
339
340 double area = 0.0;
341
342 auto cell_area = grid->cell_area();
343
344 IceModelVec::AccessList list{&geometry.cell_type, &geometry.ice_thickness};
345 for (Points p(*grid); p; p.next()) {
346 const int i = p.i(), j = p.j();
347
348 if (geometry.cell_type.ocean(i, j) and
349 geometry.ice_thickness(i, j) >= thickness_threshold) {
350 area += cell_area;
351 }
352 }
353
354 return GlobalSum(grid->com, area);
355 }
356
357
358 //! Computes the sea level rise that would result if all the ice were melted.
359 double sea_level_rise_potential(const Geometry &geometry, double thickness_threshold) {
360 auto config = geometry.ice_thickness.grid()->ctx()->config();
361
362 const double
363 water_density = config->get_number("constants.fresh_water.density"),
364 ice_density = config->get_number("constants.ice.density"),
365 ocean_area = config->get_number("constants.global_ocean_area");
366
367 const double
368 volume = ice_volume_not_displacing_seawater(geometry,
369 thickness_threshold),
370 additional_water_volume = (ice_density / water_density) * volume,
371 sea_level_change = additional_water_volume / ocean_area;
372
373 return sea_level_change;
374 }
375
376
377 /*!
378 * @brief Set no_model_mask variable to have value 1 in strip of width 'strip' m around
379 * edge of computational domain, and value 0 otherwise.
380 */
381 void set_no_model_strip(const IceGrid &grid, double width, IceModelVec2Int &result) {
382
383 if (width <= 0.0) {
384 return;
385 }
386
387 IceModelVec::AccessList list(result);
388
389 for (Points p(grid); p; p.next()) {
390 const int i = p.i(), j = p.j();
391
392 if (in_null_strip(grid, i, j, width)) {
393 result(i, j) = 1;
394 } else {
395 result(i, j) = 0;
396 }
397 }
398
399 result.update_ghosts();
400 }
401
402 } // end of namespace pism