PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
ElevationChange.cc
Go to the documentation of this file.
1 // Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023 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/ElevationChange.hh"
20 #include "pism/coupler/util/options.hh"
21 #include "pism/coupler/util/lapse_rates.hh"
22 #include "pism/geometry/Geometry.hh"
23 #include "pism/util/array/Forcing.hh"
24 
25 namespace pism {
26 namespace surface {
27 
28 ElevationChange::ElevationChange(std::shared_ptr<const Grid> g, std::shared_ptr<SurfaceModel> in)
29  : SurfaceModel(g, in) {
30 
31  {
32  m_smb_lapse_rate = m_config->get_number("surface.elevation_change.smb.lapse_rate",
33  "(m / s) / m");
34  // convert from [m s-1 / m] to [kg m-2 s-1 / m]
35  m_smb_lapse_rate *= m_config->get_number("constants.ice.density");
36 
37  m_smb_exp_factor = m_config->get_number("surface.elevation_change.smb.exp_factor");
38  }
39 
40  {
41  auto method = m_config->get_string("surface.elevation_change.smb.method");
42  m_smb_method = method == "scale" ? SCALE : SHIFT;
43  }
44 
45  m_temp_lapse_rate = m_config->get_number("surface.elevation_change.temperature_lapse_rate",
46  "K / m");
47 
48  {
49  ForcingOptions opt(*m_grid->ctx(), "surface.elevation_change");
50 
51  unsigned int buffer_size = m_config->get_number("input.forcing.buffer_size");
52 
54 
55  m_reference_surface = std::make_shared<array::Forcing>(m_grid, file, "usurf",
56  "", // no standard name
57  buffer_size, opt.periodic, LINEAR);
58  m_reference_surface->metadata(0)
59  .long_name("ice surface elevation")
60  .units("m")
61  .standard_name("surface_altitude");
62  }
63 
69 }
70 
71 void ElevationChange::init_impl(const Geometry &geometry) {
72  using units::convert;
73 
74  m_input_model->init(geometry);
75 
76  m_log->message(2,
77  " [using temperature and mass balance lapse corrections]\n");
78 
79  m_log->message(2,
80  " ice upper-surface temperature lapse rate: %3.3f K per km\n",
81  convert(m_sys, m_temp_lapse_rate, "K / m", "K / km"));
82 
83  if (m_smb_method == SHIFT) {
84  double ice_density = m_config->get_number("constants.ice.density");
85  m_log->message(2,
86  " ice-equivalent surface mass balance lapse rate: %3.3f m year-1 per km\n",
87  convert(m_sys, m_smb_lapse_rate, "kg / (m2 second)", "kg / (m2 year)") / ice_density);
88  } else {
89  m_log->message(2,
90  " surface mass balance scaling factor with temperature: %3.3f Kelvin-1\n",
92  }
93 
94  ForcingOptions opt(*m_grid->ctx(), "surface.elevation_change");
95  m_reference_surface->init(opt.filename, opt.periodic);
96 }
97 
98 void ElevationChange::update_impl(const Geometry &geometry, double t, double dt) {
99 
100  m_input_model->update(geometry, t, dt);
101 
102  m_reference_surface->update(t, dt);
103  m_reference_surface->interp(t + 0.5*dt);
104 
105  const array::Scalar &surface = geometry.ice_surface_elevation;
106 
107  m_temperature->copy_from(m_input_model->temperature());
110 
111  m_mass_flux->copy_from(m_input_model->mass_flux());
112 
113  switch (m_smb_method) {
114  case SCALE:
115  {
116  array::AccessScope list{&surface, m_reference_surface.get(), m_mass_flux.get()};
117 
118  for (auto p = m_grid->points(); p; p.next()) {
119  const int i = p.i(), j = p.j();
120 
121  double dT = -m_temp_lapse_rate * (surface(i, j) - (*m_reference_surface)(i, j));
122 
123  (*m_mass_flux)(i, j) *= exp(m_smb_exp_factor * dT);
124  }
125  }
126  break;
127  default:
128  case SHIFT:
129  {
132  }
133  break;
134  }
135 
136  // This modifier changes m_mass_flux, so we need to compute accumulation, melt, and
137  // runoff.
141 
142 }
143 
145  return *m_mass_flux;
146 }
147 
149  return *m_temperature;
150 }
151 
153  return *m_accumulation;
154 }
155 
157  return *m_melt;
158 }
159 
161  return *m_runoff;
162 }
163 
164 
165 } // end of namespace surface
166 } // end of namespace pism
const units::System::Ptr m_sys
unit system used by this component
Definition: Component.hh:160
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:56
array::Scalar2 ice_surface_elevation
Definition: Geometry.hh:57
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
Definition: Array.hh:65
virtual void init_impl(const Geometry &geometry)
virtual void update_impl(const Geometry &geometry, double t, double dt)
std::shared_ptr< array::Scalar > m_mass_flux
const array::Scalar & melt_impl() const
const array::Scalar & accumulation_impl() const
std::shared_ptr< array::Forcing > m_reference_surface
std::shared_ptr< array::Scalar > m_temperature
const array::Scalar & temperature_impl() const
const array::Scalar & runoff_impl() const
ElevationChange(std::shared_ptr< const Grid > g, std::shared_ptr< SurfaceModel > in)
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)
Definition: SurfaceModel.cc:73
void dummy_accumulation(const array::Scalar &smb, array::Scalar &result)
std::shared_ptr< array::Scalar > m_melt
static std::shared_ptr< array::Scalar > allocate_temperature(std::shared_ptr< const Grid > grid)
Definition: SurfaceModel.cc:92
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.
Definition: SurfaceModel.hh:42
@ PISM_NETCDF3
Definition: IO_Flags.hh:57
@ PISM_READONLY
open an existing file for reading only
Definition: IO_Flags.hh:72
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
void lapse_rate_correction(const array::Scalar &surface, const array::Scalar &reference_surface, double lapse_rate, array::Scalar &result)
Definition: lapse_rates.cc:26
std::string filename
Definition: options.hh:33