tch_warming.py - 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
---
tch_warming.py (6007B)
---
1 import PISM
2 from PISM.util import convert
3 import numpy as np
4 import pylab as plt
5
6 ctx = PISM.Context()
7 unit_system = ctx.unit_system
8 config = ctx.config
9
10 config.set_string("grid.ice_vertical_spacing", "equal")
11
12 k = config.get_number("constants.ice.thermal_conductivity")
13 T_melting = config.get_number("constants.fresh_water.melting_point_temperature")
14
15 EC = ctx.enthalpy_converter
16 pressure = np.vectorize(EC.pressure)
17 enthalpy = np.vectorize(EC.enthalpy)
18 temperature = np.vectorize(EC.temperature)
19 melting_temperature = np.vectorize(EC.melting_temperature)
20 water_fraction = np.vectorize(EC.water_fraction)
21
22
23 class EnthalpyColumn(object):
24 "Set up the grid and arrays needed to run column solvers"
25
26 def __init__(self, prefix, Mz, dt, Lz=1000.0):
27 self.Lz = Lz
28 self.z = np.linspace(0, self.Lz, Mz)
29
30 param = PISM.GridParameters()
31 param.Lx = 1e5
32 param.Ly = 1e5
33 param.z = PISM.DoubleVector(self.z)
34 param.Mx = 3
35 param.My = 3
36 param.Mz = Mz
37
38 param.ownership_ranges_from_options(1)
39
40 self.dt = dt
41
42 self.grid = PISM.IceGrid(ctx.ctx, param)
43 grid = self.grid
44
45 self.enthalpy = PISM.model.createEnthalpyVec(grid)
46
47 self.strain_heating = PISM.model.createStrainHeatingVec(grid)
48 self.strain_heating.set(0.0)
49
50 self.u, self.v, self.w = PISM.model.create3DVelocityVecs(grid)
51 self.u.set(0.0)
52 self.v.set(0.0)
53 self.w.set(0.0)
54
55 self.sys = PISM.enthSystemCtx(grid.z(), prefix,
56 grid.dx(), grid.dy(), self.dt,
57 config,
58 self.enthalpy,
59 self.u, self.v, self.w,
60 self.strain_heating,
61 EC)
62
63 def init(self):
64 ice_thickness = self.Lz
65 self.sys.init(1, 1, False, ice_thickness)
66
67
68 def ch_heat_flux(T_ice, T_ch, k, R):
69 "Heat flux from the CH system into the ice. See equation 1."
70 return (k / R**2) * (T_ch - T_ice)
71
72
73 def T_surface(time, mean, amplitude, summer_peak_day):
74 "Surface temperature (cosine yearly cycle)"
75 day_length = 86400
76 summer_peak = summer_peak_day * day_length
77 year_length = float(365 * day_length)
78
79 t = np.mod(time - summer_peak, year_length) / year_length
80
81 return mean + amplitude * np.cos(2 * np.pi * t)
82
83
84 def T_steady_state(H, z, T_surface, G, ice_k):
85 "Steady state ice temperature given T_surface and the geothermal flux G."
86 return T_surface + (H - z) * (G / ice_k)
87
88
89 def E_steady_state(H, z, T_surface, G, ice_k):
90 return enthalpy(T_steady_state(H, z, T_surface, G, ice_k),
91 0.0, pressure(H - z))
92
93
94 def run(T_final_years=10.0, dt_days=1, Lz=1000, Mz=101, R=20, omega=0.005):
95 """Run the one-column cryohydrologic warming setup"""
96
97 T_final = convert(T_final_years, "years", "seconds")
98 dt = convert(dt_days, "days", "seconds")
99
100 ice = EnthalpyColumn("energy.enthalpy", Mz, dt, Lz=Lz)
101 ch = EnthalpyColumn("energy.ch_warming", Mz, dt, Lz=Lz)
102
103 H = ice.Lz
104
105 z = np.array(ice.sys.z())
106 P = pressure(H - z)
107
108 z_coarse = np.array(ice.grid.z())
109 P_coarse = pressure(H - z_coarse)
110
111 T_mean_annual = 268.15 # mean annual temperature, Kelvin
112 T_amplitude = 6 # surface temperature aplitude, Kelvin
113 summer_peak_day = 365/2
114 G = 0.0 # geothermal flux, W/m^2
115
116 E_initial = E_steady_state(H, z_coarse, T_mean_annual, G, k)
117
118 with PISM.vec.Access(nocomm=[ice.enthalpy, ice.strain_heating, ice.u, ice.v, ice.w,
119 ch.enthalpy, ch.strain_heating, ch.u, ch.v, ch.w]):
120
121 # set initial state for the ice enthalpy
122 ice.enthalpy.set_column(1, 1, E_initial)
123
124 # set initial state for the CH system enthalpy
125 ch.enthalpy.set_column(1, 1, E_initial)
126
127 times = []
128 T_ice = []
129 T_ch = []
130 W_ice = []
131 W_ch = []
132 t = 0.0
133 while t < T_final:
134 times.append(t)
135
136 # compute the heat flux from the CH system into the ice
137 T_s = T_surface(t, T_mean_annual, T_amplitude, summer_peak_day)
138
139 T_column_ice = temperature(ice.enthalpy.get_column_vector(1, 1), P_coarse)
140 T_column_ch = temperature(ch.enthalpy.get_column_vector(1, 1), P_coarse)
141
142 if R > 0:
143 Q = ch_heat_flux(T_column_ice, T_column_ch, k, R)
144 else:
145 Q = np.zeros_like(P_coarse)
146
147 # add it to the strain heating term
148 ice.strain_heating.set_column(1, 1, Q)
149 ch.strain_heating.set_column(1, 1, -Q)
150
151 if T_s > T_melting:
152 E_s = EC.enthalpy(T_melting, 0.0, 0.0)
153 else:
154 E_s = EC.enthalpy(T_s, 0.0, 0.0)
155
156 # set boundary conditions and update the ice column
157 ice.init()
158 ice.sys.set_surface_dirichlet_bc(E_s)
159 ice.sys.set_basal_heat_flux(G)
160 x = ice.sys.solve()
161 T_ice.append(temperature(x, P))
162 W_ice.append(water_fraction(x, P))
163 ice.sys.fine_to_coarse(x, 1, 1, ice.enthalpy)
164
165 if T_s > T_melting:
166 E_s = EC.enthalpy(T_melting, omega, 0.0)
167 # Re-set enthalpy in the CH system during the melt season
168 ch.enthalpy.set_column(1, 1, enthalpy(melting_temperature(P), omega, P))
169 else:
170 E_s = EC.enthalpy(T_s, 0.0, 0.0)
171
172 # set boundary conditions and update the CH column
173 ch.init()
174 ch.sys.set_surface_dirichlet_bc(E_s)
175 ch.sys.set_basal_heat_flux(G)
176 x = ch.sys.solve()
177 T_ch.append(temperature(x, P))
178 W_ch.append(water_fraction(x, P))
179 ch.sys.fine_to_coarse(x, 1, 1, ch.enthalpy)
180
181 t += dt
182
183 return z, np.array(times), np.array(T_ice), np.array(W_ice), np.array(T_ch), np.array(W_ch)