tForceThickness.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
---
tForceThickness.cc (12786B)
---
1 // Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 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 "ForceThickness.hh"
20 #include "pism/util/IceGrid.hh"
21 #include "pism/util/Vars.hh"
22
23 #include "pism/util/ConfigInterface.hh"
24 #include "pism/util/Mask.hh"
25 #include "pism/util/pism_options.hh"
26 #include "pism/util/error_handling.hh"
27 #include "pism/util/io/io_helpers.hh"
28 #include "pism/util/pism_utilities.hh"
29 #include "pism/util/IceModelVec2CellType.hh"
30 #include "pism/util/MaxTimestep.hh"
31 #include "pism/util/io/File.hh"
32 #include "pism/geometry/Geometry.hh"
33
34 namespace pism {
35 namespace surface {
36
37 ///// "Force-to-thickness" mechanism
38 ForceThickness::ForceThickness(IceGrid::ConstPtr g, std::shared_ptr<SurfaceModel> input)
39 : SurfaceModel(g, input) {
40
41 m_alpha = m_config->get_number("surface.force_to_thickness.alpha", "s-1");
42 m_alpha_ice_free_factor = m_config->get_number("surface.force_to_thickness.ice_free_alpha_factor");
43 m_ice_free_thickness_threshold = m_config->get_number("surface.force_to_thickness.ice_free_thickness_threshold");
44 m_start_time = m_config->get_number("surface.force_to_thickness.start_time", "seconds");
45
46 m_target_thickness.create(m_grid, "thk", WITHOUT_GHOSTS);
47 // will set attributes in init()
48
49 m_ftt_mask.create(m_grid, "ftt_mask", WITHOUT_GHOSTS);
50 m_ftt_mask.set_attrs("diagnostic",
51 "mask specifying where to apply the force-to-thickness mechanism",
52 "", "", "", 0); // no units and no standard name
53 m_ftt_mask.set(1.0); // default: applied in whole domain
54 m_ftt_mask.metadata().set_output_type(PISM_INT);
55 m_ftt_mask.metadata().set_time_independent(true);
56
57 m_mass_flux = allocate_mass_flux(g);
58
59 m_accumulation = allocate_accumulation(g);
60 m_melt = allocate_melt(g);
61 m_runoff = allocate_runoff(g);
62 }
63
64 ForceThickness::~ForceThickness() {
65 // empty
66 }
67
68 void ForceThickness::init_impl(const Geometry &geometry) {
69
70 m_input_model->init(geometry);
71
72 m_log->message(2,
73 "* Initializing force-to-thickness mass-balance modifier...\n");
74
75 std::string input_file = m_config->get_string("surface.force_to_thickness_file");
76
77 if (input_file.empty()) {
78 throw RuntimeError(PISM_ERROR_LOCATION, "surface.force_to_thickness_file cannot be empty");
79 }
80
81 m_log->message(2,
82 " alpha = %.6f year-1 for -force_to_thickness mechanism\n"
83 " alpha = %.6f year-1 in areas with target ice thickness of less than %.3f meters\n",
84 units::convert(m_sys, m_alpha, "s-1", "year-1"),
85 m_alpha_ice_free_factor * units::convert(m_sys, m_alpha, "s-1", "year-1"),
86 m_ice_free_thickness_threshold);
87
88 // check of the input file is really there and regrid the target thickness
89 File file(m_grid->com, input_file, PISM_GUESS, PISM_READONLY);
90
91 m_log->message(2,
92 " reading target thickness 'thk' from %s ...\n"
93 " (this field will appear in output file as 'ftt_target_thk')\n",
94 input_file.c_str());
95 {
96 m_target_thickness.metadata(0).set_name("thk"); // name to read by
97 // set attributes for the read stage; see below for reset
98 m_target_thickness.set_attrs("diagnostic",
99 "target thickness for force-to-thickness mechanism (hit this at end of run)",
100 "m", "m",
101 "land_ice_thickness", 0); // standard_name *to read by*
102
103 m_target_thickness.regrid(input_file, CRITICAL);
104
105 // reset name to avoid confusion; set attributes again to overwrite "read by" choices above
106 m_target_thickness.metadata(0).set_name("ftt_target_thk");
107 m_target_thickness.set_attrs("diagnostic",
108 "target thickness for force-to-thickness mechanism (wants to hit this at end of run)",
109 "m", "m",
110 "", 0); // no CF standard_name, to put it mildly
111 }
112
113 {
114 m_log->message(2,
115 " reading force-to-thickness mask 'ftt_mask' from %s ...\n",
116 input_file.c_str());
117 m_ftt_mask.regrid(input_file, CRITICAL);
118 }
119 }
120
121 /*!
122 If `-force_to_thickness_file` `foo.nc` is in use then vthktarget will have a target ice thickness
123 map. Let \f$H_{\text{tar}}\f$ be this target thickness,
124 and let \f$H\f$ be the current model thickness. Recall that the mass continuity
125 equation solved by GeometryEvolution is
126 \f[ \frac{\partial H}{\partial t} = M - S - \nabla\cdot \mathbf{q} \f]
127 and that this procedure is supposed to produce \f$M\f$.
128 In this context, the semantics of `-force_to_thickness` are that \f$M\f$ is modified
129 by a multiple of the difference between the target thickness and the current thickness.
130 In particular, the \f$\Delta M\f$ that is produced here is
131 \f[\Delta M = \alpha (H_{\text{tar}} - H)\f]
132 where \f$\alpha>0\f$ is determined below. Note \f$\Delta M\f$ is positive in
133 areas where \f$H_{\text{tar}} > H\f$, so we are adding mass there, and we are ablating
134 in the other case.
135
136 Let \f$t_s\f$ be the start time and \f$t_e\f$ the end time for the run.
137 Without flow or basal mass balance, or any surface mass balance other than the
138 \f$\Delta M\f$ computed here, we are solving
139 \f[ \frac{\partial H}{\partial t} = \alpha (H_{\text{tar}} - H) \f]
140 Let's assume \f$H(t_s)=H_0\f$. This initial value problem has solution
141 \f$H(t) = H_{\text{tar}} + (H_0 - H_{\text{tar}}) e^{-\alpha (t-t_s)}\f$
142 and so
143 \f[ H(t_e) = H_{\text{tar}} + (H_0 - H_{\text{tar}}) e^{-\alpha (t_e-t_s)} \f]
144
145 The constant \f$\alpha\f$ has a default value `pism_config:surface.force_to_thickness.alpha`.
146
147 The next example uses files generated from the EISMINT-Greenland experiment;
148 see the corresponding chapter of the User's Manual.
149
150 Suppose we regard the SSL2 run as a spin-up to reach a better temperature field.
151 It is a spinup in which the surface was allowed to evolve. Assume the
152 early file `green20km_y1.nc` has the target thickness, because it essentially
153 has the input thickness. This script adds a 500 a run, to finalize the spinup,
154 in which the ice sheet geometry goes from the the thickness values in
155 `green_ssl2_110ka.nc` to values very close to those in `green20km_y1.nc`:
156 \code
157 #!/bin/bash
158
159 NN=8 # default number of processors
160 if [ $# -gt 0 ] ; then # if user says "test_ftt.sh 8" then NN = 8
161 NN="$1"
162 fi
163
164 # set MPIDO if using different MPI execution command, for example:
165 # $ export PISM_MPIDO="aprun -n "
166 if [ -n "${PISM_MPIDO:+1}" ] ; then # check if env var is already set
167 echo "$SCRIPTNAME PISM_MPIDO = $PISM_MPIDO (already set)"
168 else
169 PISM_MPIDO="mpiexec -n "
170 echo "$SCRIPTNAME PISM_MPIDO = $PISM_MPIDO"
171 fi
172
173 # check if env var PISM_DO was set (i.e. PISM_DO=echo for a 'dry' run)
174 if [ -n "${PISM_DO:+1}" ] ; then # check if env var DO is already set
175 echo "$SCRIPTNAME PISM_DO = $PISM_DO (already set)"
176 else
177 PISM_DO=""
178 fi
179
180 # prefix to pism (not to executables)
181 if [ -n "${PISM_PREFIX:+1}" ] ; then # check if env var is already set
182 echo "$SCRIPTNAME PISM_PREFIX = $PISM_PREFIX (already set)"
183 else
184 PISM_PREFIX="" # just a guess
185 echo "$SCRIPTNAME PISM_PREFIX = $PISM_PREFIX"
186 fi
187
188 # set PISM_EXEC if using different executables, for example:
189 # $ export PISM_EXEC="pismr -energy cold"
190 if [ -n "${PISM_EXEC:+1}" ] ; then # check if env var is already set
191 echo "$SCRIPTNAME PISM_EXEC = $PISM_EXEC (already set)"
192 else
193 PISM_EXEC="pismr"
194 echo "$SCRIPTNAME PISM_EXEC = $PISM_EXEC"
195 fi
196
197
198 PISM="${PISM_PREFIX}${PISM_EXEC}"
199
200 cmd="$PISM_MPIDO $NN $PISM -ys -1000.0 -ye 0 -skip 5 -i green_ssl2_110ka.nc -atmosphere searise_greenland \
201 -surface pdd -pdd_fausto \
202 -o no_force.nc -ts_file ts_no_force.nc -ts_times -1000:yearly:0"
203 $PISM_DO $cmd
204
205 echo
206
207 cmd="$PISM_MPIDO $NN $PISM -ys -1000.0 -ye 0 -skip 5 -i green_ssl2_110ka.nc -atmosphere searise_greenland \
208 -surface pdd,forcing -pdd_fausto -force_to_thickness_file green20km_y1.nc \
209 -o default_force.nc -ts_file ts_default_force.nc -ts_times -1000:yearly:0"
210 $PISM_DO $cmd
211
212 echo
213
214 cmd="$PISM_MPIDO $NN $PISM -ys -1000.0 -ye 0 -skip 5 -i green_ssl2_110ka.nc -atmosphere searise_greenland \
215 -surface pdd,forcing -pdd_fausto -force_to_thickness_file green20km_y1.nc -force_to_thickness_alpha 0.005 \
216 -o weak_force.nc -ts_file ts_weak_force.nc -ts_times -1000:yearly:0"
217 $PISM_DO $cmd
218
219
220 cmd="$PISM_MPIDO $NN $PISM -ys -1000.0 -ye 0 -skip 5 -i green_ssl2_110ka.nc -atmosphere searise_greenland \
221 -surface pdd,forcing -pdd_fausto -force_to_thickness_file green20km_y1.nc -force_to_thickness_alpha 0.05 \
222 -o strong_force.nc -ts_file ts_strong_force.nc -ts_times -1000:yearly:0"
223 $PISM_DO $cmd
224
225 \endcode
226 The script also has a run with no forcing, one with forcing at a lower alpha value,
227 a factor of five smaller than the default, and one with a forcing at a higher alpha value, a factor of five higher.
228 */
229 void ForceThickness::adjust_mass_flux(double time,
230 const IceModelVec2S &ice_thickness,
231 const IceModelVec2CellType &cell_type,
232 IceModelVec2S &result) const {
233
234 if (time < m_start_time) {
235 return;
236 }
237
238 m_log->message(5,
239 " updating surface mass balance using -force_to_thickness mechanism ...");
240
241 double ice_density = m_config->get_number("constants.ice.density");
242
243 IceModelVec::AccessList list{&cell_type, &ice_thickness,
244 &m_target_thickness, &m_ftt_mask, &result};
245
246 for (Points p(*m_grid); p; p.next()) {
247 const int i = p.i(), j = p.j();
248
249 if (m_ftt_mask(i,j) > 0.5 and cell_type.grounded(i, j)) {
250 if (m_target_thickness(i,j) >= m_ice_free_thickness_threshold) {
251 result(i,j) += ice_density * m_alpha * (m_target_thickness(i,j) - ice_thickness(i,j));
252 } else {
253 result(i,j) += ice_density * m_alpha * m_alpha_ice_free_factor * (m_target_thickness(i,j) - ice_thickness(i,j));
254 }
255 }
256 }
257 // no communication needed
258 }
259
260 void ForceThickness::update_impl(const Geometry &geometry, double t, double dt) {
261 m_input_model->update(geometry, t, dt);
262
263 m_mass_flux->copy_from(m_input_model->mass_flux());
264
265 adjust_mass_flux(t,
266 geometry.ice_thickness,
267 geometry.cell_type,
268 *m_mass_flux);
269
270 dummy_accumulation(*m_mass_flux, *m_accumulation);
271 dummy_melt(*m_mass_flux, *m_melt);
272 dummy_runoff(*m_mass_flux, *m_runoff);
273 }
274
275 const IceModelVec2S &ForceThickness::mass_flux_impl() const {
276 return *m_mass_flux;
277 }
278
279 const IceModelVec2S &ForceThickness::accumulation_impl() const {
280 return *m_accumulation;
281 }
282
283 const IceModelVec2S &ForceThickness::melt_impl() const {
284 return *m_melt;
285 }
286
287 const IceModelVec2S &ForceThickness::runoff_impl() const {
288 return *m_runoff;
289 }
290
291 /*!
292 The timestep restriction is, by direct analogy, the same as for
293 \f[\frac{dy}{dt} = - \alpha y\f]
294 with explicit (forward) Euler. If \f$\Delta t\f$ is the time step then Euler is
295 \f$y_{n+1} = (1-\alpha \Delta t) y_n\f$. We require for stability that
296 \f$|y_{n+1}|\le |y_n|\f$, which is to say \f$|1-\alpha \Delta t|\le 1\f$.
297 Equivalently (since \f$\alpha \Delta t>0\f$),
298 \f[\alpha \Delta t\le 2\f]
299 Therefore we set here
300 \f[\Delta t = \frac{2}{\alpha}.\f]
301 */
302 MaxTimestep ForceThickness::max_timestep_impl(double my_t) const {
303 double max_dt = 2.0 / m_alpha;
304 MaxTimestep input_max_dt = m_input_model->max_timestep(my_t);
305
306 return std::min(input_max_dt, MaxTimestep(max_dt, "surface forcing"));
307 }
308
309 void ForceThickness::define_model_state_impl(const File &output) const {
310 m_ftt_mask.define(output);
311 m_target_thickness.define(output);
312
313 if (m_input_model != NULL) {
314 m_input_model->define_model_state(output);
315 }
316 }
317
318 void ForceThickness::write_model_state_impl(const File &output) const {
319 m_ftt_mask.write(output);
320 m_target_thickness.write(output);
321
322 if (m_input_model != NULL) {
323 m_input_model->write_model_state(output);
324 }
325 }
326
327 } // end of namespace surface
328 } // end of namespace pism