tenergy_model.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
---
tenergy_model.py (4520B)
---
1 #!/usr/bin/env python3
2
3 import PISM
4 from PISM.util import convert
5
6 ctx = PISM.Context()
7
8 ctx.log.set_threshold(1)
9
10 def create_dummy_grid():
11 "Create a dummy grid"
12 params = PISM.GridParameters(ctx.config)
13 params.ownership_ranges_from_options(ctx.size)
14 return PISM.IceGrid(ctx.ctx, params)
15
16
17 def setup():
18 global grid
19
20 # "inputs" does not own its elements, so we need to make sure
21 # these don't expire when setup() returns
22 global basal_melt_rate, ice_thickness, inputs, surface_temp
23 global climatic_mass_balance, basal_heat_flux, shelf_base_temp, cell_type
24 global u, v, w, strain_heating3
25
26 grid = create_dummy_grid()
27
28 zero = PISM.IceModelVec2S(grid, "zero", PISM.WITHOUT_GHOSTS)
29 zero.set(0.0)
30
31 cell_type = PISM.IceModelVec2CellType(grid, "mask", PISM.WITHOUT_GHOSTS)
32 cell_type.set(PISM.MASK_GROUNDED)
33
34 basal_heat_flux = PISM.IceModelVec2S(grid, "bheatflx", PISM.WITHOUT_GHOSTS)
35 basal_heat_flux.set(convert(10, "mW m-2", "W m-2"))
36
37 ice_thickness = PISM.model.createIceThicknessVec(grid)
38 ice_thickness.set(4000.0)
39 # TemperatureModel needs ice_thickness to set enthalpy in restart(...)
40 grid.variables().add(ice_thickness)
41
42 shelf_base_temp = PISM.IceModelVec2S(grid, "shelfbtemp", PISM.WITHOUT_GHOSTS)
43 shelf_base_temp.set(260.0)
44
45 surface_temp = PISM.IceModelVec2S(grid, "surface_temp", PISM.WITHOUT_GHOSTS)
46 surface_temp.set(260.0)
47
48 strain_heating3 = PISM.IceModelVec3(grid, "sigma", PISM.WITHOUT_GHOSTS)
49
50 u = PISM.IceModelVec3(grid, "u", PISM.WITHOUT_GHOSTS)
51
52 v = PISM.IceModelVec3(grid, "v", PISM.WITHOUT_GHOSTS)
53
54 w = PISM.IceModelVec3(grid, "w", PISM.WITHOUT_GHOSTS)
55
56 ice_thickness.set(4000.0)
57 u.set(0.0)
58 v.set(0.0)
59 w.set(0.0)
60
61 basal_melt_rate = zero
62 climatic_mass_balance = zero
63
64 inputs = PISM.EnergyModelInputs()
65
66 inputs.cell_type = cell_type
67 inputs.basal_frictional_heating = zero
68 inputs.basal_heat_flux = basal_heat_flux
69 inputs.ice_thickness = ice_thickness
70 inputs.surface_liquid_fraction = zero
71 inputs.shelf_base_temp = shelf_base_temp
72 inputs.surface_temp = surface_temp
73 inputs.till_water_thickness = zero
74 inputs.volumetric_heating_rate = strain_heating3
75 inputs.u3 = u
76 inputs.v3 = v
77 inputs.w3 = w
78
79 dt = convert(1, "years", "seconds")
80
81 def use_model(model):
82 print("* Performing a time step...")
83 model.update(0, dt, inputs)
84
85 try:
86 model.update(0, dt)
87 raise Exception("this should fail")
88 except TypeError:
89 pass
90
91 print(model.stdout_flags())
92 stats = model.stats()
93 enthalpy = model.enthalpy()
94 bmr = model.basal_melt_rate()
95
96
97 def initialize(model):
98 model.initialize(basal_melt_rate,
99 ice_thickness,
100 surface_temp,
101 climatic_mass_balance,
102 basal_heat_flux)
103
104
105 def test_interface():
106 "Use the EnergyModel interface to make sure the code runs."
107 for M in [PISM.EnthalpyModel, PISM.DummyEnergyModel, PISM.TemperatureModel]:
108
109 model = M(grid, None)
110
111 print("Testing %s..." % M)
112
113 print("* Bootstrapping using provided basal melt rate...")
114 initialize(model)
115
116 use_model(model)
117
118 try:
119 temperature = model.get_temperature()
120 except:
121 pass
122
123 pio = PISM.util.prepare_output("energy_model_state.nc")
124
125 print("* Saving the model state...")
126 model.write_model_state(pio)
127
128 print("* Restarting from a saved model state...")
129 model.restart(pio, 0)
130
131 print("* Bootstrapping from a saved model state...")
132 model.bootstrap(pio,
133 ice_thickness,
134 surface_temp,
135 climatic_mass_balance,
136 basal_heat_flux)
137
138 pio.close()
139
140
141 def test_temp_restart_from_enth():
142 enth_model = PISM.EnthalpyModel(grid, None)
143 temp_model = PISM.TemperatureModel(grid, None)
144
145 initialize(enth_model)
146
147 pio = PISM.util.prepare_output("enth_model_state.nc")
148 enth_model.write_model_state(pio)
149
150 temp_model.restart(pio, 0)
151
152
153 def test_enth_restart_from_temp():
154 enth_model = PISM.EnthalpyModel(grid, None)
155 temp_model = PISM.TemperatureModel(grid, None)
156
157 initialize(temp_model)
158
159 pio = PISM.util.prepare_output("temp_model_state.nc")
160 temp_model.write_model_state(pio)
161
162 enth_model.restart(pio, 0)
163
164
165 setup()
166
167 test_interface()
168 test_temp_restart_from_enth()
169 test_enth_restart_from_temp()