tDischargeRouting.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch) (HTM) git clone git://src.adamsgaard.dk/pism (DIR) Log (DIR) Files (DIR) Refs (DIR) LICENSE --- tDischargeRouting.cc (6723B) --- 1 // Copyright (C) 2018, 2019 Andy Aschwanden and Constantine Khroulev 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 "DischargeRouting.hh" 20 21 #include "pism/util/IceGrid.hh" 22 #include "pism/geometry/Geometry.hh" 23 #include "pism/coupler/util/options.hh" 24 #include "FrontalMeltPhysics.hh" 25 26 namespace pism { 27 namespace frontalmelt { 28 29 DischargeRouting::DischargeRouting(IceGrid::ConstPtr grid) 30 : FrontalMelt(grid, nullptr) { 31 32 m_frontal_melt_rate = allocate_frontal_melt_rate(grid, 1); 33 34 m_log->message(2, 35 "* Initializing the frontal melt model\n" 36 " using the Rignot/Xu parameterization\n" 37 " and routing of subglacial discharge\n"); 38 39 unsigned int evaluations_per_year = m_config->get_number("input.forcing.evaluations_per_year"); 40 41 m_theta_ocean.reset(new IceModelVec2T(grid, "theta_ocean", 1, evaluations_per_year, LINEAR)); 42 m_theta_ocean->set_attrs("climate_forcing", 43 "potential temperature of the adjacent ocean", 44 "Celsius", "Celsius", "", 0); 45 46 m_theta_ocean->init_constant(0.0); 47 } 48 49 DischargeRouting::~DischargeRouting() { 50 // empty 51 } 52 53 void DischargeRouting::init_impl(const Geometry &geometry) { 54 (void) geometry; 55 56 ForcingOptions opt(*m_grid->ctx(), "frontal_melt.routing"); 57 58 { 59 unsigned int buffer_size = m_config->get_number("input.forcing.buffer_size"); 60 unsigned int evaluations_per_year = m_config->get_number("input.forcing.evaluations_per_year"); 61 bool periodic = opt.period > 0; 62 63 File file(m_grid->com, opt.filename, PISM_NETCDF3, PISM_READONLY); 64 65 m_theta_ocean = IceModelVec2T::ForcingField(m_grid, 66 file, 67 "theta_ocean", 68 "", // no standard name 69 buffer_size, 70 evaluations_per_year, 71 periodic, 72 LINEAR); 73 } 74 75 m_theta_ocean->set_attrs("climate_forcing", 76 "potential temperature of the adjacent ocean", 77 "Celsius", "Celsius", "", 0); 78 79 m_theta_ocean->init(opt.filename, opt.period, opt.reference_time); 80 } 81 82 /*! 83 * Initialize potential temperature from IceModelVecs instead of an input 84 * file (for testing). 85 */ 86 void DischargeRouting::initialize(const IceModelVec2S &theta) { 87 m_theta_ocean->copy_from(theta); 88 } 89 90 void DischargeRouting::update_impl(const FrontalMeltInputs &inputs, double t, double dt) { 91 92 m_theta_ocean->update(t, dt); 93 94 FrontalMeltPhysics physics(*m_config); 95 96 const IceModelVec2CellType &cell_type = inputs.geometry->cell_type; 97 const IceModelVec2S &bed_elevation = inputs.geometry->bed_elevation; 98 const IceModelVec2S &ice_thickness = inputs.geometry->ice_thickness; 99 const IceModelVec2S &sea_level_elevation = inputs.geometry->sea_level_elevation; 100 const IceModelVec2S &water_flux = *inputs.subglacial_water_flux; 101 102 IceModelVec::AccessList list 103 {&ice_thickness, &bed_elevation, &cell_type, &sea_level_elevation, 104 &water_flux, m_theta_ocean.get(), m_frontal_melt_rate.get()}; 105 106 double 107 seconds_per_day = 86400, 108 grid_spacing = 0.5 * (m_grid->dx() + m_grid->dy()); 109 110 for (Points p(*m_grid); p; p.next()) { 111 const int i = p.i(), j = p.j(); 112 113 if (cell_type.icy(i, j)) { 114 // Assume for now that thermal forcing is equal to theta_ocean. Also, thermal 115 // forcing is generally not available at the grounding line. 116 double TF = (*m_theta_ocean)(i, j); 117 118 double water_depth = std::max(sea_level_elevation(i, j) - bed_elevation(i, j), 0.0), 119 submerged_front_area = water_depth * grid_spacing; 120 121 // Convert subglacial water flux (m^2/s) to an "effective subglacial freshwater 122 // velocity" or flux per unit area of ice front in m/day (see Xu et al 2013, section 123 // 2, paragraph 11). 124 // 125 // [flux] = m^2 / s, so 126 // [flux * grid_spacing] = m^3 / s, so 127 // [flux * grid_spacing / submerged_front_area] = m / s, and 128 // [flux * grid_spacing * (s / day) / submerged_front_area] = m / day 129 double Q_sg = water_flux(i, j) * grid_spacing; 130 double q_sg = Q_sg / submerged_front_area * seconds_per_day; 131 132 (*m_frontal_melt_rate)(i, j) = physics.frontal_melt_from_undercutting(water_depth, q_sg, TF); 133 // convert from m / day to m / s 134 (*m_frontal_melt_rate)(i, j) /= seconds_per_day; 135 } else { 136 (*m_frontal_melt_rate)(i, j) = 0.0; 137 } 138 } // end of the loop over grid points 139 140 // Set frontal melt rate *near* grounded termini to the average of grounded icy 141 // neighbors: front retreat code uses values at these locations (the rest is for 142 // visualization). 143 144 m_frontal_melt_rate->update_ghosts(); 145 146 const Direction dirs[] = {North, East, South, West}; 147 148 for (Points p(*m_grid); p; p.next()) { 149 const int i = p.i(), j = p.j(); 150 151 if (apply(cell_type, i, j) and cell_type.ice_free(i, j)) { 152 153 auto R = m_frontal_melt_rate->star(i, j); 154 auto M = cell_type.int_star(i, j); 155 156 int N = 0; 157 double R_sum = 0.0; 158 for (int n = 0; n < 4; ++n) { 159 Direction direction = dirs[n]; 160 if (mask::grounded_ice(M[direction]) or 161 (m_include_floating_ice and mask::icy(M[direction]))) { 162 R_sum += R[direction]; 163 N++; 164 } 165 } 166 167 if (N > 0) { 168 (*m_frontal_melt_rate)(i, j) = R_sum / N; 169 } 170 } 171 } 172 } 173 174 const IceModelVec2S& DischargeRouting::frontal_melt_rate_impl() const { 175 return *m_frontal_melt_rate; 176 } 177 178 MaxTimestep DischargeRouting::max_timestep_impl(double t) const { 179 180 auto dt = m_theta_ocean->max_timestep(t); 181 182 if (dt.finite()) { 183 return MaxTimestep(dt.value(), "frontal_melt routing"); 184 } else { 185 return MaxTimestep("frontal_melt routing"); 186 } 187 } 188 189 } // end of namespace frontalmelt 190 } // end of namespace pism