texactTestH.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
---
texactTestH.c (2693B)
---
1 /*
2 Copyright (C) 2004-2006, 2014, 2016 Jed Brown and 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 <stdio.h>
22 #include <math.h>
23 #include "exactTestH.h"
24
25 static const double SperA = 31556926.0; /* seconds per year; 365.2422 days */
26
27 int exactH_old(const double f, const double tIN, const double r,
28 double *H, double *M) {
29
30 const double n = 3.0;
31 const double H0 = 3600.0, R0=750000.0;
32 /* t0 = (beta/Gamma) * pow((2n+1)/((n+1)(1-f)),n) * (pow(R0,n+1)/pow(H0,2n+1))
33 when beta=2; */
34 double t0 = (15208.0 / pow(1-f,n)) * SperA;
35 /* t0 = 40033 years; for test C with isostasy f = rho_ice/rho_rock with
36 rho_ice = 910 and rho_rock = 3300 kg/m^3 */
37 double lambda, alpha, beta, t0post, Rmargin;
38 double t;
39
40 t = tIN;
41
42 if (t < t0) { /* t <= t0: version of test C */
43 lambda = 5.0;
44 alpha = -1.0; /* alpha=(2-(n+1)*lambda)/(5*n+3) */
45 beta = 2.0; /* beta=(1+(2*n+1)*lambda)/(5*n+3) */
46 } else { /* t >= t0: version of test B */
47 t0post = (t0 / 2.0) * (1.0 / 18.0); /* reset t and t0 */
48 t = t - t0 + t0post; /* reset to Halfar w. f */
49 t0 = t0post;
50 lambda = 0.0;
51 alpha = 1.0 / 9.0; /* alpha=(2-(n+1)*lambda)/(5*n+3)=1/9 */
52 beta = 1.0 / 18.0; /* beta=(1+(2*n+1)*lambda)/(5*n+3)=1/18 */
53 }
54
55 Rmargin = R0 * pow(t / t0, beta);
56 if (r < Rmargin)
57 *H = H0 * pow(t / t0, -alpha)
58 * pow(1.0 - pow(pow(t / t0, -beta) * (r / R0), (n + 1.0) / n),
59 n / (2.0 * n + 1.0));
60 else
61 *H = 0.0;
62
63 if (t > 0.1*SperA)
64 *M = (lambda / t) * (*H);
65 else { /* when less than 0.1 year, avoid division by time */
66 Rmargin = R0 * pow(0.1*SperA / t0, beta);
67 if (r < Rmargin)
68 *M = lambda * H0 / t0; /* constant value in disc of Rmargin radius */
69 else
70 *M = 0.0;
71 }
72
73 return 0;
74 }
75
76 struct TestHParameters exactH(const double f, const double t, const double r) {
77 struct TestHParameters result;
78 result.error_code = exactH_old(f, t, r, &result.H, &result.M);
79 return result;
80 }