PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
StuffAsAnomaly.cc
Go to the documentation of this file.
1 // Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2021, 2023 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 "pism/coupler/surface/StuffAsAnomaly.hh"
20 #include "pism/util/Grid.hh"
21 #include "pism/util/Time.hh"
22 #include "pism/util/pism_utilities.hh"
23 #include "pism/util/MaxTimestep.hh"
24 
25 namespace pism {
26 namespace surface {
27 
28 StuffAsAnomaly::StuffAsAnomaly(std::shared_ptr<const Grid> g, std::shared_ptr<SurfaceModel> input)
29  : SurfaceModel(g, input),
30  m_mass_flux(m_grid, "climatic_mass_balance", WITHOUT_GHOSTS),
31  m_mass_flux_0(m_grid, "mass_flux_0", WITHOUT_GHOSTS),
32  m_mass_flux_input(m_grid, "climatic_mass_balance", WITHOUT_GHOSTS),
33  m_temp(m_grid, "ice_surface_temp", WITHOUT_GHOSTS),
34  m_temp_0(m_grid, "ice_surface_temp_0", WITHOUT_GHOSTS),
35  m_temp_input(m_grid, "ice_surface_temp", WITHOUT_GHOSTS) {
36 
37  m_mass_flux.set_attrs("climate_state",
38  "surface mass balance (accumulation/ablation) rate",
39  "kg m-2 s-1",
40  "land_ice_surface_specific_mass_balance_flux");
41  m_mass_flux.metadata()["output_units"] = "kg m-2 year-1";
42 
43  m_temp.set_attrs("climate_state", "ice temperature at the ice surface",
44  "K", "");
45 
46  // create special variables
47  m_mass_flux_0.set_attrs("internal", "surface mass flux at the beginning of a run",
48  "kg m-2 s-1", "land_ice_surface_specific_mass_balance_flux");
49 
50  m_mass_flux_input.set_attrs("model_state", "surface mass flux to apply anomalies to",
51  "kg m-2 s-1", "land_ice_surface_specific_mass_balance_flux");
52 
53  m_temp_0.set_attrs("internal", "ice-surface temperature and the beginning of a run", "K",
54  "");
55 
56  m_temp_input.set_attrs("model_state", "ice-surface temperature to apply anomalies to",
57  "K", "");
58 }
59 
60 void StuffAsAnomaly::init_impl(const Geometry &geometry) {
61  if (m_input_model != NULL) {
62  m_input_model->init(geometry);
63  }
64 
65  InputOptions opts = process_input_options(m_grid->com, m_config);
66 
67  m_log->message(2,
68  "* Initializing the 'turn_into_anomaly' modifier\n"
69  " (it applies climate data as anomalies relative to 'ice_surface_temp' and 'climatic_mass_balance'\n"
70  " read from '%s'.\n", opts.filename.c_str());
71 
72  if (opts.type == INIT_BOOTSTRAP) {
73  m_mass_flux_input.regrid(opts.filename, CRITICAL);
74  m_temp_input.regrid(opts.filename, CRITICAL);
75  } else {
76  m_mass_flux_input.read(opts.filename, opts.record);
77  m_temp_input.read(opts.filename, opts.record);
78  }
79 }
80 
81 MaxTimestep StuffAsAnomaly::max_timestep_impl(double t) const {
82  (void) t;
83  return MaxTimestep("surface turn_into_anomaly");
84 }
85 
86 void StuffAsAnomaly::update_impl(const Geometry &geometry, double t, double dt) {
87 
88  if (m_input_model != NULL) {
89  m_input_model->update(geometry, t, dt);
90  m_temp.copy_from(m_input_model->temperature());
91  m_mass_flux.copy_from(m_input_model->mass_flux());
92 
93  // if we are at the beginning of the run...
94  if (t < time().start() + 1) {
95  // this is goofy, but time-steps are usually longer than 1 second, so it should work
96  m_temp_0.copy_from(m_temp);
97  m_mass_flux_0.copy_from(m_mass_flux);
98  }
99  }
100 
101  array::AccessScope list{&m_mass_flux, &m_mass_flux_0, &m_mass_flux_input,
102  &m_temp, &m_temp_0, &m_temp_input};
103 
104  for (auto p = m_grid->points(); p; p.next()) {
105  const int i = p.i(), j = p.j();
106 
107  m_mass_flux(i, j) = m_mass_flux(i, j) - m_mass_flux_0(i, j) + m_mass_flux_input(i, j);
108  m_temp(i, j) = m_temp(i, j) - m_temp_0(i, j) + m_temp_input(i, j);
109  }
110 }
111 
112 const array::Scalar &StuffAsAnomaly::mass_flux_impl() const {
113  return m_mass_flux;
114 }
115 
116 const array::Scalar &StuffAsAnomaly::temperature_impl() const {
117  return m_temp;
118 }
119 
120 } // end of namespace surface
121 } // end of namespace pism
@ WITHOUT_GHOSTS
Definition: Array.hh:62
InputOptions process_input_options(MPI_Comm com, Config::ConstPtr config)
Definition: Component.cc:43
static const double g
Definition: exactTestP.cc:36
@ INIT_BOOTSTRAP
Definition: Component.hh:56