tAnomaly.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
---
tAnomaly.cc (5693B)
---
1 // Copyright (C) 2011, 2012, 2013, 2014, 2015, 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 #include "Anomaly.hh"
20
21 #include "pism/util/ConfigInterface.hh"
22 #include "pism/util/IceGrid.hh"
23 #include "pism/util/io/io_helpers.hh"
24 #include "pism/coupler/util/options.hh"
25
26 namespace pism {
27 namespace atmosphere {
28
29 Anomaly::Anomaly(IceGrid::ConstPtr g, std::shared_ptr<AtmosphereModel> in)
30 : AtmosphereModel(g, in) {
31
32 ForcingOptions opt(*m_grid->ctx(), "atmosphere.anomaly");
33
34 {
35 unsigned int buffer_size = m_config->get_number("input.forcing.buffer_size");
36 unsigned int evaluations_per_year = m_config->get_number("input.forcing.evaluations_per_year");
37 bool periodic = opt.period > 0;
38
39 File file(m_grid->com, opt.filename, PISM_NETCDF3, PISM_READONLY);
40
41 m_air_temp_anomaly = IceModelVec2T::ForcingField(m_grid,
42 file,
43 "air_temp_anomaly",
44 "", // no standard name
45 buffer_size,
46 evaluations_per_year,
47 periodic,
48 LINEAR);
49
50 m_precipitation_anomaly = IceModelVec2T::ForcingField(m_grid,
51 file,
52 "precipitation_anomaly",
53 "", // no standard name
54 buffer_size,
55 evaluations_per_year,
56 periodic);
57 }
58
59 m_air_temp_anomaly->set_attrs("climate_forcing",
60 "anomaly of the near-surface air temperature",
61 "Kelvin", "Kelvin", "", 0);
62
63 m_precipitation_anomaly->set_attrs("climate_forcing",
64 "anomaly of the ice-equivalent precipitation rate",
65 "kg m-2 second-1", "kg m-2 year-1", "", 0);
66
67 m_precipitation = allocate_precipitation(g);
68 m_temperature = allocate_temperature(g);
69 }
70
71 Anomaly::~Anomaly()
72 {
73 // empty
74 }
75
76 void Anomaly::init_impl(const Geometry &geometry) {
77 m_input_model->init(geometry);
78
79 ForcingOptions opt(*m_grid->ctx(), "atmosphere.anomaly");
80
81 m_log->message(2,
82 "* Initializing the -atmosphere ...,anomaly code...\n");
83
84 m_log->message(2,
85 " reading anomalies from %s ...\n",
86 opt.filename.c_str());
87
88 m_air_temp_anomaly->init(opt.filename, opt.period, opt.reference_time);
89 m_precipitation_anomaly->init(opt.filename, opt.period, opt.reference_time);
90 }
91
92 void Anomaly::update_impl(const Geometry &geometry, double t, double dt) {
93 m_input_model->update(geometry, t, dt);
94
95 m_precipitation_anomaly->update(t, dt);
96 m_air_temp_anomaly->update(t, dt);
97
98 m_precipitation_anomaly->average(t, dt);
99 m_air_temp_anomaly->average(t, dt);
100
101 // precipitation
102 {
103 m_precipitation->copy_from(m_input_model->mean_precipitation());
104 m_precipitation->add(1.0, *m_precipitation_anomaly);
105 }
106
107 // temperature
108 {
109 m_temperature->copy_from(m_input_model->mean_annual_temp());
110 m_temperature->add(1.0, *m_air_temp_anomaly);
111 }
112 }
113
114 const IceModelVec2S& Anomaly::mean_precipitation_impl() const {
115 return *m_precipitation;
116 }
117
118 const IceModelVec2S& Anomaly::mean_annual_temp_impl() const {
119 return *m_temperature;
120 }
121
122 void Anomaly::begin_pointwise_access_impl() const {
123 m_input_model->begin_pointwise_access();
124 m_air_temp_anomaly->begin_access();
125 m_precipitation_anomaly->begin_access();
126 }
127
128 void Anomaly::end_pointwise_access_impl() const {
129 m_input_model->end_pointwise_access();
130 m_precipitation_anomaly->end_access();
131 m_air_temp_anomaly->end_access();
132 }
133
134 void Anomaly::init_timeseries_impl(const std::vector<double> &ts) const {
135 AtmosphereModel::init_timeseries_impl(ts);
136
137 m_air_temp_anomaly->init_interpolation(ts);
138
139 m_precipitation_anomaly->init_interpolation(ts);
140 }
141
142 void Anomaly::temp_time_series_impl(int i, int j, std::vector<double> &result) const {
143 m_input_model->temp_time_series(i, j, result);
144
145 m_temp_anomaly.reserve(m_ts_times.size());
146 m_air_temp_anomaly->interp(i, j, m_temp_anomaly);
147
148 for (unsigned int k = 0; k < m_ts_times.size(); ++k) {
149 result[k] += m_temp_anomaly[k];
150 }
151 }
152
153 void Anomaly::precip_time_series_impl(int i, int j, std::vector<double> &result) const {
154 m_input_model->precip_time_series(i, j, result);
155
156 m_mass_flux_anomaly.reserve(m_ts_times.size());
157 m_precipitation_anomaly->interp(i, j, m_mass_flux_anomaly);
158
159 for (unsigned int k = 0; k < m_ts_times.size(); ++k) {
160 result[k] += m_mass_flux_anomaly[k];
161 }
162 }
163
164 } // end of namespace atmosphere
165 } // end of namespace pism