tYearlyCycle.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
---
tYearlyCycle.cc (6765B)
---
1 // Copyright (C) 2008-2019 Ed Bueler, Constantine Khroulev, Ricarda Winkelmann,
2 // Gudfinna Adalgeirsdottir and Andy Aschwanden
3 //
4 // This file is part of PISM.
5 //
6 // PISM is free software; you can redistribute it and/or modify it under the
7 // terms of the GNU General Public License as published by the Free Software
8 // Foundation; either version 3 of the License, or (at your option) any later
9 // version.
10 //
11 // PISM is distributed in the hope that it will be useful, but WITHOUT ANY
12 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
13 // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
14 // details.
15 //
16 // You should have received a copy of the GNU General Public License
17 // along with PISM; if not, write to the Free Software
18 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
19
20 // Implementation of the atmosphere model using constant-in-time precipitation
21 // and a cosine yearly cycle for near-surface air temperatures.
22
23 #include <gsl/gsl_math.h> // M_PI
24
25 #include "YearlyCycle.hh"
26 #include "pism/util/Time.hh"
27 #include "pism/util/IceGrid.hh"
28 #include "pism/util/ConfigInterface.hh"
29 #include "pism/util/io/io_helpers.hh"
30 #include "pism/util/pism_utilities.hh"
31
32 namespace pism {
33 namespace atmosphere {
34
35 YearlyCycle::YearlyCycle(IceGrid::ConstPtr g)
36 : AtmosphereModel(g),
37 m_air_temp_mean_annual(m_grid, "air_temp_mean_annual", WITHOUT_GHOSTS),
38 m_air_temp_mean_summer(m_grid, "air_temp_mean_summer", WITHOUT_GHOSTS),
39 m_precipitation(m_grid, "precipitation", WITHOUT_GHOSTS) {
40
41 m_snow_temp_summer_day = m_config->get_number("atmosphere.fausto_air_temp.summer_peak_day");
42
43 m_air_temp_mean_annual.set_attrs("diagnostic",
44 "mean annual near-surface air temperature (without sub-year time-dependence or forcing)",
45 "K", "K",
46 "", 0); // no CF standard_name
47 m_air_temp_mean_annual.metadata().set_string("source", m_reference);
48
49 m_air_temp_mean_summer.set_attrs("diagnostic",
50 "mean summer (NH: July/ SH: January) near-surface air temperature (without sub-year time-dependence or forcing)",
51 "Kelvin", "Kelvin",
52 "", 0); // no CF standard_name
53 m_air_temp_mean_summer.metadata().set_string("source", m_reference);
54
55 m_precipitation.set_attrs("model_state", "precipitation rate",
56 "kg m-2 second-1", "kg m-2 year-1", "precipitation_flux", 0);
57 m_precipitation.set_time_independent(true);
58 }
59
60 YearlyCycle::~YearlyCycle() {
61 // empty
62 }
63
64 //! Reads in the precipitation data from the input file.
65 void YearlyCycle::init_impl(const Geometry &geometry) {
66 (void) geometry;
67
68 InputOptions opts = process_input_options(m_grid->com, m_config);
69 init_internal(opts.filename, opts.type == INIT_BOOTSTRAP, opts.record);
70 }
71
72 //! Read precipitation data from a given file.
73 void YearlyCycle::init_internal(const std::string &input_filename, bool do_regrid,
74 unsigned int start) {
75 // read precipitation rate from file
76 m_log->message(2,
77 " reading mean annual ice-equivalent precipitation rate 'precipitation'\n"
78 " from %s ... \n",
79 input_filename.c_str());
80 if (do_regrid == true) {
81 m_precipitation.regrid(input_filename, CRITICAL); // fails if not found!
82 } else {
83 m_precipitation.read(input_filename, start); // fails if not found!
84 }
85 }
86
87 void YearlyCycle::define_model_state_impl(const File &output) const {
88 m_precipitation.define(output);
89 }
90
91 void YearlyCycle::write_model_state_impl(const File &output) const {
92 m_precipitation.write(output);
93 }
94
95 //! Copies the stored precipitation field into result.
96 const IceModelVec2S& YearlyCycle::mean_precipitation_impl() const {
97 return m_precipitation;
98 }
99
100 //! Copies the stored mean annual near-surface air temperature field into result.
101 const IceModelVec2S& YearlyCycle::mean_annual_temp_impl() const {
102 return m_air_temp_mean_annual;
103 }
104
105 //! Copies the stored mean summer near-surface air temperature field into result.
106 const IceModelVec2S& YearlyCycle::mean_summer_temp() const {
107 return m_air_temp_mean_summer;
108 }
109
110 void YearlyCycle::init_timeseries_impl(const std::vector<double> &ts) const {
111 // constants related to the standard yearly cycle
112 const double
113 summerday_fraction = m_grid->ctx()->time()->day_of_the_year_to_day_fraction(m_snow_temp_summer_day);
114
115 size_t N = ts.size();
116
117 m_ts_times.resize(N);
118 m_cosine_cycle.resize(N);
119 for (unsigned int k = 0; k < m_ts_times.size(); k++) {
120 double tk = m_grid->ctx()->time()->year_fraction(ts[k]) - summerday_fraction;
121
122 m_ts_times[k] = ts[k];
123 m_cosine_cycle[k] = cos(2.0 * M_PI * tk);
124 }
125 }
126
127 void YearlyCycle::precip_time_series_impl(int i, int j, std::vector<double> &result) const {
128 result.resize(m_ts_times.size());
129 for (unsigned int k = 0; k < m_ts_times.size(); k++) {
130 result[k] = m_precipitation(i,j);
131 }
132 }
133
134 void YearlyCycle::temp_time_series_impl(int i, int j, std::vector<double> &result) const {
135 result.resize(m_ts_times.size());
136 for (unsigned int k = 0; k < m_ts_times.size(); ++k) {
137 result[k] = m_air_temp_mean_annual(i,j) + (m_air_temp_mean_summer(i,j) - m_air_temp_mean_annual(i,j)) * m_cosine_cycle[k];
138 }
139 }
140
141 void YearlyCycle::begin_pointwise_access_impl() const {
142 m_air_temp_mean_annual.begin_access();
143 m_air_temp_mean_summer.begin_access();
144 m_precipitation.begin_access();
145 }
146
147 void YearlyCycle::end_pointwise_access_impl() const {
148 m_air_temp_mean_annual.end_access();
149 m_air_temp_mean_summer.end_access();
150 m_precipitation.end_access();
151 }
152
153 namespace diagnostics {
154
155 /*! @brief Mean summer near-surface air temperature. */
156 class MeanSummerTemperature : public Diag<YearlyCycle>
157 {
158 public:
159 MeanSummerTemperature(const YearlyCycle *m)
160 : Diag<YearlyCycle>(m) {
161
162 /* set metadata: */
163 m_vars = {SpatialVariableMetadata(m_sys, "air_temp_mean_summer")};
164
165 set_attrs("mean summer near-surface air temperature used in the cosine yearly cycle", "",
166 "Kelvin", "Kelvin", 0);
167 }
168 private:
169 IceModelVec::Ptr compute_impl() const {
170
171 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "air_temp_mean_summer", WITHOUT_GHOSTS));
172 result->metadata(0) = m_vars[0];
173
174 result->copy_from(model->mean_summer_temp());
175
176 return result;
177 }
178 };
179 } // end of namespace diagnostics
180
181 DiagnosticList YearlyCycle::diagnostics_impl() const {
182 DiagnosticList result = AtmosphereModel::diagnostics_impl();
183
184 result["air_temp_mean_summer"] = Diagnostic::Ptr(new diagnostics::MeanSummerTemperature(this));
185
186 return result;
187 }
188
189 } // end of namespace atmosphere
190 } // end of namespace pism