tIceEISModel.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
(HTM) git clone git://src.adamsgaard.dk/pism
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tIceEISModel.cc (6119B)
---
1 // Copyright (C) 2004-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 <gsl/gsl_math.h> // M_PI
20
21 #include "IceEISModel.hh"
22
23 #include "pism/util/Context.hh"
24 #include "pism/util/ConfigInterface.hh"
25 #include "pism/util/IceGrid.hh"
26
27 #include "pism/coupler/ocean/Constant.hh"
28 #include "pism/coupler/ocean/Initialization.hh"
29
30 #include "pism/coupler/SeaLevel.hh"
31 #include "pism/coupler/ocean/sea_level/Initialization.hh"
32
33 #include "pism/coupler/surface/EISMINTII.hh"
34 #include "pism/coupler/surface/Initialization.hh"
35
36 #include "pism/earth/BedDef.hh"
37
38 namespace pism {
39
40 IceEISModel::IceEISModel(IceGrid::Ptr g, Context::Ptr context, char experiment)
41 : IceModel(g, context), m_experiment(experiment) {
42
43 // the following flag must be here in constructor because
44 // IceModel::createVecs() uses it non-polythermal methods; can be
45 // overridden by the command-line option "-energy enthalpy"
46 m_config->set_flag("energy.temperature_based", true);
47
48 // see EISMINT II description; choose no ocean interaction,
49 m_config->set_flag("ocean.always_grounded", true);
50
51 // purely SIA, and E=1
52 m_config->set_number("stress_balance.sia.enhancement_factor", 1.0);
53
54 // none use bed smoothing & bed roughness parameterization
55 m_config->set_number("stress_balance.sia.bed_smoother.range", 0.0);
56
57 // basal melt does not change computation of mass continuity or vertical velocity:
58 m_config->set_flag("geometry.update.use_basal_melt_rate", false);
59
60 // Make bedrock thermal material properties into ice properties. Note that
61 // zero thickness bedrock layer is the default, but we want the ice/rock
62 // interface segment to have geothermal flux applied directly to ice without
63 // jump in material properties at base.
64 m_config->set_number("energy.bedrock_thermal.density",
65 m_config->get_number("constants.ice.density"));
66 m_config->set_number("energy.bedrock_thermal.conductivity",
67 m_config->get_number("constants.ice.thermal_conductivity"));
68 m_config->set_number("energy.bedrock_thermal.specific_heat_capacity",
69 m_config->get_number("constants.ice.specific_heat_capacity"));
70
71 // no sliding + SIA
72 m_config->set_string("stress_balance.model", "sia");
73 }
74
75 void IceEISModel::allocate_couplers() {
76
77 // Climate will always come from intercomparison formulas.
78 if (not m_surface) {
79 std::shared_ptr<surface::SurfaceModel> surface(new surface::EISMINTII(m_grid, m_experiment));
80 m_surface.reset(new surface::InitializationHelper(m_grid, surface));
81 m_submodels["surface process model"] = m_surface.get();
82 }
83
84 if (not m_ocean) {
85 std::shared_ptr<ocean::OceanModel> ocean(new ocean::Constant(m_grid));
86 m_ocean.reset(new ocean::InitializationHelper(m_grid, ocean));
87 m_submodels["ocean model"] = m_ocean.get();
88 }
89
90 if (not m_sea_level) {
91 using namespace ocean::sea_level;
92 std::shared_ptr<SeaLevel> sea_level(new SeaLevel(m_grid));
93 m_sea_level.reset(new InitializationHelper(m_grid, sea_level));
94 m_submodels["sea level forcing"] = m_sea_level.get();
95 }
96 }
97
98 void generate_trough_topography(IceModelVec2S &result) {
99 // computation based on code by Tony Payne, 6 March 1997:
100 // http://homepages.vub.ac.be/~phuybrec/eismint/topog2.f
101
102 IceGrid::ConstPtr grid = result.grid();
103
104 const double
105 b0 = 1000.0, // plateau elevation
106 L = 750.0e3, // half-width of computational domain
107 w = 200.0e3, // trough width
108 slope = b0 / L,
109 dx61 = (2.0 * L) / 60; // = 25.0e3
110
111 IceModelVec::AccessList list(result);
112 for (Points p(*grid); p; p.next()) {
113 const int i = p.i(), j = p.j();
114
115 const double nsd = i * grid->dx(), ewd = j * grid->dy();
116 if ((nsd >= (27 - 1) * dx61) && (nsd <= (35 - 1) * dx61) &&
117 (ewd >= (31 - 1) * dx61) && (ewd <= (61 - 1) * dx61)) {
118 result(i,j) = 1000.0 - std::max(0.0, slope * (ewd - L) * cos(M_PI * (nsd - L) / w));
119 } else {
120 result(i,j) = 1000.0;
121 }
122 }
123 }
124
125 void generate_mound_topography(IceModelVec2S &result) {
126 // computation based on code by Tony Payne, 6 March 1997:
127 // http://homepages.vub.ac.be/~phuybrec/eismint/topog2.f
128
129 IceGrid::ConstPtr grid = result.grid();
130
131 const double slope = 250.0;
132 const double w = 150.0e3; // mound width
133
134 IceModelVec::AccessList list(result);
135 for (Points p(*grid); p; p.next()) {
136 const int i = p.i(), j = p.j();
137
138 const double nsd = i * grid->dx(), ewd = j * grid->dy();
139 result(i,j) = fabs(slope * sin(M_PI * ewd / w) + slope * cos(M_PI * nsd / w));
140 }
141 }
142
143 void IceEISModel::initialize_2d() {
144
145 m_log->message(2,
146 "initializing variables from EISMINT II experiment %c formulas... \n",
147 m_experiment);
148
149 // set bed topography
150 if (m_experiment == 'I' or m_experiment == 'J') {
151 generate_trough_topography(m_geometry.bed_elevation);
152 } else if (m_experiment == 'K' or m_experiment == 'L') {
153 generate_mound_topography(m_geometry.bed_elevation);
154 } else {
155 m_geometry.bed_elevation.set(0.0);
156 }
157
158 m_geometry.sea_level_elevation.set(0.0);
159
160 // set uplift
161 IceModelVec2S bed_uplift(m_grid, "uplift", WITHOUT_GHOSTS);
162 bed_uplift.set(0.0);
163
164 // start with zero ice
165 m_geometry.ice_thickness.set(0.0);
166
167 m_beddef->bootstrap(m_geometry.bed_elevation, bed_uplift, m_geometry.ice_thickness,
168 m_geometry.sea_level_elevation);
169 }
170
171 } // end of namespace pism