tHayhurstCalving.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
---
tHayhurstCalving.cc (6094B)
---
1 /* Copyright (C) 2016, 2017, 2018, 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 "HayhurstCalving.hh"
21
22 #include "pism/util/IceGrid.hh"
23 #include "pism/util/error_handling.hh"
24 #include "pism/util/IceModelVec2CellType.hh"
25 #include "pism/geometry/Geometry.hh"
26
27 namespace pism {
28 namespace calving {
29
30 HayhurstCalving::HayhurstCalving(IceGrid::ConstPtr grid)
31 : Component(grid),
32 m_calving_rate(grid, "hayhurst_calving_rate", WITH_GHOSTS)
33 {
34 m_calving_rate.set_attrs("diagnostic",
35 "horizontal calving rate due to Hayhurst calving",
36 "m s-1", "m day-1", "", 0);
37
38 }
39
40 HayhurstCalving::~HayhurstCalving() {
41 // empty
42 }
43
44 void HayhurstCalving::init() {
45
46 m_log->message(2,
47 "* Initializing the 'Hayhurst calving' mechanism...\n");
48
49 m_B_tilde = m_config->get_number("calving.hayhurst_calving.B_tilde");
50 m_exponent_r = m_config->get_number("calving.hayhurst_calving.exponent_r");
51 m_sigma_threshold = m_config->get_number("calving.hayhurst_calving.sigma_threshold", "Pa");
52
53 m_log->message(2,
54 " B tilde parameter: %3.3f MPa-%3.3f yr-1.\n", m_B_tilde, m_exponent_r);
55 m_log->message(2,
56 " Hayhurst calving threshold: %3.3f MPa.\n",
57 convert(m_sys, m_sigma_threshold, "Pa", "MPa"));
58
59 if (fabs(m_grid->dx() - m_grid->dy()) / std::min(m_grid->dx(), m_grid->dy()) > 1e-2) {
60 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
61 "-calving hayhurst_calving using a non-square grid cell is not implemented (yet);\n"
62 "dx = %f, dy = %f, relative difference = %f",
63 m_grid->dx(), m_grid->dy(),
64 fabs(m_grid->dx() - m_grid->dy()) / std::max(m_grid->dx(), m_grid->dy()));
65 }
66
67 }
68
69 void HayhurstCalving::update(const IceModelVec2CellType &cell_type,
70 const IceModelVec2S &ice_thickness,
71 const IceModelVec2S &sea_level,
72 const IceModelVec2S &bed_elevation) {
73
74 using std::min;
75
76 const double
77 ice_density = m_config->get_number("constants.ice.density"),
78 water_density = m_config->get_number("constants.sea_water.density"),
79 gravity = m_config->get_number("constants.standard_gravity"),
80 // convert "Pa" to "MPa" and "m yr-1" to "m s-1"
81 unit_scaling = pow(1e-6, m_exponent_r) * convert(m_sys, 1.0, "m year-1", "m second-1");
82
83 IceModelVec::AccessList list{&ice_thickness, &cell_type, &m_calving_rate, &sea_level,
84 &bed_elevation};
85
86 for (Points pt(*m_grid); pt; pt.next()) {
87 const int i = pt.i(), j = pt.j();
88
89 double water_depth = sea_level(i, j) - bed_elevation(i, j);
90
91 if (cell_type.icy(i, j) and water_depth > 0.0) {
92 // note that ice_thickness > 0 at icy locations
93 assert(ice_thickness(i, j) > 0);
94
95 double H = ice_thickness(i, j);
96
97 // Note that for ice at floatation water_depth = H * (ice_density / water_density),
98 // so omega cannot exceed ice_density / water_density.
99 double omega = water_depth / H;
100
101 // Extend the calving parameterization to ice shelves. This tweak should (I hope)
102 // prevent a feedback in which the creation of an ice shelf reduces the calving
103 // rate, which leads to an advance of the front and an even lower calving rate, etc.
104 if (omega > ice_density / water_density) {
105 // ice at the front is floating
106 double freeboard = (1.0 - ice_density / water_density) * H;
107 H = water_depth + freeboard;
108 omega = water_depth / H;
109 }
110
111 // [\ref Mercenier2018] maximum tensile stress approximation
112 double sigma_0 = (0.4 - 0.45 * pow(omega - 0.065, 2.0)) * ice_density * gravity * H;
113
114 // ensure that sigma_0 - m_sigma_threshold >= 0
115 sigma_0 = std::max(sigma_0, m_sigma_threshold);
116
117 // [\ref Mercenier2018] equation 22
118 m_calving_rate(i, j) = (m_B_tilde * unit_scaling *
119 (1.0 - pow(omega, 2.8)) *
120 pow(sigma_0 - m_sigma_threshold, m_exponent_r) * H);
121 } else { // end of "if (ice_free_ocean and next_to_floating)"
122 m_calving_rate(i, j) = 0.0;
123 }
124 } // end of loop over grid points
125
126 // Set calving rate *near* grounded termini to the average of grounded icy
127 // neighbors: front retreat code uses values at these locations (the rest is for
128 // visualization).
129
130 m_calving_rate.update_ghosts();
131
132 const Direction dirs[] = {North, East, South, West};
133
134 for (Points p(*m_grid); p; p.next()) {
135 const int i = p.i(), j = p.j();
136
137 if (cell_type.ice_free(i, j) and cell_type.next_to_ice(i, j) ) {
138
139 auto R = m_calving_rate.star(i, j);
140 auto M = cell_type.int_star(i, j);
141
142 int N = 0;
143 double R_sum = 0.0;
144 for (int n = 0; n < 4; ++n) {
145 Direction direction = dirs[n];
146 if (mask::icy(M[direction])) {
147 R_sum += R[direction];
148 N++;
149 }
150 }
151
152 if (N > 0) {
153 m_calving_rate(i, j) = R_sum / N;
154 }
155 }
156 }
157 }
158
159 const IceModelVec2S &HayhurstCalving::calving_rate() const {
160 return m_calving_rate;
161 }
162
163 DiagnosticList HayhurstCalving::diagnostics_impl() const {
164 return {{"hayhurst_calving_rate", Diagnostic::wrap(m_calving_rate)}};
165 }
166
167 } // end of namespace calving
168 } // end of namespace pism