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