URI:
       tStuffAsAnomaly.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
       ---
       tStuffAsAnomaly.cc (4588B)
       ---
            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 "StuffAsAnomaly.hh"
           20 #include "pism/util/IceGrid.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(IceGrid::ConstPtr 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().set_string("glaciological_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 StuffAsAnomaly::~StuffAsAnomaly() {
           61   // empty
           62 }
           63 
           64 void StuffAsAnomaly::init_impl(const Geometry &geometry) {
           65   if (m_input_model != NULL) {
           66     m_input_model->init(geometry);
           67   }
           68 
           69   InputOptions opts = process_input_options(m_grid->com, m_config);
           70 
           71   m_log->message(2,
           72              "* Initializing the 'turn_into_anomaly' modifier\n"
           73              "  (it applies climate data as anomalies relative to 'ice_surface_temp' and 'climatic_mass_balance'\n"
           74              "  read from '%s'.\n", opts.filename.c_str());
           75 
           76   if (opts.type == INIT_BOOTSTRAP) {
           77     m_mass_flux_input.regrid(opts.filename, CRITICAL);
           78     m_temp_input.regrid(opts.filename, CRITICAL);
           79   } else {
           80     m_mass_flux_input.read(opts.filename, opts.record);
           81     m_temp_input.read(opts.filename, opts.record);
           82   }
           83 }
           84 
           85 MaxTimestep StuffAsAnomaly::max_timestep_impl(double t) const {
           86   (void) t;
           87   return MaxTimestep("surface turn_into_anomaly");
           88 }
           89 
           90 void StuffAsAnomaly::update_impl(const Geometry &geometry, double t, double dt) {
           91 
           92   if (m_input_model != NULL) {
           93     m_input_model->update(geometry, t, dt);
           94     m_temp.copy_from(m_input_model->temperature());
           95     m_mass_flux.copy_from(m_input_model->mass_flux());
           96 
           97     // if we are at the beginning of the run...
           98     if (t < m_grid->ctx()->time()->start() + 1) {
           99       // this is goofy, but time-steps are usually longer than 1 second, so it should work
          100       m_temp_0.copy_from(m_temp);
          101       m_mass_flux_0.copy_from(m_mass_flux);
          102     }
          103   }
          104 
          105   IceModelVec::AccessList list{&m_mass_flux, &m_mass_flux_0, &m_mass_flux_input,
          106       &m_temp, &m_temp_0, &m_temp_input};
          107 
          108   for (Points p(*m_grid); p; p.next()) {
          109     const int i = p.i(), j = p.j();
          110 
          111     m_mass_flux(i, j) = m_mass_flux(i, j) - m_mass_flux_0(i, j) + m_mass_flux_input(i, j);
          112     m_temp(i, j)      = m_temp(i, j) - m_temp_0(i, j) + m_temp_input(i, j);
          113   }
          114 }
          115 
          116 const IceModelVec2S &StuffAsAnomaly::mass_flux_impl() const {
          117   return m_mass_flux;
          118 }
          119 
          120 const IceModelVec2S &StuffAsAnomaly::temperature_impl() const {
          121   return m_temp;
          122 }
          123 
          124 } // end of namespace surface
          125 } // end of namespace pism