PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
EISMINTII.cc
Go to the documentation of this file.
1 /* Copyright (C) 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023 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 <cmath> // std::pow(), std::sqrt()
20 #include <algorithm> // std::min
21 
22 #include "pism/coupler/surface/EISMINTII.hh"
23 #include "pism/util/ConfigInterface.hh"
24 #include "pism/util/pism_options.hh"
25 #include "pism/util/Grid.hh"
26 #include "pism/util/MaxTimestep.hh"
27 
28 namespace pism {
29 namespace surface {
30 
31 EISMINTII::EISMINTII(std::shared_ptr<const Grid> g, int experiment)
32  : PSFormulas(g), m_experiment(experiment) {
33  // empty
34 }
35 
36 void EISMINTII::init_impl(const Geometry &geometry) {
37  (void) geometry;
38 
39  using units::convert;
40 
41  m_log->message(2,
42  "setting parameters for surface mass balance"
43  " and temperature in EISMINT II experiment %c ... \n",
44  m_experiment);
45 
46  // EISMINT II specified values for parameters
47  m_S_b = convert(m_sys, 1.0e-2 * 1e-3, "year-1", "second-1"); // Grad of accum rate change
48  m_S_T = 1.67e-2 * 1e-3; // K/m Temp gradient
49 
50  // these are for A,E,G,H,I,K:
51  m_M_max = convert(m_sys, 0.5, "m year-1", "m second-1"); // Max accumulation
52  m_R_el = 450.0e3; // Distance to equil line (SMB=0)
53  m_T_min = 238.15;
54 
55  switch (m_experiment) {
56  case 'B': // supposed to start from end of experiment A and:
57  m_T_min = 243.15;
58  break;
59  case 'C':
60  case 'J':
61  case 'L': // supposed to start from end of experiment A (for C;
62  // resp I and K for J and L) and:
63  m_M_max = convert(m_sys, 0.25, "m year-1", "m second-1");
64  m_R_el = 425.0e3;
65  break;
66  case 'D': // supposed to start from end of experiment A and:
67  m_R_el = 425.0e3;
68  break;
69  case 'F': // start with zero ice and:
70  m_T_min = 223.15;
71  break;
72  }
73 
74  // if user specifies Tmin, Tmax, Mmax, Sb, ST, Rel, then use that (override above)
75  m_T_min = options::Real(m_sys, "-Tmin", "T min, Kelvin", "Kelvin", m_T_min);
76 
77  options::Real Mmax(m_sys, "-Mmax", "Maximum accumulation, m year-1",
78  "m year-1",
79  convert(m_sys, m_M_max, "m second-1", "m year-1"));
80  if (Mmax.is_set()) {
81  m_M_max = convert(m_sys, Mmax, "m year-1", "m second-1");
82  }
83 
84  options::Real Sb(m_sys, "-Sb", "radial gradient of accumulation rate, (m year-1)/km",
85  "m year-1/km",
86  convert(m_sys, m_S_b, "m second-1/m", "m year-1/km"));
87  if (Sb.is_set()) {
88  m_S_b = convert(m_sys, Sb, "m year-1/km", "m second-1/m");
89  }
90 
91  options::Real ST(m_sys, "-ST", "radial gradient of surface temperature, K/km",
92  "K/km",
93  convert(m_sys, m_S_T, "K/m", "K/km"));
94  if (ST.is_set()) {
95  m_S_T = convert(m_sys, ST, "K/km", "K/m");
96  }
97 
98  options::Real Rel(m_sys, "-Rel", "radial distance to equilibrium line, km",
99  "km",
100  convert(m_sys, m_R_el, "m", "km"));
101  if (Rel.is_set()) {
102  m_R_el = convert(m_sys, Rel, "km", "m");
103  }
104 
106 }
107 
109  (void) t;
110  return MaxTimestep("surface EISMINT-II");
111 }
112 
114  using std::pow;
115  using std::sqrt;
116 
117  // center of the accumulation and surface temperature patterns
118  double cx = 0.0, cy = 0.0;
119  if (m_experiment == 'E') {
120  cx += 100.0e3;
121  cy += 100.0e3;
122  }
123 
124  array::AccessScope list{m_temperature.get(), m_mass_flux.get()};
125 
126  for (auto p = m_grid->points(); p; p.next()) {
127  const int i = p.i(), j = p.j();
128 
129  const double r = sqrt(pow(m_grid->x(i) - cx, 2) + pow(m_grid->y(j) - cy, 2));
130 
131  // accumulation (formula (7) in [Payne et al 2000])
132  (*m_mass_flux)(i,j) = std::min(m_M_max, m_S_b * (m_R_el-r));
133 
134  // surface temperature (formula (8) in [Payne et al 2000])
135  (*m_temperature)(i,j) = m_T_min + m_S_T * r;
136  }
137 
138  // convert from "m second-1" to "kg m-2 s-1"
139  m_mass_flux->scale(m_config->get_number("constants.ice.density"));
140 }
141 
142 void EISMINTII::update_impl(const Geometry &geometry, double t, double dt) {
143  (void) t;
144  (void) dt;
145  (void) geometry;
146 
150 
151 }
152 
153 } // end of namespace surface
154 } // end of namespace pism
const units::System::Ptr m_sys
unit system used by this component
Definition: Component.hh:160
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
Combines the max. time step with the flag indicating if a restriction is active. Makes is possible to...
Definition: MaxTimestep.hh:31
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
Definition: Array.hh:65
bool is_set() const
Definition: options.hh:35
EISMINTII(std::shared_ptr< const Grid > g, int experiment)
Definition: EISMINTII.cc:31
void update_impl(const Geometry &geometry, double t, double dt)
Definition: EISMINTII.cc:142
virtual MaxTimestep max_timestep_impl(double t) const
Definition: EISMINTII.cc:108
void init_impl(const Geometry &geometry)
Definition: EISMINTII.cc:36
std::shared_ptr< array::Scalar > m_temperature
Definition: Formulas.hh:50
std::shared_ptr< array::Scalar > m_mass_flux
Definition: Formulas.hh:49
void dummy_accumulation(const array::Scalar &smb, array::Scalar &result)
std::shared_ptr< array::Scalar > m_melt
std::shared_ptr< array::Scalar > m_runoff
void dummy_melt(const array::Scalar &smb, array::Scalar &result)
std::shared_ptr< array::Scalar > m_accumulation
void dummy_runoff(const array::Scalar &smb, array::Scalar &result)
double min(const array::Scalar &input)
Finds minimum over all the values in an array::Scalar object. Ignores ghosts.
Definition: Scalar.cc:193
static const double ST
double convert(System::Ptr system, double input, const std::string &spec1, const std::string &spec2)
Convert a quantity from unit1 to unit2.
Definition: Units.cc:70
static const double g
Definition: exactTestP.cc:36