PISM, A Parallel Ice Sheet Model  stable v2.0.6 committed by Constantine Khrulev on 2023-01-23 15:14:38 -0900
ISMIP6Climate.cc
Go to the documentation of this file.
1 // Copyright (C) 2019, 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 "ISMIP6Climate.hh"
20 
21 #include "pism/util/IceGrid.hh"
22 #include "pism/coupler/util/options.hh"
23 #include "pism/util/io/io_helpers.hh"
24 #include "pism/geometry/Geometry.hh"
25 
26 namespace pism {
27 namespace surface {
28 
29 ISMIP6::ISMIP6(IceGrid::ConstPtr grid, std::shared_ptr<atmosphere::AtmosphereModel> input)
30  : SurfaceModel(grid),
31  m_mass_flux_reference(m_grid, "climatic_mass_balance", WITHOUT_GHOSTS),
32  m_temperature_reference(m_grid, "ice_surface_temp", WITHOUT_GHOSTS),
33  m_surface_reference(m_grid, "usurf", WITHOUT_GHOSTS)
34 {
35  (void) input;
36 
37  // allocate model outputs
40 
41  // set metadata of reference fields
42  {
43  m_mass_flux_reference.set_attrs("climate_forcing",
44  "reference surface mass balance rate",
45  "kg m-2 s-1", "kg m-2 year-1",
46  "land_ice_surface_specific_mass_balance_flux", 0);
47 
48  auto smb_max = m_config->get_number("surface.given.smb_max", "kg m-2 second-1");
49  m_mass_flux_reference.metadata()["valid_range"] = {-smb_max, smb_max};
51 
52  m_surface_reference.set_attrs("climate_forcing",
53  "reference surface altitude",
54  "m", "m", "surface_altitude", 0);
55 
56  m_surface_reference.metadata()["valid_range"] = {0.0, m_grid->Lz()};
58 
59  m_temperature_reference.set_attrs("climate_forcing",
60  "reference temperature",
61  "Kelvin", "Kelvin", "", 0);
62 
63  m_temperature_reference.metadata()["valid_range"] = {0.0, 373.15};
65  }
66 
67  // allocate storage for time-dependent inputs
68  ForcingOptions opt(*m_grid->ctx(), "surface.ismip6");
69 
70  {
71  unsigned int buffer_size = m_config->get_number("input.forcing.buffer_size");
72 
73  File file(m_grid->com, opt.filename, PISM_NETCDF3, PISM_READONLY);
74 
75  {
77  file,
78  "climatic_mass_balance_anomaly",
79  "", // no standard name
80  buffer_size,
81  opt.periodic);
82 
83  m_mass_flux_anomaly->set_attrs("climate_forcing",
84  "surface mass balance rate anomaly",
85  "kg m-2 s-1", "kg m-2 year-1",
86  "", 0);
87 
88  }
89 
90  {
92  file,
93  "climatic_mass_balance_gradient",
94  "", // no standard name
95  buffer_size,
96  opt.periodic);
97 
98  m_mass_flux_gradient->set_attrs("climate_forcing",
99  "surface mass balance rate elevation lapse rate",
100  "kg m-2 s-1 m-1", "kg m-2 year-1 m-1",
101  "", 0);
102  }
103 
104  {
106  file,
107  "ice_surface_temp_anomaly",
108  "", // no standard name
109  buffer_size,
110  opt.periodic);
111 
112  m_temperature_anomaly->set_attrs("climate_forcing",
113  "ice surface temperature anomaly",
114  "Kelvin", "Kelvin", "", 0);
115  }
116 
117  {
119  file,
120  "ice_surface_temp_gradient",
121  "", // no standard name
122  buffer_size,
123  opt.periodic);
124 
125  m_temperature_gradient->set_attrs("climate_forcing",
126  "ice surface temperature elevation lapse rate",
127  "Kelvin m-1", "Kelvin m-1", "", 0);
128  }
129 
130  }
131 }
132 
133 void ISMIP6::init_impl(const Geometry &geometry) {
134  (void) geometry;
135 
136  m_log->message(2, "* Initializing the ISMIP6 surface model...\n");
137 
138  {
139  // File with reference surface elevation, temperature, and climatic mass balance
140  auto reference_filename = m_config->get_string("surface.ismip6.reference_file");
141  File reference_file(m_grid->com, reference_filename, PISM_GUESS, PISM_READONLY);
142 
143  m_mass_flux_reference.regrid(reference_file, CRITICAL);
144  m_surface_reference.regrid(reference_file, CRITICAL);
145  m_temperature_reference.regrid(reference_file, CRITICAL);
146  }
147 
148  {
149  ForcingOptions opt(*m_grid->ctx(), "surface.ismip6");
150 
151  m_mass_flux_anomaly->init(opt.filename, opt.periodic);
152  m_mass_flux_gradient->init(opt.filename, opt.periodic);
153 
154  m_temperature_anomaly->init(opt.filename, opt.periodic);
155  m_temperature_gradient->init(opt.filename, opt.periodic);
156  }
157 }
158 
159 void ISMIP6::update_impl(const Geometry &geometry, double t, double dt) {
160 
161  // inputs
162  const IceModelVec2S &h = geometry.ice_surface_elevation;
163  const IceModelVec2S &h_ref = m_surface_reference;
164  const IceModelVec2S &T_ref = m_temperature_reference;
165  const IceModelVec2S &SMB_ref = m_mass_flux_reference;
166 
171 
172  // outputs
175 
176  // get time-dependent input fields at the current time
177  {
178  aT.update(t, dt);
179  aSMB.update(t, dt);
180  dTdz.update(t, dt);
181  dSMBdz.update(t, dt);
182 
183  aT.average(t, dt);
184  aSMB.average(t, dt);
185  dTdz.average(t, dt);
186  dSMBdz.average(t, dt);
187  }
188 
189  // From http://www.climate-cryosphere.org/wiki/index.php?title=ISMIP6-Projections-Greenland:
190  // SMB(x,y,t) = SMB_ref(x,y) + aSMB(x,y,t) + dSMBdz(x,y,t) * [h(x,y,t) - h_ref(x,y)]
191 
192  IceModelVec::AccessList list{&h, &h_ref,
193  &SMB, &SMB_ref, &aSMB, &dSMBdz,
194  &T, &T_ref, &aT, &dTdz};
195 
196  for (Points p(*m_grid); p; p.next()) {
197  const int i = p.i(), j = p.j();
198 
199  SMB(i, j) = SMB_ref(i, j) + aSMB(i, j) + dSMBdz(i, j) * (h(i, j) - h_ref(i, j));
200  T(i, j) = T_ref(i, j) + aT(i, j) + dTdz(i, j) * (h(i, j) - h_ref(i, j));
201  }
202 
204  dummy_melt(SMB, *m_melt);
206 }
207 
209  using std::min;
210 
211  auto dt = m_temperature_anomaly->max_timestep(t);
212  dt = min(dt, m_temperature_gradient->max_timestep(t));
213  dt = min(dt, m_mass_flux_anomaly->max_timestep(t));
214  dt = min(dt, m_mass_flux_gradient->max_timestep(t));
215 
216  if (dt.finite()) {
217  return MaxTimestep(dt.value(), "surface ISMIP6");
218  } else {
219  return MaxTimestep("surface ISMIP6");
220  }
221 }
222 
224  return *m_mass_flux;
225 }
226 
228  return *m_temperature;
229 }
230 
232  return *m_accumulation;
233 }
234 
236  return *m_melt;
237 }
238 
240  return *m_runoff;
241 }
242 
243 } // end of namespace surface
244 } // end of namespace pism
Makes sure that we call begin_access() and end_access() for all accessed IceModelVecs.
Definition: iceModelVec.hh:59
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)
void average(double t, double dt)
void update(double t, double dt)
Read some data to make sure that the interval (t, t + dt) is covered.
A class for storing and accessing 2D time-series (for climate forcing)
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition: iceModelVec.cc:533
void set_attrs(const std::string &pism_intent, const std::string &long_name, const std::string &units, const std::string &glaciological_units, const std::string &standard_name, unsigned int component)
Sets NetCDF attributes of an IceModelVec object.
Definition: iceModelVec.cc:399
void regrid(const std::string &filename, RegriddingFlag flag, double default_value=0.0)
Definition: iceModelVec.cc:838
void set_time_independent(bool flag)
Set the time independent flag for all variables corresponding to this IceModelVec instance.
Definition: iceModelVec.cc:175
Combines the max. time step with the flag indicating if a restriction is active. Makes is possible to...
Definition: MaxTimestep.hh:31
const IceModelVec2S & runoff_impl() const
std::shared_ptr< IceModelVec2T > m_temperature_anomaly
IceModelVec2S m_temperature_reference
std::shared_ptr< IceModelVec2T > m_temperature_gradient
const IceModelVec2S & mass_flux_impl() const
const IceModelVec2S & temperature_impl() const
void update_impl(const Geometry &geometry, double t, double dt)
std::shared_ptr< IceModelVec2T > m_mass_flux_anomaly
std::shared_ptr< IceModelVec2T > m_mass_flux_gradient
void init_impl(const Geometry &geometry)
IceModelVec2S m_surface_reference
MaxTimestep max_timestep_impl(double t) const
const IceModelVec2S & melt_impl() const
IceModelVec2S::Ptr m_mass_flux
ISMIP6(IceGrid::ConstPtr g, std::shared_ptr< atmosphere::AtmosphereModel > input)
IceModelVec2S::Ptr m_temperature
IceModelVec2S m_mass_flux_reference
const IceModelVec2S & accumulation_impl() const
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
@ CRITICAL
Definition: IO_Flags.hh:70
@ PISM_GUESS
Definition: IO_Flags.hh:41
@ PISM_NETCDF3
Definition: IO_Flags.hh:41
double min(const IceModelVec2S &input)
Finds minimum over all the values in an IceModelVec2S object. Ignores ghosts.
@ PISM_READONLY
open an existing file for reading only
Definition: IO_Flags.hh:49
@ WITHOUT_GHOSTS
Definition: iceModelVec.hh:49
std::string filename
Definition: options.hh:33