tEISMINTII.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
---
tEISMINTII.cc (4776B)
---
1 /* Copyright (C) 2014, 2015, 2016, 2017, 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
20 #include "EISMINTII.hh"
21 #include "pism/coupler/AtmosphereModel.hh"
22 #include "pism/util/ConfigInterface.hh"
23 #include "pism/util/pism_options.hh"
24 #include "pism/util/IceGrid.hh"
25 #include "pism/util/MaxTimestep.hh"
26
27 namespace pism {
28 namespace surface {
29
30 EISMINTII::EISMINTII(IceGrid::ConstPtr g, int experiment)
31 : PSFormulas(g), m_experiment(experiment) {
32 // empty
33 }
34
35 EISMINTII::~EISMINTII() {
36 // empty
37 }
38
39 void EISMINTII::init_impl(const Geometry &geometry) {
40 (void) geometry;
41
42 using units::convert;
43
44 m_log->message(2,
45 "setting parameters for surface mass balance"
46 " and temperature in EISMINT II experiment %c ... \n",
47 m_experiment);
48
49 // EISMINT II specified values for parameters
50 m_S_b = convert(m_sys, 1.0e-2 * 1e-3, "year-1", "second-1"); // Grad of accum rate change
51 m_S_T = 1.67e-2 * 1e-3; // K/m Temp gradient
52
53 // these are for A,E,G,H,I,K:
54 m_M_max = convert(m_sys, 0.5, "m year-1", "m second-1"); // Max accumulation
55 m_R_el = 450.0e3; // Distance to equil line (SMB=0)
56 m_T_min = 238.15;
57
58 switch (m_experiment) {
59 case 'B': // supposed to start from end of experiment A and:
60 m_T_min = 243.15;
61 break;
62 case 'C':
63 case 'J':
64 case 'L': // supposed to start from end of experiment A (for C;
65 // resp I and K for J and L) and:
66 m_M_max = convert(m_sys, 0.25, "m year-1", "m second-1");
67 m_R_el = 425.0e3;
68 break;
69 case 'D': // supposed to start from end of experiment A and:
70 m_R_el = 425.0e3;
71 break;
72 case 'F': // start with zero ice and:
73 m_T_min = 223.15;
74 break;
75 }
76
77 // if user specifies Tmin, Tmax, Mmax, Sb, ST, Rel, then use that (override above)
78 m_T_min = options::Real("-Tmin", "T min, Kelvin", m_T_min);
79
80 options::Real Mmax("-Mmax", "Maximum accumulation, m year-1",
81 convert(m_sys, m_M_max, "m second-1", "m year-1"));
82 if (Mmax.is_set()) {
83 m_M_max = convert(m_sys, Mmax, "m year-1", "m second-1");
84 }
85
86 options::Real Sb("-Sb", "radial gradient of accumulation rate, (m year-1)/km",
87 convert(m_sys, m_S_b, "m second-1/m", "m year-1/km"));
88 if (Sb.is_set()) {
89 m_S_b = convert(m_sys, Sb, "m year-1/km", "m second-1/m");
90 }
91
92 options::Real ST("-ST", "radial gradient of surface temperature, K/km",
93 convert(m_sys, m_S_T, "K/m", "K/km"));
94 if (ST.is_set()) {
95 m_S_T = convert(m_sys, ST, "K/km", "K/m");
96 }
97
98 options::Real Rel("-Rel", "radial distance to equilibrium line, km",
99 convert(m_sys, m_R_el, "m", "km"));
100 if (Rel.is_set()) {
101 m_R_el = convert(m_sys, Rel, "km", "m");
102 }
103
104 initialize_using_formulas();
105 }
106
107 MaxTimestep EISMINTII::max_timestep_impl(double t) const {
108 (void) t;
109 return MaxTimestep("surface EISMINT-II");
110 }
111
112 void EISMINTII::initialize_using_formulas() {
113
114 // center of the accumulation and surface temperature patterns
115 double cx = 0.0, cy = 0.0;
116 if (m_experiment == 'E') {
117 cx += 100.0e3;
118 cy += 100.0e3;
119 }
120
121 IceModelVec::AccessList list{m_temperature.get(), m_mass_flux.get()};
122
123 for (Points p(*m_grid); p; p.next()) {
124 const int i = p.i(), j = p.j();
125
126 const double r = sqrt(PetscSqr(m_grid->x(i) - cx) + PetscSqr(m_grid->y(j) - cy));
127
128 // accumulation (formula (7) in [Payne et al 2000])
129 (*m_mass_flux)(i,j) = std::min(m_M_max, m_S_b * (m_R_el-r));
130
131 // surface temperature (formula (8) in [Payne et al 2000])
132 (*m_temperature)(i,j) = m_T_min + m_S_T * r;
133 }
134
135 // convert from "m second-1" to "kg m-2 s-1"
136 m_mass_flux->scale(m_config->get_number("constants.ice.density"));
137 }
138
139 void EISMINTII::update_impl(const Geometry &geometry, double t, double dt) {
140 (void) t;
141 (void) dt;
142 (void) geometry;
143
144 dummy_accumulation(*m_mass_flux, *m_accumulation);
145 dummy_melt(*m_mass_flux, *m_melt);
146 dummy_runoff(*m_mass_flux, *m_runoff);
147
148 }
149
150 } // end of namespace surface
151 } // end of namespace pism