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)