PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
DischargeRouting.cc
Go to the documentation of this file.
1 // Copyright (C) 2018, 2019, 2021, 2022, 2023 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 "pism/coupler/frontalmelt/DischargeRouting.hh"
20 
21 #include "pism/util/Grid.hh"
22 #include "pism/geometry/Geometry.hh"
23 #include "pism/coupler/util/options.hh"
24 #include "pism/coupler/frontalmelt/FrontalMeltPhysics.hh"
25 #include "pism/util/array/Forcing.hh"
26 
27 namespace pism {
28 namespace frontalmelt {
29 
30 DischargeRouting::DischargeRouting(std::shared_ptr<const Grid> grid)
31  : FrontalMelt(grid, nullptr),
32  m_frontal_melt_rate(grid, "frontal_melt_rate") {
33 
35  .long_name("frontal melt rate")
36  .units("m s-1")
37  .output_units("m day-1");
38 
39  m_log->message(2, "* Initializing the frontal melt model\n"
40  " using the Rignot/Xu parameterization\n"
41  " and routing of subglacial discharge\n");
42 
43  m_theta_ocean = array::Forcing::Constant(grid, "theta_ocean", 0.0);
44 }
45 
46 void DischargeRouting::init_impl(const Geometry &geometry) {
47  (void)geometry;
48 
49  ForcingOptions opt(*m_grid->ctx(), "frontal_melt.routing");
50 
51  {
52  unsigned int buffer_size = m_config->get_number("input.forcing.buffer_size");
53 
55 
56  m_theta_ocean = std::make_shared<array::Forcing>(m_grid, file, "theta_ocean",
57  "", // no standard name
58  buffer_size, opt.periodic, LINEAR);
59  }
60 
61  m_theta_ocean->metadata(0)
62  .long_name("potential temperature of the adjacent ocean")
63  .units("Celsius");
64 
65  m_theta_ocean->init(opt.filename, opt.periodic);
66 }
67 
68 /*!
69  * Initialize potential temperature from an array instead of an input
70  * file (for testing).
71  */
73  m_theta_ocean->copy_from(theta);
74 }
75 
76 void DischargeRouting::update_impl(const FrontalMeltInputs &inputs, double t, double dt) {
77 
78  m_theta_ocean->update(t, dt);
79 
80  FrontalMeltPhysics physics(*m_config);
81 
82  const auto &cell_type = inputs.geometry->cell_type;
83  const array::Scalar &bed_elevation = inputs.geometry->bed_elevation;
84  const array::Scalar &ice_thickness = inputs.geometry->ice_thickness;
85  const array::Scalar &sea_level_elevation = inputs.geometry->sea_level_elevation;
86  const array::Scalar &water_flux = *inputs.subglacial_water_flux;
87 
88  array::AccessScope list{ &ice_thickness, &bed_elevation, &cell_type,
89  &sea_level_elevation, &water_flux, m_theta_ocean.get(),
91 
92  double seconds_per_day = 86400, grid_spacing = 0.5 * (m_grid->dx() + m_grid->dy());
93 
94  for (auto p = m_grid->points(); p; p.next()) {
95  const int i = p.i(), j = p.j();
96 
97  if (cell_type.icy(i, j)) {
98  // Assume for now that thermal forcing is equal to theta_ocean. Also, thermal
99  // forcing is generally not available at the grounding line.
100  double TF = (*m_theta_ocean)(i, j);
101 
102  double water_depth = std::max(sea_level_elevation(i, j) - bed_elevation(i, j), 0.0),
103  submerged_front_area = water_depth * grid_spacing;
104 
105  // Convert subglacial water flux (m^2/s) to an "effective subglacial freshwater
106  // velocity" or flux per unit area of ice front in m/day (see Xu et al 2013, section
107  // 2, paragraph 11).
108  //
109  // [flux] = m^2 / s, so
110  // [flux * grid_spacing] = m^3 / s, so
111  // [flux * grid_spacing / submerged_front_area] = m / s, and
112  // [flux * grid_spacing * (s / day) / submerged_front_area] = m / day
113  double Q_sg = water_flux(i, j) * grid_spacing;
114  double q_sg = Q_sg / submerged_front_area * seconds_per_day;
115 
116  m_frontal_melt_rate(i, j) = physics.frontal_melt_from_undercutting(water_depth, q_sg, TF);
117  // convert from m / day to m / s
118  m_frontal_melt_rate(i, j) /= seconds_per_day;
119  } else {
120  m_frontal_melt_rate(i, j) = 0.0;
121  }
122  } // end of the loop over grid points
123 
124  // Set frontal melt rate *near* grounded termini to the average of grounded icy
125  // neighbors: front retreat code uses values at these locations (the rest is for
126  // visualization).
127 
129 
130  for (auto p = m_grid->points(); p; p.next()) {
131  const int i = p.i(), j = p.j();
132 
133  if (apply(cell_type, i, j) and cell_type.ice_free(i, j)) {
134 
135  auto R = m_frontal_melt_rate.star(i, j);
136  auto M = cell_type.star_int(i, j);
137 
138  int N = 0;
139  double R_sum = 0.0;
140  for (auto d : { North, East, South, West }) {
141  if (mask::grounded_ice(M[d]) or (m_include_floating_ice and mask::icy(M[d]))) {
142  R_sum += R[d];
143  N++;
144  }
145  }
146 
147  if (N > 0) {
148  m_frontal_melt_rate(i, j) = R_sum / N;
149  }
150  }
151  }
152 }
153 
155  return m_frontal_melt_rate;
156 }
157 
159 
160  auto dt = m_theta_ocean->max_timestep(t);
161 
162  if (dt.finite()) {
163  return {dt.value(), "frontal_melt routing"};
164  }
165 
166  return {"frontal_melt routing"};
167 }
168 
169 } // end of namespace frontalmelt
170 } // end of namespace pism
std::shared_ptr< const Grid > grid() const
Definition: Component.cc:105
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
const Geometry * geometry
Definition: FrontalMelt.hh:32
const array::Scalar * subglacial_water_flux
Definition: FrontalMelt.hh:35
array::Scalar1 sea_level_elevation
Definition: Geometry.hh:48
array::CellType2 cell_type
Definition: Geometry.hh:55
array::Scalar2 ice_thickness
Definition: Geometry.hh:51
array::Scalar2 bed_elevation
Definition: Geometry.hh:47
Combines the max. time step with the flag indicating if a restriction is active. Makes is possible to...
Definition: MaxTimestep.hh:31
VariableMetadata & long_name(const std::string &input)
VariableMetadata & output_units(const std::string &input)
VariableMetadata & units(const std::string &input)
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
Definition: Array.hh:65
stencils::Star< T > star(int i, int j) const
Definition: Array2D.hh:79
void update_ghosts()
Updates ghost points.
Definition: Array.cc:693
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition: Array.cc:553
static std::shared_ptr< Forcing > Constant(std::shared_ptr< const Grid > grid, const std::string &short_name, double value)
Definition: Forcing.cc:147
void update_impl(const FrontalMeltInputs &inputs, double t, double dt)
MaxTimestep max_timestep_impl(double t) const
DischargeRouting(std::shared_ptr< const Grid > g)
const array::Scalar & frontal_melt_rate_impl() const
void init_impl(const Geometry &geometry)
void initialize(const array::Scalar &theta)
std::shared_ptr< array::Forcing > m_theta_ocean
double frontal_melt_from_undercutting(double ice_thickness, double discharge_flux, double potential_temperature) const
bool apply(const array::CellType1 &M, int i, int j) const
Definition: FrontalMelt.cc:246
A very rudimentary PISM frontal melt model.
Definition: FrontalMelt.hh:43
double max(const array::Scalar &input)
Finds maximum over all the values in an array::Scalar object. Ignores ghosts.
Definition: Scalar.cc:165
@ PISM_NETCDF3
Definition: IO_Flags.hh:57
@ PISM_READONLY
open an existing file for reading only
Definition: IO_Flags.hh:72
bool icy(int M)
Ice-filled cell (grounded or floating).
Definition: Mask.hh:48
bool grounded_ice(int M)
Definition: Mask.hh:51
@ North
Definition: stencils.hh:24
@ East
Definition: stencils.hh:24
@ South
Definition: stencils.hh:24
@ West
Definition: stencils.hh:24
std::string filename
Definition: options.hh:33