tAtmosphereModel.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
---
tAtmosphereModel.cc (8850B)
---
1 /* Copyright (C) 2016, 2017, 2018, 2019 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 <gsl/gsl_math.h> // GSL_NAN
21
22 #include "pism/coupler/AtmosphereModel.hh"
23 #include "pism/util/Time.hh"
24 #include "pism/util/error_handling.hh"
25 #include "pism/util/MaxTimestep.hh"
26
27 namespace pism {
28 namespace atmosphere {
29
30 IceModelVec2S::Ptr AtmosphereModel::allocate_temperature(IceGrid::ConstPtr grid) {
31 IceModelVec2S::Ptr result(new IceModelVec2S(grid, "air_temp", WITHOUT_GHOSTS));
32
33 result->set_attrs("climate_forcing", "mean annual near-surface air temperature",
34 "Kelvin", "Kelvin", "", 0);
35
36 return result;
37 }
38
39 IceModelVec2S::Ptr AtmosphereModel::allocate_precipitation(IceGrid::ConstPtr grid) {
40 IceModelVec2S::Ptr result(new IceModelVec2S(grid, "precipitation", WITHOUT_GHOSTS));
41 result->set_attrs("climate_forcing", "precipitation rate",
42 "kg m-2 second-1", "kg m-2 year-1",
43 "precipitation_flux", 0);
44
45 return result;
46 }
47
48 AtmosphereModel::AtmosphereModel(IceGrid::ConstPtr g)
49 : Component(g) {
50 // empty
51 }
52
53 AtmosphereModel::AtmosphereModel(IceGrid::ConstPtr g,
54 std::shared_ptr<AtmosphereModel> input)
55 :Component(g), m_input_model(input) {
56 // empty
57 }
58
59 AtmosphereModel::~AtmosphereModel() {
60 // empty
61 }
62
63 void AtmosphereModel::init(const Geometry &geometry) {
64 this->init_impl(geometry);
65 }
66
67 void AtmosphereModel::update(const Geometry &geometry, double t, double dt) {
68 this->update_impl(geometry, t, dt);
69 }
70
71 const IceModelVec2S& AtmosphereModel::mean_precipitation() const {
72 return this->mean_precipitation_impl();
73 }
74
75 const IceModelVec2S& AtmosphereModel::mean_annual_temp() const {
76 return this->mean_annual_temp_impl();
77 }
78
79 void AtmosphereModel::begin_pointwise_access() const {
80 this->begin_pointwise_access_impl();
81 }
82
83 void AtmosphereModel::end_pointwise_access() const {
84 this->end_pointwise_access_impl();
85 }
86
87 void AtmosphereModel::init_timeseries(const std::vector<double> &ts) const {
88 this->init_timeseries_impl(ts);
89 }
90
91 void AtmosphereModel::precip_time_series(int i, int j, std::vector<double> &result) const {
92 result.resize(m_ts_times.size());
93 this->precip_time_series_impl(i, j, result);
94 }
95
96 void AtmosphereModel::temp_time_series(int i, int j, std::vector<double> &result) const {
97 result.resize(m_ts_times.size());
98 this->temp_time_series_impl(i, j, result);
99 }
100
101 namespace diagnostics {
102
103 /*! @brief Instantaneous near-surface air temperature. */
104 class AirTemperatureSnapshot : public Diag<AtmosphereModel> {
105 public:
106 AirTemperatureSnapshot(const AtmosphereModel *m)
107 : Diag<AtmosphereModel>(m) {
108
109 /* set metadata: */
110 m_vars = {SpatialVariableMetadata(m_sys, "air_temp_snapshot")};
111
112 set_attrs("instantaneous value of the near-surface air temperature",
113 "", // no standard name
114 "Kelvin", "Kelvin", 0);
115 }
116 protected:
117 IceModelVec::Ptr compute_impl() const {
118
119 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "air_temp_snapshot", WITHOUT_GHOSTS));
120 result->metadata(0) = m_vars[0];
121
122 std::vector<double> current_time(1, m_grid->ctx()->time()->current());
123 std::vector<double> temperature(1, 0.0);
124
125 model->init_timeseries(current_time);
126
127 model->begin_pointwise_access();
128
129 IceModelVec::AccessList list(*result);
130 ParallelSection loop(m_grid->com);
131 try {
132 for (Points p(*m_grid); p; p.next()) {
133 const int i = p.i(), j = p.j();
134
135 model->temp_time_series(i, j, temperature);
136
137 (*result)(i, j) = temperature[0];
138 }
139 } catch (...) {
140 loop.failed();
141 }
142 loop.check();
143
144 model->end_pointwise_access();
145
146 return result;
147 }
148 };
149
150 /*! @brief Effective near-surface mean-annual air temperature. */
151 class AirTemperature : public Diag<AtmosphereModel> {
152 public:
153 AirTemperature(const AtmosphereModel *m)
154 : Diag<AtmosphereModel>(m) {
155
156 /* set metadata: */
157 m_vars = {SpatialVariableMetadata(m_sys, "effective_air_temp")};
158
159 set_attrs("effective mean-annual near-surface air temperature", "",
160 "Kelvin", "Kelvin", 0);
161 }
162 protected:
163 IceModelVec::Ptr compute_impl() const {
164
165 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "effective_air_temp", WITHOUT_GHOSTS));
166 result->metadata(0) = m_vars[0];
167
168 result->copy_from(model->mean_annual_temp());
169
170 return result;
171 }
172 };
173
174 /*! @brief Effective precipitation rate (average over time step). */
175 class Precipitation : public Diag<AtmosphereModel> {
176 public:
177 Precipitation(const AtmosphereModel *m)
178 : Diag<AtmosphereModel>(m) {
179
180 /* set metadata: */
181 m_vars = {SpatialVariableMetadata(m_sys, "effective_precipitation")};
182
183 set_attrs("effective precipitation rate",
184 "precipitation_flux",
185 "kg m-2 second-1", "kg m-2 year-1", 0);
186 }
187 protected:
188 IceModelVec::Ptr compute_impl() const {
189
190 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "effective_precipitation", WITHOUT_GHOSTS));
191 result->metadata(0) = m_vars[0];
192
193 result->copy_from(model->mean_precipitation());
194
195 return result;
196 }
197 };
198
199 } // end of namespace diagnostics
200
201 void AtmosphereModel::update_impl(const Geometry &geometry, double t, double dt) {
202 if (m_input_model) {
203 m_input_model->update_impl(geometry, t, dt);
204 }
205 }
206
207 MaxTimestep AtmosphereModel::max_timestep_impl(double my_t) const {
208 if (m_input_model) {
209 return m_input_model->max_timestep(my_t);
210 } else {
211 return MaxTimestep("atmosphere model");
212 }
213 }
214
215 DiagnosticList AtmosphereModel::diagnostics_impl() const {
216 using namespace diagnostics;
217
218 DiagnosticList result = {
219 {"air_temp_snapshot", Diagnostic::Ptr(new AirTemperatureSnapshot(this))},
220 {"effective_air_temp", Diagnostic::Ptr(new AirTemperature(this))},
221 {"effective_precipitation", Diagnostic::Ptr(new Precipitation(this))},
222 };
223
224 if (m_input_model) {
225 result = combine(result, m_input_model->diagnostics());
226 }
227
228 return result;
229 }
230
231 TSDiagnosticList AtmosphereModel::ts_diagnostics_impl() const {
232 if (m_input_model) {
233 return m_input_model->ts_diagnostics();
234 }
235
236 return {};
237 }
238
239 void AtmosphereModel::define_model_state_impl(const File &output) const {
240 if (m_input_model) {
241 m_input_model->define_model_state(output);
242 }
243 }
244
245 void AtmosphereModel::write_model_state_impl(const File &output) const {
246 if (m_input_model) {
247 m_input_model->write_model_state(output);
248 }
249 }
250
251 const IceModelVec2S& AtmosphereModel::mean_precipitation_impl() const {
252 if (m_input_model) {
253 return m_input_model->mean_precipitation();
254 } else {
255 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "no input model");
256 }
257 }
258
259 const IceModelVec2S& AtmosphereModel::mean_annual_temp_impl() const {
260 if (m_input_model) {
261 return m_input_model->mean_annual_temp();
262 } else {
263 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "no input model");
264 }
265 }
266
267 void AtmosphereModel::begin_pointwise_access_impl() const {
268 if (m_input_model) {
269 m_input_model->begin_pointwise_access();
270 } else {
271 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "no input model");
272 }
273 }
274
275 void AtmosphereModel::end_pointwise_access_impl() const {
276 if (m_input_model) {
277 m_input_model->end_pointwise_access();
278 } else {
279 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "no input model");
280 }
281 }
282
283 void AtmosphereModel::temp_time_series_impl(int i, int j, std::vector<double> &result) const {
284 if (m_input_model) {
285 m_input_model->temp_time_series(i, j, result);
286 } else {
287 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "no input model");
288 }
289 }
290
291 void AtmosphereModel::precip_time_series_impl(int i, int j, std::vector<double> &result) const {
292 if (m_input_model) {
293 m_input_model->precip_time_series(i, j, result);
294 } else {
295 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "no input model");
296 }
297 }
298
299 void AtmosphereModel::init_timeseries_impl(const std::vector<double> &ts) const {
300 if (m_input_model) {
301 m_input_model->init_timeseries(ts);
302 m_ts_times = ts;
303 } else {
304 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "no input model");
305 }
306 }
307
308 } // end of namespace atmosphere
309 } // end of namespace pism