PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
ConstantPIK.cc
Go to the documentation of this file.
1 // Copyright (C) 2008-2019, 2021, 2022, 2023 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 "pism/coupler/ocean/ConstantPIK.hh"
21 #include "pism/util/ConfigInterface.hh"
22 #include "pism/util/Grid.hh"
23 #include "pism/util/MaxTimestep.hh"
24 #include "pism/geometry/Geometry.hh"
25 
26 namespace pism {
27 namespace ocean {
28 
29 PIK::PIK(std::shared_ptr<const Grid> g)
31  // empty
32 }
33 
34 void PIK::init_impl(const Geometry &geometry) {
35  (void) geometry;
36 
37  m_log->message(2,
38  "* Initializing the constant (PIK) ocean model...\n");
39 
40  double
41  ice_density = m_config->get_number("constants.ice.density"),
42  water_density = m_config->get_number("constants.sea_water.density"),
43  g = m_config->get_number("constants.standard_gravity");
44 
45  compute_average_water_column_pressure(geometry, ice_density, water_density, g,
47 }
48 
50  (void) t;
51  return MaxTimestep("ocean PIK");
52 }
53 
54 void PIK::update_impl(const Geometry &geometry, double t, double dt) {
55  (void) t;
56  (void) dt;
57 
58  const array::Scalar &H = geometry.ice_thickness;
59 
60  // Set shelf base temperature to the melting temperature at the base (depth within the
61  // ice equal to ice thickness).
63 
65 
66  double
67  ice_density = m_config->get_number("constants.ice.density"),
68  water_density = m_config->get_number("constants.sea_water.density"),
69  g = m_config->get_number("constants.standard_gravity");
70 
71  compute_average_water_column_pressure(geometry, ice_density, water_density, g,
73 }
74 
75 /*!
76  * Compute melting temperature at a given depth within the ice.
77  */
79  array::Scalar &result) const {
80  const double
81  T0 = m_config->get_number("constants.fresh_water.melting_point_temperature"), // K
82  beta_CC = m_config->get_number("constants.ice.beta_Clausius_Clapeyron"),
83  g = m_config->get_number("constants.standard_gravity"),
84  ice_density = m_config->get_number("constants.ice.density");
85 
86  array::AccessScope list{&depth, &result};
87 
88  for (auto p = m_grid->points(); p; p.next()) {
89  const int i = p.i(), j = p.j();
90  const double pressure = ice_density * g * depth(i,j); // FIXME task #7297
91  // result is set to melting point at depth
92  result(i,j) = T0 - beta_CC * pressure;
93  }
94 }
95 
96 //! \brief Computes mass flux in [kg m-2 s-1].
97 /*!
98  * Assumes that mass flux is proportional to the shelf-base heat flux.
99  */
100 void PIK::mass_flux(const array::Scalar &ice_thickness, array::Scalar &result) const {
101  const double
102  melt_factor = m_config->get_number("ocean.pik_melt_factor"),
103  L = m_config->get_number("constants.fresh_water.latent_heat_of_fusion"),
104  sea_water_density = m_config->get_number("constants.sea_water.density"),
105  ice_density = m_config->get_number("constants.ice.density"),
106  c_p_ocean = 3974.0, // J/(K*kg), specific heat capacity of ocean mixed layer
107  gamma_T = 1e-4, // m/s, thermal exchange velocity
108  ocean_salinity = 35.0, // g/kg
109  T_ocean = units::convert(m_sys, -1.7, "Celsius", "Kelvin"); //Default in PISM-PIK
110 
111  //FIXME: gamma_T should be a function of the friction velocity, not a const
112 
113  array::AccessScope list{&ice_thickness, &result};
114 
115  for (auto p = m_grid->points(); p; p.next()) {
116  const int i = p.i(), j = p.j();
117 
118  // compute T_f(i, j) according to beckmann_goosse03, which has the
119  // meaning of the freezing temperature of the ocean water directly
120  // under the shelf, (of salinity 35psu) [this is related to the
121  // Pressure Melting Temperature, see beckmann_goosse03 eq. 2 for
122  // details]
123  double
124  shelfbaseelev = - (ice_density / sea_water_density) * ice_thickness(i,j),
125  T_f = 273.15 + (0.0939 -0.057 * ocean_salinity + 7.64e-4 * shelfbaseelev);
126  // add 273.15 to convert from Celsius to Kelvin
127 
128  // compute ocean_heat_flux according to beckmann_goosse03
129  // positive, if T_oc > T_ice ==> heat flux FROM ocean TO ice
130  double ocean_heat_flux = melt_factor * sea_water_density * c_p_ocean * gamma_T * (T_ocean - T_f); // in W/m^2
131 
132  // TODO: T_ocean -> field!
133 
134  // shelfbmassflux is positive if ice is freezing on; here it is always negative:
135  // same sign as ocean_heat_flux (positive if massflux FROM ice TO ocean)
136  result(i,j) = ocean_heat_flux / (L * ice_density); // m s-1
137 
138  // convert from [m s-1] to [kg m-2 s-1]:
139  result(i,j) *= ice_density;
140  }
141 }
142 
143 } // end of namespace ocean
144 } // end of namespace pism
const units::System::Ptr m_sys
unit system used by this component
Definition: Component.hh:160
const Config::ConstPtr m_config
configuration database used by this component
Definition: Component.hh:158
const Logger::ConstPtr m_log
logger (for easy access)
Definition: Component.hh:162
const std::shared_ptr< const Grid > m_grid
grid used by this component
Definition: Component.hh:156
array::Scalar2 ice_thickness
Definition: Geometry.hh:51
Combines the max. time step with the flag indicating if a restriction is active. Makes is possible to...
Definition: MaxTimestep.hh:31
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
Definition: Array.hh:65
std::shared_ptr< array::Scalar > m_shelf_base_mass_flux
std::shared_ptr< array::Scalar > m_shelf_base_temperature
std::shared_ptr< array::Scalar > m_water_column_pressure
Definition: OceanModel.hh:73
MaxTimestep max_timestep_impl(double t) const
Definition: ConstantPIK.cc:49
void melting_point_temperature(const array::Scalar &depth, array::Scalar &result) const
Definition: ConstantPIK.cc:78
PIK(std::shared_ptr< const Grid > g)
Definition: ConstantPIK.cc:29
void init_impl(const Geometry &geometry)
Definition: ConstantPIK.cc:34
void mass_flux(const array::Scalar &ice_thickness, array::Scalar &result) const
Computes mass flux in [kg m-2 s-1].
Definition: ConstantPIK.cc:100
void update_impl(const Geometry &geometry, double my_t, double my_dt)
Definition: ConstantPIK.cc:54
static const double L
Definition: exactTestL.cc:40
static const double beta_CC
Definition: exactTestO.c:23
bool ocean(int M)
An ocean cell (floating ice or ice-free).
Definition: Mask.hh:40
void compute_average_water_column_pressure(const Geometry &geometry, double ice_density, double water_density, double g, array::Scalar &result)
Definition: OceanModel.cc:217
double convert(System::Ptr system, double input, const std::string &spec1, const std::string &spec2)
Convert a quantity from unit1 to unit2.
Definition: Units.cc:70
static const double g
Definition: exactTestP.cc:36