tConstantPIK.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
---
tConstantPIK.cc (5156B)
---
1 // Copyright (C) 2008-2019 Ed Bueler, Constantine Khroulev, Ricarda Winkelmann,
2 // Gudfinna Adalgeirsdottir, Andy Aschwanden and Torsten Albrecht
3 //
4 // This file is part of PISM.
5 //
6 // PISM is free software; you can redistribute it and/or modify it under the
7 // terms of the GNU General Public License as published by the Free Software
8 // Foundation; either version 3 of the License, or (at your option) any later
9 // version.
10 //
11 // PISM is distributed in the hope that it will be useful, but WITHOUT ANY
12 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
13 // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
14 // details.
15 //
16 // You should have received a copy of the GNU General Public License
17 // along with PISM; if not, write to the Free Software
18 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
19
20 #include "ConstantPIK.hh"
21 #include "pism/util/Vars.hh"
22 #include "pism/util/ConfigInterface.hh"
23 #include "pism/util/IceGrid.hh"
24 #include "pism/util/iceModelVec.hh"
25 #include "pism/util/pism_options.hh"
26 #include "pism/util/io/io_helpers.hh"
27 #include "pism/util/MaxTimestep.hh"
28 #include "pism/util/pism_utilities.hh"
29 #include "pism/geometry/Geometry.hh"
30
31 namespace pism {
32 namespace ocean {
33
34 PIK::PIK(IceGrid::ConstPtr g)
35 : CompleteOceanModel(g) {
36 // empty
37 }
38
39 PIK::~PIK() {
40 // empty
41 }
42
43 void PIK::init_impl(const Geometry &geometry) {
44 (void) geometry;
45
46 m_log->message(2,
47 "* Initializing the constant (PIK) ocean model...\n");
48 }
49
50 MaxTimestep PIK::max_timestep_impl(double t) const {
51 (void) t;
52 return MaxTimestep("ocean PIK");
53 }
54
55 void PIK::update_impl(const Geometry &geometry, double t, double dt) {
56 (void) t;
57 (void) dt;
58
59 const IceModelVec2S &H = geometry.ice_thickness;
60
61 // Set shelf base temperature to the melting temperature at the base (depth within the
62 // ice equal to ice thickness).
63 melting_point_temperature(H, *m_shelf_base_temperature);
64
65 mass_flux(H, *m_shelf_base_mass_flux);
66
67 m_melange_back_pressure_fraction->set(m_config->get_number("ocean.melange_back_pressure_fraction"));
68 }
69
70 /*!
71 * Compute melting temperature at a given depth within the ice.
72 */
73 void PIK::melting_point_temperature(const IceModelVec2S &depth,
74 IceModelVec2S &result) const {
75 const double
76 T0 = m_config->get_number("constants.fresh_water.melting_point_temperature"), // K
77 beta_CC = m_config->get_number("constants.ice.beta_Clausius_Clapeyron"),
78 g = m_config->get_number("constants.standard_gravity"),
79 ice_density = m_config->get_number("constants.ice.density");
80
81 IceModelVec::AccessList list{&depth, &result};
82
83 for (Points p(*m_grid); p; p.next()) {
84 const int i = p.i(), j = p.j();
85 const double pressure = ice_density * g * depth(i,j); // FIXME task #7297
86 // result is set to melting point at depth
87 result(i,j) = T0 - beta_CC * pressure;
88 }
89 }
90
91 //! \brief Computes mass flux in [kg m-2 s-1].
92 /*!
93 * Assumes that mass flux is proportional to the shelf-base heat flux.
94 */
95 void PIK::mass_flux(const IceModelVec2S &ice_thickness, IceModelVec2S &result) const {
96 const double
97 melt_factor = m_config->get_number("ocean.pik_melt_factor"),
98 L = m_config->get_number("constants.fresh_water.latent_heat_of_fusion"),
99 sea_water_density = m_config->get_number("constants.sea_water.density"),
100 ice_density = m_config->get_number("constants.ice.density"),
101 c_p_ocean = 3974.0, // J/(K*kg), specific heat capacity of ocean mixed layer
102 gamma_T = 1e-4, // m/s, thermal exchange velocity
103 ocean_salinity = 35.0, // g/kg
104 T_ocean = units::convert(m_sys, -1.7, "Celsius", "Kelvin"); //Default in PISM-PIK
105
106 //FIXME: gamma_T should be a function of the friction velocity, not a const
107
108 IceModelVec::AccessList list{&ice_thickness, &result};
109
110 for (Points p(*m_grid); p; p.next()) {
111 const int i = p.i(), j = p.j();
112
113 // compute T_f(i, j) according to beckmann_goosse03, which has the
114 // meaning of the freezing temperature of the ocean water directly
115 // under the shelf, (of salinity 35psu) [this is related to the
116 // Pressure Melting Temperature, see beckmann_goosse03 eq. 2 for
117 // details]
118 double
119 shelfbaseelev = - (ice_density / sea_water_density) * ice_thickness(i,j),
120 T_f = 273.15 + (0.0939 -0.057 * ocean_salinity + 7.64e-4 * shelfbaseelev);
121 // add 273.15 to convert from Celsius to Kelvin
122
123 // compute ocean_heat_flux according to beckmann_goosse03
124 // positive, if T_oc > T_ice ==> heat flux FROM ocean TO ice
125 double ocean_heat_flux = melt_factor * sea_water_density * c_p_ocean * gamma_T * (T_ocean - T_f); // in W/m^2
126
127 // TODO: T_ocean -> field!
128
129 // shelfbmassflux is positive if ice is freezing on; here it is always negative:
130 // same sign as ocean_heat_flux (positive if massflux FROM ice TO ocean)
131 result(i,j) = ocean_heat_flux / (L * ice_density); // m s-1
132
133 // convert from [m s-1] to [kg m-2 s-1]:
134 result(i,j) *= ice_density;
135 }
136 }
137
138 } // end of namespace ocean
139 } // end of namespace pism