Constant.cc (3148B)
1 // Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018 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 "Constant.hh" 20 21 #include "pism/util/Vars.hh" 22 #include "pism/util/ConfigInterface.hh" 23 #include "pism/util/IceGrid.hh" 24 #include "pism/util/iceModelVec.hh" 25 #include "pism/util/MaxTimestep.hh" 26 #include "pism/util/pism_utilities.hh" 27 #include "pism/geometry/Geometry.hh" 28 29 namespace pism { 30 namespace ocean { 31 32 Constant::Constant(IceGrid::ConstPtr g) 33 : CompleteOceanModel(g) { 34 // empty 35 } 36 37 Constant::~Constant() { 38 // empty 39 } 40 41 void Constant::update_impl(const Geometry &geometry, double t, double dt) { 42 (void) t; 43 (void) dt; 44 45 const IceModelVec2S &ice_thickness = geometry.ice_thickness; 46 47 const double 48 melt_rate = m_config->get_number("ocean.constant.melt_rate", "m second-1"), 49 ice_density = m_config->get_number("constants.ice.density"), 50 mass_flux = melt_rate * ice_density; 51 52 melting_point_temperature(ice_thickness, *m_shelf_base_temperature); 53 54 m_shelf_base_mass_flux->set(mass_flux); 55 56 m_melange_back_pressure_fraction->set(m_config->get_number("ocean.melange_back_pressure_fraction")); 57 } 58 59 void Constant::init_impl(const Geometry &geometry) { 60 (void) geometry; 61 62 if (not m_config->get_flag("ocean.always_grounded")) { 63 m_log->message(2, "* Initializing the constant ocean model...\n"); 64 m_log->message(2, " Sub-shelf melt rate set to %f m/year.\n", 65 m_config->get_number("ocean.constant.melt_rate", "m year-1")); 66 67 } 68 } 69 70 MaxTimestep Constant::max_timestep_impl(double t) const { 71 (void) t; 72 return MaxTimestep("ocean constant"); 73 } 74 75 /*! 76 * Compute melting point temperature of ice at the depth `depth`. 77 */ 78 void Constant::melting_point_temperature(const IceModelVec2S& depth, 79 IceModelVec2S &result) const { 80 const double 81 T0 = m_config->get_number("constants.fresh_water.melting_point_temperature"), 82 beta_CC = m_config->get_number("constants.ice.beta_Clausius_Clapeyron"), 83 g = m_config->get_number("constants.standard_gravity"), 84 ice_density = m_config->get_number("constants.ice.density"); 85 86 IceModelVec::AccessList list{&depth, &result}; 87 88 for (Points p(*m_grid); p; p.next()) { 89 const int i = p.i(), j = p.j(); 90 const double pressure = ice_density * g * depth(i, j); // FIXME issue #15 91 92 result(i, j) = T0 - beta_CC * pressure; 93 } 94 } 95 96 } // end of namespape ocean 97 } // end of namespace pism