tenergy.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
---
tenergy.cc (6648B)
---
1 // Copyright (C) 2004-2011, 2013, 2014, 2015, 2016, 2017, 2018 Jed Brown, Ed Bueler and Constantine Khroulev
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 <cassert>
20
21 #include "IceModel.hh"
22
23 #include "pism/energy/BedThermalUnit.hh"
24 #include "pism/util/IceGrid.hh"
25 #include "pism/util/Mask.hh"
26 #include "pism/util/ConfigInterface.hh"
27 #include "pism/util/error_handling.hh"
28 #include "pism/util/pism_utilities.hh"
29 #include "pism/coupler/SurfaceModel.hh"
30 #include "pism/util/EnthalpyConverter.hh"
31 #include "pism/util/Profiling.hh"
32
33 #include "pism/hydrology/Hydrology.hh"
34 #include "pism/stressbalance/StressBalance.hh"
35 #include "pism/energy/EnergyModel.hh"
36 #include "pism/energy/utilities.hh"
37
38 namespace pism {
39
40 //! \file energy.cc Methods of IceModel which address conservation of energy.
41 //! Common to enthalpy (polythermal) and temperature (cold-ice) methods.
42
43 void IceModel::bedrock_thermal_model_step() {
44
45 const Profiling &profiling = m_ctx->profiling();
46
47 IceModelVec2S &basal_enthalpy = m_work2d[2];
48
49 m_energy_model->enthalpy().getHorSlice(basal_enthalpy, 0.0);
50
51 bedrock_surface_temperature(m_geometry.sea_level_elevation,
52 m_geometry.cell_type,
53 m_geometry.bed_elevation,
54 m_geometry.ice_thickness,
55 basal_enthalpy,
56 m_surface->temperature(),
57 m_bedtoptemp);
58
59 profiling.begin("btu");
60 m_btu->update(m_bedtoptemp, t_TempAge, dt_TempAge);
61 profiling.end("btu");
62 }
63
64 //! Manage the solution of the energy equation, and related parallel communication.
65 void IceModel::energy_step() {
66
67 // operator-splitting occurs here (ice and bedrock energy updates are split):
68 // tell BedThermalUnit* btu that we have an ice base temp; it will return
69 // the z=0 value of geothermal flux when called inside temperatureStep() or
70 // enthalpyStep()
71 bedrock_thermal_model_step();
72
73 m_energy_model->update(t_TempAge, dt_TempAge, energy_model_inputs());
74
75 m_stdout_flags = m_energy_model->stdout_flags() + m_stdout_flags;
76 }
77
78 //! @brief Combine basal melt rate in grounded and floating areas.
79 /**
80 * Grounded basal melt rate is computed as a part of the energy
81 * (enthalpy or temperature) step; floating basal melt rate is
82 * provided by an ocean model.
83 *
84 * This method updates IceModel::basal_melt_rate (in meters per second
85 * ice-equivalent).
86 *
87 * The sub shelf mass flux provided by an ocean model is in [kg m-2
88 * s-1], so we divide by the ice density to convert to [m second-1].
89 */
90 void IceModel::combine_basal_melt_rate(const Geometry &geometry,
91 const IceModelVec2S &shelf_base_mass_flux,
92 const IceModelVec2S &grounded_basal_melt_rate,
93 IceModelVec2S &result) {
94
95 const bool sub_gl = (m_config->get_flag("geometry.grounded_cell_fraction") and
96 m_config->get_flag("energy.basal_melt.use_grounded_cell_fraction"));
97
98 IceModelVec::AccessList list{&geometry.cell_type,
99 &grounded_basal_melt_rate, &shelf_base_mass_flux, &result};
100 if (sub_gl) {
101 list.add(geometry.cell_grounded_fraction);
102 }
103
104 double ice_density = m_config->get_number("constants.ice.density");
105
106 for (Points p(*m_grid); p; p.next()) {
107 const int i = p.i(), j = p.j();
108
109 double lambda = 1.0; // 1.0 corresponds to the grounded case
110 // Note: here we convert shelf base mass flux from [kg m-2 s-1] to [m s-1]:
111 const double
112 M_shelf_base = shelf_base_mass_flux(i,j) / ice_density;
113
114 // Use the fractional floatation mask to adjust the basal melt
115 // rate near the grounding line:
116 if (sub_gl) {
117 lambda = geometry.cell_grounded_fraction(i,j);
118 } else if (geometry.cell_type.ocean(i,j)) {
119 lambda = 0.0;
120 }
121 result(i,j) = lambda * grounded_basal_melt_rate(i, j) + (1.0 - lambda) * M_shelf_base;
122 }
123 }
124
125 //! @brief Compute the temperature seen by the top of the bedrock thermal layer.
126 void bedrock_surface_temperature(const IceModelVec2S &sea_level,
127 const IceModelVec2CellType &cell_type,
128 const IceModelVec2S &bed_topography,
129 const IceModelVec2S &ice_thickness,
130 const IceModelVec2S &basal_enthalpy,
131 const IceModelVec2S &ice_surface_temperature,
132 IceModelVec2S &result) {
133
134 IceGrid::ConstPtr grid = result.grid();
135 Config::ConstPtr config = grid->ctx()->config();
136
137 const double
138 T0 = config->get_number("constants.fresh_water.melting_point_temperature"),
139 beta_CC_grad_sea_water = (config->get_number("constants.ice.beta_Clausius_Clapeyron") *
140 config->get_number("constants.sea_water.density") *
141 config->get_number("constants.standard_gravity")); // K m-1
142
143 EnthalpyConverter::Ptr EC = grid->ctx()->enthalpy_converter();
144
145 IceModelVec::AccessList list{&cell_type, &bed_topography, &sea_level, &ice_thickness,
146 &ice_surface_temperature, &basal_enthalpy, &result};
147 ParallelSection loop(grid->com);
148 try {
149 for (Points p(*grid); p; p.next()) {
150 const int i = p.i(), j = p.j();
151
152 if (cell_type.grounded(i,j)) {
153 if (cell_type.ice_free(i,j)) { // no ice: sees air temp
154 result(i,j) = ice_surface_temperature(i,j);
155 } else { // ice: sees temp of base of ice
156 const double pressure = EC->pressure(ice_thickness(i,j));
157 result(i,j) = EC->temperature(basal_enthalpy(i,j), pressure);
158 }
159 } else { // floating: apply pressure melting temp as top of bedrock temp
160 result(i,j) = T0 - (sea_level(i, j) - bed_topography(i,j)) * beta_CC_grad_sea_water;
161 }
162 }
163 } catch (...) {
164 loop.failed();
165 }
166 loop.check();
167 }
168
169 } // end of namespace pism