tPIK.cc - 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
---
tPIK.cc (10589B)
---
1 // Copyright (C) 2009-2018 Ricarda Winkelmann, Torsten Albrecht, Constantine Khrulev
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 2 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 // Implementation of the atmosphere model using constant-in-time precipitation
20 // and a cosine yearly cycle for near-surface air temperatures.
21
22 // This includes the PIK temperature parameterization.
23
24 #include "PIK.hh"
25
26 #include "pism/geometry/Geometry.hh"
27 #include "pism/util/ConfigInterface.hh"
28 #include "pism/util/IceGrid.hh"
29 #include "pism/util/MaxTimestep.hh"
30 #include "pism/util/error_handling.hh"
31
32 namespace pism {
33 namespace atmosphere {
34
35 PIK::PIK(IceGrid::ConstPtr g)
36 : YearlyCycle(g) {
37
38 auto parameterization = m_config->get_string("atmosphere.pik.parameterization");
39
40 std::map<std::string, Parameterization>
41 models = {{"martin", MARTIN},
42 {"huybrechts_dewolde", HUYBRECHTS_DEWOLDE},
43 {"martin_huybrechts_dewolde", MARTIN_HUYBRECHTS_DEWOLDE},
44 {"era_interim", ERA_INTERIM},
45 {"era_interim_sin", ERA_INTERIM_SIN},
46 {"era_interim_lon", ERA_INTERIM_LON}};
47
48 if (models.find(parameterization) == models.end()) {
49 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
50 "invalid pik parameterization: %s",
51 parameterization.c_str());
52 }
53
54 m_parameterization = models[parameterization];
55 }
56
57 PIK::~PIK() {
58 // empty
59 }
60
61 void PIK::init_impl(const Geometry &geometry) {
62
63 m_log->message(2,
64 "* Initializing PIK atmosphere model with air temperature parameterization based on \n"
65 " Huybrechts & De Wolde (1999) or multiple regression analysis of ERA INTERIM data...\n"
66 " Uses a time-independent precipitation field read from a file.");
67
68 m_reference = "Winkelmann et al.";
69
70 auto precip_file = m_config->get_string("atmosphere.pik.file");
71
72 if (not precip_file.empty()) {
73 YearlyCycle::init_internal(precip_file,
74 true, /* do regrid */
75 0 /* start (irrelevant) */);
76 } else {
77 // try to read precipitation from the input (-i) file
78 YearlyCycle::init_impl(geometry);
79 }
80
81 switch (m_parameterization) {
82 case HUYBRECHTS_DEWOLDE:
83 m_log->message(2,
84 " Parameterization based on: Huybrechts & De Wolde (1999).\n");
85 break;
86 case ERA_INTERIM:
87 m_log->message(2,
88 " Parameterization based on: multiple regression analysis of ERA INTERIM data.\n");
89 break;
90 case ERA_INTERIM_SIN:
91 m_log->message(2,
92 " Parameterization based on: multiple regression analysis of ERA INTERIM data with a sin(lat) dependence.\n");
93 break;
94 case ERA_INTERIM_LON:
95 m_log->message(2,
96 " Parameterization based on: multiple regression analysis of ERA INTERIM data with a cos(lon) dependence.\n");
97 break;
98 case MARTIN_HUYBRECHTS_DEWOLDE:
99 m_log->message(2,
100 " Mean annual temperature: as in Martin et al (2011).\n"
101 " Mean summer temperature: anomaly to the parameterization used by Huybrechts & De Wolde (1999).\n");
102 break;
103 default:
104 case MARTIN:
105 m_log->message(2,
106 " Mean annual temperature: as in Martin et al (2011).\n"
107 " No seasonal variation in air temperature.\n");
108 }
109 }
110
111 MaxTimestep PIK::max_timestep_impl(double t) const {
112 (void) t;
113 return MaxTimestep("atmosphere pik");
114 }
115
116 /*!
117 * See equation C1 in HuybrechtsdeWolde.
118 */
119 static double huybrechts_dewolde_mean_annual(double surface_elevation, double latitude) {
120 double gamma_a = surface_elevation < 1500.0 ? -0.005102 : -0.014285;
121 return 273.15 + 34.46 + gamma_a * surface_elevation - 0.68775 * latitude * (-1.0);
122 }
123
124 /*!
125 * See equation C2 in HuybrechtsdeWolde.
126 */
127 static double huybrechts_dewolde_mean_summer(double surface_elevation, double latitude) {
128 return 273.15 + 16.81 - 0.00692 * surface_elevation - 0.27937 * latitude * (-1.0);
129 }
130
131 /*!
132 * Parameterization of mean annual and mean summer near-surface temperature as in
133 * Huybrechts & DeWolde (1999)
134 */
135 static void huybrechts_dewolde(const Geometry &geometry, IceModelVec2S &T_ma, IceModelVec2S &T_ms) {
136 IceGrid::ConstPtr grid = T_ma.grid();
137
138 const IceModelVec2S
139 &h = geometry.ice_surface_elevation,
140 &lat = geometry.latitude;
141
142 IceModelVec::AccessList list{&h, &lat, &T_ma, &T_ms};
143
144 for (Points p(*grid); p; p.next()) {
145 const int i = p.i(), j = p.j();
146
147 T_ma(i, j) = huybrechts_dewolde_mean_annual(h(i, j), lat(i, j));
148 T_ms(i, j) = huybrechts_dewolde_mean_summer(h(i, j), lat(i, j));
149 }
150 }
151
152 /*!
153 * Parametrization based on multiple regression analysis of ERA INTERIM data
154 */
155 static void era_interim(const Geometry &geometry, IceModelVec2S &T_ma, IceModelVec2S &T_ms) {
156 IceGrid::ConstPtr grid = T_ma.grid();
157
158 const IceModelVec2S
159 &h = geometry.ice_surface_elevation,
160 &lat = geometry.latitude;
161
162 IceModelVec::AccessList list{&h, &lat, &T_ma, &T_ms};
163
164 for (Points p(*grid); p; p.next()) {
165 const int i = p.i(), j = p.j();
166
167 T_ma(i, j) = 273.15 + 29.2 - 0.0082 * h(i, j) - 0.576 * lat(i, j) * (-1.0);
168 T_ms(i, j) = 273.15 + 16.5 - 0.0068 * h(i, j) - 0.248 * lat(i, j) * (-1.0);
169 }
170 }
171
172 /*!
173 * Parametrization based on multiple regression analysis of ERA INTERIM data with sin(lat)
174 */
175 static void era_interim_sin(const Geometry &geometry, IceModelVec2S &T_ma, IceModelVec2S &T_ms) {
176 IceGrid::ConstPtr grid = T_ma.grid();
177
178 const IceModelVec2S
179 &h = geometry.ice_surface_elevation,
180 &lat = geometry.latitude;
181
182 IceModelVec::AccessList list{&h, &lat, &T_ma, &T_ms};
183
184 for (Points p(*grid); p; p.next()) {
185 const int i = p.i(), j = p.j();
186
187 T_ma(i, j) = 273.15 - 2.0 - 0.0082 * h(i, j) + 18.4 * (sin(3.1415 * lat(i, j) / 180.0) + 0.8910) / (1 - 0.8910);
188 T_ms(i, j) = 273.15 + 3.2 - 0.0067 * h(i, j) + 8.3 * (sin(3.1415 * lat(i, j) / 180.0) + 0.8910) / (1 - 0.8910);
189 }
190 }
191
192 /*!
193 * Parametrization based on multiple regression analysis of ERA INTERIM data with cos(lon)
194 */
195 static void era_interim_lon(const Geometry &geometry, IceModelVec2S &T_ma, IceModelVec2S &T_ms) {
196 IceGrid::ConstPtr grid = T_ma.grid();
197
198 const IceModelVec2S
199 &h = geometry.ice_surface_elevation,
200 &lat = geometry.latitude,
201 &lon = geometry.longitude;
202
203 IceModelVec::AccessList list{&h, &lat, &lon, &T_ma, &T_ms};
204
205 for (Points p(*grid); p; p.next()) {
206 const int i = p.i(), j = p.j();
207
208 double hmod = std::max(1000.0, h(i, j));
209 T_ma(i, j) = 273.15 + 37.49 - 0.0095 * hmod - 0.644 * lat(i, j) * (-1.0) + 2.146 * cos(3.1415 * (lon(i, j) + 110.0) / 180.0);
210 T_ms(i, j) = 273.15 + 15.74 - 0.0083 * hmod - 0.196 * lat(i, j) * (-1.0) + 0.225 * cos(3.1415 * (lon(i, j) + 110.0) / 180.0);
211
212 }
213 }
214
215 /*!
216 * Parameterization of the mean annual near-surface air temperature, see equation (1) in
217 * Martin et al, 2011.
218 */
219 static double martin2011_mean_annual(double elevation, double latitude) {
220 return 273.15 + 30 - 0.0075 * elevation - 0.68775 * latitude * (-1.0);
221 }
222
223 /*!
224 * - annual mean temperature as in Martin et al. (2011)
225 * - no seasonal variation of air temperature
226 */
227 static void martin2011(const Geometry &geometry, IceModelVec2S &T_ma, IceModelVec2S &T_ms) {
228 IceGrid::ConstPtr grid = T_ma.grid();
229
230 const IceModelVec2S
231 &h = geometry.ice_surface_elevation,
232 &lat = geometry.latitude;
233
234 IceModelVec::AccessList list{&h, &lat, &T_ma, &T_ms};
235
236 for (Points p(*grid); p; p.next()) {
237 const int i = p.i(), j = p.j();
238
239 T_ma(i, j) = martin2011_mean_annual(h(i, j), lat(i, j));
240 T_ms(i, j) = T_ma(i, j);
241 }
242 }
243
244 /*!
245 * - annual mean temperature as in Martin et al. (2011)
246 * - summer mean temperature computed as an anomaly to Huybrechts & DeWolde (1999)
247 */
248 static void martin_huybrechts_dewolde(const Geometry &geometry, IceModelVec2S &T_ma, IceModelVec2S &T_ms) {
249 IceGrid::ConstPtr grid = T_ma.grid();
250
251 const IceModelVec2S
252 &h = geometry.ice_surface_elevation,
253 &lat = geometry.latitude;
254
255 IceModelVec::AccessList list{&h, &lat, &T_ma, &T_ms};
256
257 for (Points p(*grid); p; p.next()) {
258 const int i = p.i(), j = p.j();
259
260 // mean annual surface temperature as in Martin et al. 2011, equation (1)
261 T_ma(i, j) = 273.15 + 30 - 0.0075 * h(i, j) - 0.68775 * lat(i, j) * (-1.0);
262
263 double
264 TMA = huybrechts_dewolde_mean_annual(h(i, j), lat(i, j)),
265 TMS = huybrechts_dewolde_mean_summer(h(i, j), lat(i, j));
266
267 T_ms(i, j) = T_ma(i, j) + (TMS - TMA);
268 }
269 }
270
271
272 /*!
273 * Updates mean annual and mean summer (January) near-surface air temperatures. Note that
274 * the precipitation rate is time-independent and does not need to be updated.
275 */
276 void PIK::update_impl(const Geometry &geometry, double t, double dt) {
277 (void) t;
278 (void) dt;
279
280 if (geometry.latitude.metadata().has_attribute("missing_at_bootstrap")) {
281 throw RuntimeError(PISM_ERROR_LOCATION,
282 "latitude variable was missing at bootstrap;\n"
283 "PIK atmosphere model depends on latitude and would return nonsense!");
284 }
285
286 switch (m_parameterization) {
287 case HUYBRECHTS_DEWOLDE:
288 huybrechts_dewolde(geometry, m_air_temp_mean_annual, m_air_temp_mean_summer);
289 break;
290 case ERA_INTERIM:
291 era_interim(geometry, m_air_temp_mean_annual, m_air_temp_mean_summer);
292 break;
293 case ERA_INTERIM_SIN:
294 era_interim_sin(geometry, m_air_temp_mean_annual, m_air_temp_mean_summer);
295 break;
296 case ERA_INTERIM_LON:
297 era_interim_lon(geometry, m_air_temp_mean_annual, m_air_temp_mean_summer);
298 break;
299 case MARTIN_HUYBRECHTS_DEWOLDE:
300 martin_huybrechts_dewolde(geometry, m_air_temp_mean_annual, m_air_temp_mean_summer);
301 break;
302 default:
303 case MARTIN:
304 martin2011(geometry, m_air_temp_mean_annual, m_air_temp_mean_summer);
305 break;
306 }
307 }
308
309 } // end of namespace atmosphere
310 } // end of namespace pism