PISM, A Parallel Ice Sheet Model  stable v2.0.5 committed by Constantine Khrulev on 2022-10-14 09:56:26 -0800
NoGLRetreat.cc
Go to the documentation of this file.
1 /* Copyright (C) 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 "NoGLRetreat.hh"
21 
22 #include "pism/geometry/Geometry.hh"
23 #include "pism/util/Diagnostic.hh"
24 #include "pism/util/pism_utilities.hh" // combine()
25 
26 namespace pism {
27 namespace surface {
28 
30  std::shared_ptr<SurfaceModel> input)
31  : SurfaceModel(grid, input),
32  m_smb_adjustment(grid, "smb_adjustment", WITHOUT_GHOSTS),
33  m_min_ice_thickness(grid, "minimum_ice_thickness", WITHOUT_GHOSTS) {
34 
35  m_smb_adjustment.metadata()["units"] = "kg m-2 s-1";
36 
41 }
42 
43 void NoGLRetreat::init_impl(const Geometry &geometry) {
44  m_input_model->init(geometry);
45 
46  m_log->message(2,
47  "* Initializing a SMB adjustment preventing grounding line retreat...\n");
48 
49  const auto &ice_thickness = geometry.ice_thickness;
50  const auto &sea_level = geometry.sea_level_elevation;
51  const auto &bed = geometry.bed_elevation;
52 
53  double rho_i = m_config->get_number("constants.ice.density");
54  double rho_w = m_config->get_number("constants.sea_water.density");
55  double eps = m_config->get_number("geometry.ice_free_thickness_standard");
56 
57  IceModelVec::AccessList list{&sea_level, &bed, &ice_thickness,
59 
60  for (Points p(*m_grid); p; p.next()) {
61  const int i = p.i(), j = p.j();
62 
63  double H_min = 0.0;
64  if (sea_level(i, j) > bed(i, j)) {
65  // compute the minimum grounded ice thickess
66  H_min = (sea_level(i, j) - bed(i, j)) * (rho_w / rho_i) + eps;
67 
68  // exclude areas where the ice thickness is below this threshold (i.e. the area is
69  // ice free or covered by floating ice)
70  if (ice_thickness(i, j) < H_min) {
71  H_min = 0.0;
72  }
73  }
74 
75  m_min_ice_thickness(i, j) = H_min;
76 
77  }
78 }
79 
80 void NoGLRetreat::update_impl(const Geometry &geometry, double t, double dt) {
81 
82  m_input_model->update(geometry, t, dt);
83 
84  const auto &mass_flux = m_input_model->mass_flux();
85  const auto &ice_thickness = geometry.ice_thickness;
86 
87  double rho_i = m_config->get_number("constants.ice.density");
88 
90  &ice_thickness,
91  m_mass_flux.get(),
94 
95  for (Points p(*m_grid); p; p.next()) {
96  const int i = p.i(), j = p.j();
97 
98  double SMB_old = mass_flux(i, j);
99  double SMB_new = SMB_old;
100 
101  // convert from kg/(m^2 s) to m:
102  double H_min = m_min_ice_thickness(i, j);
103  double H = ice_thickness(i, j);
104  double dH = mass_flux(i, j) * (dt / rho_i);
105  double H_new = H + dH;
106 
107  if (H_min > 0.0 and H_new < H_min) {
108  SMB_new = (H_min - H) * (rho_i / dt);
109  }
110 
111  (*m_mass_flux)(i, j) = SMB_new;
112  m_smb_adjustment(i, j) = SMB_new - SMB_old;
113  }
114 
118 }
119 
121  return *m_mass_flux;
122 }
123 
125  return *m_accumulation;
126 }
127 
129  return *m_melt;
130 }
131 
133  return *m_runoff;
134 }
135 
137  return m_smb_adjustment;
138 }
139 
140 namespace diagnostics {
141 
142 class SMBAdjustment : public DiagAverageRate<NoGLRetreat>
143 {
144 public:
147  "no_gl_retreat_smb_adjustment",
148  RATE)
149  {
150 
151  m_vars = {{m_sys, "no_gl_retreat_smb_adjustment"}};
152  m_accumulator.metadata()["units"] = "kg m-2";
153 
154  set_attrs("SMB adjustment needed to maintain grounded ice extent",
155  "", // no standard name
156  "kg m-2 s-1",
157  "kg m-2 year-1",
158  0);
159  m_vars[0]["cell_methods"] = "time: mean";
160 
161  m_vars[0]["_FillValue"] = {to_internal(m_fill_value)};
162  }
163 
164 protected:
166  return model->smb_adjustment();
167  }
168 };
169 
170 } // end of namespace diagnostics
171 
173  return combine({{"no_gl_retreat_smb_adjustment",
175  m_input_model->diagnostics());
176 }
177 
178 } // end of namespace surface
179 } // 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
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
const Model * model
Definition: Diagnostic.hh:166
double m_fill_value
fill value (used often enough to justify storing it)
Definition: Diagnostic.hh:114
const units::System::Ptr m_sys
the unit system
Definition: Diagnostic.hh:108
double to_internal(double x) const
Definition: Diagnostic.cc:58
std::vector< SpatialVariableMetadata > m_vars
metadata corresponding to NetCDF variables
Definition: Diagnostic.hh:112
std::shared_ptr< Diagnostic > Ptr
Definition: Diagnostic.hh:64
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
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
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition: iceModelVec.cc:533
DiagnosticList diagnostics_impl() const
Definition: NoGLRetreat.cc:172
const IceModelVec2S & runoff_impl() const
Definition: NoGLRetreat.cc:132
const IceModelVec2S & accumulation_impl() const
Definition: NoGLRetreat.cc:124
IceModelVec2S m_smb_adjustment
Definition: NoGLRetreat.hh:46
const IceModelVec2S & mass_flux_impl() const
Definition: NoGLRetreat.cc:120
NoGLRetreat(IceGrid::ConstPtr g, std::shared_ptr< SurfaceModel > input)
Definition: NoGLRetreat.cc:29
void update_impl(const Geometry &geometry, double t, double dt)
Definition: NoGLRetreat.cc:80
const IceModelVec2S & melt_impl() const
Definition: NoGLRetreat.cc:128
IceModelVec2S m_min_ice_thickness
Definition: NoGLRetreat.hh:47
const IceModelVec2S & smb_adjustment() const
Definition: NoGLRetreat.cc:136
void init_impl(const Geometry &geometry)
Definition: NoGLRetreat.cc:43
IceModelVec2S::Ptr m_mass_flux
Definition: NoGLRetreat.hh:45
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
std::shared_ptr< SurfaceModel > m_input_model
void dummy_accumulation(const IceModelVec2S &smb, IceModelVec2S &result)
const IceModelVec2S & mass_flux() const
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
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