tBedThermalUnit.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
---
tBedThermalUnit.cc (7986B)
---
1 // Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 Ed Bueler 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 <gsl/gsl_math.h> // GSL_NAN
20
21 #include "BedThermalUnit.hh"
22 #include "pism/util/io/File.hh"
23 #include "pism/util/Vars.hh"
24 #include "pism/util/IceGrid.hh"
25 #include "pism/util/pism_options.hh"
26 #include "pism/util/ConfigInterface.hh"
27 #include "pism/util/error_handling.hh"
28 #include "pism/util/MaxTimestep.hh"
29 #include "pism/util/pism_utilities.hh"
30
31 #include "BTU_Full.hh"
32 #include "BTU_Minimal.hh"
33
34 namespace pism {
35 namespace energy {
36
37 BTUGrid::BTUGrid(Context::ConstPtr ctx) {
38 Mbz = (unsigned int) ctx->config()->get_number("grid.Mbz");
39 Lbz = ctx->config()->get_number("grid.Lbz");
40 }
41
42
43 BTUGrid BTUGrid::FromOptions(Context::ConstPtr ctx) {
44 BTUGrid result(ctx);
45
46 Config::ConstPtr config = ctx->config();
47 InputOptions opts = process_input_options(ctx->com(), config);
48
49 if (opts.type == INIT_RESTART) {
50 // If we're initializing from a file we need to get the number of bedrock
51 // levels and the depth of the bed thermal layer from it:
52 File input_file(ctx->com(), opts.filename, PISM_NETCDF3, PISM_READONLY);
53
54 if (input_file.find_variable("litho_temp")) {
55 grid_info info(input_file, "litho_temp", ctx->unit_system(),
56 CELL_CENTER); // grid registration is irrelevant
57
58 result.Mbz = info.z_len;
59 result.Lbz = -info.z_min;
60 } else {
61 // override values we got using config.get_number() in the constructor
62 result.Mbz = 1;
63 result.Lbz = 0;
64 }
65
66 input_file.close();
67 } else {
68 // Bootstrapping or initializing without an input file.
69 result.Mbz = config->get_number("grid.Mbz");
70 result.Lbz = config->get_number("grid.Lbz");
71
72 if (result.Mbz == 1) {
73 result.Lbz = 0;
74 result.Mbz = 1;
75 }
76 }
77
78 return result;
79 }
80
81 /*! Allocate a complete or minimal bedrock thermal unit depending on the number of bedrock levels.
82 *
83 */
84 BedThermalUnit* BedThermalUnit::FromOptions(IceGrid::ConstPtr grid,
85 Context::ConstPtr ctx) {
86
87 BTUGrid bedrock_grid = BTUGrid::FromOptions(ctx);
88
89 if (bedrock_grid.Mbz > 1) {
90 return new energy::BTU_Full(grid, bedrock_grid);
91 } else {
92 return new energy::BTU_Minimal(grid);
93 }
94 }
95
96
97 BedThermalUnit::BedThermalUnit(IceGrid::ConstPtr g)
98 : Component(g) {
99
100 {
101 m_top_surface_flux.create(m_grid, "heat_flux_from_bedrock", WITHOUT_GHOSTS);
102 m_top_surface_flux.set_attrs("diagnostic", "upward geothermal flux at the top bedrock surface",
103 "W m-2", "mW m-2",
104 "upward_geothermal_heat_flux_at_ground_level_in_land_ice", 0);
105 m_top_surface_flux.metadata().set_string("comment", "positive values correspond to an upward flux");
106 }
107 {
108 m_bottom_surface_flux.create(m_grid, "bheatflx", WITHOUT_GHOSTS);
109 // PROPOSED standard_name = lithosphere_upward_heat_flux
110 m_bottom_surface_flux.set_attrs("model_state",
111 "upward geothermal flux at the bottom bedrock surface",
112 "W m-2", "mW m-2", "", 0);
113
114 m_bottom_surface_flux.metadata().set_string("comment", "positive values correspond to an upward flux");
115 m_bottom_surface_flux.set_time_independent(true);
116 }
117 }
118
119 BedThermalUnit::~BedThermalUnit() {
120 // empty
121 }
122
123 void BedThermalUnit::init(const InputOptions &opts) {
124 this->init_impl(opts);
125 }
126
127 //! \brief Initialize the bedrock thermal unit.
128 void BedThermalUnit::init_impl(const InputOptions &opts) {
129 auto input_file = m_config->get_string("energy.bedrock_thermal.file");
130
131 if (not input_file.empty()) {
132 m_log->message(2, " - Reading geothermal flux from '%s' ...\n",
133 input_file.c_str());
134
135 m_bottom_surface_flux.regrid(input_file, CRITICAL);
136 } else {
137 m_log->message(2,
138 " - Parameter %s is not set. Reading geothermal flux from '%s'...\n",
139 "energy.bedrock_thermal.file",
140 opts.filename.c_str());
141
142 switch (opts.type) {
143 case INIT_RESTART:
144 m_bottom_surface_flux.read(opts.filename, opts.record);
145 break;
146 case INIT_BOOTSTRAP:
147 m_bottom_surface_flux.regrid(opts.filename, OPTIONAL,
148 m_config->get_number("bootstrapping.defaults.geothermal_flux"));
149 break;
150 case INIT_OTHER:
151 default:
152 initialize_bottom_surface_flux();
153 }
154 }
155
156 regrid("bedrock thermal layer", m_bottom_surface_flux, REGRID_WITHOUT_REGRID_VARS);
157 }
158
159 void BedThermalUnit::initialize_bottom_surface_flux() {
160 const double heat_flux = m_config->get_number("bootstrapping.defaults.geothermal_flux");
161
162 m_log->message(2, " using constant geothermal flux %f W m-2 ...\n",
163 heat_flux);
164
165 m_bottom_surface_flux.set(heat_flux);
166 }
167
168 /** Returns the vertical spacing used by the bedrock grid.
169 *
170 */
171 double BedThermalUnit::vertical_spacing() const {
172 return this->vertical_spacing_impl();
173 }
174
175 /*!
176 * Return the depth of the bedrock thermal layer.
177 */
178 double BedThermalUnit::depth() const {
179 return this->depth_impl();
180 }
181
182 /*!
183 * Return the number of levels in the bedrock thermal layer.
184 */
185 unsigned int BedThermalUnit::Mz() const {
186 return this->Mz_impl();
187 }
188
189 void BedThermalUnit::define_model_state_impl(const File &output) const {
190 m_bottom_surface_flux.define(output);
191 }
192
193 void BedThermalUnit::write_model_state_impl(const File &output) const {
194 m_bottom_surface_flux.write(output);
195 }
196
197 DiagnosticList BedThermalUnit::diagnostics_impl() const {
198 DiagnosticList result = {
199 {"bheatflx", Diagnostic::wrap(m_bottom_surface_flux)},
200 {"heat_flux_from_bedrock", Diagnostic::Ptr(new BTU_geothermal_flux_at_ground_level(this))}};
201
202 if (m_config->get_flag("output.ISMIP6")) {
203 result["hfgeoubed"] = Diagnostic::Ptr(new BTU_geothermal_flux_at_ground_level(this));
204 }
205 return result;
206 }
207
208 void BedThermalUnit::update(const IceModelVec2S &bedrock_top_temperature,
209 double t, double dt) {
210 this->update_impl(bedrock_top_temperature, t, dt);
211 }
212
213 const IceModelVec2S& BedThermalUnit::flux_through_top_surface() const {
214 return m_top_surface_flux;
215 }
216
217 const IceModelVec2S& BedThermalUnit::flux_through_bottom_surface() const {
218 return m_bottom_surface_flux;
219 }
220
221 BTU_geothermal_flux_at_ground_level::BTU_geothermal_flux_at_ground_level(const BedThermalUnit *m)
222 : Diag<BedThermalUnit>(m) {
223
224 auto ismip6 = m_config->get_flag("output.ISMIP6");
225
226 // set metadata:
227 m_vars = {SpatialVariableMetadata(m_sys, ismip6 ? "hfgeoubed" : "heat_flux_from_bedrock")};
228
229 set_attrs("upward geothermal flux at the top bedrock surface",
230 (ismip6 ?
231 "upward_geothermal_heat_flux_in_land_ice" :
232 "upward_geothermal_heat_flux_at_ground_level_in_land_ice"),
233 "W m-2", "mW m-2", 0);
234 m_vars[0].set_string("comment",
235 "positive values correspond to an upward flux");
236 }
237
238 IceModelVec::Ptr BTU_geothermal_flux_at_ground_level::compute_impl() const {
239 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "hfgeoubed", WITHOUT_GHOSTS));
240 result->metadata() = m_vars[0];
241
242 result->copy_from(model->flux_through_top_surface());
243
244 return result;
245 }
246
247 } // end of namespace energy
248 } // end of namespace pism