tBedDef.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
---
tBedDef.cc (7644B)
---
1 // Copyright (C) 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 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 "BedDef.hh"
20 #include "pism/util/pism_utilities.hh"
21 #include "pism/util/Time.hh"
22 #include "pism/util/Vars.hh"
23 #include "pism/util/IceGrid.hh"
24 #include "pism/util/ConfigInterface.hh"
25
26 namespace pism {
27 namespace bed {
28
29 BedDef::BedDef(IceGrid::ConstPtr g)
30 : Component(g) {
31
32 const unsigned int WIDE_STENCIL = m_config->get_number("grid.max_stencil_width");
33
34 m_topg.create(m_grid, "topg", WITH_GHOSTS, WIDE_STENCIL);
35 m_topg.set_attrs("model_state", "bedrock surface elevation",
36 "m", "m", "bedrock_altitude", 0);
37
38 m_topg_last.create(m_grid, "topg", WITH_GHOSTS, WIDE_STENCIL);
39 m_topg_last.set_attrs("model_state", "bedrock surface elevation",
40 "m", "m", "bedrock_altitude", 0);
41
42 m_uplift.create(m_grid, "dbdt", WITHOUT_GHOSTS);
43 m_uplift.set_attrs("model_state", "bedrock uplift rate",
44 "m s-1", "mm year-1", "tendency_of_bedrock_altitude", 0);
45 }
46
47 BedDef::~BedDef() {
48 // empty
49 }
50
51 const IceModelVec2S& BedDef::bed_elevation() const {
52 return m_topg;
53 }
54
55 const IceModelVec2S& BedDef::uplift() const {
56 return m_uplift;
57 }
58
59 void BedDef::define_model_state_impl(const File &output) const {
60 m_uplift.define(output);
61 m_topg.define(output);
62 }
63
64 void BedDef::write_model_state_impl(const File &output) const {
65 m_uplift.write(output);
66 m_topg.write(output);
67 }
68
69 DiagnosticList BedDef::diagnostics_impl() const {
70 DiagnosticList result;
71 result = {
72 {"dbdt", Diagnostic::wrap(m_uplift)},
73 {"topg", Diagnostic::wrap(m_topg)}
74 };
75
76 return result;
77 }
78
79 void BedDef::init(const InputOptions &opts, const IceModelVec2S &ice_thickness,
80 const IceModelVec2S &sea_level_elevation) {
81 this->init_impl(opts, ice_thickness, sea_level_elevation);
82 }
83
84 //! Initialize using provided bed elevation and uplift.
85 void BedDef::bootstrap(const IceModelVec2S &bed_elevation,
86 const IceModelVec2S &bed_uplift,
87 const IceModelVec2S &ice_thickness,
88 const IceModelVec2S &sea_level_elevation) {
89 this->bootstrap_impl(bed_elevation, bed_uplift, ice_thickness, sea_level_elevation);
90 }
91
92 void BedDef::bootstrap_impl(const IceModelVec2S &bed_elevation,
93 const IceModelVec2S &bed_uplift,
94 const IceModelVec2S &ice_thickness,
95 const IceModelVec2S &sea_level_elevation) {
96 m_topg.copy_from(bed_elevation);
97 m_uplift.copy_from(bed_uplift);
98
99 // suppress a compiler warning:
100 (void) ice_thickness;
101 (void) sea_level_elevation;
102 }
103
104 void BedDef::update(const IceModelVec2S &ice_thickness,
105 const IceModelVec2S &sea_level_elevation,
106 double t, double dt) {
107 this->update_impl(ice_thickness, sea_level_elevation, t, dt);
108 }
109
110 //! Initialize from the context (input file and the "variables" database).
111 void BedDef::init_impl(const InputOptions &opts, const IceModelVec2S &ice_thickness,
112 const IceModelVec2S &sea_level_elevation) {
113 (void) ice_thickness;
114 (void) sea_level_elevation;
115
116 switch (opts.type) {
117 case INIT_RESTART:
118 // read bed elevation and uplift rate from file
119 m_log->message(2,
120 " reading bed topography and uplift from %s ... \n",
121 opts.filename.c_str());
122 // re-starting
123 m_topg.read(opts.filename, opts.record); // fails if not found!
124 m_uplift.read(opts.filename, opts.record); // fails if not found!
125 break;
126 case INIT_BOOTSTRAP:
127 // bootstrapping
128 m_topg.regrid(opts.filename, OPTIONAL,
129 m_config->get_number("bootstrapping.defaults.bed"));
130 m_uplift.regrid(opts.filename, OPTIONAL,
131 m_config->get_number("bootstrapping.defaults.uplift"));
132 break;
133 case INIT_OTHER:
134 default:
135 {
136 // do nothing
137 }
138 }
139
140 // process -regrid_file and -regrid_vars
141 regrid("bed deformation", m_topg);
142 // uplift is not a part of the model state, but the user may want to take it from a -regrid_file
143 // during bootstrapping
144 regrid("bed deformation", m_uplift);
145
146 std::string uplift_file = m_config->get_string("bed_deformation.bed_uplift_file");
147 if (not uplift_file.empty()) {
148 m_log->message(2,
149 " reading bed uplift from %s ... \n",
150 uplift_file.c_str());
151 m_uplift.regrid(uplift_file, CRITICAL);
152 }
153
154 std::string correction_file = m_config->get_string("bed_deformation.bed_topography_delta_file");
155 if (not correction_file.empty()) {
156 apply_topg_offset(correction_file);
157 }
158
159 // this should be the last thing we do here
160 m_topg_last.copy_from(m_topg);
161 }
162
163 /*!
164 * Apply a correction to the bed topography by reading topg_delta from filename.
165 */
166 void BedDef::apply_topg_offset(const std::string &filename) {
167 m_log->message(2, " Adding a bed topography correction read in from %s...\n",
168 filename.c_str());
169
170 IceModelVec2S topg_delta;
171 topg_delta.create(m_grid, "topg_delta", WITHOUT_GHOSTS);
172 topg_delta.set_attrs("internal", "bed topography correction",
173 "meters", "meters", "", 0);
174
175 topg_delta.regrid(filename, CRITICAL);
176
177 m_topg.add(1.0, topg_delta);
178 }
179
180 //! Compute bed uplift (dt is in seconds).
181 void BedDef::compute_uplift(const IceModelVec2S &bed, const IceModelVec2S &bed_last,
182 double dt, IceModelVec2S &result) {
183 bed.add(-1, bed_last, result);
184 //! uplift = (topg - topg_last) / dt
185 result.scale(1.0 / dt);
186 }
187
188 double compute_load(double bed, double ice_thickness, double sea_level,
189 double ice_density, double ocean_density) {
190
191 double
192 ice_load = ice_thickness,
193 ocean_depth = std::max(sea_level - bed, 0.0),
194 ocean_load = (ocean_density / ice_density) * ocean_depth;
195
196 // this excludes the load of ice shelves
197 return ice_load > ocean_load ? ice_load : 0.0;
198 }
199
200 /*! Compute the load on the bedrock in units of ice-equivalent thickness.
201 *
202 */
203 void compute_load(const IceModelVec2S &bed_elevation,
204 const IceModelVec2S &ice_thickness,
205 const IceModelVec2S &sea_level_elevation,
206 IceModelVec2S &result) {
207
208 Config::ConstPtr config = result.grid()->ctx()->config();
209
210 const double
211 ice_density = config->get_number("constants.ice.density"),
212 ocean_density = config->get_number("constants.sea_water.density");
213
214 IceModelVec::AccessList list{&bed_elevation, &ice_thickness, &sea_level_elevation, &result};
215
216 for (Points p(*result.grid()); p; p.next()) {
217 const int i = p.i(), j = p.j();
218
219 result(i, j) = compute_load(bed_elevation(i, j),
220 ice_thickness(i, j),
221 sea_level_elevation(i, j),
222 ice_density, ocean_density);
223 }
224 }
225
226 } // end of namespace bed
227 } // end of namespace pism