tpart_grid_threshold_thickness.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
       ---
       tpart_grid_threshold_thickness.cc (2317B)
       ---
            1 // Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018 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'.
           31 double part_grid_threshold_thickness(StarStencil<int> M,
           32                                      StarStencil<double> H,
           33                                      StarStencil<double> h,
           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