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