PISM, A Parallel Ice Sheet Model  stable v2.0.4 committed by Constantine Khrulev on 2022-05-25 12:02:27 -0800
OceanModel.cc
Go to the documentation of this file.
1 /* Copyright (C) 2013, 2014, 2015, 2016, 2017, 2018, 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 
20 #include <gsl/gsl_math.h> // GSL_NAN
21 
22 #include "pism/coupler/OceanModel.hh"
23 #include "pism/util/iceModelVec.hh"
24 #include "pism/util/MaxTimestep.hh"
25 #include "pism/util/pism_utilities.hh"
26 #include "pism/geometry/Geometry.hh"
27 
28 namespace pism {
29 
30 namespace ocean {
31 
33  IceModelVec2S::Ptr result(new IceModelVec2S(g, "shelfbtemp", WITHOUT_GHOSTS));
34  result->set_attrs("diagnostic",
35  "ice temperature at the bottom of floating ice",
36  "Kelvin", "Kelvin", "", 0);
37  return result;
38 }
39 
41  IceModelVec2S::Ptr result(new IceModelVec2S(g, "shelfbmassflux", WITHOUT_GHOSTS));
42 
43  result->set_attrs("diagnostic", "shelf base mass flux",
44  "kg m-2 s-1", "kg m-2 year-1", "", 0);
45 
46  return result;
47 }
48 
51  "average_water_column_pressure",
53  result->set_attrs("diagnostic",
54  "vertically-averaged water column pressure",
55  "Pa", "Pa", "", 0);
56 
57  return result;
58 }
59 
60 // "modifier" constructor
61 OceanModel::OceanModel(IceGrid::ConstPtr g, std::shared_ptr<OceanModel> input)
62  : Component(g), m_input_model(input) {
63 
64  if (not input) {
66  }
67 }
68 
69 // "model" constructor
71  : OceanModel(g, std::shared_ptr<OceanModel>()) {
72  // empty
73 }
74 
75 void OceanModel::init(const Geometry &geometry) {
76  this->init_impl(geometry);
77 }
78 
79 void OceanModel::init_impl(const Geometry &geometry) {
80  if (m_input_model) {
81  m_input_model->init(geometry);
82  } else {
83  double
84  ice_density = m_config->get_number("constants.ice.density"),
85  water_density = m_config->get_number("constants.sea_water.density"),
86  g = m_config->get_number("constants.standard_gravity");
87 
88  compute_average_water_column_pressure(geometry, ice_density, water_density, g,
90  }
91 }
92 
93 void OceanModel::update(const Geometry &geometry, double t, double dt) {
94  this->update_impl(geometry, t, dt);
95 }
96 
97 
100 }
101 
104 }
105 
108 }
109 
110 // pass-through default implementations for "modifiers"
111 
112 void OceanModel::update_impl(const Geometry &geometry, double t, double dt) {
113  if (m_input_model) {
114  m_input_model->update(geometry, t, dt);
115  } else {
116  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "no input model");
117  }
118 }
119 
121  if (m_input_model) {
122  return m_input_model->max_timestep(t);
123  } else {
124  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "no input model");
125  }
126 }
127 
128 void OceanModel::define_model_state_impl(const File &output) const {
129  if (m_input_model) {
130  return m_input_model->define_model_state(output);
131  } else {
132  // no state to define
133  }
134 }
135 
136 void OceanModel::write_model_state_impl(const File &output) const {
137  if (m_input_model) {
138  return m_input_model->write_model_state(output);
139  } else {
140  // no state to write
141  }
142 }
143 
145  if (m_input_model) {
146  return m_input_model->shelf_base_temperature();
147  } else {
148  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "no input model");
149  }
150 }
151 
153  if (m_input_model) {
154  return m_input_model->shelf_base_mass_flux();
155  } else {
156  throw RuntimeError::formatted(PISM_ERROR_LOCATION, "no input model");
157  }
158 }
159 
161  if (m_input_model) {
162  return m_input_model->average_water_column_pressure();
163  } else {
164  return *m_water_column_pressure;
165  }
166 }
167 
168 namespace diagnostics {
169 
170 /*! @brief Shelf base temperature. */
171 class PO_shelf_base_temperature : public Diag<OceanModel>
172 {
173 public:
175  : Diag<OceanModel>(m) {
176 
177  /* set metadata: */
178  m_vars = {SpatialVariableMetadata(m_sys, "shelfbtemp")};
179 
180  set_attrs("ice temperature at the basal surface of ice shelves", "",
181  "Kelvin", "Kelvin", 0);
182  }
183 protected:
185 
186  IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "shelfbtemp", WITHOUT_GHOSTS));
187  result->metadata(0) = m_vars[0];
188 
189  result->copy_from(model->shelf_base_temperature());
190 
191  return result;
192  }
193 };
194 
195 
196 
197 /*! @brief Shelf base mass flux. */
198 class PO_shelf_base_mass_flux : public Diag<OceanModel>
199 {
200 public:
202  : Diag<OceanModel>(m) {
203 
204  /* set metadata: */
205  m_vars = {SpatialVariableMetadata(m_sys, "shelfbmassflux")};
206 
207  set_attrs("mass flux at the basal surface of ice shelves", "",
208  "kg m-2 s-1", "kg m-2 s-1", 0);
209  }
210 protected:
212 
213  IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "shelfbmassflux", WITHOUT_GHOSTS));
214  result->metadata(0) = m_vars[0];
215 
216  result->copy_from(model->shelf_base_mass_flux());
217 
218  return result;
219  }
220 };
221 
222 } // end of namespace diagnostics
223 
225  using namespace diagnostics;
226  DiagnosticList result = {
227  {"shelfbtemp", Diagnostic::Ptr(new PO_shelf_base_temperature(this))},
228  {"shelfbmassflux", Diagnostic::Ptr(new PO_shelf_base_mass_flux(this))}
229  };
230 
231  if (m_input_model) {
232  return combine(m_input_model->diagnostics(), result);
233  } else {
234  return result;
235  }
236 }
237 
239  if (m_input_model) {
240  return m_input_model->ts_diagnostics();
241  } else {
242  return {};
243  }
244 }
245 
247  double ice_density,
248  double water_density,
249  double g,
250  IceModelVec2S &result) {
251 
252  auto grid = result.grid();
253 
254  const IceModelVec2S
255  &bed = geometry.bed_elevation,
256  &H = geometry.ice_thickness,
257  &z_s = geometry.sea_level_elevation;
258 
259  IceModelVec::AccessList l{&bed, &H, &z_s, &result};
260 
261  ParallelSection loop(grid->com);
262  try {
263  for (Points p(*grid); p; p.next()) {
264  const int i = p.i(), j = p.j();
265 
266  result(i, j) = pism::average_water_column_pressure(H(i, j), bed(i, j), z_s(i, j),
267  ice_density, water_density, g);
268  }
269  } catch (...) {
270  loop.failed();
271  }
272  loop.check();
273 
274 
275 }
276 
277 
278 } // end of namespace ocean
279 } // 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 Config::ConstPtr m_config
configuration database used by this component
Definition: Component.hh:138
DiagnosticList diagnostics() const
Definition: Component.cc:89
A class defining a common interface for most PISM sub-models.
Definition: Component.hh:101
const OceanModel * model
Definition: Diagnostic.hh:166
A template derived from Diagnostic, adding a "Model".
Definition: Diagnostic.hh:161
const units::System::Ptr m_sys
the unit system
Definition: Diagnostic.hh:108
std::vector< SpatialVariableMetadata > m_vars
metadata corresponding to NetCDF variables
Definition: Diagnostic.hh:112
std::shared_ptr< Diagnostic > Ptr
Definition: Diagnostic.hh:64
IceGrid::ConstPtr m_grid
the grid
Definition: Diagnostic.hh:106
void set_attrs(const std::string &long_name, const std::string &standard_name, const std::string &units, const std::string &glaciological_units, unsigned int N=0)
A method for setting common variable attributes.
Definition: Diagnostic.cc:132
High-level PISM I/O class.
Definition: File.hh:51
IceModelVec2S bed_elevation
Definition: Geometry.hh:49
IceModelVec2S sea_level_elevation
Definition: Geometry.hh:50
IceModelVec2S ice_thickness
Definition: Geometry.hh:53
std::shared_ptr< const IceGrid > ConstPtr
Definition: IceGrid.hh:233
std::shared_ptr< IceModelVec2S > Ptr
Definition: iceModelVec.hh:341
std::shared_ptr< IceModelVec > Ptr
Definition: iceModelVec.hh:206
IceGrid::ConstPtr grid() const
Definition: iceModelVec.cc:128
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.
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
Spatial NetCDF variable (corresponding to a 2D or 3D scalar field).
std::shared_ptr< OceanModel > m_input_model
Definition: OceanModel.hh:69
static IceModelVec2S::Ptr allocate_shelf_base_mass_flux(IceGrid::ConstPtr g)
Definition: OceanModel.cc:40
virtual MaxTimestep max_timestep_impl(double t) const
Definition: OceanModel.cc:120
virtual void define_model_state_impl(const File &output) const
The default (empty implementation).
Definition: OceanModel.cc:128
IceModelVec2S::Ptr m_water_column_pressure
Definition: OceanModel.hh:70
virtual const IceModelVec2S & shelf_base_temperature_impl() const
Definition: OceanModel.cc:144
const IceModelVec2S & shelf_base_temperature() const
Definition: OceanModel.cc:102
virtual const IceModelVec2S & shelf_base_mass_flux_impl() const
Definition: OceanModel.cc:152
const IceModelVec2S & shelf_base_mass_flux() const
Definition: OceanModel.cc:98
virtual const IceModelVec2S & average_water_column_pressure_impl() const
Definition: OceanModel.cc:160
OceanModel(IceGrid::ConstPtr g, std::shared_ptr< OceanModel > input)
Definition: OceanModel.cc:61
virtual DiagnosticList diagnostics_impl() const
Definition: OceanModel.cc:224
virtual TSDiagnosticList ts_diagnostics_impl() const
Definition: OceanModel.cc:238
void update(const Geometry &geometry, double t, double dt)
Definition: OceanModel.cc:93
static void compute_average_water_column_pressure(const Geometry &geometry, double ice_density, double water_density, double g, IceModelVec2S &result)
Definition: OceanModel.cc:246
virtual void init_impl(const Geometry &geometry)
Definition: OceanModel.cc:79
void init(const Geometry &geometry)
Definition: OceanModel.cc:75
virtual void write_model_state_impl(const File &output) const
The default (empty implementation).
Definition: OceanModel.cc:136
virtual void update_impl(const Geometry &geometry, double t, double dt)
Definition: OceanModel.cc:112
const IceModelVec2S & average_water_column_pressure() const
Definition: OceanModel.cc:106
static IceModelVec2S::Ptr allocate_water_column_pressure(IceGrid::ConstPtr g)
Definition: OceanModel.cc:49
static IceModelVec2S::Ptr allocate_shelf_base_temperature(IceGrid::ConstPtr g)
Definition: OceanModel.cc:32
A very rudimentary PISM ocean model.
Definition: OceanModel.hh:36
#define PISM_ERROR_LOCATION
bool ocean(int M)
An ocean cell (floating ice or ice-free).
Definition: Mask.hh:39
double average_water_column_pressure(double ice_thickness, double bed, double floatation_level, double rho_ice, double rho_water, double g)
static const double g
Definition: exactTestP.cc:39
std::map< std::string, TSDiagnostic::Ptr > TSDiagnosticList
Definition: Diagnostic.hh:346
std::map< std::string, Diagnostic::Ptr > DiagnosticList
Definition: Diagnostic.hh:117
T combine(const T &a, const T &b)
@ WITHOUT_GHOSTS
Definition: iceModelVec.hh:49