URI: 
       tPrescribedRetreat.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
       ---
       tPrescribedRetreat.cc (3493B)
       ---
            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 
           20 #include "PrescribedRetreat.hh"
           21 #include "pism/coupler/util/options.hh"
           22 
           23 namespace pism {
           24 
           25 PrescribedRetreat::PrescribedRetreat(IceGrid::ConstPtr grid)
           26   : Component(grid) {
           27   ForcingOptions opt(*m_grid->ctx(), "geometry.front_retreat.prescribed");
           28   {
           29     unsigned int buffer_size = m_config->get_number("input.forcing.buffer_size");
           30     unsigned int evaluations_per_year = m_config->get_number("input.forcing.evaluations_per_year");
           31     bool periodic = opt.period > 0;
           32 
           33     File file(m_grid->com, opt.filename, PISM_NETCDF3, PISM_READONLY);
           34 
           35     m_retreat_mask = IceModelVec2T::ForcingField(m_grid,
           36                                                  file,
           37                                                  "land_ice_area_fraction_retreat",
           38                                                  "", // no standard name
           39                                                  buffer_size,
           40                                                  evaluations_per_year,
           41                                                  periodic);
           42     m_retreat_mask->set_attrs("forcing", "maximum ice extent mask",
           43                               "1", "1", "", 0);
           44   }
           45 }
           46 
           47 PrescribedRetreat::~PrescribedRetreat() {
           48   // empty
           49 }
           50 
           51 void PrescribedRetreat::init() {
           52 
           53   ForcingOptions opt(*m_grid->ctx(), "geometry.front_retreat.prescribed");
           54 
           55   m_log->message(2,
           56                  "* Initializing the prescribed front retreat mechanism\n"
           57                  "  using a time-dependent ice extent mask '%s' in '%s'...\n",
           58                  m_retreat_mask->get_name().c_str(),
           59                  opt.filename.c_str());
           60 
           61   m_retreat_mask->init(opt.filename, opt.period, opt.reference_time);
           62 }
           63 
           64 void PrescribedRetreat::update(double t,
           65                                double dt,
           66                                IceModelVec2S& ice_thickness,
           67                                IceModelVec2S& ice_area_specific_volume) {
           68   m_retreat_mask->update(t, dt);
           69   m_retreat_mask->average(t, dt);
           70 
           71   IceModelVec::AccessList list{m_retreat_mask.get(), &ice_thickness, &ice_area_specific_volume};
           72 
           73   for (Points p(*m_grid); p; p.next()) {
           74     const int i = p.i(), j = p.j();
           75 
           76     double f = (*m_retreat_mask)(i, j);
           77 
           78     if (f <= 0.0) {
           79       ice_area_specific_volume(i, j) = 0.0;
           80       ice_thickness(i, j)            = 0.0;
           81     } else if (f < 1.0) {
           82       ice_area_specific_volume(i, j) = ice_thickness(i, j) * f;
           83       ice_thickness(i, j)            = 0.0;
           84     } else {
           85       // M == 1.0: do nothing
           86     }
           87   }
           88 }
           89 
           90 MaxTimestep PrescribedRetreat::max_timestep_impl(double t) const {
           91   auto dt = m_retreat_mask->max_timestep(t);
           92 
           93   if (dt.finite()) {
           94     return MaxTimestep(dt.value(), "prescribed ice retreat");
           95   } else {
           96     return MaxTimestep("prescribed ice retreat");
           97   }
           98 }
           99 
          100 } // end of namespace pism