tlocalMassBalance.hh - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tlocalMassBalance.hh (8796B)
---
1 // Copyright (C) 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2017, 2018 Ed Bueler 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 #ifndef __localMassBalance_hh
20 #define __localMassBalance_hh
21
22
23 #include <gsl/gsl_rng.h>
24
25 #include "pism/util/iceModelVec.hh" // only needed for FaustoGrevePDDObject
26
27 namespace pism {
28 namespace surface {
29
30 //! \brief Base class for a model which computes surface mass flux rate (ice
31 //! thickness per time) from precipitation and temperature.
32 /*!
33 This is a process model. At each spatial location, it uses a 1D array, with a
34 time dimension, for the temperature used in melting snow or ice. At each spatial
35 location it assumes the precipitation is time-independent.
36
37 This process model does not know its location on the ice sheet, but
38 simply computes the surface mass balance from three quantities:
39 - the time interval \f$[t,t+\Delta t]\f$,
40 - time series of values of surface temperature at \f$N\f$ equally-spaced
41 times in the time interval
42 - a scalar precipation rate which is taken to apply in the time interval.
43
44 This model also uses degree day factors passed-in in DegreeDayFactors `ddf`,
45 and the standard deviation `pddStdDev`. The latter is the standard deviation of the
46 modeled temperature away from the input temperature time series which contains
47 the part of location-dependent temperature cycle on the time interval.
48
49 @note
50 - Please avoid using `config.get...("...")` calls
51 inside those methods of this class which are called inside loops over
52 spatial grids. Doing otherwise increases computational costs.
53 - This base class should be more general. For instance, it could allow as
54 input a time series for precipation rate.
55 */
56 class LocalMassBalance {
57 public:
58
59 //! A struct which holds degree day factors.
60 /*!
61 Degree day factors convert positive degree days (=PDDs) into amount of melt.
62 */
63 struct DegreeDayFactors {
64 //! m day^-1 K^-1; ice-equivalent amount of snow melted, per PDD
65 double snow;
66 //! m day^-1 K^-1; ice-equivalent amount of ice melted, per PDD
67 double ice;
68 //! fraction of melted snow which refreezes as ice
69 double refreeze_fraction;
70 };
71
72 LocalMassBalance(Config::ConstPtr config, units::System::Ptr system);
73 virtual ~LocalMassBalance();
74
75 std::string method() const;
76
77 virtual unsigned int get_timeseries_length(double dt) = 0;
78
79 //! Count positive degree days (PDDs). Returned value in units of K day.
80 /*! Inputs T[0],...,T[N-1] are temperatures (K) at times t, t+dt_series, ..., t+(N-1)dt_series.
81 Inputs `t`, `dt_series` are in seconds. */
82 virtual void get_PDDs(double dt_series,
83 const std::vector<double> &S,
84 const std::vector<double> &T,
85 std::vector<double> &PDDs) = 0;
86
87 /*! Remove rain from precipitation. */
88 virtual void get_snow_accumulation(const std::vector<double> &T,
89 std::vector<double> &precip_rate) = 0;
90
91 class Changes {
92 public:
93 Changes();
94
95 double firn_depth;
96 double snow_depth;
97 double melt;
98 double runoff;
99 double smb;
100 };
101
102 /**
103 * Take a step of the PDD model.
104 *
105 * @param[in] ddf degree day factors
106 * @param[in] PDDs number of positive degree days during the time step [K day]
107 * @param[in] old_firn_depth firn depth [ice equivalent meters]
108 * @param[in] old_snow_depth snow depth [ice equivalent meters]
109 * @param[in] accumulation total solid (snow) accumulation during the time-step [ice equivalent meters]
110 */
111 virtual Changes step(const DegreeDayFactors &ddf,
112 double PDDs,
113 double ice_thickness,
114 double old_firn_depth,
115 double old_snow_depth,
116 double accumulation) = 0;
117
118 protected:
119 std::string m_method;
120
121 const Config::ConstPtr m_config;
122 const units::System::Ptr m_unit_system;
123 const double m_seconds_per_day;
124 };
125
126
127 //! A PDD implementation which computes the local mass balance based on an expectation integral.
128 /*!
129 The expected number of positive degree days is computed by an integral in \ref CalovGreve05.
130
131 Uses degree day factors which are location-independent.
132 */
133 class PDDMassBalance : public LocalMassBalance {
134
135 public:
136 PDDMassBalance(Config::ConstPtr config, units::System::Ptr system);
137 virtual ~PDDMassBalance() {}
138
139 virtual unsigned int get_timeseries_length(double dt);
140 virtual void get_PDDs(double dt_series,
141 const std::vector<double> &S,
142 const std::vector<double> &T,
143 std::vector<double> &PDDs);
144
145 void get_snow_accumulation(const std::vector<double> &T,
146 std::vector<double> &precip_rate);
147
148 Changes step(const DegreeDayFactors &ddf,
149 double PDDs,
150 double ice_thickness,
151 double firn_depth,
152 double snow_depth,
153 double accumulation);
154
155 protected:
156 double CalovGreveIntegrand(double sigma, double TacC);
157
158 bool precip_as_snow, //!< interpret all the precipitation as snow (no rain)
159 refreeze_ice_melt; //!< refreeze melted ice
160 double Tmin, //!< the temperature below which all precipitation is snow
161 Tmax; //!< the temperature above which all precipitation is rain
162 double pdd_threshold_temp; //!< threshold temperature for the PDD computation
163 };
164
165
166 //! An alternative PDD implementation which simulates a random process to get the number of PDDs.
167 /*!
168 Uses a GSL random number generator. Significantly slower because new random numbers are
169 generated for each grid point.
170
171 The way the number of positive degree-days are used to produce a surface mass balance
172 is identical to the base class PDDMassBalance.
173
174 \note
175 A more realistic pattern for the variability of surface melting might have correlation
176 with appropriate spatial and temporal ranges.
177 */
178 class PDDrandMassBalance : public PDDMassBalance {
179 public:
180
181 enum Kind {NOT_REPEATABLE = 0, REPEATABLE = 1};
182
183 PDDrandMassBalance(Config::ConstPtr config,
184 units::System::Ptr system,
185 Kind repeatable);
186 virtual ~PDDrandMassBalance();
187
188 virtual unsigned int get_timeseries_length(double dt);
189
190 virtual void get_PDDs(double dt_series,
191 const std::vector<double> &S,
192 const std::vector<double> &T,
193 std::vector<double> &PDDs);
194 protected:
195 gsl_rng *pddRandGen;
196 };
197
198
199 /*!
200 The PDD scheme described by Formula (6) in [\ref Faustoetal2009] requires
201 special knowledge of latitude and mean July temp to set degree day factors
202 for Greenland.
203
204 These formulas are inherited by [\ref Faustoetal2009] from [\ref Greve2005geothermal].
205 There was, apparently, tuning in [\ref Greve2005geothermal] which mixed ice
206 dynamical ideas and surface process ideas. That is, these formulas and parameter
207 choices arise from looking at margin shape. This may not be a good source of
208 PDD parameters.
209
210 This may become a derived class of a LocationDependentPDDObject, if the idea
211 is needed more in the future.
212 */
213 class FaustoGrevePDDObject {
214
215 public:
216 FaustoGrevePDDObject(IceGrid::ConstPtr g);
217 virtual ~FaustoGrevePDDObject();
218
219 void update_temp_mj(const IceModelVec2S &surfelev,
220 const IceModelVec2S &lat,
221 const IceModelVec2S &lon);
222
223 /*! If this method is called, it is assumed that i,j is in the ownership range
224 for IceModelVec2S temp_mj. */
225 LocalMassBalance::DegreeDayFactors degree_day_factors(int i, int j, double latitude);
226
227 protected:
228 IceGrid::ConstPtr m_grid;
229 const Config::ConstPtr m_config;
230
231 double m_beta_ice_w;
232 double m_beta_snow_w;
233 double m_T_c;
234 double m_T_w;
235 double m_beta_ice_c;
236 double m_beta_snow_c;
237 double m_fresh_water_density;
238 double m_ice_density;
239 double m_pdd_fausto_latitude_beta_w;
240 double m_refreeze_fraction;
241
242 IceModelVec2S m_temp_mj;
243 };
244
245
246 } // end of namespace surface
247 } // end of namespace pism
248
249 #endif