PISM, A Parallel Ice Sheet Model  stable v2.0.5 committed by Constantine Khrulev on 2022-10-14 09:56:26 -0800
Elevation.cc
Go to the documentation of this file.
1 // Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018 Andy Aschwanden 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 "Elevation.hh"
20 
21 #include "pism/util/iceModelVec.hh"
22 #include "pism/util/io/File.hh"
23 #include "pism/util/Vars.hh"
24 #include "pism/util/IceGrid.hh"
25 #include "pism/util/ConfigInterface.hh"
26 #include "pism/util/error_handling.hh"
27 #include "pism/util/pism_options.hh"
28 #include "pism/util/io/io_helpers.hh"
29 #include "pism/util/MaxTimestep.hh"
30 #include "pism/util/pism_utilities.hh"
31 #include "pism/geometry/Geometry.hh"
32 
33 namespace pism {
34 namespace surface {
35 
36 ///// Elevation-dependent temperature and surface mass balance.
37 Elevation::Elevation(IceGrid::ConstPtr grid, std::shared_ptr<atmosphere::AtmosphereModel> input)
38  : SurfaceModel(grid) {
39  (void) input;
40 
43 }
44 
45 void Elevation::init_impl(const Geometry &geometry) {
46  (void) geometry;
47 
48  bool limits_set = false;
49 
50  m_log->message(2,
51  "* Initializing the constant-in-time surface processes model Elevation. Setting...\n");
52 
53  units::Converter K(m_sys, "Celsius", "Kelvin");
54  // options
55  {
56  // ice surface temperature
57  {
58  options::RealList IST("-ice_surface_temp", "ice surface temperature parameterization",
59  {-5.0, 0.0, 1325.0, 1350.0});
60 
61  if (IST->size() != 4) {
62  throw RuntimeError(PISM_ERROR_LOCATION, "option -ice_surface_temp requires an argument"
63  " (comma-separated list of 4 numbers)");
64  }
65  m_T_min = K(IST[0]);
66  m_T_max = K(IST[1]);
67  m_z_T_min = IST[2];
68  m_z_T_max = IST[3];
69  }
70 
71  // climatic mass balance
72  units::Converter meter_per_second(m_sys, "m year-1", "m second-1");
73  {
74  options::RealList CMB("-climatic_mass_balance",
75  "climatic mass balance parameterization",
76  {-3.0, 4.0, 1100.0, 1450.0, 1700.0});
77  if (CMB->size() != 5) {
78  throw RuntimeError(PISM_ERROR_LOCATION, "-climatic_mass_balance requires an argument"
79  " (comma-separated list of 5 numbers)");
80  }
81  m_M_min = meter_per_second(CMB[0]);
82  m_M_max = meter_per_second(CMB[1]);
83  m_z_M_min = CMB[2];
84  m_z_ELA = CMB[3];
85  m_z_M_max = CMB[4];
86  }
87 
88  // limits of the climatic mass balance
89  {
90  options::RealList limits("-climatic_mass_balance_limits",
91  "lower and upper limits of the climatic mass balance", {});
92  limits_set = limits.is_set();
93  if (limits.is_set()) {
94  if (limits->size() != 2) {
95  throw RuntimeError(PISM_ERROR_LOCATION, "-climatic_mass_balance_limits requires an argument"
96  " (a comma-separated list of 2 numbers)");
97  }
98  m_M_limit_min = meter_per_second(limits[0]);
99  m_M_limit_max = meter_per_second(limits[1]);
100  } else {
103  }
104  }
105  }
106 
107  units::Converter meter_per_year(m_sys, "m second-1", "m year-1");
108  m_log->message(3,
109  " temperature at %.0f m a.s.l. = %.2f deg C\n"
110  " temperature at %.0f m a.s.l. = %.2f deg C\n"
111  " mass balance below %.0f m a.s.l. = %.2f m year-1\n"
112  " mass balance at %.0f m a.s.l. = %.2f m year-1\n"
113  " mass balance at %.0f m a.s.l. = %.2f m year-1\n"
114  " mass balance above %.0f m a.s.l. = %.2f m year-1\n"
115  " equilibrium line altitude z_ELA = %.2f m a.s.l.\n",
118  m_z_M_min, meter_per_year(m_M_limit_min),
120  m_z_M_max, meter_per_year(m_M_max),
121  m_z_M_max, meter_per_year(m_M_limit_max),
122  m_z_ELA);
123 
124  // parameterizing the ice surface temperature 'ice_surface_temp'
125  m_log->message(2, " - parameterizing the ice surface temperature 'ice_surface_temp' ... \n");
126  m_log->message(2,
127  " ice temperature at the ice surface (T = ice_surface_temp) is piecewise-linear function\n"
128  " of surface altitude (usurf):\n"
129  " / %2.2f K for usurf < %.f m\n"
130  " T = | %5.2f K + %5.3f * (usurf - %.f m) for %.f m < usurf < %.f m\n"
131  " \\ %5.2f K for %.f m < usurf\n",
134  m_T_max, m_z_T_max);
135 
136  // parameterizing the ice surface mass balance 'climatic_mass_balance'
137  m_log->message(2,
138 
139  " - parameterizing the ice surface mass balance 'climatic_mass_balance' ... \n");
140 
141  if (limits_set) {
142  m_log->message(2,
143  " - option '-climatic_mass_balance_limits' seen, limiting upper and lower bounds ... \n");
144  }
145 
146  m_log->message(2,
147  " surface mass balance (M = climatic_mass_balance) is piecewise-linear function\n"
148  " of surface altitue (usurf):\n"
149  " / %5.2f m year-1 for usurf < %3.f m\n"
150  " M = | %5.3f 1/a * (usurf-%.0f m) for %3.f m < usurf < %3.f m\n"
151  " \\ %5.3f 1/a * (usurf-%.0f m) for %3.f m < usurf < %3.f m\n"
152  " \\ %5.2f m year-1 for %3.f m < usurf\n",
153  meter_per_year(m_M_limit_min), m_z_M_min,
154  meter_per_year(-m_M_min)/(m_z_ELA - m_z_M_min), m_z_ELA, m_z_M_min, m_z_ELA,
155  meter_per_year(m_M_max)/(m_z_M_max - m_z_ELA), m_z_ELA, m_z_ELA, m_z_M_max,
156  meter_per_year(m_M_limit_max), m_z_M_max);
157 }
158 
160  (void) t;
161  return MaxTimestep("surface 'elevation'");
162 }
163 
164 void Elevation::update_impl(const Geometry &geometry, double t, double dt) {
165  (void) t;
166  (void) dt;
167 
170 
174 
175 }
176 
178  return *m_mass_flux;
179 }
180 
182  return *m_temperature;
183 }
184 
186  return *m_accumulation;
187 }
188 
190  return *m_melt;
191 }
192 
194  return *m_runoff;
195 }
196 
197 void Elevation::compute_mass_flux(const IceModelVec2S &surface, IceModelVec2S &result) const {
198  double dabdz = -m_M_min/(m_z_ELA - m_z_M_min);
199  double dacdz = m_M_max/(m_z_M_max - m_z_ELA);
200 
201  IceModelVec::AccessList list{&result, &surface};
202 
203  ParallelSection loop(m_grid->com);
204  try {
205  for (Points p(*m_grid); p; p.next()) {
206  const int i = p.i(), j = p.j();
207 
208  double z = surface(i, j);
209 
210  if (z < m_z_M_min) {
211  result(i, j) = m_M_limit_min;
212  }
213  else if ((z >= m_z_M_min) and (z < m_z_ELA)) {
214  result(i, j) = dabdz * (z - m_z_ELA);
215  }
216  else if ((z >= m_z_ELA) and (z <= m_z_M_max)) {
217  result(i, j) = dacdz * (z - m_z_ELA);
218  }
219  else if (z > m_z_M_max) {
220  result(i, j) = m_M_limit_max;
221  }
222  else {
224  "Elevation::compute_mass_flux: HOW DID I GET HERE?");
225  }
226  }
227  } catch (...) {
228  loop.failed();
229  }
230  loop.check();
231 
232  // convert from m second-1 ice equivalent to kg m-2 s-1:
233  result.scale(m_config->get_number("constants.ice.density"));
234 }
235 
236 void Elevation::compute_temperature(const IceModelVec2S &surface, IceModelVec2S &result) const {
237 
238  IceModelVec::AccessList list{&result, &surface};
239 
240  double dTdz = (m_T_max - m_T_min)/(m_z_T_max - m_z_T_min);
241  ParallelSection loop(m_grid->com);
242  try {
243  for (Points p(*m_grid); p; p.next()) {
244  const int i = p.i(), j = p.j();
245 
246  double z = surface(i, j);
247 
248  if (z <= m_z_T_min) {
249  result(i, j) = m_T_min;
250  }
251  else if ((z > m_z_T_min) and (z < m_z_T_max)) {
252  result(i, j) = m_T_min + dTdz * (z - m_z_T_min);
253  }
254  else if (z >= m_z_T_max) {
255  result(i, j) = m_T_max;
256  }
257  else {
259  "Elevation::compute_temperature: HOW DID I GET HERE?");
260  }
261  }
262  } catch (...) {
263  loop.failed();
264  }
265  loop.check();
266 }
267 
268 } // end of namespace surface
269 } // end of namespace pism
Makes sure that we call begin_access() and end_access() for all accessed IceModelVecs.
Definition: iceModelVec.hh:59
IceGrid::ConstPtr grid() const
Definition: Component.cc:105
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
IceModelVec2S ice_surface_elevation
Definition: Geometry.hh:59
std::shared_ptr< const IceGrid > ConstPtr
Definition: IceGrid.hh:233
void scale(double alpha)
Result: v <- v * alpha. Calls VecScale.
Definition: iceModelVec.cc:252
Combines the max. time step with the flag indicating if a restriction is active. Makes is possible to...
Definition: MaxTimestep.hh:31
void failed()
Indicates a failure of a parallel section.
bool is_set() const
Definition: options.hh:35
const IceModelVec2S & mass_flux_impl() const
Definition: Elevation.cc:177
virtual const IceModelVec2S & accumulation_impl() const
Definition: Elevation.cc:185
Elevation(IceGrid::ConstPtr grid, std::shared_ptr< atmosphere::AtmosphereModel > input)
Definition: Elevation.cc:37
virtual const IceModelVec2S & runoff_impl() const
Definition: Elevation.cc:193
void init_impl(const Geometry &geometry)
Definition: Elevation.cc:45
void update_impl(const Geometry &geometry, double t, double dt)
Definition: Elevation.cc:164
void compute_temperature(const IceModelVec2S &surface, IceModelVec2S &result) const
Definition: Elevation.cc:236
void compute_mass_flux(const IceModelVec2S &surface, IceModelVec2S &result) const
Definition: Elevation.cc:197
IceModelVec2S::Ptr m_mass_flux
Definition: Elevation.hh:53
const IceModelVec2S & temperature_impl() const
Definition: Elevation.cc:181
virtual const IceModelVec2S & melt_impl() const
Definition: Elevation.cc:189
MaxTimestep max_timestep_impl(double t) const
Definition: Elevation.cc:159
IceModelVec2S::Ptr m_temperature
Definition: Elevation.hh:54
static IceModelVec2S::Ptr allocate_mass_flux(IceGrid::ConstPtr grid)
Definition: SurfaceModel.cc:76
IceModelVec2S::Ptr m_melt
IceModelVec2S::Ptr m_accumulation
static IceModelVec2S::Ptr allocate_temperature(IceGrid::ConstPtr grid)
Definition: SurfaceModel.cc:93
void dummy_accumulation(const IceModelVec2S &smb, IceModelVec2S &result)
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
#define PISM_ERROR_LOCATION
static double K(double psi_x, double psi_y, double speed, double epsilon)