tElevationChange.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
---
tElevationChange.cc (7417B)
---
1 // Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020 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 "ElevationChange.hh"
20
21 #include <cmath> // std::exp()
22
23 #include "pism/coupler/util/options.hh"
24 #include "pism/coupler/util/lapse_rates.hh"
25 #include "pism/util/io/io_helpers.hh"
26 #include "pism/util/pism_utilities.hh"
27 #include "pism/util/pism_options.hh"
28 #include "pism/geometry/Geometry.hh"
29
30 namespace pism {
31 namespace atmosphere {
32
33 ElevationChange::ElevationChange(IceGrid::ConstPtr grid, std::shared_ptr<AtmosphereModel> in)
34 : AtmosphereModel(grid, in),
35 m_surface(grid, "ice_surface_elevation", WITHOUT_GHOSTS) {
36
37 m_precip_lapse_rate = m_config->get_number("atmosphere.elevation_change.precipitation.lapse_rate",
38 "(kg m-2 / s) / m");
39
40 m_precip_exp_factor = m_config->get_number("atmosphere.precip_exponential_factor_for_temperature");
41
42 m_temp_lapse_rate = m_config->get_number("atmosphere.elevation_change.temperature_lapse_rate",
43 "K / m");
44
45 {
46 auto method = m_config->get_string("atmosphere.elevation_change.precipitation.method");
47 m_precip_method = method == "scale" ? SCALE : SHIFT;
48 }
49
50 {
51 ForcingOptions opt(*m_grid->ctx(), "atmosphere.elevation_change");
52
53 unsigned int buffer_size = m_config->get_number("input.forcing.buffer_size");
54 unsigned int evaluations_per_year = m_config->get_number("input.forcing.evaluations_per_year");
55 bool periodic = opt.period > 0;
56
57 File file(m_grid->com, opt.filename, PISM_NETCDF3, PISM_READONLY);
58
59 m_reference_surface = IceModelVec2T::ForcingField(m_grid,
60 file,
61 "usurf",
62 "", // no standard name
63 buffer_size,
64 evaluations_per_year,
65 periodic,
66 LINEAR);
67 m_reference_surface->set_attrs("climate_forcing", "ice surface elevation",
68 "m", "m", "surface_altitude", 0);
69 }
70
71 m_precipitation = allocate_precipitation(grid);
72 m_temperature = allocate_temperature(grid);
73 }
74
75 ElevationChange::~ElevationChange() {
76 // empty
77 }
78
79 void ElevationChange::init_impl(const Geometry &geometry) {
80 using units::convert;
81
82 m_input_model->init(geometry);
83
84 m_log->message(2,
85 " [using elevation-change-dependent adjustments of air temperature and precipitation]\n");
86
87 m_log->message(2,
88 " air temperature lapse rate: %3.3f K per km\n",
89 convert(m_sys, m_temp_lapse_rate, "K / m", "K / km"));
90
91 if (m_precip_method == SHIFT) {
92 m_log->message(2,
93 " precipitation lapse rate: %3.3f (kg m-2 year-1) per km\n",
94 convert(m_sys, m_precip_lapse_rate, "(kg m-2 / s) / m", "(kg m-2 / year) / km"));
95 } else {
96 m_log->message(2,
97 " precipitation scaling factor with temperature: %3.3f Kelvin-1\n",
98 m_precip_exp_factor);
99 }
100
101 ForcingOptions opt(*m_grid->ctx(), "atmosphere.elevation_change");
102
103 m_reference_surface->init(opt.filename, opt.period, opt.reference_time);
104 }
105
106 void ElevationChange::update_impl(const Geometry &geometry, double t, double dt) {
107
108 m_input_model->update(geometry, t, dt);
109
110 m_reference_surface->update(t, dt);
111 m_reference_surface->interp(t + 0.5*dt);
112
113 // make a copy of the surface elevation so that it is available in methods computing
114 // temperature and precipitation time series
115 m_surface.copy_from(geometry.ice_surface_elevation);
116
117 // temperature
118 {
119 m_temperature->copy_from(m_input_model->mean_annual_temp());
120
121 lapse_rate_correction(m_surface, *m_reference_surface,
122 m_temp_lapse_rate, *m_temperature);
123 }
124
125 // precipitation
126 {
127 m_precipitation->copy_from(m_input_model->mean_precipitation());
128
129 switch (m_precip_method) {
130 case SCALE:
131 {
132 const IceModelVec2S &input_temperature = m_input_model->mean_annual_temp();
133
134 IceModelVec::AccessList list{m_temperature.get(),
135 &input_temperature,
136 m_precipitation.get()};
137
138 for (Points p(*m_grid); p; p.next()) {
139 const int i = p.i(), j = p.j();
140
141 double dT = (*m_temperature)(i, j) - input_temperature(i, j);
142
143 (*m_precipitation)(i, j) *= std::exp(m_precip_exp_factor * dT);
144 }
145
146 }
147 break;
148 case SHIFT:
149 default:
150 {
151 lapse_rate_correction(m_surface, *m_reference_surface,
152 m_precip_lapse_rate, *m_precipitation);
153 }
154 break;
155 }
156 }
157 }
158
159 const IceModelVec2S& ElevationChange::mean_annual_temp_impl() const {
160 return *m_temperature;
161 }
162
163 const IceModelVec2S& ElevationChange::mean_precipitation_impl() const {
164 return *m_precipitation;
165 }
166
167 void ElevationChange::begin_pointwise_access_impl() const {
168 m_input_model->begin_pointwise_access();
169
170 m_reference_surface->begin_access();
171 m_surface.begin_access();
172 }
173
174 void ElevationChange::end_pointwise_access_impl() const {
175 m_input_model->end_pointwise_access();
176
177 m_reference_surface->end_access();
178 m_surface.end_access();
179 }
180
181 void ElevationChange::init_timeseries_impl(const std::vector<double> &ts) const {
182 AtmosphereModel::init_timeseries_impl(ts);
183
184 m_reference_surface->init_interpolation(ts);
185 }
186
187 void ElevationChange::temp_time_series_impl(int i, int j, std::vector<double> &result) const {
188 std::vector<double> usurf(m_ts_times.size());
189
190 m_input_model->temp_time_series(i, j, result);
191
192 m_reference_surface->interp(i, j, usurf);
193
194 for (unsigned int m = 0; m < m_ts_times.size(); ++m) {
195 result[m] -= m_temp_lapse_rate * (m_surface(i, j) - usurf[m]);
196 }
197 }
198
199 void ElevationChange::precip_time_series_impl(int i, int j, std::vector<double> &result) const {
200 auto N = m_ts_times.size();
201 std::vector<double> usurf(N);
202
203 m_input_model->precip_time_series(i, j, result);
204
205 m_reference_surface->interp(i, j, usurf);
206
207 switch (m_precip_method) {
208 case SCALE:
209 {
210 for (unsigned int m = 0; m < N; ++m) {
211 double dT = -m_temp_lapse_rate * (m_surface(i, j) - usurf[m]);
212 result[m] *= std::exp(m_precip_exp_factor * dT);
213 }
214 }
215 break;
216 case SHIFT:
217 for (unsigned int m = 0; m < N; ++m) {
218 result[m] -= m_precip_lapse_rate * (m_surface(i, j) - usurf[m]);
219 }
220 break;
221 }
222 }
223
224 } // end of namespace atmosphere
225 } // end of namespace pism