tDischargeGiven.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
---
tDischargeGiven.cc (7676B)
---
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 "DischargeGiven.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 DischargeGiven::DischargeGiven(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 " UAF-UT\n");
37
38 unsigned int evaluations_per_year = m_config->get_number("input.forcing.evaluations_per_year");
39
40 m_theta_ocean.reset(new IceModelVec2T(grid, "theta_ocean", 1, evaluations_per_year));
41 m_theta_ocean->set_attrs("climate_forcing",
42 "potential temperature of the adjacent ocean",
43 "Celsius", "Celsius", "", 0);
44
45 m_theta_ocean->init_constant(0.0);
46
47 m_subglacial_discharge.reset(new IceModelVec2T(grid,
48 "subglacial_discharge", 1,
49 evaluations_per_year));
50 m_subglacial_discharge->set_attrs("climate_forcing",
51 "subglacial discharge",
52 "kg m-2 s-1", "kg m-2 year-1", "", 0);
53
54 m_subglacial_discharge->init_constant(0.0);
55 }
56
57 DischargeGiven::~DischargeGiven() {
58 // empty
59 }
60
61 void DischargeGiven::init_impl(const Geometry &geometry) {
62 (void) geometry;
63
64 ForcingOptions opt(*m_grid->ctx(), "frontal_melt.discharge_given");
65
66 {
67 unsigned int buffer_size = m_config->get_number("input.forcing.buffer_size");
68 unsigned int evaluations_per_year = m_config->get_number("input.forcing.evaluations_per_year");
69 bool periodic = opt.period > 0;
70
71 File file(m_grid->com, opt.filename, PISM_NETCDF3, PISM_READONLY);
72
73 m_theta_ocean = IceModelVec2T::ForcingField(m_grid,
74 file,
75 "theta_ocean",
76 "", // no standard name
77 buffer_size,
78 evaluations_per_year,
79 periodic);
80
81 m_subglacial_discharge = IceModelVec2T::ForcingField(m_grid,
82 file,
83 "subglacial_discharge",
84 "", // no standard name
85 buffer_size,
86 evaluations_per_year,
87 periodic);
88 }
89
90 m_theta_ocean->set_attrs("climate_forcing",
91 "potential temperature of the adjacent ocean",
92 "Celsius", "Celsius", "", 0);
93
94 m_theta_ocean->init(opt.filename, opt.period, opt.reference_time);
95
96 m_subglacial_discharge->set_attrs("climate_forcing",
97 "subglacial discharge",
98 "kg m-2 s-1", "kg m-2 year-1", "", 0);
99
100 m_subglacial_discharge->init(opt.filename, opt.period, opt.reference_time);
101 }
102
103 /*!
104 * Initialize potential temperature from IceModelVecs instead of an input
105 * file (for testing).
106 */
107 void DischargeGiven::initialize(const IceModelVec2S &theta, const IceModelVec2S &sgl) {
108 m_theta_ocean->copy_from(theta);
109 m_subglacial_discharge->copy_from(sgl);
110 }
111
112 void DischargeGiven::update_impl(const FrontalMeltInputs &inputs, double t, double dt) {
113
114 m_theta_ocean->update(t, dt);
115 m_subglacial_discharge->update(t, dt);
116
117 m_theta_ocean->average(t, dt);
118 m_subglacial_discharge->average(t, dt);
119
120 FrontalMeltPhysics physics(*m_config);
121
122 const IceModelVec2CellType &cell_type = inputs.geometry->cell_type;
123 const IceModelVec2S &bed_elevation = inputs.geometry->bed_elevation;
124 const IceModelVec2S &sea_level_elevation = inputs.geometry->sea_level_elevation;
125
126 IceModelVec::AccessList list
127 {&bed_elevation, &cell_type, &sea_level_elevation,
128 m_theta_ocean.get(), m_subglacial_discharge.get(), m_frontal_melt_rate.get()};
129
130 double
131 water_density = m_config->get_number("constants.fresh_water.density"),
132 seconds_per_day = 86400;
133
134 for (Points p(*m_grid); p; p.next()) {
135 const int i = p.i(), j = p.j();
136
137 if (cell_type.icy(i, j)) {
138 // Assume for now that thermal forcing is equal to theta_ocean. Also, thermal
139 // forcing is generally not available at the grounding line.
140 double TF = (*m_theta_ocean)(i, j);
141
142 // Convert subglacial discharge (kg /(m^2 * s)) to an "effective subglacial freshwater
143 // velocity" or volume flux per unit area of ice front in m/day (see Xu et al 2013,
144 // section 2, paragraph 11).
145 //
146 // [flux] = kg / (m^2 * s), so
147 // [flux / water_density] = m / s, so
148 // [flux / water_density * (s / day) ] = m / day
149 double q_sg = ((*m_subglacial_discharge)(i, j) / water_density) * seconds_per_day;
150
151 double water_depth = sea_level_elevation(i, j) - bed_elevation(i, j);
152
153 (*m_frontal_melt_rate)(i, j) = physics.frontal_melt_from_undercutting(water_depth, q_sg, TF);
154 // convert from m / day to m / s
155 (*m_frontal_melt_rate)(i, j) /= seconds_per_day;
156 } else {
157 (*m_frontal_melt_rate)(i, j) = 0.0;
158 }
159 } // end of the loop over grid points
160
161 // Set frontal melt rate *near* grounded termini to the average of grounded icy
162 // neighbors: front retreat code uses values at these locations (the rest is for
163 // visualization).
164
165 m_frontal_melt_rate->update_ghosts();
166
167 const Direction dirs[] = {North, East, South, West};
168
169 for (Points p(*m_grid); p; p.next()) {
170 const int i = p.i(), j = p.j();
171
172 if (apply(cell_type, i, j) and cell_type.ice_free(i, j)) {
173
174 auto R = m_frontal_melt_rate->star(i, j);
175 auto M = cell_type.int_star(i, j);
176
177 int N = 0;
178 double R_sum = 0.0;
179 for (int n = 0; n < 4; ++n) {
180 Direction direction = dirs[n];
181 if (mask::grounded_ice(M[direction]) or
182 (m_include_floating_ice and mask::icy(M[direction]))) {
183 R_sum += R[direction];
184 N++;
185 }
186 }
187
188 if (N > 0) {
189 (*m_frontal_melt_rate)(i, j) = R_sum / N;
190 }
191 }
192 }
193 }
194
195 const IceModelVec2S& DischargeGiven::frontal_melt_rate_impl() const {
196 return *m_frontal_melt_rate;
197 }
198
199 MaxTimestep DischargeGiven::max_timestep_impl(double t) const {
200
201 auto dt = std::min(m_theta_ocean->max_timestep(t),
202 m_subglacial_discharge->max_timestep(t));
203
204 if (dt.finite()) {
205 return MaxTimestep(dt.value(), "frontal_melt discharge_given");
206 } else {
207 return MaxTimestep("frontal_melt discharge_given");
208 }
209 }
210
211 } // end of namespace frontalmelt
212 } // end of namespace pism