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 (6052B)
---
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 #include "pism/coupler/util/options.hh"
21 #include "pism/coupler/util/lapse_rates.hh"
22 #include "pism/util/io/io_helpers.hh"
23 #include "pism/util/pism_utilities.hh"
24 #include "pism/util/pism_options.hh"
25 #include "pism/geometry/Geometry.hh"
26
27 namespace pism {
28 namespace surface {
29
30 ElevationChange::ElevationChange(IceGrid::ConstPtr g, std::shared_ptr<SurfaceModel> in)
31 : SurfaceModel(g, in) {
32
33 {
34 m_smb_lapse_rate = m_config->get_number("surface.elevation_change.smb.lapse_rate",
35 "(m / s) / m");
36 // convert from [m s-1 / m] to [kg m-2 s-1 / m]
37 m_smb_lapse_rate *= m_config->get_number("constants.ice.density");
38
39 m_smb_exp_factor = m_config->get_number("surface.elevation_change.smb.exp_factor");
40 }
41
42 {
43 auto method = m_config->get_string("surface.elevation_change.smb.method");
44 m_smb_method = method == "scale" ? SCALE : SHIFT;
45 }
46
47 m_temp_lapse_rate = m_config->get_number("surface.elevation_change.temperature_lapse_rate",
48 "K / m");
49
50 {
51 ForcingOptions opt(*m_grid->ctx(), "surface.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_mass_flux = allocate_mass_flux(g);
72 m_temperature = allocate_temperature(g);
73 m_accumulation = allocate_accumulation(g);
74 m_melt = allocate_melt(g);
75 m_runoff = allocate_runoff(g);
76 }
77
78 ElevationChange::~ElevationChange() {
79 // empty
80 }
81
82 void ElevationChange::init_impl(const Geometry &geometry) {
83 using units::convert;
84
85 m_input_model->init(geometry);
86
87 m_log->message(2,
88 " [using temperature and mass balance lapse corrections]\n");
89
90 m_log->message(2,
91 " ice upper-surface temperature lapse rate: %3.3f K per km\n",
92 convert(m_sys, m_temp_lapse_rate, "K / m", "K / km"));
93
94 if (m_smb_method == SHIFT) {
95 double ice_density = m_config->get_number("constants.ice.density");
96 m_log->message(2,
97 " ice-equivalent surface mass balance lapse rate: %3.3f m year-1 per km\n",
98 convert(m_sys, m_smb_lapse_rate, "kg / (m2 second)", "kg / (m2 year)") / ice_density);
99 } else {
100 m_log->message(2,
101 " surface mass balance scaling factor with temperature: %3.3f Kelvin-1\n",
102 m_smb_exp_factor);
103 }
104
105 ForcingOptions opt(*m_grid->ctx(), "surface.elevation_change");
106 m_reference_surface->init(opt.filename, opt.period, opt.reference_time);
107 }
108
109 void ElevationChange::update_impl(const Geometry &geometry, double t, double dt) {
110
111 m_input_model->update(geometry, t, dt);
112
113 m_reference_surface->update(t, dt);
114 m_reference_surface->interp(t + 0.5*dt);
115
116 const IceModelVec2S &surface = geometry.ice_surface_elevation;
117
118 m_temperature->copy_from(m_input_model->temperature());
119 lapse_rate_correction(surface, *m_reference_surface,
120 m_temp_lapse_rate, *m_temperature);
121
122 m_mass_flux->copy_from(m_input_model->mass_flux());
123
124 switch (m_smb_method) {
125 case SCALE:
126 {
127 IceModelVec::AccessList list{&surface, m_reference_surface.get(), m_mass_flux.get()};
128
129 for (Points p(*m_grid); p; p.next()) {
130 const int i = p.i(), j = p.j();
131
132 double dT = -m_temp_lapse_rate * (surface(i, j) - (*m_reference_surface)(i, j));
133
134 (*m_mass_flux)(i, j) *= exp(m_smb_exp_factor * dT);
135 }
136 }
137 break;
138 default:
139 case SHIFT:
140 {
141 lapse_rate_correction(surface, *m_reference_surface,
142 m_smb_lapse_rate, *m_mass_flux);
143 }
144 break;
145 }
146
147 // This modifier changes m_mass_flux, so we need to compute accumulation, melt, and
148 // runoff.
149 dummy_accumulation(*m_mass_flux, *m_accumulation);
150 dummy_melt(*m_mass_flux, *m_melt);
151 dummy_runoff(*m_mass_flux, *m_runoff);
152
153 }
154
155 const IceModelVec2S &ElevationChange::mass_flux_impl() const {
156 return *m_mass_flux;
157 }
158
159 const IceModelVec2S &ElevationChange::temperature_impl() const {
160 return *m_temperature;
161 }
162
163 const IceModelVec2S &ElevationChange::accumulation_impl() const {
164 return *m_accumulation;
165 }
166
167 const IceModelVec2S &ElevationChange::melt_impl() const {
168 return *m_melt;
169 }
170
171 const IceModelVec2S &ElevationChange::runoff_impl() const {
172 return *m_runoff;
173 }
174
175
176 } // end of namespace surface
177 } // end of namespace pism