tfrontretreat.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
---
tfrontretreat.cc (8592B)
---
1 // Copyright (C) 2004--2020 Torsten Albrecht and Constantine Khroulev
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 "IceModel.hh"
20
21 #include "pism/util/IceGrid.hh"
22 #include "pism/util/Mask.hh"
23 #include "pism/util/ConfigInterface.hh"
24 #include "pism/util/pism_utilities.hh"
25
26 #include "pism/frontretreat/FrontRetreat.hh"
27 #include "pism/frontretreat/util/IcebergRemover.hh"
28 #include "pism/frontretreat/calving/CalvingAtThickness.hh"
29 #include "pism/frontretreat/calving/EigenCalving.hh"
30 #include "pism/frontretreat/calving/FloatKill.hh"
31 #include "pism/frontretreat/calving/HayhurstCalving.hh"
32 #include "pism/frontretreat/calving/vonMisesCalving.hh"
33
34 #include "pism/energy/EnergyModel.hh"
35 #include "pism/coupler/FrontalMelt.hh"
36 #include "pism/stressbalance/ShallowStressBalance.hh"
37 #include "pism/hydrology/Hydrology.hh"
38 #include "pism/frontretreat/util/remove_narrow_tongues.hh"
39 #include "pism/frontretreat/PrescribedRetreat.hh"
40
41 namespace pism {
42
43 void IceModel::front_retreat_step() {
44 const bool
45 add_values = true,
46 insert_values = false;
47
48 // compute retreat rates due to eigencalving, von Mises calving, Hayhurst calving,
49 // and frontal melt.
50 // We do this first to make sure that all three mechanisms use the same ice geometry.
51 {
52 if (m_eigen_calving) {
53 m_eigen_calving->update(m_geometry.cell_type,
54 m_stress_balance->shallow()->velocity());
55 }
56
57 if (m_hayhurst_calving) {
58 m_hayhurst_calving->update(m_geometry.cell_type,
59 m_geometry.ice_thickness,
60 m_geometry.sea_level_elevation,
61 m_geometry.bed_elevation);
62 }
63
64 if (m_vonmises_calving) {
65 // FIXME: consider computing vertically-averaged hardness here and providing that
66 // instead of using ice thickness and enthalpy.
67 m_vonmises_calving->update(m_geometry.cell_type,
68 m_geometry.ice_thickness,
69 m_stress_balance->shallow()->velocity(),
70 m_energy_model->enthalpy());
71 }
72
73 if (m_frontal_melt) {
74 IceModelVec2S &flux_magnitude = m_work2d[0];
75
76 flux_magnitude.set_to_magnitude(m_subglacial_hydrology->flux());
77
78 FrontalMeltInputs inputs;
79
80 inputs.geometry = &m_geometry;
81 inputs.subglacial_water_flux = &flux_magnitude;
82
83 m_frontal_melt->update(inputs, m_time->current(), m_dt);
84 }
85 }
86
87 IceModelVec2S
88 &old_H = m_work2d[0],
89 &old_Href = m_work2d[1];
90
91 // frontal melt
92 if (m_frontal_melt) {
93 assert(m_front_retreat);
94
95 old_H.copy_from(m_geometry.ice_thickness);
96 old_Href.copy_from(m_geometry.ice_area_specific_volume);
97
98 // apply frontal melt rate
99 m_front_retreat->update_geometry(m_dt, m_geometry, m_ssa_dirichlet_bc_mask,
100 m_frontal_melt->retreat_rate(),
101 m_geometry.ice_area_specific_volume,
102 m_geometry.ice_thickness);
103
104 compute_geometry_change(m_geometry.ice_thickness,
105 m_geometry.ice_area_specific_volume,
106 old_H, old_Href,
107 insert_values,
108 m_thickness_change.frontal_melt);
109 } else {
110 m_thickness_change.frontal_melt.set(0.0);
111 }
112
113 // calving
114 if (m_eigen_calving or m_vonmises_calving or m_hayhurst_calving or
115 m_float_kill_calving or m_thickness_threshold_calving) {
116
117 old_H.copy_from(m_geometry.ice_thickness);
118 old_Href.copy_from(m_geometry.ice_area_specific_volume);
119
120 if (m_eigen_calving or m_vonmises_calving or m_hayhurst_calving) {
121 assert(m_front_retreat);
122
123 IceModelVec2S &retreat_rate = m_work2d[2];
124 retreat_rate.set(0.0);
125
126 if (m_eigen_calving) {
127 retreat_rate.add(1.0, m_eigen_calving->calving_rate());
128 }
129
130 if (m_hayhurst_calving) {
131 retreat_rate.add(1.0, m_hayhurst_calving->calving_rate());
132 }
133
134 if (m_vonmises_calving) {
135 retreat_rate.add(1.0, m_vonmises_calving->calving_rate());
136 }
137
138 m_front_retreat->update_geometry(m_dt, m_geometry, m_ssa_dirichlet_bc_mask,
139 retreat_rate,
140 m_geometry.ice_area_specific_volume,
141 m_geometry.ice_thickness);
142
143 auto thickness_threshold = m_config->get_number("stress_balance.ice_free_thickness_standard");
144
145 m_geometry.ensure_consistency(thickness_threshold);
146
147 if (m_eigen_calving or m_vonmises_calving or m_hayhurst_calving) {
148 remove_narrow_tongues(m_geometry, m_geometry.ice_thickness);
149
150 m_geometry.ensure_consistency(thickness_threshold);
151 }
152 }
153
154 if (m_float_kill_calving) {
155 m_float_kill_calving->update(m_geometry.cell_type, m_geometry.ice_thickness);
156 }
157
158 if (m_thickness_threshold_calving) {
159 m_thickness_threshold_calving->update(m_geometry.cell_type, m_geometry.ice_thickness);
160 }
161
162 compute_geometry_change(m_geometry.ice_thickness,
163 m_geometry.ice_area_specific_volume,
164 old_H, old_Href,
165 insert_values,
166 m_thickness_change.calving);
167 } else {
168 m_thickness_change.calving.set(0.0);
169 }
170
171 // prescribed retreat
172
173 if (m_prescribed_retreat) {
174 old_H.copy_from(m_geometry.ice_thickness);
175 old_Href.copy_from(m_geometry.ice_area_specific_volume);
176
177 m_prescribed_retreat->update(m_time->current(), m_dt,
178 m_geometry.ice_thickness,
179 m_geometry.ice_area_specific_volume);
180
181 compute_geometry_change(m_geometry.ice_thickness,
182 m_geometry.ice_area_specific_volume,
183 old_H, old_Href,
184 insert_values,
185 m_thickness_change.forced_retreat);
186
187 } else {
188 m_thickness_change.forced_retreat.set(0.0);
189 }
190
191
192 // Changes above may create icebergs; here we remove them and account for additional
193 // mass losses.
194 {
195 old_H.copy_from(m_geometry.ice_thickness);
196 old_Href.copy_from(m_geometry.ice_area_specific_volume);
197
198 enforce_consistency_of_geometry(REMOVE_ICEBERGS);
199
200 compute_geometry_change(m_geometry.ice_thickness,
201 m_geometry.ice_area_specific_volume,
202 old_H, old_Href,
203 add_values,
204 m_thickness_change.calving);
205 }
206 }
207
208 /**
209 * Compute the change in ice geometry from "old" to "current".
210 *
211 * Units: ice equivalent meters.
212 *
213 * @param thickness current ice thickness
214 * @param Href current "reference ice thickness"
215 * @param thickness_old old ice thickness
216 * @param Href_old old "reference ice thickness"
217 * @param[in] flag if `flag == ADD_VALUES`, add computed values to `output`, otherwise
218 * overwrite them
219 * @param[in,out] output computed change
220 */
221 void IceModel::compute_geometry_change(const IceModelVec2S &thickness,
222 const IceModelVec2S &Href,
223 const IceModelVec2S &thickness_old,
224 const IceModelVec2S &Href_old,
225 bool add_values,
226 IceModelVec2S &output) {
227
228 IceModelVec::AccessList list{&thickness, &thickness_old,
229 &Href, &Href_old, &output};
230
231 for (Points p(*m_grid); p; p.next()) {
232 const int i = p.i(), j = p.j();
233
234 const double
235 H_old = thickness_old(i, j) + Href_old(i, j),
236 H_new = thickness(i, j) + Href(i, j),
237 change = H_new - H_old;
238
239 if (add_values) {
240 output(i, j) += change;
241 } else {
242 output(i, j) = change;
243 }
244 }
245 }
246
247 } // end of namespace pism