PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
FrontalMelt.cc
Go to the documentation of this file.
1 /* Copyright (C) 2018, 2019, 2020, 2022, 2023 Constantine Khroulev and Andy Aschwanden
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 "pism/coupler/FrontalMelt.hh"
21 #include "pism/util/MaxTimestep.hh"
22 #include "pism/util/pism_utilities.hh" // combine()
23 #include "pism/geometry/Geometry.hh"
24 #include "pism/geometry/part_grid_threshold_thickness.hh"
25 #include "pism/util/Mask.hh" // GeometryCalculator
26 
27 namespace pism {
28 
30  geometry = nullptr;
31 
32  subglacial_water_flux = nullptr;
33 }
34 
35 namespace frontalmelt {
36 
37 /*!
38  * Compute retreat rate corresponding to a given frontal melt rate.
39  *
40  * This code computes the fraction of the front submerged in ocean water and uses it to
41  * scale the provided melt rate.
42  */
44  const array::Scalar &frontal_melt_rate,
45  array::Scalar &result) const {
46 
48 
49  const array::Scalar2
50  &bed_elevation = geometry.bed_elevation;
51  const array::Scalar1
52  &sea_level_elevation = geometry.sea_level_elevation,
53  &surface_elevation = geometry.ice_surface_elevation,
54  &ice_thickness = geometry.ice_thickness;
55  const array::CellType1 &cell_type = geometry.cell_type;
56 
57  const double
58  ice_density = m_config->get_number("constants.ice.density"),
59  alpha = ice_density / m_config->get_number("constants.sea_water.density");
60 
61  array::AccessScope list{&cell_type, &frontal_melt_rate, &sea_level_elevation,
62  &bed_elevation, &surface_elevation, &ice_thickness, &result};
63 
64  ParallelSection loop(m_grid->com);
65  try {
66  for (auto p = m_grid->points(); p; p.next()) {
67  const int i = p.i(), j = p.j();
68 
69  if (cell_type.ice_free_ocean(i, j) and cell_type.next_to_ice(i, j)) {
70  const double
71  bed = bed_elevation(i, j),
72  sea_level = sea_level_elevation(i, j);
73 
74  auto H = ice_thickness.star(i, j);
75  auto h = surface_elevation.star(i, j);
76  auto M = cell_type.star_int(i, j);
77 
78  double H_threshold = part_grid_threshold_thickness(M, H, h, bed);
79 
80  int m = gc.mask(sea_level, bed, H_threshold);
81 
82  double H_submerged = (mask::grounded(m) ?
83  std::max(sea_level - bed, 0.0) :
84  alpha * H_threshold);
85 
86  result(i, j) = (H_submerged / H_threshold) * frontal_melt_rate(i, j);
87  } else {
88  result(i, j) = 0.0;
89  }
90  }
91  } catch (...) {
92  loop.failed();
93  }
94  loop.check();
95 }
96 
97 // "modifier" constructor
98 FrontalMelt::FrontalMelt(std::shared_ptr<const Grid> g, std::shared_ptr<FrontalMelt> input)
99  : Component(g),
100  m_input_model(input),
101  m_retreat_rate(m_grid, "retreat_rate_due_to_frontal_melt") {
102 
104  .long_name("retreat rate due to frontal melt")
105  .units("m s-1")
106  .output_units("m day-1");
107 
108  m_include_floating_ice = m_config->get_flag("frontal_melt.include_floating_ice");
109 }
110 
111 // "model" constructor
112 FrontalMelt::FrontalMelt(std::shared_ptr<const Grid> g)
113  : FrontalMelt(g, nullptr) {
114  // empty
115 }
116 
117 void FrontalMelt::init(const Geometry &geometry) {
118  this->init_impl(geometry);
119 }
120 
121 void FrontalMelt::init_impl(const Geometry &geometry) {
122  if (m_input_model) {
123  m_input_model->init(geometry);
124  }
125 }
126 
127 void FrontalMelt::update(const FrontalMeltInputs &inputs, double t, double dt) {
128  this->update_impl(inputs, t, dt);
129 
131 }
132 
134  return frontal_melt_rate_impl();
135 }
136 
138  return m_retreat_rate;
139 }
140 
141 // pass-through default implementations for "modifiers"
142 
143 void FrontalMelt::update_impl(const FrontalMeltInputs &inputs, double t, double dt) {
144  if (m_input_model) {
145  m_input_model->update(inputs, t, dt);
146  } else {
147  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "no input model");
148  }
149 }
150 
152  if (m_input_model) {
153  return m_input_model->max_timestep(t);
154  } else {
155  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "no input model");
156  }
157 }
158 
159 void FrontalMelt::define_model_state_impl(const File &output) const {
160  if (m_input_model) {
161  return m_input_model->define_model_state(output);
162  } else {
163  // no state to define
164  }
165 }
166 
167 void FrontalMelt::write_model_state_impl(const File &output) const {
168  if (m_input_model) {
169  return m_input_model->write_model_state(output);
170  }
171  // no state to write
172 }
173 
174 namespace diagnostics {
175 
176 /*! @brief Report frontal melt rate. */
177 class FrontalMeltRate : public DiagAverageRate<FrontalMelt>
178 {
179 public:
181  : DiagAverageRate<FrontalMelt>(m, "frontal_melt_rate", RATE) {
182 
183  m_accumulator.metadata()["units"] = "m";
184 
185  m_vars = { { m_sys, "frontal_melt_rate" } };
186  m_vars[0]
187  .long_name("frontal melt rate")
188  .units("m second-1")
189  .output_units("m day-1");
190  m_vars[0]["cell_methods"] = "time: mean";
191  m_vars[0]["_FillValue"] = {to_internal(m_fill_value)};
192  }
193 
194 protected:
196  return model->frontal_melt_rate();
197  }
198 };
199 
200 /*! @brief Report retreat rate due to frontal melt. */
201 class FrontalMeltRetreatRate : public DiagAverageRate<FrontalMelt>
202 {
203 public:
205  : DiagAverageRate<FrontalMelt>(m, "frontal_melt_retreat_rate", RATE) {
206  m_accumulator.metadata()["units"] = "m";
207 
208  m_vars = { { m_sys, "frontal_melt_retreat_rate" } };
209  m_vars[0]
210  .long_name("retreat rate due to frontal melt")
211  .units("m second-1")
212  .output_units("m year-1");
213  m_vars[0]["cell_methods"] = "time: mean";
214  m_vars[0]["_FillValue"] = {to_internal(m_fill_value)};
215  m_vars[0]["comment"] = "takes into account what part of the front is submerged";
216  }
217 
218 protected:
220  return model->retreat_rate();
221  }
222 };
223 
224 } // end of namespace diagnostics
225 
227  using namespace diagnostics;
228  DiagnosticList result = {
229  {"frontal_melt_rate", Diagnostic::Ptr(new FrontalMeltRate(this))},
230  {"frontal_melt_retreat_rate", Diagnostic::Ptr(new FrontalMeltRetreatRate(this))}
231  };
232 
233  if (m_input_model) {
234  return combine(m_input_model->diagnostics(), result);
235  }
236  return result;
237 }
238 
240  if (m_input_model) {
241  return m_input_model->ts_diagnostics();
242  }
243  return {};
244 }
245 
246 bool FrontalMelt::apply(const array::CellType1 &M, int i, int j) const {
247  // icy and grounded_ice cells are included for visualization only (values at these
248  // locations have no effect)
250  return (M.ice_free_ocean(i, j) and M.next_to_ice(i, j)) or M.icy(i, j);
251  }
252 
253  return (M.ice_free_ocean(i, j) and M.next_to_grounded_ice(i, j)) or M.grounded_ice(i, j);
254 }
255 
256 
257 } // end of namespace frontalmelt
258 } // end of namespace pism
const Config::ConstPtr m_config
configuration database used by this component
Definition: Component.hh:158
const std::shared_ptr< const Grid > m_grid
grid used by this component
Definition: Component.hh:156
DiagnosticList diagnostics() const
Definition: Component.cc:89
A class defining a common interface for most PISM sub-models.
Definition: Component.hh:118
const Model * model
Definition: Diagnostic.hh:172
double m_fill_value
fill value (used often enough to justify storing it)
Definition: Diagnostic.hh:122
const units::System::Ptr m_sys
the unit system
Definition: Diagnostic.hh:116
double to_internal(double x) const
Definition: Diagnostic.cc:59
std::vector< SpatialVariableMetadata > m_vars
metadata corresponding to NetCDF variables
Definition: Diagnostic.hh:120
std::shared_ptr< Diagnostic > Ptr
Definition: Diagnostic.hh:65
High-level PISM I/O class.
Definition: File.hh:56
const Geometry * geometry
Definition: FrontalMelt.hh:32
const array::Scalar * subglacial_water_flux
Definition: FrontalMelt.hh:35
int mask(double sea_level, double bed, double thickness) const
Definition: Mask.hh:138
array::Scalar1 sea_level_elevation
Definition: Geometry.hh:48
array::Scalar2 ice_surface_elevation
Definition: Geometry.hh:57
array::CellType2 cell_type
Definition: Geometry.hh:55
array::Scalar2 ice_thickness
Definition: Geometry.hh:51
array::Scalar2 bed_elevation
Definition: Geometry.hh:47
Combines the max. time step with the flag indicating if a restriction is active. Makes is possible to...
Definition: MaxTimestep.hh:31
void failed()
Indicates a failure of a parallel section.
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
VariableMetadata & long_name(const std::string &input)
VariableMetadata & output_units(const std::string &input)
VariableMetadata & units(const std::string &input)
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
Definition: Array.hh:65
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition: Array.cc:553
bool next_to_ice(int i, int j) const
Ice-free margin (at least one of four neighbors has ice).
Definition: CellType.hh:87
stencils::Star< int > star_int(int i, int j) const
Definition: Scalar.hh:72
bool next_to_grounded_ice(int i, int j) const
Definition: CellType.hh:96
bool ice_free_ocean(int i, int j) const
Definition: CellType.hh:58
bool grounded_ice(int i, int j) const
Definition: CellType.hh:46
bool icy(int i, int j) const
Definition: CellType.hh:42
virtual const array::Scalar & frontal_melt_rate_impl() const =0
virtual MaxTimestep max_timestep_impl(double t) const
Definition: FrontalMelt.cc:151
void init(const Geometry &geometry)
Definition: FrontalMelt.cc:117
void update(const FrontalMeltInputs &inputs, double t, double dt)
Definition: FrontalMelt.cc:127
void compute_retreat_rate(const Geometry &geometry, const array::Scalar &frontal_melt_rate, array::Scalar &result) const
Definition: FrontalMelt.cc:43
std::shared_ptr< FrontalMelt > m_input_model
Definition: FrontalMelt.hh:77
virtual DiagnosticList diagnostics_impl() const
Definition: FrontalMelt.cc:226
virtual void update_impl(const FrontalMeltInputs &inputs, double t, double dt)
Definition: FrontalMelt.cc:143
const array::Scalar & retreat_rate() const
Definition: FrontalMelt.cc:137
bool apply(const array::CellType1 &M, int i, int j) const
Definition: FrontalMelt.cc:246
FrontalMelt(std::shared_ptr< const Grid > g, std::shared_ptr< FrontalMelt > input)
Definition: FrontalMelt.cc:98
virtual void define_model_state_impl(const File &output) const
The default (empty implementation).
Definition: FrontalMelt.cc:159
const array::Scalar & frontal_melt_rate() const
Definition: FrontalMelt.cc:133
virtual TSDiagnosticList ts_diagnostics_impl() const
Definition: FrontalMelt.cc:239
virtual void init_impl(const Geometry &geometry)
Definition: FrontalMelt.cc:121
virtual void write_model_state_impl(const File &output) const
The default (empty implementation).
Definition: FrontalMelt.cc:167
A very rudimentary PISM frontal melt model.
Definition: FrontalMelt.hh:43
Report retreat rate due to frontal melt.
Definition: FrontalMelt.cc:202
#define PISM_ERROR_LOCATION
double max(const array::Scalar &input)
Finds maximum over all the values in an array::Scalar object. Ignores ghosts.
Definition: Scalar.cc:165
bool grounded(int M)
Grounded cell (grounded ice or ice-free).
Definition: Mask.hh:44
static const double g
Definition: exactTestP.cc:36
std::map< std::string, TSDiagnostic::Ptr > TSDiagnosticList
Definition: Diagnostic.hh:343
std::map< std::string, Diagnostic::Ptr > DiagnosticList
Definition: Diagnostic.hh:125
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'...
T combine(const T &a, const T &b)