Loading [MathJax]/extensions/tex2jax.js
PISM, A Parallel Ice Sheet Model 2.2.1-cd005eec8 committed by Constantine Khrulev on 2025-03-07
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
iCMthermo.cc
Go to the documentation of this file.
1// Copyright (C) 2004-2018, 2020, 2021, 2022, 2023 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 <cmath>
20
21#include "pism/verification/tests/exactTestsFG.hh"
22#include "pism/verification/tests/exactTestK.h"
23#include "pism/verification/tests/exactTestO.h"
24#include "pism/verification/iceCompModel.hh"
25
26#include "pism/stressbalance/StressBalance.hh"
27#include "pism/util/Time.hh"
28#include "pism/util/Grid.hh"
29#include "pism/util/error_handling.hh"
30#include "pism/earth/BedDef.hh"
31#include "pism/util/ConfigInterface.hh"
32#include "pism/util/pism_utilities.hh"
33#include "pism/verification/BTU_Verification.hh"
34#include "pism/energy/TemperatureModel.hh"
35#include "pism/coupler/SurfaceModel.hh"
36#include "pism/coupler/OceanModel.hh"
37#include "pism/hydrology/Hydrology.hh"
38
39namespace pism {
40
41// boundary conditions for tests F, G (same as EISMINT II Experiment F)
42const double IceCompModel::m_ST = 1.67e-5;
43const double IceCompModel::m_Tmin = 223.15; // K
44const double IceCompModel::m_LforFG = 750000; // m
45const double IceCompModel::m_ApforG = 200; // m
46
48
50
51 array::Scalar &bedtoptemp = *m_work2d[1];
52 array::Scalar &basal_enthalpy = *m_work2d[2];
53 extract_surface(m_energy_model->enthalpy(), 0.0, basal_enthalpy);
54
59 basal_enthalpy,
60 m_surface->temperature(),
61 bedtoptemp);
62
63 m_btu->update(bedtoptemp, t_TempAge, dt_TempAge);
64
65 energy::Inputs inputs;
66 {
67 inputs.basal_frictional_heating = &m_stress_balance->basal_frictional_heating();
68 inputs.basal_heat_flux = &m_btu->flux_through_top_surface(); // bedrock thermal layer
70 inputs.ice_thickness = &m_geometry.ice_thickness; // geometry
71 inputs.shelf_base_temp = &m_ocean->shelf_base_temperature(); // ocean model
72 inputs.surface_liquid_fraction = &m_surface->liquid_water_fraction(); // surface model
73 inputs.surface_temp = &m_surface->temperature(); // surface model
74 inputs.till_water_thickness = &m_subglacial_hydrology->till_water_thickness();
75
76 inputs.volumetric_heating_rate = &m_stress_balance->volumetric_strain_heating();
77 inputs.u3 = &m_stress_balance->velocity_u();
78 inputs.v3 = &m_stress_balance->velocity_v();
79 inputs.w3 = &m_stress_balance->velocity_w();
80
81 inputs.check(); // make sure all data members were set
82 }
83
84 if ((m_testname == 'F') || (m_testname == 'G')) {
85 // Compute compensatory strain heating (fills strain_heating3_comp).
87
88 // Add computed strain heating to the compensatory part.
90
91 // Use the result.
93 }
94
95 m_energy_model->update(t_TempAge, dt_TempAge, inputs);
96
98}
99
101
103
104 const double time = m_testname == 'F' ? 0.0 : m_time->current();
105 const double A = m_testname == 'F' ? 0.0 : m_ApforG;
106
107 for (auto p = m_grid->points(); p; p.next()) {
108 const int i = p.i(), j = p.j();
109
110 // avoid singularity at origin
111 const double r = std::max(grid::radius(*m_grid, i, j), 1.0);
112
113 if (r > m_LforFG - 1.0) { // if (essentially) outside of sheet
114 m_geometry.ice_thickness(i, j) = 0.0;
115 } else {
116 m_geometry.ice_thickness(i, j) = exactFG(time, r, m_grid->z(), A).H;
117 }
118 }
119
121
122 {
123 array::Scalar bed_topography(m_grid, "topg");
124 bed_topography.set(0.0);
125
126 array::Scalar bed_uplift(m_grid, "uplift");
127 bed_uplift.set(0.0);
128
129 array::Scalar sea_level(m_grid, "sea_level");
130 sea_level.set(0.0);
131
132 m_beddef->bootstrap(bed_topography, bed_uplift, m_geometry.ice_thickness,
133 sea_level);
134 }
135}
136
138
139 array::Scalar bed_topography(m_grid, "topg");
140 bed_topography.set(0.0);
141
142 array::Scalar bed_uplift(m_grid, "uplift");
143 bed_uplift.set(0.0);
144
145 array::Scalar sea_level(m_grid, "sea_level");
146 sea_level.set(0.0);
147
149
150 m_beddef->bootstrap(bed_topography, bed_uplift, m_geometry.ice_thickness,
151 sea_level);
152}
153
155
156 const double
157 ice_rho = m_config->get_number("constants.ice.density"),
158 ice_c = m_config->get_number("constants.ice.specific_heat_capacity");
159
160 // before temperature and flow step, set strain_heating_c from exact values
161
163
164 const double time = m_testname == 'F' ? 0.0 : m_time->current();
165 const double A = m_testname == 'F' ? 0.0 : m_ApforG;
166
167 for (auto p = m_grid->points(); p; p.next()) {
168 const int i = p.i(), j = p.j();
169
170 double r = std::max(grid::radius(*m_grid, i, j), 1.0); // avoid singularity at origin
171
172 if (r > m_LforFG - 1.0) { // outside of sheet
174 } else {
175 TestFGParameters P = exactFG(time, r, m_grid->z(), A);
176
177 m_strain_heating3_comp.set_column(i, j, (P.Sigc).data());
178 }
179 }
180
181 // scale strain_heating to J/(s m^3)
182 m_strain_heating3_comp.scale(ice_rho * ice_c);
183}
184
186 double &gavTerr) {
187 double maxTerr = 0.0, avTerr = 0.0, avcount = 0.0;
188
189 if (m_testname != 'F' and m_testname != 'G') {
190 throw RuntimeError(PISM_ERROR_LOCATION, "temperature errors only computable for tests F and G");
191 }
192
193 const double time = m_testname == 'F' ? 0.0 : m_time->current();
194 const double A = m_testname == 'F' ? 0.0 : m_ApforG;
195
196 auto *m = dynamic_cast<energy::TemperatureModel*>(m_energy_model.get());
197 const auto &ice_temperature = m->temperature();
198
199 array::AccessScope list{&m_geometry.ice_thickness, &ice_temperature};
200
201 ParallelSection loop(m_grid->com);
202 try {
203 for (auto p = m_grid->points(); p; p.next()) {
204 const int i = p.i(), j = p.j();
205
206 const double r = grid::radius(*m_grid, i, j);
207 const double *T = ice_temperature.get_column(i, j);
208
209 // only evaluate error if inside sheet and not at central
210 // singularity
211 if ((r >= 1.0) and (r <= m_LforFG - 1.0)) {
212 TestFGParameters P = exactFG(time, r, m_grid->z(), A);
213
214 // only evaluate error if below ice surface
215 const int ks = m_grid->kBelowHeight(m_geometry.ice_thickness(i, j));
216 for (int k = 0; k < ks; k++) {
217 const double Terr = fabs(T[k] - P.T[k]);
218 maxTerr = std::max(maxTerr, Terr);
219 avcount += 1.0;
220 avTerr += Terr;
221 }
222 }
223 }
224 } catch (...) {
225 loop.failed();
226 }
227 loop.check();
228
229 gmaxTerr = GlobalMax(m_grid->com, maxTerr);
230 gavTerr = GlobalSum(m_grid->com, avTerr);
231 double gavcount = GlobalSum(m_grid->com, avcount);
232 gavTerr = gavTerr / std::max(gavcount, 1.0); // avoid division by zero
233}
234
235
236void IceCompModel::computeIceBedrockTemperatureErrors(double &gmaxTerr, double &gavTerr,
237 double &gmaxTberr, double &gavTberr) {
238
239 if (m_testname != 'K' and m_testname != 'O') {
240 throw RuntimeError(PISM_ERROR_LOCATION, "ice and bedrock temperature errors only computable for tests K and O");
241 }
242
243 double maxTerr = 0.0, avTerr = 0.0, avcount = 0.0;
244 double maxTberr = 0.0, avTberr = 0.0, avbcount = 0.0;
245
246 const double *Tb, *T;
247 std::vector<double> Tex(m_grid->Mz());
248
249 auto *my_btu = dynamic_cast<energy::BTU_Verification*>(m_btu.get());
250 if (my_btu == NULL) {
251 throw RuntimeError(PISM_ERROR_LOCATION, "BTU_Verification is required");
252 }
253
254 const auto &bedrock_temp = my_btu->temperature();
255
256 auto zblevels = bedrock_temp.levels();
257 auto Mbz = zblevels.size();
258 std::vector<double> Tbex(Mbz);
259
260 switch (m_testname) {
261 case 'K':
262 for (unsigned int k = 0; k < m_grid->Mz(); k++) {
264 Tex[k] = K.T;
265 }
266 for (unsigned int k = 0; k < Mbz; k++) {
267 TestKParameters K = exactK(m_time->current(), zblevels[k], m_bedrock_is_ice_forK);
268 Tbex[k] = K.T;
269 }
270 break;
271 case 'O':
272 for (unsigned int k = 0; k < m_grid->Mz(); k++) {
273 Tex[k] = exactO(m_grid->z(k)).TT;
274 }
275 for (unsigned int k = 0; k < Mbz; k++) {
276 Tbex[k] = exactO(zblevels[k]).TT;
277 }
278 break;
279 default:
280 throw RuntimeError(PISM_ERROR_LOCATION, "ice and bedrock temperature errors only for tests K and O");
281 }
282
283 auto *m = dynamic_cast<energy::TemperatureModel*>(m_energy_model.get());
284 const auto &ice_temperature = m->temperature();
285
286 array::AccessScope list{&ice_temperature, &bedrock_temp};
287
288 for (auto p = m_grid->points(); p; p.next()) {
289 const int i = p.i(), j = p.j();
290
291 Tb = bedrock_temp.get_column(i, j);
292 for (unsigned int kb = 0; kb < Mbz; kb++) {
293 const double Tberr = fabs(Tb[kb] - Tbex[kb]);
294 maxTberr = std::max(maxTberr, Tberr);
295 avbcount += 1.0;
296 avTberr += Tberr;
297 }
298 T = ice_temperature.get_column(i, j);
299 for (unsigned int k = 0; k < m_grid->Mz(); k++) {
300 const double Terr = fabs(T[k] - Tex[k]);
301 maxTerr = std::max(maxTerr, Terr);
302 avcount += 1.0;
303 avTerr += Terr;
304 }
305 }
306
307 gmaxTerr = GlobalMax(m_grid->com, maxTerr);
308 gavTerr = GlobalSum(m_grid->com, avTerr);
309 double gavcount = GlobalSum(m_grid->com, avcount);
310 gavTerr = gavTerr/std::max(gavcount, 1.0); // avoid division by zero
311
312 gmaxTberr = GlobalMax(m_grid->com, maxTberr);
313 gavTberr = GlobalSum(m_grid->com, avTberr);
314 double gavbcount = GlobalSum(m_grid->com, avbcount);
315 gavTberr = gavTberr/std::max(gavbcount, 1.0); // avoid division by zero
316}
317
318
319void IceCompModel::computeBasalTemperatureErrors(double &gmaxTerr, double &gavTerr, double &centerTerr) {
320
321 if (m_testname != 'F' and m_testname != 'G') {
322 throw RuntimeError(PISM_ERROR_LOCATION, "temperature errors only computable for tests F and G");
323 }
324
325 double
326 Texact = 0.0,
327 domeT = 0.0,
328 domeTexact = 0.0,
329 Terr = 0.0,
330 avTerr = 0.0;
331
332 const double time = m_testname == 'F' ? 0.0 : m_time->current();
333 const double A = m_testname == 'F' ? 0.0 : m_ApforG;
334 std::vector<double> z(1, 0.0);
335
336 auto *m = dynamic_cast<energy::TemperatureModel*>(m_energy_model.get());
337 const auto &ice_temperature = m->temperature();
338
339 array::AccessScope list(ice_temperature);
340
341 ParallelSection loop(m_grid->com);
342 try {
343 for (auto p = m_grid->points(); p; p.next()) {
344 const int i = p.i(), j = p.j();
345
346 double r = std::max(grid::radius(*m_grid, i, j), 1.0);
347
348 if (r > m_LforFG - 1.0) { // outside of sheet
349 Texact = m_Tmin + m_ST * r; // = Ts
350 } else {
351 Texact = exactFG(time, r, z, A).T[0];
352 }
353
354 const double Tbase = ice_temperature.get_column(i,j)[0];
355 if (i == ((int)m_grid->Mx() - 1) / 2 and
356 j == ((int)m_grid->My() - 1) / 2) {
357 domeT = Tbase;
358 domeTexact = Texact;
359 }
360 // compute maximum errors
361 Terr = std::max(Terr, fabs(Tbase - Texact));
362 // add to sums for average errors
363 avTerr += fabs(Tbase - Texact);
364 }
365 } catch (...) {
366 loop.failed();
367 }
368 loop.check();
369
370 double gdomeT, gdomeTexact;
371
372 gmaxTerr = GlobalMax(m_grid->com, Terr);
373 gavTerr = GlobalSum(m_grid->com, avTerr);
374 gavTerr = gavTerr/(m_grid->Mx()*m_grid->My());
375 gdomeT = GlobalMax(m_grid->com, domeT);
376 gdomeTexact = GlobalMax(m_grid->com, domeTexact);
377 centerTerr = fabs(gdomeT - gdomeTexact);
378}
379
380
381void IceCompModel::compute_strain_heating_errors(double &gmax_strain_heating_err, double &gav_strain_heating_err) {
382 double max_strain_heating_err = 0.0, av_strain_heating_err = 0.0, avcount = 0.0;
383
384 if (m_testname != 'F' and m_testname != 'G') {
385 throw RuntimeError(PISM_ERROR_LOCATION, "strain-heating (strain_heating) errors only computable for tests F and G");
386 }
387
388 const double
389 ice_rho = m_config->get_number("constants.ice.density"),
390 ice_c = m_config->get_number("constants.ice.specific_heat_capacity");
391
392 const auto &strain_heating3 = m_stress_balance->volumetric_strain_heating();
393
394 array::AccessScope list{&m_geometry.ice_thickness, &strain_heating3};
395
396 const double time = m_testname == 'F' ? 0.0 : m_time->current();
397 const double A = m_testname == 'F' ? 0.0 : m_ApforG;
398
399 ParallelSection loop(m_grid->com);
400 try {
401 for (auto p = m_grid->points(); p; p.next()) {
402 const int i = p.i(), j = p.j();
403
404 double r = grid::radius(*m_grid, i, j);
405 if ((r >= 1.0) && (r <= m_LforFG - 1.0)) {
406 // only evaluate error if inside sheet and not at central singularity
407
408 TestFGParameters P = exactFG(time, r, m_grid->z(), A);
409
410 for (unsigned int k = 0; k < m_grid->Mz(); k++) {
411 // scale exact strain_heating to J/(s m^3)
412 P.Sig[k] *= ice_rho * ice_c;
413 }
414
415 const unsigned int ks = m_grid->kBelowHeight(m_geometry.ice_thickness(i, j));
416 const double *strain_heating = strain_heating3.get_column(i, j);
417
418 for (unsigned int k = 0; k < ks; k++) { // only evaluate error if below ice surface
419 const double strain_heating_err = fabs(strain_heating[k] - P.Sig[k]);
420 max_strain_heating_err = std::max(max_strain_heating_err, strain_heating_err);
421 avcount += 1.0;
422 av_strain_heating_err += strain_heating_err;
423 }
424 }
425 }
426 } catch (...) {
427 loop.failed();
428 }
429 loop.check();
430
431 gmax_strain_heating_err = GlobalMax(m_grid->com, max_strain_heating_err);
432 gav_strain_heating_err = GlobalSum(m_grid->com, av_strain_heating_err);
433 double gavcount = GlobalSum(m_grid->com, avcount);
434 gav_strain_heating_err = gav_strain_heating_err/std::max(gavcount, 1.0); // avoid div by zero
435}
436
437
438void IceCompModel::computeSurfaceVelocityErrors(double &gmaxUerr, double &gavUerr,
439 double &gmaxWerr, double &gavWerr) {
440 double
441 maxUerr = 0.0,
442 maxWerr = 0.0,
443 avUerr = 0.0,
444 avWerr = 0.0;
445
446 const array::Array3D
447 &u3 = m_stress_balance->velocity_u(),
448 &v3 = m_stress_balance->velocity_v(),
449 &w3 = m_stress_balance->velocity_w();
450
451 array::AccessScope list{&m_geometry.ice_thickness, &u3, &v3, &w3};
452
453 const double time = m_testname == 'F' ? 0.0 : m_time->current();
454 const double A = m_testname == 'F' ? 0.0 : m_ApforG;
455
456 ParallelSection loop(m_grid->com);
457 try {
458 for (auto p = m_grid->points(); p; p.next()) {
459 const int i = p.i(), j = p.j();
460
461 const double H = m_geometry.ice_thickness(i, j);
462 std::vector<double> z(1, H);
463
464 const double
465 x = m_grid->x(i),
466 y = m_grid->y(j),
467 r = grid::radius(*m_grid, i, j);
468
469 if ((r >= 1.0) and (r <= m_LforFG - 1.0)) {
470 // only evaluate error if inside sheet and not at central singularity
471
472 TestFGParameters P = exactFG(time, r, z, A);
473
474 const double
475 uex = (x/r) * P.U[0],
476 vex = (y/r) * P.U[0];
477 // note that because interpolate does linear interpolation and H(i, j) is not exactly at
478 // a grid point, this causes nonzero errors
479 double
480 u_err = u3.interpolate(i, j, H) - uex,
481 v_err = v3.interpolate(i, j, H) - vex,
482 Uerr = sqrt(u_err * u_err + v_err * v_err);
483 maxUerr = std::max(maxUerr, Uerr);
484 avUerr += Uerr;
485 const double Werr = fabs(w3.interpolate(i, j, H) - P.w[0]);
486 maxWerr = std::max(maxWerr, Werr);
487 avWerr += Werr;
488 }
489 }
490 } catch (...) {
491 loop.failed();
492 }
493 loop.check();
494
495 gmaxUerr = GlobalMax(m_grid->com, maxUerr);
496 gmaxWerr = GlobalMax(m_grid->com, maxWerr);
497 gavUerr = GlobalSum(m_grid->com, avUerr);
498 gavUerr = gavUerr/(m_grid->Mx()*m_grid->My());
499 gavWerr = GlobalSum(m_grid->com, avWerr);
500 gavWerr = gavWerr/(m_grid->Mx()*m_grid->My());
501}
502
503
504void IceCompModel::computeBasalMeltRateErrors(double &gmaxbmelterr, double &gminbmelterr) {
505 double
506 maxbmelterr = -9.99e40,
507 minbmelterr = 9.99e40;
508
509 if (m_testname != 'O') {
510 throw RuntimeError(PISM_ERROR_LOCATION, "basal melt rate errors are only computable for test O");
511 }
512
513 // we just need one constant from exact solution:
514 double bmelt = exactO(0.0).bmelt;
515
516 const array::Scalar &basal_melt_rate = m_energy_model->basal_melt_rate();
517
518 array::AccessScope list(basal_melt_rate);
519
520 for (auto p = m_grid->points(); p; p.next()) {
521 const int i = p.i(), j = p.j();
522
523 double err = fabs(basal_melt_rate(i, j) - bmelt);
524 maxbmelterr = std::max(maxbmelterr, err);
525 minbmelterr = std::min(minbmelterr, err);
526 }
527
528 gmaxbmelterr = GlobalMax(m_grid->com, maxbmelterr);
529 gminbmelterr = GlobalMin(m_grid->com, minbmelterr);
530}
531
532} // end of namespace pism
array::Scalar1 sea_level_elevation
Definition Geometry.hh:48
array::CellType2 cell_type
Definition Geometry.hh:55
array::Scalar2 ice_thickness
Definition Geometry.hh:51
array::Scalar2 bed_elevation
Definition Geometry.hh:47
void compute_strain_heating_errors(double &gmax_strain_heating_err, double &gav_strain_heating_err)
Definition iCMthermo.cc:381
array::Array3D m_strain_heating3_comp
void computeSurfaceVelocityErrors(double &gmaxUerr, double &gavUerr, double &gmaxWerr, double &gavWerr)
Definition iCMthermo.cc:438
void getCompSourcesTestFG()
Definition iCMthermo.cc:154
static const double m_ST
static const double m_Tmin
void computeIceBedrockTemperatureErrors(double &gmaxTerr, double &gavTerr, double &gmaxTberr, double &gavTberr)
Definition iCMthermo.cc:236
void computeTemperatureErrors(double &gmaxTerr, double &gavTerr)
Definition iCMthermo.cc:185
static const double m_LforFG
void computeBasalMeltRateErrors(double &gmaxbmelterr, double &gminbmelterr)
Definition iCMthermo.cc:504
virtual void energy_step()
Manage the solution of the energy equation, and related parallel communication.
Definition iCMthermo.cc:47
static const double m_ApforG
void computeBasalTemperatureErrors(double &gmaxTerr, double &gavTerr, double &centerTerr)
Definition iCMthermo.cc:319
std::shared_ptr< stressbalance::StressBalance > m_stress_balance
Definition IceModel.hh:395
std::shared_ptr< ocean::OceanModel > m_ocean
Definition IceModel.hh:280
std::shared_ptr< surface::SurfaceModel > m_surface
Definition IceModel.hh:279
const Time::Ptr m_time
Time manager.
Definition IceModel.hh:246
Geometry m_geometry
Definition IceModel.hh:288
std::vector< std::shared_ptr< array::Scalar2 > > m_work2d
Definition IceModel.hh:393
std::string m_stdout_flags
Definition IceModel.hh:321
Config::Ptr m_config
Configuration flags and parameters.
Definition IceModel.hh:238
std::shared_ptr< energy::BedThermalUnit > m_btu
Definition IceModel.hh:259
double t_TempAge
time of last update for enthalpy/temperature
Definition IceModel.hh:313
std::unique_ptr< hydrology::Hydrology > m_subglacial_hydrology
Definition IceModel.hh:254
double dt_TempAge
enthalpy/temperature and age time-steps
Definition IceModel.hh:315
std::shared_ptr< energy::EnergyModel > m_energy_model
Definition IceModel.hh:260
const std::shared_ptr< Grid > m_grid
Computational grid.
Definition IceModel.hh:236
std::shared_ptr< bed::BedDef > m_beddef
Definition IceModel.hh:284
void failed()
Indicates a failure of a parallel section.
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
Definition Array.hh:64
double interpolate(int i, int j, double z) const
Return value of scalar quantity at level z (m) above base of ice (by linear interpolation).
Definition Array3D.cc:89
void set_column(int i, int j, double c)
Set all values of scalar quantity to given a single value in a particular column.
Definition Array3D.cc:50
A virtual class collecting methods common to ice and bedrock 3D fields.
Definition Array3D.hh:33
void scale(double alpha)
Result: v <- v * alpha. Calls VecScale.
Definition Array.cc:224
void set(double c)
Result: v[j] <- c for all j.
Definition Array.cc:629
void add(double alpha, const Array &x)
Result: v <- v + alpha * x. Calls VecAXPY.
Definition Array.cc:206
void update_ghosts()
Updates ghost points.
Definition Array.cc:615
const array::Scalar * surface_liquid_fraction
const array::Scalar1 * ice_thickness
const array::Array3D * w3
const array::Array3D * v3
const array::Scalar * shelf_base_temp
const array::Scalar * basal_heat_flux
const array::Scalar * basal_frictional_heating
const array::Scalar * till_water_thickness
const array::Array3D * u3
const array::Scalar * surface_temp
const array::CellType * cell_type
const array::Array3D * volumetric_heating_rate
const array::Array3D & temperature() const
#define PISM_ERROR_LOCATION
struct TestKParameters exactK(double t, double z, int bedrock_is_ice)
Definition exactTestK.c:129
struct TestOParameters exactO(double z)
Definition exactTestO.c:112
double radius(const Grid &grid, int i, int j)
Returns the distance from the point (i,j) to the origin.
Definition Grid.cc:1377
void GlobalMax(MPI_Comm comm, double *local, double *result, int count)
static const double k
Definition exactTestP.cc:42
TestFGParameters exactFG(double t, double r, const std::vector< double > &z, double Cp)
void bedrock_surface_temperature(const array::Scalar &sea_level, const array::CellType &cell_type, const array::Scalar &bed_topography, const array::Scalar &ice_thickness, const array::Scalar &basal_enthalpy, const array::Scalar &ice_surface_temperature, array::Scalar &result)
Compute the temperature seen by the top of the bedrock thermal layer.
Definition energy.cc:121
void GlobalMin(MPI_Comm comm, double *local, double *result, int count)
void GlobalSum(MPI_Comm comm, double *local, double *result, int count)
std::vector< double > U
std::vector< double > Sig
std::vector< double > Sigc
std::vector< double > w
std::vector< double > T