ttestsFG.f90 - 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
---
ttestsFG.f90 (7395B)
---
1 module testsFG
2 ! Copyright (C) 2005-2007 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 ! TESTSFG is a Fortran 90 implementation of two exact solutions for a
22 ! thermocoupled ice sheet. Reference:
23 !
24 ! E. Bueler, J. Brown, and C. Lingle (2007). "Exact solutions to the
25 ! thermomechanically coupled shallow ice approximation: effective tools
26 ! for verification", J. Glaciol., J. Glaciol., vol. 53 no. 182, 499--516.
27 !
28 ! ELB 3/29/05; 7/27/07; 7/29/08
29 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
30
31 public :: testF, testG
32 private :: bothexact, p3, p4
33
34 ! DOUBLE PRECISION DESIRABLE:
35 integer, parameter, public :: kind=8
36
37 real(kind), parameter, public :: SperA = 31556926.0 ! 365.2422 days
38 real(kind), parameter, public :: g=9.81 ! m/s^2; accel of gravity
39 real(kind), parameter, public :: Rgas=8.314 ! J/(mol K)
40
41 ! ice properties; parameters which appear in constitutive relation
42 real(kind), parameter, public :: rho=910.0 ! kg/m^3; density
43 real(kind), parameter, public :: k=2.1 ! J/m K s; thermal conductivity
44 real(kind), parameter, public :: cpheat=2009.0! J/kg K; specific heat capacity
45 real(kind), parameter, public :: n=3 ! Glen exponent
46 ! next two are EISMINT II values; Paterson-Budd for T<263
47 real(kind), parameter, public :: A=3.615e-13 ! Pa^-3 s^-1
48 real(kind), parameter, public :: Q=6.0e4 ! J/mol
49
50 ! EISMINT II temperature boundary condition (Experiment F):
51 real(kind), parameter, public :: Ggeo=.042 ! J/m^2 s; geo. heat flux
52 real(kind), parameter, public :: ST=1.67e-5 ! K m^-1
53 real(kind), parameter, public :: Tmin=223.15 ! K
54
55 ! parameters describing extent of sheet
56 real(kind), parameter, public :: H0=3000.0 ! m
57 real(kind), parameter, public :: L=750000.0 ! m
58
59 ! period and magnitude of perturbation; inactive in Test F:
60 real(kind), parameter, public :: Tp=2000.0*SperA ! s
61 real(kind), parameter, public :: Cp=200.0 ! m
62
63 contains
64
65 subroutine testF(r,z,H,T,U,w,Sig,M,Sigc)
66 real(kind), intent(in) :: r, z
67 real(kind), intent(out) :: H, T, U, w, Sig, M, Sigc
68 call bothexact(0.0_kind,r,z,0.0_kind,H,T,U,w,Sig,M,Sigc)
69 end subroutine testF
70
71 subroutine testG(t,r,z,H,TT,U,w,Sig,M,Sigc)
72 real(kind), intent(in) :: t, r, z
73 real(kind), intent(out) :: H, TT, U, w, Sig, M, Sigc
74 call bothexact(t,r,z,Cp,H,TT,U,w,Sig,M,Sigc)
75 end subroutine testG
76
77 subroutine bothexact(t,r,z,Cp,H,TT,U,w,Sig,M,Sigc)
78 real(kind), intent(in) :: t, r, z, Cp
79 real(kind), intent(out) :: H, TT, U, w, Sig, M, Sigc
80
81 real(kind), parameter :: pi = 3.14159265358979
82 real(kind), parameter :: Kcond=k/(rho*cpheat) ! constant in temp eqn
83
84 ! declare all temporary quantities real(kind); computed in blocks below
85 real(kind) pow, Hconst, s, lamhat, f, goft, Ts, nusqrt, nu
86 real(kind) lamhatr, fr, Hr, mu, I3, surfArr, Uconst, omega
87 real(kind) Sigmu, lamhatrr, frr, Hrr, Tsr, nur, mur, phi, gamma, I4
88 real(kind) I4H, divQ, Ht, nut, dTt, Tr, Tz, Tzz
89
90 if (r<=0 .or. r>=L) then
91 print *,'code and derivation assume 0<r<L'
92 stop
93 end if
94
95 ! compute H from analytical steady state Hs (Test D) plus perturbation
96 pow = n/(2*n+2)
97 Hconst = H0/((1-1/n)**pow)
98 s = r/L
99 lamhat = (1+1/n)*s - (1/n) + (1-s)**(1+1/n) - s**(1+1/n);
100 if (r>0.3*L .and. r<0.9*L) then
101 f = ( cos(pi*(r-0.6*L)/(0.6*L)) )**2
102 else
103 f = 0.0
104 end if
105 goft = Cp*sin(2.0*pi*t/Tp)
106 H = Hconst*(lamhat)**pow + goft*f
107
108 ! compute TT = temperature
109 Ts = Tmin+ST*r
110 nusqrt = sqrt( 1 + (4.0*H*Ggeo)/(k*Ts) )
111 nu = ( k*Ts/(2.0*Ggeo) )*( 1 + nusqrt )
112 TT = Ts * (nu+H) / (nu+z)
113
114 ! compute surface slope and horizontal velocity
115 lamhatr = ((1+1/n)/L)*( 1 - (1-s)**(1/n) - s**(1/n) )
116 if (r>0.3*L .and. r<0.9*L) then
117 fr = -(pi/(0.6*L)) * sin(2.0*pi*(r-0.6*L)/(0.6*L))
118 else
119 fr = 0.0
120 end if
121 Hr = Hconst * pow * lamhat**(pow-1) * lamhatr + goft*fr ! chain rule
122 if (Hr>0) then
123 print *,'code and derivation assume H_r negative for all 0<r<L'
124 stop
125 end if
126 mu = Q/(Rgas*Ts*(nu+H))
127 I3 = p3(mu*H) * exp(mu*H) - p3(mu*(H-z)) * exp(mu*(H-z))
128 surfArr = exp(-Q/(Rgas*Ts))
129 Uconst = 2.0 * (rho*g)**n * A
130 omega = Uconst * (-Hr)**n * surfArr * mu**(-n-1)
131 U = omega * I3
132
133 ! compute strain heating
134 Sigmu = -(Q*(nu+z)) / (Rgas*Ts*(nu+H))
135 Sig = (Uconst*g/cpheat) * exp(Sigmu) * ( abs(Hr)*(H-z) )**(n+1)
136
137 ! compute vertical velocity
138 lamhatrr = ((1+1/n) / (n*L*L)) * ( (1-s)**((1/n)-1) - s**((1/n)-1) )
139 if (r>0.3*L .and. r<0.9*L) then
140 frr = -(2.0*pi*pi/(0.36*L*L)) * cos(2.0*pi*(r-0.6*L)/(0.6*L))
141 else
142 frr = 0.0
143 end if
144 Hrr = Hconst*pow*(pow-1)*(lamhat)**(pow-2) * lamhatr**2 + &
145 Hconst*pow*(lamhat)**(pow-1)*lamhatrr + goft*frr
146 Tsr = ST
147 nur = (k*Tsr/(2.0*Ggeo)) * (1 + nusqrt) + &
148 (1/Ts) * (Hr*Ts-H*Tsr) / nusqrt
149 mur = (-Q/(Rgas*Ts*Ts*(nu+H)**2)) * (Tsr*(nu+H)+Ts*(nur+Hr))
150 phi = 1/r + n*Hrr/Hr + Q*Tsr/(Rgas*Ts*Ts) - (n+1)*mur/mu ! division by r
151 gamma = mu**n * exp(mu*H) * (mur*H+mu*Hr) * H**n
152 I4 = p4(mu*H) * exp(mu*H) - p4(mu*(H-z)) * exp(mu*(H-z))
153 w = omega * ((mur/mu - phi)*I4/mu + (phi*(H-z)+Hr)*I3 - gamma*z)
154
155 ! compute compensatory accumulation M
156 I4H = p4(mu*H) * exp(mu*H) - 24
157 divQ = - omega * (mur/mu - phi) * I4H / mu + omega * gamma * H
158 Ht = (Cp*2.0*pi/Tp) * cos(2.0*pi*t/Tp) * f
159 M = Ht + divQ
160
161 ! compute compensatory heating
162 nut = Ht/nusqrt
163 dTt = Ts * ((nut+Ht)*(nu+z)-(nu+H)*nut) * (nu+z)**(-2)
164 Tr = Tsr*(nu+H)/(nu+z) + Ts * ((nur+Hr)*(nu+z)-(nu+H)*nur) * (nu+z)**(-2)
165 Tz = -Ts * (nu+H) * (nu+z)**(-2)
166 Tzz = 2.0 * Ts * (nu+H) * (nu+z)**(-3)
167 Sigc = dTt + U*Tr + w*Tz - Kcond*Tzz - Sig
168 end subroutine bothexact
169
170 function p3(x)
171 real(kind), intent(in) :: x
172 real(kind) :: p3
173 ! p_3=x^3-3*x^2+6*x-6, using Horner's
174 p3 = -6 + x*(6 + x*(-3 + x))
175 end function p3
176
177 function p4(x)
178 real(kind), intent(in) :: x
179 real(kind) :: p4
180 ! p_4=x^4-4*x^3+12*x^2-24*x+24, using Horner's
181 p4 = 24 + x*(-24 + x*(12 + x*(-4 + x)))
182 end function p4
183
184 end module testsFG