tISMIP6Climate.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
---
tISMIP6Climate.cc (9571B)
---
1 // Copyright (C) 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 #include "ISMIP6Climate.hh"
20
21 #include "pism/util/IceGrid.hh"
22 #include "pism/coupler/util/options.hh"
23 #include "pism/util/io/io_helpers.hh"
24 #include "pism/geometry/Geometry.hh"
25
26 namespace pism {
27 namespace surface {
28
29 ISMIP6::ISMIP6(IceGrid::ConstPtr grid, std::shared_ptr<atmosphere::AtmosphereModel> input)
30 : SurfaceModel(grid),
31 m_mass_flux_reference(m_grid, "climatic_mass_balance", WITHOUT_GHOSTS),
32 m_temperature_reference(m_grid, "ice_surface_temp", WITHOUT_GHOSTS),
33 m_surface_reference(m_grid, "usurf", WITHOUT_GHOSTS)
34 {
35 (void) input;
36
37 // allocate model outputs
38 m_temperature = allocate_temperature(m_grid);
39 m_mass_flux = allocate_mass_flux(m_grid);
40
41 // set metadata of reference fields
42 {
43 m_mass_flux_reference.set_attrs("climate_forcing",
44 "reference surface mass balance rate",
45 "kg m-2 s-1", "kg m-2 year-1",
46 "land_ice_surface_specific_mass_balance_flux", 0);
47
48 auto smb_max = m_config->get_number("surface.given.smb_max", "kg m-2 second-1");
49 m_mass_flux_reference.metadata().set_numbers("valid_range", {-smb_max, smb_max});
50 m_mass_flux_reference.set_time_independent(true);
51
52 m_surface_reference.set_attrs("climate_forcing",
53 "reference surface altitude",
54 "m", "m", "surface_altitude", 0);
55
56 m_surface_reference.metadata().set_numbers("valid_range", {0.0, m_grid->Lz()});
57 m_surface_reference.set_time_independent(true);
58
59 m_temperature_reference.set_attrs("climate_forcing",
60 "reference temperature",
61 "Kelvin", "Kelvin", "", 0);
62
63 m_temperature_reference.metadata().set_numbers("valid_range", {0.0, 373.15});
64 m_temperature_reference.set_time_independent(true);
65 }
66
67 // allocate storage for time-dependent inputs
68 ForcingOptions opt(*m_grid->ctx(), "surface.ismip6");
69
70 {
71 unsigned int buffer_size = m_config->get_number("input.forcing.buffer_size");
72 unsigned int evaluations_per_year = m_config->get_number("input.forcing.evaluations_per_year");
73 bool periodic = opt.period > 0;
74
75 File file(m_grid->com, opt.filename, PISM_NETCDF3, PISM_READONLY);
76
77 {
78 m_mass_flux_anomaly = IceModelVec2T::ForcingField(m_grid,
79 file,
80 "climatic_mass_balance_anomaly",
81 "", // no standard name
82 buffer_size,
83 evaluations_per_year,
84 periodic);
85
86 m_mass_flux_anomaly->set_attrs("climate_forcing",
87 "surface mass balance rate anomaly",
88 "kg m-2 s-1", "kg m-2 year-1",
89 "", 0);
90
91 }
92
93 {
94 m_mass_flux_gradient = IceModelVec2T::ForcingField(m_grid,
95 file,
96 "climatic_mass_balance_gradient",
97 "", // no standard name
98 buffer_size,
99 evaluations_per_year,
100 periodic);
101
102 m_mass_flux_gradient->set_attrs("climate_forcing",
103 "surface mass balance rate elevation lapse rate",
104 "kg m-2 s-1 m-1", "kg m-2 year-1 m-1",
105 "", 0);
106 }
107
108 {
109 m_temperature_anomaly = IceModelVec2T::ForcingField(m_grid,
110 file,
111 "ice_surface_temp_anomaly",
112 "", // no standard name
113 buffer_size,
114 evaluations_per_year,
115 periodic);
116
117 m_temperature_anomaly->set_attrs("climate_forcing",
118 "ice surface temperature anomaly",
119 "Kelvin", "Kelvin", "", 0);
120 }
121
122 {
123 m_temperature_gradient = IceModelVec2T::ForcingField(m_grid,
124 file,
125 "ice_surface_temp_gradient",
126 "", // no standard name
127 buffer_size,
128 evaluations_per_year,
129 periodic);
130
131 m_temperature_gradient->set_attrs("climate_forcing",
132 "ice surface temperature elevation lapse rate",
133 "Kelvin m-1", "Kelvin m-1", "", 0);
134 }
135
136 }
137 }
138
139 ISMIP6::~ISMIP6() {
140 // empty
141 }
142
143 void ISMIP6::init_impl(const Geometry &geometry) {
144 (void) geometry;
145
146 m_log->message(2, "* Initializing the ISMIP6 surface model...\n");
147
148 {
149 // File with reference surface elevation, temperature, and climatic mass balance
150 auto reference_filename = m_config->get_string("surface.ismip6.reference_file");
151 File reference_file(m_grid->com, reference_filename, PISM_GUESS, PISM_READONLY);
152
153 m_mass_flux_reference.regrid(reference_file, CRITICAL);
154 m_surface_reference.regrid(reference_file, CRITICAL);
155 m_temperature_reference.regrid(reference_file, CRITICAL);
156 }
157
158 {
159 ForcingOptions opt(*m_grid->ctx(), "surface.ismip6");
160
161 m_mass_flux_anomaly->init(opt.filename, opt.period, opt.reference_time);
162 m_mass_flux_gradient->init(opt.filename, opt.period, opt.reference_time);
163
164 m_temperature_anomaly->init(opt.filename, opt.period, opt.reference_time);
165 m_temperature_gradient->init(opt.filename, opt.period, opt.reference_time);
166 }
167 }
168
169 void ISMIP6::update_impl(const Geometry &geometry, double t, double dt) {
170
171 // inputs
172 const IceModelVec2S &h = geometry.ice_surface_elevation;
173 const IceModelVec2S &h_ref = m_surface_reference;
174 const IceModelVec2S &T_ref = m_temperature_reference;
175 const IceModelVec2S &SMB_ref = m_mass_flux_reference;
176
177 IceModelVec2T &dTdz = *m_temperature_gradient;
178 IceModelVec2T &dSMBdz = *m_mass_flux_gradient;
179 IceModelVec2T &aT = *m_temperature_anomaly;
180 IceModelVec2T &aSMB = *m_mass_flux_anomaly;
181
182 // outputs
183 IceModelVec2S &T = *m_temperature;
184 IceModelVec2S &SMB = *m_mass_flux;
185
186 // get time-dependent input fields at the current time
187 {
188 aT.update(t, dt);
189 aSMB.update(t, dt);
190 dTdz.update(t, dt);
191 dSMBdz.update(t, dt);
192
193 aT.average(t, dt);
194 aSMB.average(t, dt);
195 dTdz.average(t, dt);
196 dSMBdz.average(t, dt);
197 }
198
199 // From http://www.climate-cryosphere.org/wiki/index.php?title=ISMIP6-Projections-Greenland:
200 // SMB(x,y,t) = SMB_ref(x,y) + aSMB(x,y,t) + dSMBdz(x,y,t) * [h(x,y,t) - h_ref(x,y)]
201
202 IceModelVec::AccessList list{&h, &h_ref,
203 &SMB, &SMB_ref, &aSMB, &dSMBdz,
204 &T, &T_ref, &aT, &dTdz};
205
206 for (Points p(*m_grid); p; p.next()) {
207 const int i = p.i(), j = p.j();
208
209 SMB(i, j) = SMB_ref(i, j) + aSMB(i, j) + dSMBdz(i, j) * (h(i, j) - h_ref(i, j));
210 T(i, j) = T_ref(i, j) + aT(i, j) + dTdz(i, j) * (h(i, j) - h_ref(i, j));
211 }
212
213 dummy_accumulation(SMB, *m_accumulation);
214 dummy_melt(SMB, *m_melt);
215 dummy_runoff(SMB, *m_runoff);
216 }
217
218 MaxTimestep ISMIP6::max_timestep_impl(double t) const {
219 using std::min;
220
221 auto dt = m_temperature_anomaly->max_timestep(t);
222 dt = min(dt, m_temperature_gradient->max_timestep(t));
223 dt = min(dt, m_mass_flux_anomaly->max_timestep(t));
224 dt = min(dt, m_mass_flux_gradient->max_timestep(t));
225
226 if (dt.finite()) {
227 return MaxTimestep(dt.value(), "surface ISMIP6");
228 } else {
229 return MaxTimestep("surface ISMIP6");
230 }
231 }
232
233 const IceModelVec2S &ISMIP6::mass_flux_impl() const {
234 return *m_mass_flux;
235 }
236
237 const IceModelVec2S &ISMIP6::temperature_impl() const {
238 return *m_temperature;
239 }
240
241 const IceModelVec2S &ISMIP6::accumulation_impl() const {
242 return *m_accumulation;
243 }
244
245 const IceModelVec2S &ISMIP6::melt_impl() const {
246 return *m_melt;
247 }
248
249 const IceModelVec2S &ISMIP6::runoff_impl() const {
250 return *m_runoff;
251 }
252
253 } // end of namespace surface
254 } // end of namespace pism