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