PISM, A Parallel Ice Sheet Model  stable v2.0.4 committed by Constantine Khrulev on 2022-05-25 12:02: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 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'.
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  const Direction dirs[] = {North, East, South, West};
44  for (int n = 0; n < 4; ++n) {
45  Direction direction = dirs[n];
46  if (mask::icy(M[direction])) {
47  H_average += H[direction];
48  h_average += h[direction];
49  N++;
50  }
51  }
52 
53  if (N == 0) {
54  // If there are no "icy" neighbors, return the threshold thickness
55  // of zero, forcing Href to be converted to H immediately.
56  return 0.0;
57  }
58 
59  H_average = H_average / N;
60  h_average = h_average / N;
61 
62  if (bed_elevation + H_average > h_average) {
63  H_threshold = h_average - bed_elevation;
64  } else {
65  H_threshold = H_average;
66  }
67 
68  // make sure that the returned threshold thickness is non-negative:
69  return std::max(H_threshold, 0.0);
70 }
71 
72 } // end of namespace pism
#define n
Definition: exactTestM.c:37
bool icy(int M)
Ice-filled cell (grounded or floating).
Definition: Mask.hh:47
double max(const IceModelVec2S &input)
Finds maximum over all the values in an IceModelVec2S object. Ignores ghosts.
Direction
Definition: stencils.hh:24
@ North
Definition: stencils.hh:24
@ East
Definition: stencils.hh:24
@ South
Definition: stencils.hh:24
@ West
Definition: stencils.hh:24
double part_grid_threshold_thickness(stencils::Star< int > M, stencils::Star< double > H, stencils::Star< double > h, double bed_elevation)
Compute threshold thickness used when deciding if a partially-filled cell should be considered 'full'...
Star stencil points (in the map-plane).
Definition: stencils.hh:30