PISM, A Parallel Ice Sheet Model  stable v2.0.5 committed by Constantine Khrulev on 2022-10-14 09:56:26 -0800
ElevationChange.cc
Go to the documentation of this file.
1 // Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021 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 "ElevationChange.hh"
20 #include "pism/coupler/util/options.hh"
21 #include "pism/coupler/util/lapse_rates.hh"
22 #include "pism/util/io/io_helpers.hh"
23 #include "pism/util/pism_utilities.hh"
24 #include "pism/util/pism_options.hh"
25 #include "pism/geometry/Geometry.hh"
26 
27 namespace pism {
28 namespace surface {
29 
30 ElevationChange::ElevationChange(IceGrid::ConstPtr g, std::shared_ptr<SurfaceModel> in)
31  : SurfaceModel(g, in) {
32 
33  {
34  m_smb_lapse_rate = m_config->get_number("surface.elevation_change.smb.lapse_rate",
35  "(m / s) / m");
36  // convert from [m s-1 / m] to [kg m-2 s-1 / m]
37  m_smb_lapse_rate *= m_config->get_number("constants.ice.density");
38 
39  m_smb_exp_factor = m_config->get_number("surface.elevation_change.smb.exp_factor");
40  }
41 
42  {
43  auto method = m_config->get_string("surface.elevation_change.smb.method");
44  m_smb_method = method == "scale" ? SCALE : SHIFT;
45  }
46 
47  m_temp_lapse_rate = m_config->get_number("surface.elevation_change.temperature_lapse_rate",
48  "K / m");
49 
50  {
51  ForcingOptions opt(*m_grid->ctx(), "surface.elevation_change");
52 
53  unsigned int buffer_size = m_config->get_number("input.forcing.buffer_size");
54 
55  File file(m_grid->com, opt.filename, PISM_NETCDF3, PISM_READONLY);
56 
58  file,
59  "usurf",
60  "", // no standard name
61  buffer_size,
62  opt.periodic,
63  LINEAR);
64  m_reference_surface->set_attrs("climate_forcing", "ice surface elevation",
65  "m", "m", "surface_altitude", 0);
66  }
67 
73 }
74 
75 void ElevationChange::init_impl(const Geometry &geometry) {
76  using units::convert;
77 
78  m_input_model->init(geometry);
79 
80  m_log->message(2,
81  " [using temperature and mass balance lapse corrections]\n");
82 
83  m_log->message(2,
84  " ice upper-surface temperature lapse rate: %3.3f K per km\n",
85  convert(m_sys, m_temp_lapse_rate, "K / m", "K / km"));
86 
87  if (m_smb_method == SHIFT) {
88  double ice_density = m_config->get_number("constants.ice.density");
89  m_log->message(2,
90  " ice-equivalent surface mass balance lapse rate: %3.3f m year-1 per km\n",
91  convert(m_sys, m_smb_lapse_rate, "kg / (m2 second)", "kg / (m2 year)") / ice_density);
92  } else {
93  m_log->message(2,
94  " surface mass balance scaling factor with temperature: %3.3f Kelvin-1\n",
96  }
97 
98  ForcingOptions opt(*m_grid->ctx(), "surface.elevation_change");
99  m_reference_surface->init(opt.filename, opt.periodic);
100 }
101 
102 void ElevationChange::update_impl(const Geometry &geometry, double t, double dt) {
103 
104  m_input_model->update(geometry, t, dt);
105 
106  m_reference_surface->update(t, dt);
107  m_reference_surface->interp(t + 0.5*dt);
108 
109  const IceModelVec2S &surface = geometry.ice_surface_elevation;
110 
111  m_temperature->copy_from(m_input_model->temperature());
114 
115  m_mass_flux->copy_from(m_input_model->mass_flux());
116 
117  switch (m_smb_method) {
118  case SCALE:
119  {
120  IceModelVec::AccessList list{&surface, m_reference_surface.get(), m_mass_flux.get()};
121 
122  for (Points p(*m_grid); p; p.next()) {
123  const int i = p.i(), j = p.j();
124 
125  double dT = -m_temp_lapse_rate * (surface(i, j) - (*m_reference_surface)(i, j));
126 
127  (*m_mass_flux)(i, j) *= exp(m_smb_exp_factor * dT);
128  }
129  }
130  break;
131  default:
132  case SHIFT:
133  {
136  }
137  break;
138  }
139 
140  // This modifier changes m_mass_flux, so we need to compute accumulation, melt, and
141  // runoff.
145 
146 }
147 
149  return *m_mass_flux;
150 }
151 
153  return *m_temperature;
154 }
155 
157  return *m_accumulation;
158 }
159 
161  return *m_melt;
162 }
163 
165  return *m_runoff;
166 }
167 
168 
169 } // end of namespace surface
170 } // end of namespace pism
Makes sure that we call begin_access() and end_access() for all accessed IceModelVecs.
Definition: iceModelVec.hh:59
const units::System::Ptr m_sys
unit system used by this component
Definition: Component.hh:140
const Config::ConstPtr m_config
configuration database used by this component
Definition: Component.hh:138
const Logger::ConstPtr m_log
logger (for easy access)
Definition: Component.hh:142
const IceGrid::ConstPtr m_grid
grid used by this component
Definition: Component.hh:136
High-level PISM I/O class.
Definition: File.hh:51
IceModelVec2S ice_surface_elevation
Definition: Geometry.hh:59
std::shared_ptr< const IceGrid > ConstPtr
Definition: IceGrid.hh:233
static std::shared_ptr< IceModelVec2T > ForcingField(IceGrid::ConstPtr grid, const File &file, const std::string &short_name, const std::string &standard_name, int max_buffer_size, bool periodic, InterpolationType interpolation_type=PIECEWISE_CONSTANT)
const IceModelVec2S & accumulation_impl() const
virtual void init_impl(const Geometry &geometry)
virtual void update_impl(const Geometry &geometry, double t, double dt)
const IceModelVec2S & mass_flux_impl() const
ElevationChange(IceGrid::ConstPtr g, std::shared_ptr< SurfaceModel > in)
const IceModelVec2S & runoff_impl() const
const IceModelVec2S & temperature_impl() const
IceModelVec2S::Ptr m_temperature
const IceModelVec2S & melt_impl() const
std::shared_ptr< IceModelVec2T > m_reference_surface
static IceModelVec2S::Ptr allocate_mass_flux(IceGrid::ConstPtr grid)
Definition: SurfaceModel.cc:76
IceModelVec2S::Ptr m_melt
static IceModelVec2S::Ptr allocate_melt(IceGrid::ConstPtr grid)
IceModelVec2S::Ptr m_accumulation
static IceModelVec2S::Ptr allocate_temperature(IceGrid::ConstPtr grid)
Definition: SurfaceModel.cc:93
std::shared_ptr< SurfaceModel > m_input_model
void dummy_accumulation(const IceModelVec2S &smb, IceModelVec2S &result)
static IceModelVec2S::Ptr allocate_accumulation(IceGrid::ConstPtr grid)
static IceModelVec2S::Ptr allocate_runoff(IceGrid::ConstPtr grid)
IceModelVec2S::Ptr m_runoff
void dummy_runoff(const IceModelVec2S &smb, IceModelVec2S &result)
void dummy_melt(const IceModelVec2S &smb, IceModelVec2S &result)
The interface of PISM's surface models.
Definition: SurfaceModel.hh:45
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:72
void lapse_rate_correction(const IceModelVec2S &surface, const IceModelVec2S &reference_surface, double lapse_rate, IceModelVec2S &result)
Definition: lapse_rates.cc:26
static const double g
Definition: exactTestP.cc:39
@ PISM_NETCDF3
Definition: IO_Flags.hh:41
@ PISM_READONLY
open an existing file for reading only
Definition: IO_Flags.hh:49
std::string filename
Definition: options.hh:33