tPicoPhysics.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
---
tPicoPhysics.cc (7705B)
---
1 /* Copyright (C) 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 #include <cmath> // sqrt
20 #include <cassert> // assert
21 #include <algorithm> // std::min
22
23 #include "PicoPhysics.hh"
24
25 #include "pism/util/ConfigInterface.hh"
26
27 namespace pism {
28 namespace ocean {
29
30 PicoPhysics::PicoPhysics(const Config &config) {
31
32 // threshold between deep ocean and continental shelf
33 m_continental_shelf_depth = config.get_number("ocean.pico.continental_shelf_depth");
34
35 // value for ocean temperature around Antarctica if no other data available (cold
36 // conditions)
37 m_T_dummy = -1.5 + config.get_number("constants.fresh_water.melting_point_temperature");
38
39 // value for ocean salinity around Antarctica if no other data available (cold
40 // conditions)
41 m_S_dummy = 34.7;
42
43 m_earth_grav = config.get_number("constants.standard_gravity");
44 m_ice_density = config.get_number("constants.ice.density");
45 m_sea_water_density = config.get_number("constants.sea_water.density");
46 m_rho_star = 1033; // kg/m^3
47 m_nu = m_ice_density / m_sea_water_density; // no unit
48
49 //Joule / kg
50 m_latentHeat = config.get_number("constants.fresh_water.latent_heat_of_fusion");
51
52 // J/(K*kg), specific heat capacity of ocean mixed layer
53 m_c_p_ocean = 3974.0;
54
55 m_lambda = m_latentHeat / m_c_p_ocean; // °C, NOTE K vs °C
56
57 // Values for linearized potential freezing point (from Xylar Asay-Davis, should be in
58 // Asay-Davis et al 2016, but not correct in there )
59 m_a_pot = -0.0572; // K/psu
60 m_b_pot = 0.0788 + 273.15; // K
61 m_c_pot = 7.77e-4; // K/dbar
62
63 // in-situ pressure melting point from Jenkins et al. 2010 paper
64 m_a_in_situ = -0.0573; // K/p_in_situu
65 m_b_in_situ = 0.0832 + 273.15; // K
66 m_c_in_situ = 7.53e-4; // K/dbar
67
68 // in-situ pressure melting point from Olbers & Hellmer 2010 paper
69 // as = -0.057; // K/psu
70 // bs = 0.0832 + 273.15; // K
71 // cs = 7.64e-4; // K/dbar
72
73 m_alpha = 7.5e-5; // 1/K
74 m_beta = 7.7e-4; // 1/psu
75
76 // m s-1, best fit value in paper
77 m_gamma_T = config.get_number("ocean.pico.heat_exchange_coefficent");
78
79 // m6 kg−1 s−1, best fit value in paper
80 m_overturning_coeff = config.get_number("ocean.pico.overturning_coefficent");
81
82 // for shelf cells where normal box model is not calculated, used in
83 // calculate_basal_melt_missing_cells(), compare POConstantPIK m/s, thermal exchange
84 // velocity for Beckmann-Goosse parameterization this is the same meltFactor as in
85 // POConstantPIK
86 m_meltFactor = config.get_number("ocean.pik_melt_factor");
87 }
88
89
90 double PicoPhysics::pressure(double ice_thickness) const {
91 // pressure in dbar, 1dbar = 10000 Pa = 1e4 kg m-1 s-2
92 return m_ice_density * m_earth_grav * ice_thickness * 1e-4;
93 }
94
95 TocBox1 PicoPhysics::Toc_box1(double area, double T_star, double Soc_box0, double Toc_box0) const {
96
97 TocBox1 result = { false, 0.0 };
98
99 const double
100 g1 = area * m_gamma_T,
101 s1 = Soc_box0 / (m_nu * m_lambda),
102 p = p_coeff(g1, s1),
103 q = p * T_star;
104
105 // This can only happen if T_star > 0.25*p, in particular T_star > 0 which can only
106 // happen for values of Toc_box0 close to the local pressure melting point
107 double D = 0.25 * p * p - q;
108 if (D < 0.0) {
109 D = 0.0;
110 result.failed = true;
111 }
112
113 // temperature for box 1, p-q formula
114 // equation A12 in the PICO paper.
115 result.value = Toc_box0 - (-0.5 * p + sqrt(D));
116
117 return result;
118 }
119
120 double PicoPhysics::Toc(double area, double temperature, double T_star, double overturning, double salinity) const {
121
122 double g1 = area * m_gamma_T;
123 double g2 = g1 / (m_nu * m_lambda);
124
125 // temperature for Box i > 1
126 return temperature + g1 * T_star / (overturning + g1 - g2 * m_a_pot * salinity); // K
127 }
128
129 double PicoPhysics::Soc_box1(double Toc_box0, double Soc_box0, double Toc) const {
130
131 return Soc_box0 - (Soc_box0 / (m_nu * m_lambda)) * (Toc_box0 - Toc);
132 }
133
134 double PicoPhysics::Soc(double salinity, double temperature_in_boundary, double Toc) const {
135
136 return salinity - salinity * (temperature_in_boundary - Toc) / (m_nu * m_lambda); // psu;
137 }
138
139
140 //! equation 5 in the PICO paper.
141 //! calculate pressure melting point from potential temperature
142 double PicoPhysics::theta_pm(double salinity, double pressure) const {
143 // using coefficients for potential temperature
144 return m_a_pot * salinity + m_b_pot - m_c_pot * pressure;
145 }
146
147 //! equation 5 in the PICO paper.
148 //! calculate pressure melting point from in-situ temperature
149 double PicoPhysics::T_pm(double salinity, double pressure) const {
150 // using coefficients for in-situ temperature
151 return m_a_in_situ * salinity + m_b_in_situ - m_c_in_situ * pressure;
152 }
153
154 //! equation 8 in the PICO paper.
155 double PicoPhysics::melt_rate(double pm_point, double Toc) const {
156 // in m/s
157 return m_gamma_T / (m_nu * m_lambda) * (Toc - pm_point);
158 }
159
160 //! Beckmann & Goosse meltrate
161 double PicoPhysics::melt_rate_beckmann_goosse(double pot_pm_point, double Toc) const {
162 // in W/m^2
163 double heat_flux = m_meltFactor * m_sea_water_density * m_c_p_ocean * m_gamma_T * (Toc - pot_pm_point);
164 // in m s-1
165 return heat_flux / (m_latentHeat * m_ice_density);
166 }
167
168 //! equation 3 in the PICO paper. See also equation 4.
169 double PicoPhysics::overturning(double Soc_box0, double Soc, double Toc_box0, double Toc) const {
170 // in m^3/s
171 return m_overturning_coeff * m_rho_star * (m_beta * (Soc_box0 - Soc) - m_alpha * (Toc_box0 - Toc));
172 }
173
174 //! See equation A6 and lines before in PICO paper
175 double PicoPhysics::T_star(double salinity, double temperature, double pressure) const {
176 double T_s = theta_pm(salinity, pressure) - temperature; // in Kelvin
177
178 // Positive values are unphysical as temperatures below the pressure melting point would
179 // be ice, but we are in the ocean. This should not occur because
180 // set_ocean_input_fields(...) sets too cold temperatures to pressure melting point +
181 // 0.001
182 return std::min(T_s, 0.0);
183 }
184
185 //! calculate p coefficent for solving the quadratic temperature equation
186 //! trough the p-q formula. See equation A12 in the PICO paper.
187 //! is only used once in Toc_box1(...)
188 double PicoPhysics::p_coeff(double g1, double s1) const {
189 // in 1 / (1/K) = K
190 // inputs g1 and overturning coefficicient are always positive
191 // so output is positive if beta*s1 > alpha
192 // which is shown in the text following equation A12
193 return g1 / (m_overturning_coeff * m_rho_star * (m_beta * s1 - m_alpha));
194 }
195
196 double PicoPhysics::gamma_T() const {
197 return m_gamma_T;
198 }
199
200 double PicoPhysics::overturning_coeff() const {
201 return m_overturning_coeff;
202 }
203
204 double PicoPhysics::T_dummy() const {
205 return m_T_dummy;
206 }
207
208 double PicoPhysics::S_dummy() const {
209 return m_S_dummy;
210 }
211
212 double PicoPhysics::ice_density() const {
213 return m_ice_density;
214 }
215
216 double PicoPhysics::continental_shelf_depth() const {
217 return m_continental_shelf_depth;
218 }
219
220 } // end of namespace ocean
221 } // end of namespace pism