tGivenTH.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
---
tGivenTH.cc (24661B)
---
1 // Copyright (C) 2011, 2012, 2013, 2014, 2015, 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 2 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 <gsl/gsl_poly.h>
20 #include <cassert>
21
22 #include "GivenTH.hh"
23 #include "pism/util/IceGrid.hh"
24 #include "pism/util/Vars.hh"
25 #include "pism/util/Time.hh"
26 #include "pism/util/ConfigInterface.hh"
27 #include "pism/geometry/Geometry.hh"
28 #include "pism/coupler/util/options.hh"
29
30 namespace pism {
31 namespace ocean {
32
33 GivenTH::Constants::Constants(const Config &config) {
34 // coefficients of the in situ melting point temperature
35 // parameterization:
36 a[0] = -0.0575;
37 a[1] = 0.0901;
38 a[2] = -7.61e-4;
39 // coefficients of the in situ melting point potential temperature
40 // parameterization:
41 b[0] = -0.0575;
42 b[1] = 0.0921;
43 b[2] = -7.85e-4;
44
45 // turbulent heat transfer coefficient
46 gamma_T = 1.00e-4; // [m/s] RG3417 Default value from Hellmer and Olbers 89
47 // turbulent salt transfer coefficient
48 gamma_S = 5.05e-7; // [m/s] RG3417 Default value from Hellmer and Olbers 89
49
50 // FIXME: this should not be hard-wired. Eventually we should be able
51 // to use the spatially-variable top-of-the-ice temperature.
52 shelf_top_surface_temperature = -20.0; // degrees Celsius
53
54 water_latent_heat_fusion = config.get_number("constants.fresh_water.latent_heat_of_fusion");
55 sea_water_density = config.get_number("constants.sea_water.density");
56 sea_water_specific_heat_capacity = config.get_number("constants.sea_water.specific_heat_capacity");
57 ice_density = config.get_number("constants.ice.density");
58 ice_specific_heat_capacity = config.get_number("constants.ice.specific_heat_capacity");
59 ice_thermal_diffusivity = config.get_number("constants.ice.thermal_conductivity") / (ice_density * ice_specific_heat_capacity);
60 limit_salinity_range = config.get_flag("ocean.three_equation_model_clip_salinity");
61 }
62
63 GivenTH::GivenTH(IceGrid::ConstPtr g)
64 : CompleteOceanModel(g, std::shared_ptr<OceanModel>()) {
65
66 ForcingOptions opt(*m_grid->ctx(), "ocean.th");
67
68 {
69 unsigned int buffer_size = m_config->get_number("input.forcing.buffer_size");
70 unsigned int evaluations_per_year = m_config->get_number("input.forcing.evaluations_per_year");
71 bool periodic = opt.period > 0;
72
73 File file(m_grid->com, opt.filename, PISM_NETCDF3, PISM_READONLY);
74
75 m_theta_ocean = IceModelVec2T::ForcingField(m_grid,
76 file,
77 "theta_ocean",
78 "", // no standard name
79 buffer_size,
80 evaluations_per_year,
81 periodic,
82 LINEAR);
83
84 m_salinity_ocean = IceModelVec2T::ForcingField(m_grid,
85 file,
86 "salinity_ocean",
87 "", // no standard name
88 buffer_size,
89 evaluations_per_year,
90 periodic,
91 LINEAR);
92 }
93
94 m_theta_ocean->set_attrs("climate_forcing",
95 "potential temperature of the adjacent ocean",
96 "Kelvin", "Kelvin", "", 0);
97
98 m_salinity_ocean->set_attrs("climate_forcing",
99 "salinity of the adjacent ocean",
100 "g/kg", "g/kg", "", 0);
101 }
102
103 GivenTH::~GivenTH() {
104 // empty
105 }
106
107 void GivenTH::init_impl(const Geometry &geometry) {
108
109 m_log->message(2,
110 "* Initializing the 3eqn melting parameterization ocean model\n"
111 " reading ocean temperature and salinity from a file...\n");
112
113 ForcingOptions opt(*m_grid->ctx(), "ocean.th");
114
115 m_theta_ocean->init(opt.filename, opt.period, opt.reference_time);
116 m_salinity_ocean->init(opt.filename, opt.period, opt.reference_time);
117
118 // read time-independent data right away:
119 if (m_theta_ocean->n_records() == 1 && m_salinity_ocean->n_records() == 1) {
120 update(geometry, m_grid->ctx()->time()->current(), 0); // dt is irrelevant
121 }
122 }
123
124 void GivenTH::update_impl(const Geometry &geometry, double t, double dt) {
125 m_theta_ocean->update(t, dt);
126 m_salinity_ocean->update(t, dt);
127
128 m_theta_ocean->average(t, dt);
129 m_salinity_ocean->average(t, dt);
130
131 Constants c(*m_config);
132
133 const IceModelVec2S &ice_thickness = geometry.ice_thickness;
134
135 IceModelVec2S &temperature = *m_shelf_base_temperature;
136 IceModelVec2S &mass_flux = *m_shelf_base_mass_flux;
137
138 IceModelVec::AccessList list{ &ice_thickness, m_theta_ocean.get(), m_salinity_ocean.get(),
139 &temperature, &mass_flux};
140
141 for (Points p(*m_grid); p; p.next()) {
142 const int i = p.i(), j = p.j();
143
144 double potential_temperature_celsius = (*m_theta_ocean)(i,j) - 273.15;
145
146 double
147 shelf_base_temp_celsius = 0.0,
148 shelf_base_massflux = 0.0;
149
150 pointwise_update(c,
151 (*m_salinity_ocean)(i,j),
152 potential_temperature_celsius,
153 ice_thickness(i,j),
154 &shelf_base_temp_celsius,
155 &shelf_base_massflux);
156
157 // Convert from Celsius to Kelvin:
158 temperature(i,j) = shelf_base_temp_celsius + 273.15;
159 mass_flux(i,j) = shelf_base_massflux;
160 }
161
162 // convert mass flux from [m s-1] to [kg m-2 s-1]:
163 m_shelf_base_mass_flux->scale(m_config->get_number("constants.ice.density"));
164 }
165
166 MaxTimestep GivenTH::max_timestep_impl(double t) const {
167 (void) t;
168
169 return MaxTimestep("ocean th");
170 }
171
172 //* Evaluate the parameterization of the melting point temperature.
173 /** The value returned is in degrees Celsius.
174 */
175 static double melting_point_temperature(GivenTH::Constants c,
176 double salinity, double ice_thickness) {
177 return c.a[0] * salinity + c.a[1] + c.a[2] * ice_thickness;
178 }
179
180 /** Melt rate, obtained by solving the salt flux balance equation.
181 *
182 * @param c model constants
183 * @param sea_water_salinity sea water salinity
184 * @param basal_salinity shelf base salinity
185 *
186 * @return shelf base melt rate, in [m/s]
187 */
188 static double shelf_base_melt_rate(GivenTH::Constants c,
189 double sea_water_salinity, double basal_salinity) {
190
191 return c.gamma_S * c.sea_water_density * (sea_water_salinity - basal_salinity) / (c.ice_density * basal_salinity);
192 }
193
194 /** @brief Compute temperature and melt rate at the base of the shelf.
195 * Based on [@ref HellmerOlbers1989] and [@ref HollandJenkins1999].
196 *
197 * We use equations for the heat and salt flux balance at the base of
198 * the shelf to compute the temperature at the base of the shelf and
199 * the sub-shelf melt rate.
200 *
201 * @note The linearized equation for the freezing point of seawater as
202 * a function of salinity and pressure (ice thickness) is only valid
203 * for salinity ranges from 4 to 40 psu (see [@ref
204 * HollandJenkins1999]).
205 *
206 * Following [@ref HellmerOlbers1989], let @f$ Q_T @f$ be the total heat
207 * flux crossing the interface between the shelf base and the ocean,
208 * @f$ Q_T^B @f$ be the amount of heat lost by the ocean due to
209 * melting of glacial ice, and @f$ Q_T^I @f$ be the conductive flux
210 * into the ice column.
211 *
212 * @f$ Q_{T} @f$ is parameterized by (see [@ref HellmerOlbers1989], equation
213 * 10):
214 *
215 * @f[ Q_{T} = \rho_W\, c_{pW}\, \gamma_{T}\, (T^B - T^W),@f]
216 *
217 * where @f$ \rho_{W} @f$ is the sea water density, @f$ c_{pW} @f$ is
218 * the heat capacity of sea water, and @f$ \gamma_{T} @f$ is a
219 * turbulent heat exchange coefficient.
220 *
221 * We assume that the difference between the basal temperature and
222 * adjacent ocean temperature @f$ T^B - T^W @f$ is well approximated
223 * by @f$ \Theta_B - \Theta_W, @f$ where @f$ \Theta_{\cdot} @f$ is the
224 * corresponding potential temperature.
225 *
226 * @f$ Q_T^B @f$ is (see [@ref HellmerOlbers1989], equation 11):
227 *
228 * @f[ Q_T^B = \rho_I\, L\, \frac{\partial h}{\partial t}, @f]
229 *
230 * where @f$ \rho_I @f$ is the ice density, @f$ L @f$ is the latent
231 * heat of fusion, and @f$ \frac{\partial h}{\partial t} @f$ is the ice thickening rate
232 * (equal to minus the melt rate).
233 *
234 * The conductive flux into the ice column is ([@ref Hellmeretal1998],
235 * equation 7):
236 *
237 * @f[ Q_T^I = \rho_I\, c_{pI}\, \kappa\, T_{\text{grad}}, @f]
238 *
239 * where @f$ \rho_I @f$ is the ice density, @f$ c_{pI} @f$ is the heat
240 * capacity of ice, @f$ \kappa @f$ is the ice thermal diffusivity, and
241 * @f$ T_{\text{grad}} @f$ is the vertical temperature gradient at the
242 * base of a column of ice.
243 *
244 * Now, the heat flux balance implies
245 *
246 * @f[ Q_T = Q_T^B + Q_T^I. @f]
247 *
248 * For the salt flux balance, we have
249 *
250 * @f[ Q_S = Q_S^B + Q_S^I, @f]
251 *
252 * where @f$ Q_S @f$ is the total salt flux across the interface, @f$
253 * Q_S^B @f$ is the basal salt flux (negative for melting), @f$ Q_S^I = 0
254 * @f$ is the salt flux due to molecular diffusion of salt through
255 * ice.
256 *
257 * @f$ Q_S @f$ is parameterized by ([@ref Hellmeretal1998], equation 13)
258 *
259 * @f[ Q_S = \rho_W\, \gamma_S\, (S^B - S^W), @f]
260 *
261 * where @f$ \gamma_S @f$ is a turbulent salt exchange coefficient,
262 * @f$ S^B @f$ is salinity at the shelf base, and @f$ S^W @f$ is the
263 * salinity of adjacent ocean.
264 *
265 * The basal salt flux @f$ Q_S^B @f$ is ([@ref
266 * Hellmeretal1998], equation 10)
267 *
268 * @f[ Q_S^B = \rho_I\, S^B\, {\frac{\partial h}{\partial t}}. @f]
269 *
270 * To avoid converting shelf base temperature to shelf base potential
271 * temperature and back, we use two linearizations of the freezing point equation
272 * for sea water for in-situ and for potential temperature, respectively:
273 *
274 * @f[ T^{B}(S,h) = a_0\cdot S + a_1 + a_2\cdot h, @f]
275 *
276 * @f[ \Theta^{B}(S,h) = b_0\cdot S + b_1 + b_2\cdot h, @f]
277 *
278 * where @f$ S @f$ is salinity and @f$ h @f$ is ice shelf thickness.
279 *
280 * The linearization coefficients for the basal temperature @f$ T^B(S,h) @f$ are
281 * taken from [@ref Hellmeretal1998], going back to [@ref FoldvikKvinge1974].
282 *
283 * Given @f$ T^B(S,h) @f$ and a function @f$ \Theta_T^B(T)
284 * @f$ one can define @f$ \Theta^B_{*}(S,h) = \Theta_T^B\left(T^B(S,h)\right) @f$.
285 *
286 * The parameterization @f$ \Theta^B(S,h) @f$ used here was produced
287 * by linearizing @f$ \Theta^B_{*}(S,h) @f$ near the melting point.
288 * (The definition of @f$ \Theta_T^B(T) @f$, converting in situ
289 * temperature into potential temperature, was adopted from FESOM
290 * [@ref Wangetal2013]).
291 *
292 * Treating ice thickness, sea water salinity, and sea water potential
293 * temperature as "known" and choosing an approximation of the
294 * temperature gradient at the base @f$ T_{\text{grad}} @f$ (see
295 * subshelf_salinity_melt(), subshelf_salinity_freeze_on(),
296 * subshelf_salinity_diffusion_only()) we can write down a system of
297 * equations
298 *
299 * @f{align*}{
300 * Q_T &= Q_T^B + Q_T^I,\\
301 * Q_S &= Q_S^B + Q_S^I,\\
302 * T^{B}(S,h) &= a_0\cdot S + a_1 + a_2\cdot h,\\
303 * \Theta^{B}(S,h) &= b_0\cdot S + b_1 + b_2\cdot h\\
304 * @f}
305 *
306 * and simplify it to produce a quadratic equation for the salinity at the shelf base, @f$ S^B @f$:
307 *
308 * @f[ A\cdot (S^B)^2 + B\cdot S^B + C = 0 @f]
309 *
310 * The coefficients @f$ A, @f$ @f$ B, @f$ and @f$ C @f$ depend on the
311 * basal temperature gradient approximation for the sub-shelf melt,
312 * sub-shelf freeze-on, and diffusion-only cases.
313 *
314 * One remaining problem is that we cannot compute the basal melt rate
315 * without making an assumption about whether there is basal melt or
316 * not, and cannot pick one of the three cases without computing the
317 * basal melt rate first.
318 *
319 * This method tries to compute basal salinity that is consistent with
320 * the corresponding basal melt rate. See the code for details.
321 *
322 * Once @f$ S_B @f$ is found by solving this quadratic equation, we can
323 * compute the basal temperature using the parameterization for @f$
324 * T^{B}(S,h) @f$.
325 *
326 * To find the basal melt rate, we solve the salt flux balance
327 * equation for @f$ {\frac{\partial h}{\partial t}}, @f$ obtaining
328 *
329 * @f[ w_b = -\frac{\partial h}{\partial t} = \frac{\gamma_S\, \rho_W\, (S^W - S^B)}{\rho_I\, S^B}. @f]
330 *
331 *
332 * @param[in] constants model constants
333 * @param[in] sea_water_salinity sea water salinity
334 * @param[in] sea_water_potential_temperature sea water potential temperature
335 * @param[in] thickness ice shelf thickness
336 * @param[out] shelf_base_temperature_out resulting basal temperature
337 * @param[out] shelf_base_melt_rate_out resulting basal melt rate
338 *
339 * @return 0 on success
340 */
341 void GivenTH::pointwise_update(const Constants &constants,
342 double sea_water_salinity,
343 double sea_water_potential_temperature,
344 double thickness,
345 double *shelf_base_temperature_out,
346 double *shelf_base_melt_rate_out) {
347
348 assert(thickness >= 0.0);
349
350 // This model works for sea water salinity in the range of [4, 40]
351 // psu. Ensure that input salinity is in this range.
352 const double
353 min_salinity = 4.0,
354 max_salinity = 40.0;
355
356 if (constants.limit_salinity_range == true) {
357 if (sea_water_salinity < min_salinity) {
358 sea_water_salinity = min_salinity;
359 } else if (sea_water_salinity > max_salinity) {
360 sea_water_salinity = max_salinity;
361 }
362 }
363
364 double basal_salinity = sea_water_salinity;
365 subshelf_salinity(constants, sea_water_salinity, sea_water_potential_temperature,
366 thickness, &basal_salinity);
367
368 // Clip basal salinity so that we can use the freezing point
369 // temperature parameterization to recover shelf base temperature.
370 if (constants.limit_salinity_range == true) {
371 if (basal_salinity <= min_salinity) {
372 basal_salinity = min_salinity;
373 } else if (basal_salinity >= max_salinity) {
374 basal_salinity = max_salinity;
375 }
376 }
377
378 *shelf_base_temperature_out = melting_point_temperature(constants, basal_salinity, thickness);
379
380 *shelf_base_melt_rate_out = shelf_base_melt_rate(constants, sea_water_salinity, basal_salinity);
381
382 // no melt if there is no ice
383 if (thickness == 0.0) {
384 *shelf_base_melt_rate_out = 0.0;
385 }
386 }
387
388
389 /** @brief Compute the basal salinity and make sure that it is
390 * consistent with the basal melt rate.
391 *
392 * @param[in] c constants
393 * @param[in] sea_water_salinity sea water salinity
394 * @param[in] sea_water_potential_temperature sea water potential temperature
395 * @param[in] thickness ice shelf thickness
396 * @param[out] shelf_base_salinity resulting shelf base salinity
397 *
398 * @return 0 on success
399 */
400 void GivenTH::subshelf_salinity(const Constants &c,
401 double sea_water_salinity,
402 double sea_water_potential_temperature,
403 double thickness,
404 double *shelf_base_salinity) {
405 double basal_salinity = sea_water_salinity;
406
407 // first, assume that there is melt at the shelf base:
408 {
409 subshelf_salinity_melt(c, sea_water_salinity, sea_water_potential_temperature,
410 thickness, &basal_salinity);
411
412 double basal_melt_rate = shelf_base_melt_rate(c, sea_water_salinity, basal_salinity);
413
414 if (basal_melt_rate > 0.0) {
415 // computed basal melt rate is consistent with the assumption used
416 // to compute basal salinity
417 *shelf_base_salinity = basal_salinity;
418 return;
419 }
420 }
421
422 // Assuming that there is melt resulted in an inconsistent
423 // (salinity, melt_rate) pair. Assume that there is freeze-on at the base.
424 {
425 subshelf_salinity_freeze_on(c, sea_water_salinity, sea_water_potential_temperature,
426 thickness, &basal_salinity);
427
428 double basal_melt_rate = shelf_base_melt_rate(c, sea_water_salinity, basal_salinity);
429
430 if (basal_melt_rate < 0.0) {
431 // computed basal melt rate is consistent with the assumption
432 // used to compute basal salinity
433 *shelf_base_salinity = basal_salinity;
434 return;
435 }
436 }
437
438 // Both assumptions (above) resulted in inconsistencies. Revert to
439 // the "diffusion-only" case, which may be less accurate, but is
440 // generic and is always consistent.
441 {
442 subshelf_salinity_diffusion_only(c, sea_water_salinity, sea_water_potential_temperature,
443 thickness, &basal_salinity);
444
445 *shelf_base_salinity = basal_salinity;
446 }
447 }
448
449 /** Compute basal salinity in the basal melt case.
450 *
451 * We use the parameterization of the temperature gradient from [@ref
452 * Hellmeretal1998], equation 13:
453 *
454 * @f[ T_{\text{grad}} = -\Delta T\, \frac{\frac{\partial h}{\partial t}}{\kappa}, @f]
455 *
456 * where @f$ \Delta T @f$ is the difference between the ice
457 * temperature at the top of the ice column and its bottom:
458 * @f$ \Delta T = T^S - T^B. @f$ With this parameterization, we have
459 *
460 * @f[ Q_T^I = \rho_I\, c_{pI}\, {\frac{\partial h}{\partial t}}\, (T^S - T^B). @f]
461 *
462 * Then the coefficients of the quadratic equation for basal salinity
463 * (see pointwise_update()) are
464 *
465 * @f{align*}{
466 * A &= a_{0}\,\gamma_S\,c_{pI}-b_{0}\,\gamma_T\,c_{pW}\\
467 * B &= \gamma_S\,\left(L-c_{pI}\,\left(T^S+a_{0}\,S^W-a_{2}\,h-a_{1}\right)\right)+
468 * \gamma_T\,c_{pW}\,\left(\Theta^W-b_{2}\,h-b_{1}\right)\\
469 * C &= -\gamma_S\,S^W\,\left(L-c_{pI}\,\left(T^S-a_{2}\,h-a_{1}\right)\right)
470 * @f}
471 *
472 * @param[in] c physical constants, stored here to avoid looking them up in a double for loop
473 * @param[in] sea_water_salinity salinity of the ocean immediately adjacent to the shelf, [g/kg]
474 * @param[in] sea_water_potential_temperature potential temperature of the sea water, [degrees Celsius]
475 * @param[in] thickness thickness of the ice shelf, [meters]
476 * @param[out] shelf_base_salinity resulting shelf base salinity
477 *
478 * @return 0 on success
479 */
480 void GivenTH::subshelf_salinity_melt(const Constants &c,
481 double sea_water_salinity,
482 double sea_water_potential_temperature,
483 double thickness,
484 double *shelf_base_salinity) {
485
486 const double
487 c_pI = c.ice_specific_heat_capacity,
488 c_pW = c.sea_water_specific_heat_capacity,
489 L = c.water_latent_heat_fusion,
490 T_S = c.shelf_top_surface_temperature,
491 S_W = sea_water_salinity,
492 Theta_W = sea_water_potential_temperature;
493
494 // We solve a quadratic equation for Sb, the salinity at the shelf
495 // base.
496 //
497 // A*Sb^2 + B*Sb + C = 0
498 const double A = c.a[0] * c.gamma_S * c_pI - c.b[0] * c.gamma_T * c_pW;
499 const double B = (c.gamma_S * (L - c_pI * (T_S + c.a[0] * S_W - c.a[2] * thickness - c.a[1])) +
500 c.gamma_T * c_pW * (Theta_W - c.b[2] * thickness - c.b[1]));
501 const double C = -c.gamma_S * S_W * (L - c_pI * (T_S - c.a[2] * thickness - c.a[1]));
502
503 double S1 = 0.0, S2 = 0.0;
504 const int n_roots = gsl_poly_solve_quadratic(A, B, C, &S1, &S2);
505
506 assert(n_roots > 0);
507 assert(S2 > 0.0); // The bigger root should be positive.
508
509 *shelf_base_salinity = S2;
510 }
511
512 /** Compute basal salinity in the basal freeze-on case.
513 *
514 * In this case we assume that the temperature gradient at the shelf base is zero:
515 *
516 * @f[ T_{\text{grad}} = 0. @f]
517 *
518 * Please see pointwise_update() for details.
519 *
520 * In this case the coefficients of the quadratic equation for the
521 * basal salinity are:
522 *
523 * @f{align*}{
524 * A &= -b_{0}\,\gamma_T\,c_{pW} \\
525 * B &= \gamma_S\,L+\gamma_T\,c_{pW}\,\left(\Theta^W-b_{2}\,h-b_{1}\right) \\
526 * C &= -\gamma_S\,S^W\,L\\
527 * @f}
528 *
529 * @param[in] c model constants
530 * @param[in] sea_water_salinity sea water salinity
531 * @param[in] sea_water_potential_temperature sea water temperature
532 * @param[in] thickness ice shelf thickness
533 * @param[out] shelf_base_salinity resulting basal salinity
534 *
535 * @return 0 on success
536 */
537 void GivenTH::subshelf_salinity_freeze_on(const Constants &c,
538 double sea_water_salinity,
539 double sea_water_potential_temperature,
540 double thickness,
541 double *shelf_base_salinity) {
542
543 const double
544 c_pW = c.sea_water_specific_heat_capacity,
545 L = c.water_latent_heat_fusion,
546 S_W = sea_water_salinity,
547 Theta_W = sea_water_potential_temperature,
548 h = thickness;
549
550 // We solve a quadratic equation for Sb, the salinity at the shelf
551 // base.
552 //
553 // A*Sb^2 + B*Sb + C = 0
554 const double A = -c.b[0] * c.gamma_T * c_pW;
555 const double B = c.gamma_S * L + c.gamma_T * c_pW * (Theta_W - c.b[2] * h - c.b[1]);
556 const double C = -c.gamma_S * S_W * L;
557
558 double S1 = 0.0, S2 = 0.0;
559 const int n_roots = gsl_poly_solve_quadratic(A, B, C, &S1, &S2);
560
561 assert(n_roots > 0);
562 assert(S2 > 0.0); // The bigger root should be positive.
563
564 *shelf_base_salinity = S2;
565 }
566
567 /** @brief Compute basal salinity in the case of no basal melt and no
568 * freeze-on, with the diffusion-only temperature distribution in the
569 * ice column.
570 *
571 * In this case the temperature gradient at the base ([@ref
572 * HollandJenkins1999], equation 21) is
573 *
574 * @f[ T_{\text{grad}} = \frac{\Delta T}{h}, @f]
575 *
576 * where @f$ h @f$ is the ice shelf thickness and @f$ \Delta T = T^S -
577 * T^B @f$ is the difference between the temperature at the top and
578 * the bottom of the shelf.
579 *
580 * In this case the coefficients of the quadratic equation for the basal salinity are:
581 *
582 * @f{align*}{
583 * A &= - \frac{b_{0}\,\gamma_T\,h\,\rho_W\,c_{pW}-a_{0}\,\rho_I\,c_{pI}\,\kappa}{h\,\rho_W}\\
584 * B &= \frac{\rho_I\,c_{pI}\,\kappa\,\left(T^S-a_{2}\,h-a_{1}\right)}{h\,\rho_W}
585 +\gamma_S\,L+\gamma_T\,c_{pW}\,\left(\Theta^W-b_{2}\,h-b_{1}\right)\\
586 * C &= -\gamma_S\,S^W\,L\\
587 * @f}
588 *
589 * @param[in] c model constants
590 * @param[in] sea_water_salinity sea water salinity
591 * @param[in] sea_water_potential_temperature sea water potential temperature
592 * @param[in] thickness ice shelf thickness
593 * @param[out] shelf_base_salinity resulting basal salinity
594 *
595 * @return 0 on success
596 */
597 void GivenTH::subshelf_salinity_diffusion_only(const Constants &c,
598 double sea_water_salinity,
599 double sea_water_potential_temperature,
600 double thickness,
601 double *shelf_base_salinity) {
602 const double
603 c_pI = c.ice_specific_heat_capacity,
604 c_pW = c.sea_water_specific_heat_capacity,
605 L = c.water_latent_heat_fusion,
606 T_S = c.shelf_top_surface_temperature,
607 S_W = sea_water_salinity,
608 Theta_W = sea_water_potential_temperature,
609 h = thickness,
610 rho_W = c.sea_water_density,
611 rho_I = c.ice_density,
612 kappa = c.ice_thermal_diffusivity;
613
614 // We solve a quadratic equation for Sb, the salinity at the shelf
615 // base.
616 //
617 // A*Sb^2 + B*Sb + C = 0
618 const double A = -(c.b[0] * c.gamma_T * h * rho_W * c_pW - c.a[0] * rho_I * c_pI * kappa) / (h * rho_W);
619 const double B = ((rho_I * c_pI * kappa * (T_S - c.a[2] * h - c.a[1])) / (h * rho_W) +
620 c.gamma_S * L + c.gamma_T * c_pW * (Theta_W - c.b[2] * h - c.b[1]));
621 const double C = -c.gamma_S * S_W * L;
622
623 double S1 = 0.0, S2 = 0.0;
624 const int n_roots = gsl_poly_solve_quadratic(A, B, C, &S1, &S2);
625
626 assert(n_roots > 0);
627 assert(S2 > 0.0); // The bigger root should be positive.
628
629 *shelf_base_salinity = S2;
630 }
631
632 } // end of namespace ocean
633 } // end of namespace pism