tatmosphere_models.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
---
tatmosphere_models.py (17706B)
---
1 #!/usr/bin/env python3
2 """
3 Tests of PISM's atmosphere models and modifiers.
4 """
5
6 import PISM
7 from PISM.testing import *
8 import os
9 import numpy as np
10 from unittest import TestCase
11
12 config = PISM.Context().config
13
14 # reduce the grid size to speed this up
15 config.set_number("grid.Mx", 3)
16 config.set_number("grid.My", 5) # non-square grid
17 config.set_number("grid.Mz", 2)
18
19 seconds_per_year = 365 * 86400
20 # ensure that this is the correct year length
21 config.set_string("time.calendar", "365_day")
22
23 # silence models' initialization messages
24 PISM.Context().log.set_threshold(1)
25
26 def write_state(model):
27 "Test writing of the model state"
28
29 o_filename = "atmosphere_model_state.nc"
30 o_diagnostics = "atmosphere_diagnostics.nc"
31
32 try:
33 output = PISM.util.prepare_output(o_filename)
34 model.define_model_state(output)
35 model.write_model_state(output)
36 output.close()
37
38 ds = model.diagnostics()
39 output = PISM.util.prepare_output(o_diagnostics)
40
41 for d in ds:
42 ds[d].define(output, PISM.PISM_DOUBLE)
43
44 for d in ds:
45 ds[d].compute().write(output)
46
47 output.close()
48
49 finally:
50 os.remove(o_filename)
51 os.remove(o_diagnostics)
52
53 def check_model(model, T, P, ts=None, Ts=None, Ps=None):
54 check(model.mean_annual_temp(), T)
55 check(model.mean_precipitation(), P)
56
57 model.init_timeseries(ts)
58
59 try:
60 model.begin_pointwise_access()
61 np.testing.assert_almost_equal(model.temp_time_series(0, 0), Ts)
62 np.testing.assert_almost_equal(model.precip_time_series(0, 0), Ps)
63 finally:
64 model.end_pointwise_access()
65
66 write_state(model)
67
68 model.max_timestep(ts[0])
69
70 def check_modifier(model, modifier, T=0.0, P=0.0, ts=None, Ts=None, Ps=None):
71 check_difference(modifier.mean_annual_temp(),
72 model.mean_annual_temp(),
73 T)
74 check_difference(modifier.mean_precipitation(),
75 model.mean_precipitation(),
76 P)
77
78 model.init_timeseries(ts)
79 modifier.init_timeseries(ts)
80
81 try:
82 model.begin_pointwise_access()
83 modifier.begin_pointwise_access()
84
85 Ts_model = np.array(model.temp_time_series(0, 0))
86 Ts_modifier = np.array(modifier.temp_time_series(0, 0))
87
88 Ps_model = np.array(model.precip_time_series(0, 0))
89 Ps_modifier = np.array(modifier.precip_time_series(0, 0))
90
91 np.testing.assert_almost_equal(Ts_modifier - Ts_model, Ts)
92 np.testing.assert_almost_equal(Ps_modifier - Ps_model, Ps)
93 finally:
94 modifier.end_pointwise_access()
95 model.end_pointwise_access()
96
97 write_state(modifier)
98
99 modifier.max_timestep(ts[0])
100
101 def precipitation(grid, value):
102 precip = PISM.IceModelVec2S(grid, "precipitation", PISM.WITHOUT_GHOSTS)
103 precip.set_attrs("climate", "precipitation", "kg m-2 s-1", "kg m-2 s-1",
104 "precipitation_flux", 0)
105 precip.set(value)
106 return precip
107
108 def air_temperature(grid, value):
109 temperature = PISM.IceModelVec2S(grid, "air_temp", PISM.WITHOUT_GHOSTS)
110 temperature.set_attrs("climate", "near-surface air temperature", "Kelvin", "Kelvin", "", 0)
111 temperature.set(value)
112 return temperature
113
114 class PIK(TestCase):
115 def setUp(self):
116 self.filename = "atmosphere_pik_input.nc"
117 self.grid = shallow_grid()
118 self.geometry = PISM.Geometry(self.grid)
119
120 self.geometry.latitude.set(-80.0)
121
122 self.P = 10.0 # this is very high, but that's fine
123
124 precipitation(self.grid, self.P).dump(self.filename)
125
126 config.set_string("atmosphere.pik.file", self.filename)
127
128 def test_atmosphere_pik(self):
129 "Model 'pik'"
130
131 parameterizations = {"martin" : (248.13, [248.13], [self.P]),
132 "huybrechts_dewolde" : (252.59, [237.973373], [self.P]),
133 "martin_huybrechts_dewolde" : (248.13, [233.51337298], [self.P]),
134 "era_interim" : (256.27, [243.0939774], [self.P]),
135 "era_interim_sin" : (255.31577, [241.7975841], [self.P]),
136 "era_interim_lon" : (248.886139, [233.3678998], [self.P])}
137
138 for p, (T, Ts, Ps) in parameterizations.items():
139 config.set_string("atmosphere.pik.parameterization", p)
140
141 model = PISM.AtmospherePIK(self.grid)
142 model.init(self.geometry)
143
144 # t and dt are irrelevant here
145 model.update(self.geometry, 0, 1)
146
147 check_model(model, T=T, P=self.P, ts=[0.5], Ts=Ts, Ps=Ps)
148
149 assert model.max_timestep(0).infinite()
150
151 try:
152 config.set_string("atmosphere.pik.parameterization", "invalid")
153 model = PISM.AtmospherePIK(self.grid)
154 assert False, "failed to catch an invalid parameterization"
155 except RuntimeError:
156 pass
157
158 def tearDown(self):
159 os.remove(self.filename)
160
161 class DeltaT(TestCase):
162 def setUp(self):
163 self.filename = "atmosphere_delta_T_input.nc"
164 self.grid = shallow_grid()
165 self.geometry = PISM.Geometry(self.grid)
166 self.geometry.ice_thickness.set(1000.0)
167 self.model = PISM.AtmosphereUniform(self.grid)
168 self.dT = -5.0
169
170 create_scalar_forcing(self.filename, "delta_T", "Kelvin", [self.dT], [0])
171
172 def tearDown(self):
173 os.remove(self.filename)
174
175 def test_atmosphere_delta_t(self):
176 "Modifier Delta_T"
177
178 modifier = PISM.AtmosphereDeltaT(self.grid, self.model)
179
180 config.set_string("atmosphere.delta_T.file", self.filename)
181
182 modifier.init(self.geometry)
183 modifier.update(self.geometry, 0, 1)
184
185 check_modifier(self.model, modifier, T=self.dT, ts=[0.5], Ts=[self.dT], Ps=[0])
186
187 class DeltaP(TestCase):
188 def setUp(self):
189 self.filename = "atmosphere_delta_P_input.nc"
190 self.grid = shallow_grid()
191 self.geometry = PISM.Geometry(self.grid)
192 self.geometry.ice_thickness.set(1000.0)
193 self.model = PISM.AtmosphereUniform(self.grid)
194 self.dP = 5.0
195
196 create_scalar_forcing(self.filename, "delta_P", "kg m-2 s-1", [self.dP], [0])
197
198 config.set_string("atmosphere.delta_P.file", self.filename)
199
200 def tearDown(self):
201 os.remove(self.filename)
202
203 def test_atmosphere_delta_p(self):
204 "Modifier 'delta_P'"
205
206 modifier = PISM.AtmosphereDeltaP(self.grid, self.model)
207
208 modifier.init(self.geometry)
209 modifier.update(self.geometry, 0, 1)
210
211 check_modifier(self.model, modifier, P=self.dP, ts=[0.5], Ts=[0], Ps=[self.dP])
212
213 class Given(TestCase):
214 def setUp(self):
215 self.filename = "atmosphere_given_input.nc"
216 self.grid = shallow_grid()
217 self.geometry = PISM.Geometry(self.grid)
218
219 self.P = 10.0
220 self.T = 250.0
221
222 output = PISM.util.prepare_output(self.filename)
223 precipitation(self.grid, self.P).write(output)
224 air_temperature(self.grid, self.T).write(output)
225
226 config.set_string("atmosphere.given.file", self.filename)
227
228 def tearDown(self):
229 os.remove(self.filename)
230
231 def test_atmosphere_given(self):
232 "Model 'given'"
233
234 model = PISM.AtmosphereFactory(self.grid).create("given")
235
236 model.init(self.geometry)
237 model.update(self.geometry, 0, 1)
238
239 check_model(model, T=self.T, P=self.P, ts=[0.5], Ts=[self.T], Ps=[self.P])
240
241 class SeaRISE(TestCase):
242 def setUp(self):
243 self.filename = "atmosphere_searise_input.nc"
244 self.grid = shallow_grid()
245
246 self.geometry = PISM.Geometry(self.grid)
247 self.geometry.latitude.set(70.0)
248 self.geometry.longitude.set(-45.0)
249 self.geometry.ice_thickness.set(2500.0)
250 self.geometry.ensure_consistency(0.0)
251
252 self.P = 10.0
253
254 output = PISM.util.prepare_output(self.filename)
255 precipitation(self.grid, self.P).write(output)
256
257 config.set_string("atmosphere.searise_greenland.file", self.filename)
258
259 def tearDown(self):
260 os.remove(self.filename)
261
262 def test_atmosphere_searise_greenland(self):
263 "Model 'searise_greenland'"
264
265 model = PISM.AtmosphereSeaRISEGreenland(self.grid)
266
267 model.init(self.geometry)
268
269 model.update(self.geometry, 0, 1)
270
271 check_model(model, P=self.P, T=251.9085, ts=[0.5], Ts=[238.66192632], Ps=[self.P])
272
273 class YearlyCycle(TestCase):
274 def setUp(self):
275 self.filename = "atmosphere_yearly_cycle.nc"
276 self.grid = shallow_grid()
277 self.geometry = PISM.Geometry(self.grid)
278
279 self.T_mean = 250.0
280 self.T_summer = 270.0
281 self.P = 15.0
282
283 output = PISM.util.prepare_output(self.filename)
284 precipitation(self.grid, self.P).write(output)
285
286 T_mean = PISM.IceModelVec2S(self.grid, "air_temp_mean_annual", PISM.WITHOUT_GHOSTS)
287 T_mean.set_attrs("climate", "mean annual near-surface air temperature", "K", "K", "", 0)
288 T_mean.set(self.T_mean)
289 T_mean.write(output)
290
291 T_summer = PISM.IceModelVec2S(self.grid, "air_temp_mean_summer", PISM.WITHOUT_GHOSTS)
292 T_summer.set_attrs("climate", "mean summer near-surface air temperature", "K", "K", "", 0)
293 T_summer.set(self.T_summer)
294 T_summer.write(output)
295
296 config.set_string("atmosphere.yearly_cycle.file", self.filename)
297
298 # FIXME: test "-atmosphere_yearly_cycle_scaling_file", too
299
300 def tearDown(self):
301 os.remove(self.filename)
302
303 def test_atmosphere_yearly_cycle(self):
304 "Model 'yearly_cycle'"
305
306 model = PISM.AtmosphereCosineYearlyCycle(self.grid)
307
308 model.init(self.geometry)
309
310 one_year = 365 * 86400.0
311 model.update(self.geometry, 0, one_year)
312
313 summer_peak = config.get_number("atmosphere.fausto_air_temp.summer_peak_day") / 365.0
314
315 ts = np.linspace(0, one_year, 13)
316 cycle = np.cos(2.0 * np.pi * (ts / one_year - summer_peak))
317 T = (self.T_summer - self.T_mean) * cycle + self.T_mean
318 P = np.zeros_like(T) + self.P
319
320 check_model(model, T=self.T_mean, P=self.P, ts=ts, Ts=T, Ps=P)
321
322 class OneStation(TestCase):
323 def setUp(self):
324 self.filename = "atmosphere_one_station.nc"
325 self.grid = shallow_grid()
326 self.geometry = PISM.Geometry(self.grid)
327 self.T = 263.15
328 self.P = 10.0
329
330 time_name = config.get_string("time.dimension_name")
331
332 output = PISM.util.prepare_output(self.filename, append_time=True)
333
334 output.redef()
335 output.define_variable("precipitation", PISM.PISM_DOUBLE, [time_name])
336 output.write_attribute("precipitation", "units", "kg m-2 s-1")
337
338 output.define_variable("air_temp", PISM.PISM_DOUBLE, [time_name])
339 output.write_attribute("air_temp", "units", "Kelvin")
340
341 output.write_variable("precipitation", [0], [1], [self.P])
342 output.write_variable("air_temp", [0], [1], [self.T])
343
344 output.close()
345
346 config.set_string("atmosphere.one_station.file", self.filename)
347
348 def tearDown(self):
349 os.remove(self.filename)
350
351 def test_atmosphere_one_station(self):
352 "Model 'weather_station'"
353
354 model = PISM.AtmosphereWeatherStation(self.grid)
355
356 model.init(self.geometry)
357
358 model.update(self.geometry, 0, 1)
359
360 check_model(model, P=self.P, T=self.T, ts=[0.5], Ts=[self.T], Ps=[self.P])
361
362 class Uniform(TestCase):
363 def setUp(self):
364 self.grid = shallow_grid()
365 self.geometry = PISM.Geometry(self.grid)
366 self.P = 5.0
367 self.T = 250.0
368
369 config.set_number("atmosphere.uniform.temperature", self.T)
370 config.set_number("atmosphere.uniform.precipitation", self.P)
371
372 def test_atmosphere_uniform(self):
373 "Model 'uniform'"
374 model = PISM.AtmosphereUniform(self.grid)
375
376 model.init(self.geometry)
377
378 model.update(self.geometry, 0, 1)
379
380 P = PISM.util.convert(self.P, "kg m-2 year-1", "kg m-2 s-1")
381 check_model(model, T=self.T, P=P, ts=[0.5], Ts=[self.T], Ps=[P])
382
383 class Anomaly(TestCase):
384 def setUp(self):
385 self.filename = "atmosphere_anomaly_input.nc"
386 self.grid = shallow_grid()
387 self.geometry = PISM.Geometry(self.grid)
388 self.geometry.ice_thickness.set(1000.0)
389 self.model = PISM.AtmosphereUniform(self.grid)
390 self.dT = -5.0
391 self.dP = 20.0
392
393 dT = PISM.IceModelVec2S(self.grid, "air_temp_anomaly", PISM.WITHOUT_GHOSTS)
394 dT.set_attrs("climate", "air temperature anomaly", "Kelvin", "Kelvin", "", 0)
395 dT.set(self.dT)
396
397 dP = PISM.IceModelVec2S(self.grid, "precipitation_anomaly", PISM.WITHOUT_GHOSTS)
398 dP.set_attrs("climate", "precipitation anomaly", "kg m-2 s-1", "kg m-2 s-1", "", 0)
399 dP.set(self.dP)
400
401 output = PISM.util.prepare_output(self.filename)
402 dT.write(output)
403 dP.write(output)
404
405 config.set_string("atmosphere.anomaly.file", self.filename)
406
407 def tearDown(self):
408 os.remove(self.filename)
409
410 def test_atmosphere_anomaly(self):
411 "Modifier 'anomaly'"
412
413 modifier = PISM.AtmosphereAnomaly(self.grid, self.model)
414
415 modifier.init(self.geometry)
416
417 modifier.update(self.geometry, 0, 1)
418
419 check_modifier(self.model, modifier, T=self.dT, P=self.dP,
420 ts=[0.5], Ts=[self.dT], Ps=[self.dP])
421
422 class PrecipScaling(TestCase):
423 def setUp(self):
424 self.filename = "atmosphere_precip_scaling_input.nc"
425 self.grid = shallow_grid()
426 self.geometry = PISM.Geometry(self.grid)
427 self.geometry.ice_thickness.set(1000.0)
428 self.model = PISM.AtmosphereUniform(self.grid)
429 self.dT = 5.0
430
431 create_scalar_forcing(self.filename, "delta_T", "Kelvin", [self.dT], [0])
432
433 config.set_string("atmosphere.precip_scaling.file", self.filename)
434
435 def tearDown(self):
436 os.remove(self.filename)
437
438 def test_atmosphere_precip_scaling(self):
439 "Modifier 'precip_scaling'"
440
441 modifier = PISM.AtmospherePrecipScaling(self.grid, self.model)
442
443 modifier.init(self.geometry)
444
445 modifier.update(self.geometry, 0, 1)
446
447 check_modifier(self.model, modifier, P=1.3373514942327523e-05,
448 ts=[0.5], Ts=[0], Ps=[1.33735149e-05])
449
450 class FracP(TestCase):
451 def setUp(self):
452 self.filename = "atmosphere_frac_P_input.nc"
453 self.grid = shallow_grid()
454 self.geometry = PISM.Geometry(self.grid)
455 self.geometry.ice_thickness.set(1000.0)
456 self.model = PISM.AtmosphereUniform(self.grid)
457 self.P_ratio = 5.0
458
459 create_scalar_forcing(self.filename, "frac_P", "1", [self.P_ratio], [0])
460
461 config.set_string("atmosphere.frac_P.file", self.filename)
462
463 def tearDown(self):
464 os.remove(self.filename)
465
466 def test_atmosphere_frac_p(self):
467 "Modifier 'frac_P'"
468
469 modifier = PISM.AtmosphereFracP(self.grid, self.model)
470
471 modifier.init(self.geometry)
472
473 modifier.update(self.geometry, 0, 1)
474
475 check_ratio(modifier.mean_precipitation(), self.model.mean_precipitation(),
476 self.P_ratio)
477
478 check_modifier(self.model, modifier, T=0, P=0.00012675505856327396,
479 ts=[0.5], Ts=[0], Ps=[0.00012676])
480
481
482 class ElevationChange(TestCase):
483 def setUp(self):
484 self.filename = "atmosphere_reference_surface.nc"
485 self.grid = shallow_grid()
486 self.dTdz = 1.0 # Kelvin per km
487 self.dPdz = 1000.0 # (kg/m^2)/year per km
488 self.dz = 1000.0 # m
489 self.dT = -self.dTdz * self.dz / 1000.0
490 self.dP = -PISM.util.convert(self.dPdz * self.dz / 1000.0, "kg m-2 year-1", "kg m-2 s-1")
491
492 self.geometry = PISM.Geometry(self.grid)
493
494 # save current surface elevation to use it as a "reference" surface elevation
495 self.geometry.ice_surface_elevation.dump(self.filename)
496
497 config.set_string("atmosphere.elevation_change.file", self.filename)
498
499 config.set_number("atmosphere.elevation_change.precipitation.lapse_rate", self.dPdz)
500
501 config.set_number("atmosphere.elevation_change.temperature_lapse_rate", self.dTdz)
502
503 def tearDown(self):
504 os.remove(self.filename)
505
506 def test_atmosphere_elevation_change_shift(self):
507 "Modifier 'elevation_change': lapse rate"
508
509 config.set_string("atmosphere.elevation_change.precipitation.method", "shift")
510
511 model = PISM.AtmosphereUniform(self.grid)
512 modifier = PISM.AtmosphereElevationChange(self.grid, model)
513
514 modifier.init(self.geometry)
515
516 # change surface elevation
517 self.geometry.ice_surface_elevation.shift(self.dz)
518
519 # check that the temperature changed accordingly
520 modifier.update(self.geometry, 0, 1)
521 check_modifier(model, modifier, T=self.dT, P=self.dP,
522 ts=[0.5], Ts=[self.dT], Ps=[self.dP])
523
524 def test_atmosphere_elevation_change_scale(self):
525 "Modifier 'elevation_change': scaling"
526
527 config.set_string("atmosphere.elevation_change.precipitation.method", "scale")
528
529 model = PISM.AtmosphereUniform(self.grid)
530 modifier = PISM.AtmosphereElevationChange(self.grid, model)
531
532 modifier.init(self.geometry)
533
534 # change surface elevation
535 self.geometry.ice_surface_elevation.shift(self.dz)
536
537 # check that the temperature changed accordingly
538 modifier.update(self.geometry, 0, 1)
539
540 C = config.get_number("atmosphere.precip_exponential_factor_for_temperature")
541 P = sample(model.mean_precipitation())
542 dP = np.exp(C * self.dT) * P - P
543
544 check_modifier(model, modifier, T=self.dT, P=dP,
545 ts=[0.5], Ts=[self.dT], Ps=[dP])