texactTestN.c - 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
---
texactTestN.c (2990B)
---
1 /*
2 Copyright (C) 2010, 2014, 2016 Ed Bueler
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
21 #include <math.h>
22 #include "exactTestN.h"
23
24 #define secpera 31556926.0 /* seconds per year; 365.2422 days */
25 #define g 9.81
26 #define rho 910.0 /* ice density; kg m-3 */
27 #define rhow 1028.0 /* sea water density; kg m-3 */
28 #define n 3.0 /* Glen power */
29
30 struct TestNConstants exactNConstants() {
31 double s = 0.0;
32 struct TestNConstants result;
33
34 /* geometry */
35 result.H0 = 3000.0;
36 result.L0 = 500.0e3;
37 result.xc = 0.9 * (result.L0);
38
39 /* mass balance */
40 result.a = 0.003 / secpera; /* s-1; mass balance gradient with elevation */
41 result.H_ela = (result.H0) / 1.5; /* m; H0 = 1.5 H_ela exactly */
42
43 /* sliding */
44 /* s m-1; choose k so that eqn (24) gives our L0 */
45 result.k = 9.0 * (result.H_ela) / ((result.a) * (result.L0) * (result.L0));
46
47 /* grounded calving front boundary condition, imposed at xc = .9 L0, determines
48 constant vertically-integrated longitudinal stress T; see (2.12) in Schoof (2006);
49 treats Hc = H(xc) as exactly at flotation */
50 s = (result.xc) / (result.L0);
51 result.H_xc = (result.H0) * (1.0 - s * s);
52 result.T_xc = 0.5 * (1.0 - rho / rhow) * rho * g * (result.H_xc) * (result.H_xc);
53
54 return result;
55 }
56
57 struct TestNParameters exactN(double x) {
58
59 double q = 0.0, hxx = 0.0, ux = 0.0;
60 const struct TestNConstants c = exactNConstants();
61 struct TestNParameters result;
62 result.error_code = 0;
63
64 if (x < 0.0) {
65 result.error_code = 1;
66 return result;
67 }
68
69 if (x > c.L0) {
70 result.error_code = 2;
71 return result;
72 }
73
74 q = (1.0 / n) - 1.0; /* a useful power */
75 hxx = - 2.0 * c.H0 / (c.L0 * c.L0); /* constant concavity of h(x) */
76 ux = - hxx / c.k; /* constant strain rate */
77
78 result.H = c.H0 * (1.0 - (x / c.L0) * (x / c.L0)); /* eqn (23) in Bodvardsson */
79
80 result.h_x = hxx * x;
81
82 result.u = - (result.h_x) / c.k; /* eqn (10) in Bodvardson, once SIA is dropped */
83
84 result.M = c.a * ((result.H) - c.H_ela); /* page 6 in Bodvardsson, just before eqn (23) */
85
86 result.B = c.T_xc / (2.0 * (result.H) * pow(fabs(ux),q) * ux); /* Bueler interpretation */
87
88 result.beta = c.k * rho * g * (result.H);
89
90 return result;
91 }
92
93