PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
BedThermalUnit.cc
Go to the documentation of this file.
1 // Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2023 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 <memory>
20 
21 #include "pism/energy/BedThermalUnit.hh"
22 #include "pism/util/ConfigInterface.hh"
23 #include "pism/util/Grid.hh"
24 #include "pism/util/io/File.hh"
25 
26 #include "pism/energy/BTU_Full.hh"
27 #include "pism/energy/BTU_Minimal.hh"
28 
29 namespace pism {
30 namespace energy {
31 
32 BTUGrid::BTUGrid(std::shared_ptr<const Context> ctx) {
33  Mbz = (unsigned int) ctx->config()->get_number("grid.Mbz");
34  Lbz = ctx->config()->get_number("grid.Lbz");
35 }
36 
37 
38 BTUGrid BTUGrid::FromOptions(std::shared_ptr<const Context> ctx) {
39  BTUGrid result(ctx);
40 
41  Config::ConstPtr config = ctx->config();
42  InputOptions opts = process_input_options(ctx->com(), config);
43 
44  if (opts.type == INIT_RESTART) {
45  // If we're initializing from a file we need to get the number of bedrock
46  // levels and the depth of the bed thermal layer from it:
47  File input_file(ctx->com(), opts.filename, io::PISM_NETCDF3, io::PISM_READONLY);
48 
49  if (input_file.find_variable("litho_temp")) {
50  grid::InputGridInfo info(input_file, "litho_temp", ctx->unit_system(),
51  grid::CELL_CENTER); // grid registration is irrelevant
52 
53  result.Mbz = info.z.size();
54  result.Lbz = -info.z_min;
55  } else {
56  // override values we got using config.get_number() in the constructor
57  result.Mbz = 1;
58  result.Lbz = 0;
59  }
60 
61  input_file.close();
62  } else {
63  // Bootstrapping or initializing without an input file.
64  result.Mbz = config->get_number("grid.Mbz");
65  result.Lbz = config->get_number("grid.Lbz");
66 
67  if (result.Mbz == 1) {
68  result.Lbz = 0;
69  result.Mbz = 1;
70  }
71  }
72 
73  return result;
74 }
75 
76 /*! Allocate a complete or minimal bedrock thermal unit depending on the number of bedrock levels.
77  *
78  */
79 std::shared_ptr<BedThermalUnit> BedThermalUnit::FromOptions(std::shared_ptr<const Grid> grid,
80  std::shared_ptr<const Context> ctx) {
81 
82  auto bedrock_grid = BTUGrid::FromOptions(ctx);
83 
84  if (bedrock_grid.Mbz > 1) {
85  return std::make_shared<energy::BTU_Full>(grid, bedrock_grid);
86  }
87 
88  return std::make_shared<energy::BTU_Minimal>(grid);
89 }
90 
91 
92 BedThermalUnit::BedThermalUnit(std::shared_ptr<const Grid> g)
93  : Component(g),
94  m_bottom_surface_flux(m_grid, "bheatflx"),
95  m_top_surface_flux(m_grid, "heat_flux_from_bedrock") {
96 
97  {
99  .long_name("upward geothermal flux at the top bedrock surface")
100  .units("W m-2")
101  .output_units("mW m-2")
102  .standard_name("upward_geothermal_heat_flux_at_ground_level_in_land_ice");
103  m_top_surface_flux.metadata()["comment"] = "positive values correspond to an upward flux";
104  }
105  {
106  // PROPOSED standard_name = lithosphere_upward_heat_flux
108  .long_name("upward geothermal flux at the bottom bedrock surface")
109  .units("W m-2") // note: don't convert to "mW m-2" when saving
110  .set_time_independent(true);
111 
112  m_bottom_surface_flux.metadata()["comment"] = "positive values correspond to an upward flux";
113  }
114 }
115 
117  this->init_impl(opts);
118 }
119 
120 //! \brief Initialize the bedrock thermal unit.
122  auto input_file = m_config->get_string("energy.bedrock_thermal.file");
123 
124  if (not input_file.empty()) {
125  m_log->message(2, " - Reading geothermal flux from '%s' ...\n",
126  input_file.c_str());
127 
129  } else {
130  m_log->message(2,
131  " - Parameter %s is not set. Reading geothermal flux from '%s'...\n",
132  "energy.bedrock_thermal.file",
133  opts.filename.c_str());
134 
135  switch (opts.type) {
136  case INIT_RESTART:
138  break;
139  case INIT_BOOTSTRAP:
141  "bootstrapping.defaults.geothermal_flux")));
142  break;
143  case INIT_OTHER:
144  default:
146  }
147  }
148 
149  regrid("bedrock thermal layer", m_bottom_surface_flux, REGRID_WITHOUT_REGRID_VARS);
150 }
151 
153  const double heat_flux = m_config->get_number("bootstrapping.defaults.geothermal_flux");
154 
155  m_log->message(2, " using constant geothermal flux %f W m-2 ...\n",
156  heat_flux);
157 
158  m_bottom_surface_flux.set(heat_flux);
159 }
160 
161 /** Returns the vertical spacing used by the bedrock grid.
162  *
163  */
165  return this->vertical_spacing_impl();
166 }
167 
168 /*!
169  * Return the depth of the bedrock thermal layer.
170  */
171 double BedThermalUnit::depth() const {
172  return this->depth_impl();
173 }
174 
175 /*!
176  * Return the number of levels in the bedrock thermal layer.
177  */
178 unsigned int BedThermalUnit::Mz() const {
179  return this->Mz_impl();
180 }
181 
184 }
185 
186 void BedThermalUnit::write_model_state_impl(const File &output) const {
188 }
189 
191  DiagnosticList result = {
193  {"heat_flux_from_bedrock", Diagnostic::Ptr(new BTU_geothermal_flux_at_ground_level(this))}};
194 
195  if (m_config->get_flag("output.ISMIP6")) {
196  result["hfgeoubed"] = Diagnostic::Ptr(new BTU_geothermal_flux_at_ground_level(this));
197  }
198  return result;
199 }
200 
201 void BedThermalUnit::update(const array::Scalar &bedrock_top_temperature,
202  double t, double dt) {
203  this->update_impl(bedrock_top_temperature, t, dt);
204 }
205 
207  return m_top_surface_flux;
208 }
209 
211  return m_bottom_surface_flux;
212 }
213 
215  : Diag<BedThermalUnit>(m) {
216 
217  auto ismip6 = m_config->get_flag("output.ISMIP6");
218  m_vars = { { m_sys, ismip6 ? "hfgeoubed" : "heat_flux_from_bedrock" } };
219  m_vars[0]
220  .long_name("upward geothermal flux at the top bedrock surface")
221  .standard_name((ismip6 ? "upward_geothermal_heat_flux_in_land_ice" :
222  "upward_geothermal_heat_flux_at_ground_level_in_land_ice"))
223  .units("W m-2")
224  .output_units("mW m-2");
225  m_vars[0]["comment"] = "positive values correspond to an upward flux";
226 }
227 
228 std::shared_ptr<array::Array> BTU_geothermal_flux_at_ground_level::compute_impl() const {
229  auto result = std::make_shared<array::Scalar>(m_grid, "hfgeoubed");
230  result->metadata() = m_vars[0];
231 
232  result->copy_from(model->flux_through_top_surface());
233 
234  return result;
235 }
236 
237 } // end of namespace energy
238 } // end of namespace pism
std::shared_ptr< const Grid > grid() const
Definition: Component.cc:105
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
@ REGRID_WITHOUT_REGRID_VARS
Definition: Component.hh:151
void regrid(const std::string &module_name, array::Array &variable, RegriddingFlag flag=NO_REGRID_WITHOUT_REGRID_VARS)
Definition: Component.cc:159
A class defining a common interface for most PISM sub-models.
Definition: Component.hh:118
std::shared_ptr< const Config > ConstPtr
const BedThermalUnit * model
Definition: Diagnostic.hh:172
A template derived from Diagnostic, adding a "Model".
Definition: Diagnostic.hh:166
static Ptr wrap(const T &input)
Definition: Diagnostic.hh:160
const units::System::Ptr m_sys
the unit system
Definition: Diagnostic.hh:116
std::vector< SpatialVariableMetadata > m_vars
metadata corresponding to NetCDF variables
Definition: Diagnostic.hh:120
std::shared_ptr< Diagnostic > Ptr
Definition: Diagnostic.hh:65
std::shared_ptr< const Grid > m_grid
the grid
Definition: Diagnostic.hh:114
const Config::ConstPtr m_config
Configuration flags and parameters.
Definition: Diagnostic.hh:118
VariableLookupData find_variable(const std::string &short_name, const std::string &std_name) const
Find a variable using its standard name and/or short name.
Definition: File.cc:361
void close()
Definition: File.cc:270
High-level PISM I/O class.
Definition: File.hh:56
VariableMetadata & standard_name(const std::string &input)
VariableMetadata & long_name(const std::string &input)
VariableMetadata & output_units(const std::string &input)
VariableMetadata & units(const std::string &input)
VariableMetadata & set_time_independent(bool flag)
void read(const std::string &filename, unsigned int time)
Definition: Array.cc:809
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:540
void write(const std::string &filename) const
Definition: Array.cc:800
void set(double c)
Result: v[j] <- c for all j.
Definition: Array.cc:707
void regrid(const std::string &filename, io::Default default_value)
Definition: Array.cc:814
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition: Array.cc:553
BTU_geothermal_flux_at_ground_level(const BedThermalUnit *m)
virtual std::shared_ptr< array::Array > compute_impl() const
BedThermalUnit(std::shared_ptr< const Grid > g)
virtual void initialize_bottom_surface_flux()
void update(const array::Scalar &bedrock_top_temperature, double t, double dt)
const array::Scalar & flux_through_top_surface() const
Return the upward heat flux through the top surface of the bedrock thermal layer.
array::Scalar m_bottom_surface_flux
upward heat flux through the bottom surface of the bed thermal layer
virtual void init_impl(const InputOptions &opts)
Initialize the bedrock thermal unit.
void init(const InputOptions &opts)
virtual double depth_impl() const =0
static std::shared_ptr< BedThermalUnit > FromOptions(std::shared_ptr< const Grid > g, std::shared_ptr< const Context > ctx)
virtual void define_model_state_impl(const File &output) const
The default (empty implementation).
virtual void write_model_state_impl(const File &output) const
The default (empty implementation).
virtual double vertical_spacing_impl() const =0
virtual DiagnosticList diagnostics_impl() const
virtual unsigned int Mz_impl() const =0
const array::Scalar & flux_through_bottom_surface() const
Return the upward heat flux through the bottom surface of the bedrock thermal layer.
array::Scalar m_top_surface_flux
upward heat flux through the top surface of the bed thermal layer
virtual void update_impl(const array::Scalar &bedrock_top_temperature, double t, double dt)=0
Given the temperature of the top of the bedrock, for the duration of one time-step,...
double z_min
minimal value of the z dimension
Definition: Grid.hh:81
std::vector< double > z
z coordinates
Definition: Grid.hh:90
Contains parameters of an input file grid.
Definition: Grid.hh:65
static Default Nil()
Definition: IO_Flags.hh:97
@ CELL_CENTER
Definition: Grid.hh:53
@ PISM_NETCDF3
Definition: IO_Flags.hh:57
@ PISM_READONLY
open an existing file for reading only
Definition: IO_Flags.hh:72
@ PISM_DOUBLE
Definition: IO_Flags.hh:52
InputOptions process_input_options(MPI_Comm com, Config::ConstPtr config)
Definition: Component.cc:43
static const double g
Definition: exactTestP.cc:36
@ INIT_BOOTSTRAP
Definition: Component.hh:56
@ INIT_OTHER
Definition: Component.hh:56
@ INIT_RESTART
Definition: Component.hh:56
std::map< std::string, Diagnostic::Ptr > DiagnosticList
Definition: Diagnostic.hh:125
InitializationType type
initialization type
Definition: Component.hh:61
std::string filename
name of the input file (if applicable)
Definition: Component.hh:63
unsigned int record
index of the record to re-start from
Definition: Component.hh:65
static BTUGrid FromOptions(std::shared_ptr< const Context > ctx)
BTUGrid(std::shared_ptr< const Context > ctx)