tCHSystem.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
(HTM) git clone git://src.adamsgaard.dk/pism
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tCHSystem.cc (10400B)
---
1 /* Copyright (C) 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 <algorithm> // std::max
20
21 #include "CHSystem.hh"
22
23 #include "DrainageCalculator.hh"
24 #include "pism/util/EnthalpyConverter.hh"
25 #include "pism/energy/enthSystem.hh"
26 #include "pism/util/IceModelVec2CellType.hh"
27 #include "pism/util/io/File.hh"
28 #include "utilities.hh"
29 #include "pism/util/pism_utilities.hh"
30
31 namespace pism {
32 namespace energy {
33
34 /*!
35 * Model of the energy content in the cryo-hydrologic system.
36 *
37 * This model can be used to include the influence of cryo-hydrologic warming on the
38 * internal energy of the ice.
39 *
40 * This model is based on
41 *
42 * @article{PhillipsRajaramSteffan2010,
43 * author = {Phillips, T. and Rajaram, H. and Steffen, K.},
44 * doi = {10.1029/2010GL044397},
45 * journal = {Geophy. Res. Lett.},
46 * number = {20},
47 * pages = {1--5},
48 * title = {{Cryo-hydrologic warming: A potential mechanism for rapid thermal response of ice sheets}},
49 * volume = {37},
50 * year = {2010}
51 * }
52 *
53 * We assume that during the melt season the enthalpy in the cryo-hydrologic system is
54 * equal to the enthalpy of the ice with a fixed water fraction corresponding to the
55 * amount of water remaining in the system once all of the moving water drains out (0.5%
56 * by default). During the winter the CH system is allowed to cool.
57 */
58
59 CHSystem::CHSystem(IceGrid::ConstPtr grid,
60 stressbalance::StressBalance *stress_balance)
61 : EnergyModel(grid, stress_balance) {
62
63 m_ice_enthalpy.set_name("ch_enthalpy");
64 m_ice_enthalpy.metadata().set_name("ch_enthalpy");
65 m_ice_enthalpy.metadata().set_string("long_name",
66 "enthalpy of the cryo-hydrologic system");
67 }
68
69 CHSystem::~CHSystem() {
70 // empty
71 }
72
73 void CHSystem::restart_impl(const File &input_file, int record) {
74
75 m_log->message(2, "* Restarting the cryo-hydrologic system from %s...\n",
76 input_file.filename().c_str());
77
78 init_enthalpy(input_file, false, record);
79
80 regrid_enthalpy();
81 }
82
83 void CHSystem::bootstrap_impl(const File &input_file,
84 const IceModelVec2S &ice_thickness,
85 const IceModelVec2S &surface_temperature,
86 const IceModelVec2S &climatic_mass_balance,
87 const IceModelVec2S &basal_heat_flux) {
88
89 m_log->message(2, "* Bootstrapping the cryo-hydrologic warming model from %s...\n",
90 input_file.filename().c_str());
91
92 int enthalpy_revision = m_ice_enthalpy.state_counter();
93 regrid_enthalpy();
94
95 if (enthalpy_revision == m_ice_enthalpy.state_counter()) {
96 bootstrap_ice_enthalpy(ice_thickness, surface_temperature, climatic_mass_balance,
97 basal_heat_flux, m_ice_enthalpy);
98 }
99 }
100
101 void CHSystem::initialize_impl(const IceModelVec2S &basal_melt_rate,
102 const IceModelVec2S &ice_thickness,
103 const IceModelVec2S &surface_temperature,
104 const IceModelVec2S &climatic_mass_balance,
105 const IceModelVec2S &basal_heat_flux) {
106 (void) basal_melt_rate;
107
108 m_log->message(2, "* Bootstrapping the cryo-hydrologic warming model...\n");
109
110 int enthalpy_revision = m_ice_enthalpy.state_counter();
111 regrid_enthalpy();
112
113 if (enthalpy_revision == m_ice_enthalpy.state_counter()) {
114 bootstrap_ice_enthalpy(ice_thickness, surface_temperature, climatic_mass_balance,
115 basal_heat_flux, m_ice_enthalpy);
116 }
117 }
118
119 //! Update the enthalpy of the cryo-hydrologic system.
120 /*!
121 This method updates IceModelVec3 m_work. No communication of ghosts is done.
122 */
123 void CHSystem::update_impl(double t, double dt, const Inputs &inputs) {
124 // current time does not matter here
125 (void) t;
126
127 EnthalpyConverter::Ptr EC = m_grid->ctx()->enthalpy_converter();
128
129 inputs.check();
130
131 // give them names that are a bit shorter...
132 const IceModelVec3
133 &volumetric_heat = *inputs.volumetric_heating_rate,
134 &u3 = *inputs.u3,
135 &v3 = *inputs.v3,
136 &w3 = *inputs.w3;
137
138 const IceModelVec2CellType &cell_type = *inputs.cell_type;
139
140 const IceModelVec2S
141 &basal_frictional_heating = *inputs.basal_frictional_heating,
142 &basal_heat_flux = *inputs.basal_heat_flux,
143 &ice_thickness = *inputs.ice_thickness,
144 &surface_liquid_fraction = *inputs.surface_liquid_fraction,
145 &shelf_base_temp = *inputs.shelf_base_temp,
146 &ice_surface_temp = *inputs.surface_temp;
147
148 energy::enthSystemCtx system(m_grid->z(), "energy.ch_warming", m_grid->dx(), m_grid->dy(), dt,
149 *m_config, m_ice_enthalpy, u3, v3, w3, volumetric_heat, EC);
150
151 const size_t Mz_fine = system.z().size();
152 const double dz = system.dz();
153 std::vector<double> Enthnew(Mz_fine); // new enthalpy in column
154
155 IceModelVec::AccessList list{&ice_surface_temp, &shelf_base_temp, &surface_liquid_fraction,
156 &ice_thickness, &basal_frictional_heating, &basal_heat_flux,
157 &cell_type, &u3, &v3, &w3, &volumetric_heat, &m_ice_enthalpy,
158 &m_work};
159
160 double
161 margin_threshold = m_config->get_number("energy.margin_ice_thickness_limit"),
162 T_pm = m_config->get_number("constants.fresh_water.melting_point_temperature"),
163 residual_water_fraction = m_config->get_number("energy.ch_warming.residual_water_fraction");
164
165 const std::vector<double> &z = m_grid->z();
166 const unsigned int Mz = m_grid->Mz();
167
168 ParallelSection loop(m_grid->com);
169 try {
170 for (Points pt(*m_grid); pt; pt.next()) {
171 const int i = pt.i(), j = pt.j();
172
173 const double H = ice_thickness(i, j);
174
175 if (ice_surface_temp(i, j) >= T_pm) {
176 // We use surface temperature to determine if we're in a melt season or not. It
177 // probably makes sense to use the surface mass balance instead.
178
179 double *column = m_work.get_column(i, j);
180 for (unsigned int k = 0; k < Mz; ++k) {
181 double
182 depth = std::max(H - z[k], 0.0),
183 P = EC->pressure(depth);
184 column[k] = EC->enthalpy(EC->melting_temperature(P),
185 residual_water_fraction,
186 P);
187 }
188 continue;
189 }
190
191 // enthalpy and pressures at top of ice
192 const double
193 depth_ks = H - system.ks() * dz,
194 p_ks = EC->pressure(depth_ks); // FIXME issue #15
195
196 const double Enth_ks = EC->enthalpy_permissive(ice_surface_temp(i, j),
197 surface_liquid_fraction(i, j), p_ks);
198
199 system.init(i, j,
200 marginal(ice_thickness, i, j, margin_threshold),
201 H);
202
203 const bool ice_free_column = (system.ks() == 0);
204
205 // deal completely with columns with no ice
206 if (ice_free_column) {
207 m_work.set_column(i, j, Enth_ks);
208 continue;
209 } // end of if (ice_free_column)
210
211 if (system.lambda() < 1.0) {
212 m_stats.reduced_accuracy_counter += 1; // count columns with lambda < 1
213 }
214
215 // set boundary conditions and update enthalpy
216 {
217 system.set_surface_dirichlet_bc(Enth_ks);
218
219 if (cell_type.ocean(i, j)) {
220 // floating base: Dirichlet application of known temperature from ocean coupler;
221 // assumes base of ice shelf has zero liquid fraction
222 double Enth0 = EC->enthalpy_permissive(shelf_base_temp(i, j), 0.0, EC->pressure(H));
223
224 system.set_basal_dirichlet_bc(Enth0);
225 } else {
226 // grounded
227 system.set_basal_heat_flux(basal_heat_flux(i, j) + basal_frictional_heating(i, j));
228 }
229 // solve the system
230 system.solve(Enthnew);
231 }
232
233 system.fine_to_coarse(Enthnew, i, j, m_work);
234 }
235 } catch (...) {
236 loop.failed();
237 }
238 loop.check();
239 }
240
241 void CHSystem::define_model_state_impl(const File &output) const {
242 m_ice_enthalpy.define(output);
243 }
244
245 void CHSystem::write_model_state_impl(const File &output) const {
246 m_ice_enthalpy.write(output);
247 }
248
249 /*!
250 * Compute the heat flux corresponding to the cryo-hydrologic warming.
251 *
252 * `Q = (k / R**2) * (T_ch - T_ice),`
253 *
254 * where `k` is the thermal conductivity of ice and `R` us the average spacing of
255 * channels in the cryo-hydrologic system.
256 */
257 void cryo_hydrologic_warming_flux(double k,
258 double R,
259 const IceModelVec2S &ice_thickness,
260 const IceModelVec3 &ice_enthalpy,
261 const IceModelVec3 &ch_enthalpy,
262 IceModelVec3 &result) {
263
264 auto grid = result.grid();
265
266 const auto &z = grid->z();
267 auto Mz = grid->Mz();
268
269 auto EC = grid->ctx()->enthalpy_converter();
270
271 IceModelVec::AccessList access{&ice_thickness, &ice_enthalpy, &ch_enthalpy, &result};
272
273 double C = k / (R * R);
274
275 ParallelSection loop(grid->com);
276 try {
277 for (Points p(*grid); p; p.next()) {
278 const int i = p.i(), j = p.j();
279
280 const double
281 *E_ch = ch_enthalpy.get_column(i, j),
282 *E_ice = ice_enthalpy.get_column(i, j);
283 double *Q = result.get_column(i, j);
284
285 for (unsigned int m = 0; m < Mz; ++m) {
286 double
287 depth = ice_thickness(i, j) - z[m];
288
289 if (depth > 0.0) {
290 double P = EC->pressure(depth);
291 Q[m] = std::max(C * (EC->temperature(E_ch[m], P) - EC->temperature(E_ice[m], P)),
292 0.0);
293 } else {
294 Q[m] = 0.0;
295 }
296 }
297 }
298 } catch (...) {
299 loop.failed();
300 }
301 loop.check();
302 }
303
304
305 } // end of namespace energy
306 } // end of namespace pism