PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
part_grid_threshold_thickness.cc
Go to the documentation of this file.
1 // Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2021, 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 #include <cmath>
20 
21 #include "pism/geometry/part_grid_threshold_thickness.hh"
22 #include "pism/util/Mask.hh"
23 
24 namespace pism {
25 
26 
27 //! \file part_grid_threshold_thickness.cc Methods implementing PIK option -part_grid [\ref Albrechtetal2011].
28 
29 //! @brief Compute threshold thickness used when deciding if a
30 //! partially-filled cell should be considered 'full'.
32  stencils::Star<double> ice_thickness,
33  stencils::Star<double> surface_elevation,
34  double bed_elevation) {
35  // get mean ice thickness and surface elevation over adjacent
36  // icy cells
37  double
38  H_average = 0.0,
39  h_average = 0.0,
40  H_threshold = 0.0;
41  int N = 0;
42 
43  for (auto d : {North, East, South, West}) {
44  if (mask::icy(cell_type[d])) {
45  H_average += ice_thickness[d];
46  h_average += surface_elevation[d];
47  N++;
48  }
49  }
50 
51  if (N == 0) {
52  // If there are no "icy" neighbors, return the threshold thickness
53  // of zero, forcing Href to be converted to ice_thickness immediately.
54  return 0.0;
55  }
56 
57  H_average = H_average / N;
58  h_average = h_average / N;
59 
60  if (bed_elevation + H_average > h_average) {
61  H_threshold = h_average - bed_elevation;
62  } else {
63  H_threshold = H_average;
64  }
65 
66  // make sure that the returned threshold thickness is non-negative:
67  return std::max(H_threshold, 0.0);
68 }
69 
70 } // end of namespace pism
double max(const array::Scalar &input)
Finds maximum over all the values in an array::Scalar object. Ignores ghosts.
Definition: Scalar.cc:165
bool icy(int M)
Ice-filled cell (grounded or floating).
Definition: Mask.hh:48
double part_grid_threshold_thickness(stencils::Star< int > cell_type, stencils::Star< double > ice_thickness, stencils::Star< double > surface_elevation, double bed_elevation)
Compute threshold thickness used when deciding if a partially-filled cell should be considered 'full'...
@ North
Definition: stencils.hh:24
@ East
Definition: stencils.hh:24
@ South
Definition: stencils.hh:24
@ West
Definition: stencils.hh:24
Star stencil points (in the map-plane).
Definition: stencils.hh:30