Loading [MathJax]/extensions/tex2jax.js
PISM, A Parallel Ice Sheet Model 2.2.1-cd005eec8 committed by Constantine Khrulev on 2025-03-07
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ForceThickness.cc
Go to the documentation of this file.
1// Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2022, 2023, 2024 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 "pism/coupler/surface/ForceThickness.hh"
20#include "pism/util/Grid.hh"
21
22#include "pism/util/ConfigInterface.hh"
23#include "pism/util/error_handling.hh"
24#include "pism/util/array/CellType.hh"
25#include "pism/util/MaxTimestep.hh"
26#include "pism/util/io/File.hh"
27#include "pism/geometry/Geometry.hh"
28#include "pism/util/Interpolation1D.hh"
29
30namespace pism {
31namespace surface {
32
33///// "Force-to-thickness" mechanism
34ForceThickness::ForceThickness(std::shared_ptr<const Grid> g, std::shared_ptr<SurfaceModel> input)
35 : SurfaceModel(g, input),
36 m_target_thickness(m_grid, "thk"),
37 m_ftt_mask(m_grid, "ftt_mask")
38{
39
41
42 m_alpha = m_config->get_number("surface.force_to_thickness.alpha", "s-1");
43 m_alpha_ice_free_factor = m_config->get_number("surface.force_to_thickness.ice_free_alpha_factor");
44 m_ice_free_thickness_threshold = m_config->get_number("surface.force_to_thickness.ice_free_thickness_threshold");
45 m_start_time = m_config->get_number("surface.force_to_thickness.start_time", "seconds");
46
48 .long_name("mask specifying where to apply the force-to-thickness mechanism")
51
52 m_ftt_mask.set(1.0); // default: applied in whole domain
53
55
59}
60
61void ForceThickness::init_impl(const Geometry &geometry) {
62
63 m_input_model->init(geometry);
64
65 m_log->message(2, "* Initializing force-to-thickness mass-balance modifier...\n");
66
67 std::string input_file = m_config->get_string("surface.force_to_thickness.file");
68
69 if (input_file.empty()) {
70 throw RuntimeError(PISM_ERROR_LOCATION, "surface.force_to_thickness.file cannot be empty");
71 }
72
73 m_log->message(
74 2,
75 " alpha = %.6f year-1 for -force_to_thickness mechanism\n"
76 " alpha = %.6f year-1 in areas with target ice thickness of less than %.3f meters\n",
77 units::convert(m_sys, m_alpha, "s^-1", "year^-1"),
80
81 // check of the input file is really there and regrid the target thickness
82 File file(m_grid->com, input_file, io::PISM_GUESS, io::PISM_READONLY);
83
84 m_log->message(2,
85 " reading target thickness 'thk' from %s ...\n"
86 " (this field will appear in output file as 'ftt_target_thk')\n",
87 input_file.c_str());
88 {
89 m_target_thickness.metadata(0).set_name("thk"); // name to read by
90 // set attributes for the read stage; see below for reset
92 .long_name("target thickness for force-to-thickness mechanism (hit this at end of run)")
93 .units("m")
94 .standard_name("land_ice_thickness"); // standard_name *to read by*
95
97
98 // reset name to avoid confusion; set attributes again to overwrite "read by" choices above
99 m_target_thickness.metadata(0).set_name("ftt_target_thk");
100 m_target_thickness.metadata(0).standard_name(""); // no CF standard_name, to put it mildly
101 }
102
103 {
104 m_log->message(2,
105 " reading force-to-thickness mask 'ftt_mask' from %s ...\n",
106 input_file.c_str());
107 m_ftt_mask.regrid(input_file, io::Default::Nil());
108 }
109}
110
111/*!
112If `-force_to_thickness_file` `foo.nc` is in use then vthktarget will have a target ice thickness
113map. Let \f$H_{\text{tar}}\f$ be this target thickness,
114and let \f$H\f$ be the current model thickness. Recall that the mass continuity
115equation solved by GeometryEvolution is
116 \f[ \frac{\partial H}{\partial t} = M - S - \nabla\cdot \mathbf{q} \f]
117and that this procedure is supposed to produce \f$M\f$.
118In this context, the semantics of `-force_to_thickness` are that \f$M\f$ is modified
119by a multiple of the difference between the target thickness and the current thickness.
120In particular, the \f$\Delta M\f$ that is produced here is
121 \f[\Delta M = \alpha (H_{\text{tar}} - H)\f]
122where \f$\alpha>0\f$ is determined below. Note \f$\Delta M\f$ is positive in
123areas where \f$H_{\text{tar}} > H\f$, so we are adding mass there, and we are ablating
124in the other case.
125
126Let \f$t_s\f$ be the start time and \f$t_e\f$ the end time for the run.
127Without flow or basal mass balance, or any surface mass balance other than the
128\f$\Delta M\f$ computed here, we are solving
129 \f[ \frac{\partial H}{\partial t} = \alpha (H_{\text{tar}} - H) \f]
130Let's assume \f$H(t_s)=H_0\f$. This initial value problem has solution
131\f$H(t) = H_{\text{tar}} + (H_0 - H_{\text{tar}}) e^{-\alpha (t-t_s)}\f$
132and so
133 \f[ H(t_e) = H_{\text{tar}} + (H_0 - H_{\text{tar}}) e^{-\alpha (t_e-t_s)} \f]
134
135The constant \f$\alpha\f$ has a default value `pism_config:surface.force_to_thickness.alpha`.
136
137The next example uses files generated from the EISMINT-Greenland experiment;
138see the corresponding chapter of the User's Manual.
139
140Suppose we regard the SSL2 run as a spin-up to reach a better temperature field.
141It is a spinup in which the surface was allowed to evolve. Assume the
142early file `green20km_y1.nc` has the target thickness, because it essentially
143has the input thickness. This script adds a 500 a run, to finalize the spinup,
144in which the ice sheet geometry goes from the the thickness values in
145`green_ssl2_110ka.nc` to values very close to those in `green20km_y1.nc`:
146\code
147#!/bin/bash
148
149NN=8 # default number of processors
150if [ $# -gt 0 ] ; then # if user says "test_ftt.sh 8" then NN = 8
151 NN="$1"
152fi
153
154# set MPIDO if using different MPI execution command, for example:
155# $ export PISM_MPIDO="aprun -n "
156if [ -n "${PISM_MPIDO:+1}" ] ; then # check if env var is already set
157 echo "$SCRIPTNAME PISM_MPIDO = $PISM_MPIDO (already set)"
158else
159 PISM_MPIDO="mpiexec -n "
160 echo "$SCRIPTNAME PISM_MPIDO = $PISM_MPIDO"
161fi
162
163# check if env var PISM_DO was set (i.e. PISM_DO=echo for a 'dry' run)
164if [ -n "${PISM_DO:+1}" ] ; then # check if env var DO is already set
165 echo "$SCRIPTNAME PISM_DO = $PISM_DO (already set)"
166else
167 PISM_DO=""
168fi
169
170# prefix to pism (not to executables)
171if [ -n "${PISM_PREFIX:+1}" ] ; then # check if env var is already set
172 echo "$SCRIPTNAME PISM_PREFIX = $PISM_PREFIX (already set)"
173else
174 PISM_PREFIX="" # just a guess
175 echo "$SCRIPTNAME PISM_PREFIX = $PISM_PREFIX"
176fi
177
178# set PISM_EXEC if using different executables, for example:
179# $ export PISM_EXEC="pism -energy cold"
180if [ -n "${PISM_EXEC:+1}" ] ; then # check if env var is already set
181 echo "$SCRIPTNAME PISM_EXEC = $PISM_EXEC (already set)"
182else
183 PISM_EXEC="pism"
184 echo "$SCRIPTNAME PISM_EXEC = $PISM_EXEC"
185fi
186
187
188PISM="${PISM_PREFIX}${PISM_EXEC}"
189
190cmd="$PISM_MPIDO $NN $PISM -ys -1000.0 -ye 0 -skip 5 -i green_ssl2_110ka.nc -atmosphere searise_greenland \
191 -surface pdd -pdd_fausto \
192 -o no_force.nc -ts_file ts_no_force.nc -ts_times -1000:yearly:0"
193$PISM_DO $cmd
194
195echo
196
197cmd="$PISM_MPIDO $NN $PISM -ys -1000.0 -ye 0 -skip 5 -i green_ssl2_110ka.nc -atmosphere searise_greenland \
198 -surface pdd,forcing -pdd_fausto -force_to_thickness_file green20km_y1.nc \
199 -o default_force.nc -ts_file ts_default_force.nc -ts_times -1000:yearly:0"
200$PISM_DO $cmd
201
202echo
203
204cmd="$PISM_MPIDO $NN $PISM -ys -1000.0 -ye 0 -skip 5 -i green_ssl2_110ka.nc -atmosphere searise_greenland \
205 -surface pdd,forcing -pdd_fausto -force_to_thickness_file green20km_y1.nc -force_to_thickness_alpha 0.005 \
206 -o weak_force.nc -ts_file ts_weak_force.nc -ts_times -1000:yearly:0"
207$PISM_DO $cmd
208
209
210cmd="$PISM_MPIDO $NN $PISM -ys -1000.0 -ye 0 -skip 5 -i green_ssl2_110ka.nc -atmosphere searise_greenland \
211 -surface pdd,forcing -pdd_fausto -force_to_thickness_file green20km_y1.nc -force_to_thickness_alpha 0.05 \
212 -o strong_force.nc -ts_file ts_strong_force.nc -ts_times -1000:yearly:0"
213$PISM_DO $cmd
214
215\endcode
216The script also has a run with no forcing, one with forcing at a lower alpha value,
217a factor of five smaller than the default, and one with a forcing at a higher alpha value, a factor of five higher.
218 */
220 const array::Scalar &ice_thickness,
221 const array::CellType &cell_type,
222 array::Scalar &result) const {
223
224 if (time < m_start_time) {
225 return;
226 }
227
228 m_log->message(5,
229 " updating surface mass balance using -force_to_thickness mechanism ...");
230
231 double ice_density = m_config->get_number("constants.ice.density");
232
233 array::AccessScope list{&cell_type, &ice_thickness,
234 &m_target_thickness, &m_ftt_mask, &result};
235
236 for (auto p = m_grid->points(); p; p.next()) {
237 const int i = p.i(), j = p.j();
238
239 if (m_ftt_mask(i,j) > 0.5 and cell_type.grounded(i, j)) {
241 result(i,j) += ice_density * m_alpha * (m_target_thickness(i,j) - ice_thickness(i,j));
242 } else {
243 result(i,j) += ice_density * m_alpha * m_alpha_ice_free_factor * (m_target_thickness(i,j) - ice_thickness(i,j));
244 }
245 }
246 }
247 // no communication needed
248}
249
250void ForceThickness::update_impl(const Geometry &geometry, double t, double dt) {
251 m_input_model->update(geometry, t, dt);
252
253 m_mass_flux->copy_from(m_input_model->mass_flux());
254
256 geometry.ice_thickness,
257 geometry.cell_type,
258 *m_mass_flux);
259
263}
264
268
272
274 return *m_melt;
275}
276
278 return *m_runoff;
279}
280
281/*!
282The timestep restriction is, by direct analogy, the same as for
283 \f[\frac{dy}{dt} = - \alpha y\f]
284with explicit (forward) Euler. If \f$\Delta t\f$ is the time step then Euler is
285\f$y_{n+1} = (1-\alpha \Delta t) y_n\f$. We require for stability that
286\f$|y_{n+1}|\le |y_n|\f$, which is to say \f$|1-\alpha \Delta t|\le 1\f$.
287Equivalently (since \f$\alpha \Delta t>0\f$),
288 \f[\alpha \Delta t\le 2\f]
289Therefore we set here
290 \f[\Delta t = \frac{2}{\alpha}.\f]
291 */
293 double max_dt = 2.0 / m_alpha;
294 MaxTimestep input_max_dt = m_input_model->max_timestep(my_t);
295
296 return std::min(input_max_dt, MaxTimestep(max_dt, "surface forcing"));
297}
298
302
303 if (m_input_model != NULL) {
304 m_input_model->define_model_state(output);
305 }
306}
307
309 m_ftt_mask.write(output);
311
312 if (m_input_model != NULL) {
313 m_input_model->write_model_state(output);
314 }
315}
316
317} // end of namespace surface
318} // end of namespace pism
const units::System::Ptr m_sys
unit system used by this component
Definition Component.hh:160
const Time & time() const
Definition Component.cc:109
const Config::ConstPtr m_config
configuration database used by this component
Definition Component.hh:158
const Logger::ConstPtr m_log
logger (for easy access)
Definition Component.hh:162
const std::shared_ptr< const Grid > m_grid
grid used by this component
Definition Component.hh:156
High-level PISM I/O class.
Definition File.hh:55
array::CellType2 cell_type
Definition Geometry.hh:55
array::Scalar2 ice_thickness
Definition Geometry.hh:51
Combines the max. time step with the flag indicating if a restriction is active. Makes is possible to...
VariableMetadata & long_name(const std::string &input)
VariableMetadata & units(const std::string &input)
VariableMetadata & set_name(const std::string &name)
VariableMetadata & standard_name(const std::string &input)
VariableMetadata & set_output_type(io::Type type)
VariableMetadata & set_time_independent(bool flag)
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
Definition Array.hh:64
void set_interpolation_type(InterpolationType type)
Definition Array.cc:178
void define(const File &file, io::Type default_type) const
Define variables corresponding to an Array in a file opened using file.
Definition Array.cc:463
void write(const std::string &filename) const
Definition Array.cc:722
void set(double c)
Result: v[j] <- c for all j.
Definition Array.cc:629
void regrid(const std::string &filename, io::Default default_value)
Definition Array.cc:736
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition Array.cc:476
bool grounded(int i, int j) const
Definition CellType.hh:38
"Cell type" mask. Adds convenience methods to array::Scalar.
Definition CellType.hh:30
static Default Nil()
Definition IO_Flags.hh:93
void adjust_mass_flux(double time, const array::Scalar &ice_thickness, const array::CellType &cell_type, array::Scalar &result) const
const array::Scalar & melt_impl() const
const array::Scalar & accumulation_impl() const
std::shared_ptr< array::Scalar > m_mass_flux
void write_model_state_impl(const File &output) const
The default (empty implementation).
MaxTimestep max_timestep_impl(double t) const
void define_model_state_impl(const File &output) const
The default (empty implementation).
void init_impl(const Geometry &geometry)
const array::Scalar & runoff_impl() const
void update_impl(const Geometry &geometry, double t, double dt)
ForceThickness(std::shared_ptr< const Grid > g, std::shared_ptr< SurfaceModel > input)
const array::Scalar & mass_flux_impl() const
static std::shared_ptr< array::Scalar > allocate_runoff(std::shared_ptr< const Grid > grid)
static std::shared_ptr< array::Scalar > allocate_mass_flux(std::shared_ptr< const Grid > grid)
void dummy_accumulation(const array::Scalar &smb, array::Scalar &result)
std::shared_ptr< array::Scalar > m_melt
static std::shared_ptr< array::Scalar > allocate_accumulation(std::shared_ptr< const Grid > grid)
static std::shared_ptr< array::Scalar > allocate_melt(std::shared_ptr< const Grid > grid)
std::shared_ptr< array::Scalar > m_runoff
void dummy_melt(const array::Scalar &smb, array::Scalar &result)
std::shared_ptr< SurfaceModel > m_input_model
std::shared_ptr< array::Scalar > m_accumulation
void dummy_runoff(const array::Scalar &smb, array::Scalar &result)
The interface of PISM's surface models.
#define PISM_ERROR_LOCATION
@ PISM_GUESS
Definition IO_Flags.hh:56
@ PISM_READONLY
open an existing file for reading only
Definition IO_Flags.hh:68
@ PISM_DOUBLE
Definition IO_Flags.hh:52
double convert(System::Ptr system, double input, const std::string &spec1, const std::string &spec2)
Convert a quantity from unit1 to unit2.
Definition Units.cc:70
static const double g
Definition exactTestP.cc:36