tFrontalMelt.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
---
tFrontalMelt.cc (8461B)
---
1 /* Copyright (C) 2018, 2019 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
30 FrontalMeltInputs::FrontalMeltInputs() {
31 geometry = nullptr;
32
33 subglacial_water_flux = nullptr;
34 }
35
36 namespace frontalmelt {
37
38 IceModelVec2S::Ptr FrontalMelt::allocate_frontal_melt_rate(IceGrid::ConstPtr g,
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 */
60 void FrontalMelt::compute_retreat_rate(const Geometry &geometry,
61 const IceModelVec2S &frontal_melt_rate,
62 IceModelVec2S &result) const {
63
64 GeometryCalculator gc(*m_config);
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.int_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), m_input_model(input) {
116
117 m_retreat_rate.create(m_grid, "retreat_rate_due_to_frontal_melt", WITHOUT_GHOSTS);
118 m_retreat_rate.set_attrs("diagnostic", "retreat rate due to frontal melt",
119 "m s-1", "m day-1", "", 0);
120
121 m_include_floating_ice = m_config->get_flag("frontal_melt.include_floating_ice");
122 }
123
124 // "model" constructor
125 FrontalMelt::FrontalMelt(IceGrid::ConstPtr g)
126 : FrontalMelt(g, nullptr) {
127 // empty
128 }
129
130 FrontalMelt::~FrontalMelt() {
131 // empty
132 }
133
134 void FrontalMelt::init(const Geometry &geometry) {
135 this->init_impl(geometry);
136 }
137
138 void FrontalMelt::init_impl(const Geometry &geometry) {
139 if (m_input_model) {
140 m_input_model->init(geometry);
141 }
142 }
143
144 void FrontalMelt::update(const FrontalMeltInputs &inputs, double t, double dt) {
145 this->update_impl(inputs, t, dt);
146
147 compute_retreat_rate(*inputs.geometry, frontal_melt_rate(), m_retreat_rate);
148 }
149
150 const IceModelVec2S& FrontalMelt::frontal_melt_rate() const {
151 return frontal_melt_rate_impl();
152 }
153
154 const IceModelVec2S& FrontalMelt::retreat_rate() const {
155 return m_retreat_rate;
156 }
157
158 // pass-through default implementations for "modifiers"
159
160 void FrontalMelt::update_impl(const FrontalMeltInputs &inputs, double t, double dt) {
161 if (m_input_model) {
162 m_input_model->update(inputs, t, dt);
163 } else {
164 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "no input model");
165 }
166 }
167
168 MaxTimestep FrontalMelt::max_timestep_impl(double t) const {
169 if (m_input_model) {
170 return m_input_model->max_timestep(t);
171 } else {
172 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "no input model");
173 }
174 }
175
176 void FrontalMelt::define_model_state_impl(const File &output) const {
177 if (m_input_model) {
178 return m_input_model->define_model_state(output);
179 } else {
180 // no state to define
181 }
182 }
183
184 void FrontalMelt::write_model_state_impl(const File &output) const {
185 if (m_input_model) {
186 return m_input_model->write_model_state(output);
187 } else {
188 // no state to write
189 }
190 }
191
192 namespace diagnostics {
193
194 /*! @brief Report frontal melt rate. */
195 class FrontalMeltRate : public DiagAverageRate<FrontalMelt>
196 {
197 public:
198 FrontalMeltRate(const FrontalMelt *m)
199 : DiagAverageRate<FrontalMelt>(m, "frontal_melt_rate", RATE) {
200
201 m_vars = {SpatialVariableMetadata(m_sys, "frontal_melt_rate")};
202 m_accumulator.metadata().set_string("units", "m");
203
204 set_attrs("frontal melt rate", "",
205 "m second-1", "m day-1", 0);
206 m_vars[0].set_string("cell_methods", "time: mean");
207
208 m_vars[0].set_number("_FillValue", to_internal(m_fill_value));
209 }
210
211 protected:
212 const IceModelVec2S& model_input() {
213 return model->frontal_melt_rate();
214 }
215 };
216
217 /*! @brief Report retreat rate due to frontal melt. */
218 class FrontalMeltRetreatRate : public DiagAverageRate<FrontalMelt>
219 {
220 public:
221 FrontalMeltRetreatRate(const FrontalMelt *m)
222 : DiagAverageRate<FrontalMelt>(m, "frontal_melt_retreat_rate", RATE) {
223
224 m_vars = {SpatialVariableMetadata(m_sys, "frontal_melt_retreat_rate")};
225 m_accumulator.metadata().set_string("units", "m");
226
227 set_attrs("retreat rate due to frontal melt", "",
228 "m second-1", "m year-1", 0);
229 m_vars[0].set_string("cell_methods", "time: mean");
230
231 m_vars[0].set_number("_FillValue", to_internal(m_fill_value));
232 m_vars[0].set_string("comment", "takes into account what part of the front is submerged");
233 }
234
235 protected:
236 const IceModelVec2S& model_input() {
237 return model->retreat_rate();
238 }
239 };
240
241 } // end of namespace diagnostics
242
243 DiagnosticList FrontalMelt::diagnostics_impl() const {
244 using namespace diagnostics;
245 DiagnosticList result = {
246 {"frontal_melt_rate", Diagnostic::Ptr(new FrontalMeltRate(this))},
247 {"frontal_melt_retreat_rate", Diagnostic::Ptr(new FrontalMeltRetreatRate(this))}
248 };
249
250 if (m_input_model) {
251 return combine(m_input_model->diagnostics(), result);
252 } else {
253 return result;
254 }
255 }
256
257 TSDiagnosticList FrontalMelt::ts_diagnostics_impl() const {
258 if (m_input_model) {
259 return m_input_model->ts_diagnostics();
260 } else {
261 return {};
262 }
263 }
264
265 bool FrontalMelt::apply(const IceModelVec2CellType &M, int i, int j) {
266 // icy and grounded_ice cells are included for visualization only (values at these
267 // locations have no effect)
268 if (m_include_floating_ice) {
269 return (M.ice_free_ocean(i, j) and M.next_to_ice(i, j)) or M.icy(i, j);
270 } else {
271 return (M.ice_free_ocean(i, j) and M.next_to_grounded_ice(i, j)) or M.grounded_ice(i, j);
272 }
273 }
274
275
276 } // end of namespace frontalmelt
277 } // end of namespace pism