tsurface_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
---
tsurface_models.py (26279B)
---
1 #!/usr/bin/env python3
2 """
3 Tests of PISM's surface models and modifiers.
4 """
5
6 import PISM
7 from PISM.testing import *
8 import sys
9 import os
10 import numpy as np
11 from unittest import TestCase, SkipTest
12
13 from PISM.util import convert
14
15 config = PISM.Context().config
16
17 seconds_per_year = 365 * 86400
18 # ensure that this is the correct year length
19 config.set_string("time.calendar", "365_day")
20
21 log = PISM.Context().log
22 # silence models' initialization messages
23 log.set_threshold(1)
24
25 options = PISM.PETSc.Options()
26
27 def write_state(model, filename):
28 "Write the state of the model to a file"
29
30 PISM.util.prepare_output(filename)
31 f = PISM.File(model.grid().ctx().com(),
32 filename,
33 PISM.PISM_NETCDF3,
34 PISM.PISM_READWRITE)
35 model.define_model_state(f)
36 model.write_model_state(f)
37
38 diags = model.diagnostics()
39 for k in diags.keys():
40 diags[k].compute().write(f)
41
42 f.close()
43
44 def probe_interface(model):
45 """Prove the interface of a surface model to check if its public methods run successfully."""
46 model.accumulation()
47 model.layer_mass()
48 model.layer_thickness()
49 model.liquid_water_fraction()
50 model.mass_flux()
51 model.melt()
52 model.runoff()
53 model.temperature()
54
55 model.max_timestep(0)
56
57 model.diagnostics()
58
59 # FIXME: this causes a memory leak
60 # model.ts_diagnostics()
61
62 model.grid()
63
64 def check_model(model, T, omega, SMB,
65 mass=0.0, thickness=0.0, accumulation=0.0, melt=0.0, runoff=0.0):
66 check(model.mass_flux(), SMB)
67 check(model.temperature(), T)
68 check(model.liquid_water_fraction(), omega)
69 check(model.layer_mass(), mass)
70 check(model.layer_thickness(), thickness)
71 check(model.accumulation(), accumulation)
72 check(model.melt(), melt)
73 check(model.runoff(), runoff)
74
75 def check_modifier(model, modifier,
76 T=0.0, omega=0.0, SMB=0.0, mass=0.0, thickness=0.0,
77 accumulation=0.0, melt=0.0, runoff=0.0):
78 check_difference(modifier.mass_flux(),
79 model.mass_flux(),
80 SMB)
81
82 check_difference(modifier.temperature(),
83 model.temperature(),
84 T)
85
86 check_difference(modifier.liquid_water_fraction(),
87 model.liquid_water_fraction(),
88 omega)
89
90 check_difference(modifier.layer_mass(),
91 model.layer_mass(),
92 mass)
93
94 check_difference(modifier.layer_thickness(),
95 model.layer_thickness(),
96 thickness)
97
98 check_difference(modifier.accumulation(),
99 model.accumulation(),
100 accumulation)
101
102 check_difference(modifier.melt(),
103 model.melt(),
104 melt)
105
106 check_difference(modifier.runoff(),
107 model.runoff(),
108 runoff)
109
110 def surface_simple(grid):
111 return PISM.SurfaceSimple(grid, PISM.AtmosphereUniform(grid))
112
113 def climatic_mass_balance(grid, value):
114 SMB = PISM.IceModelVec2S(grid, "climatic_mass_balance", PISM.WITHOUT_GHOSTS)
115 SMB.set_attrs("climate", "surface mass balance", "kg m-2 s-1", "kg m-2 s-1",
116 "land_ice_surface_specific_mass_balance_flux", 0)
117 SMB.set(value)
118 return SMB
119
120 def ice_surface_temp(grid, value):
121 temperature = PISM.IceModelVec2S(grid, "ice_surface_temp", PISM.WITHOUT_GHOSTS)
122 temperature.set_attrs("climate", "ice temperature at the top surface", "Kelvin", "Kelvin", "", 0)
123 temperature.set(value)
124 return temperature
125
126 class Given(TestCase):
127 def setUp(self):
128 self.filename = "surface_given_input.nc"
129 self.output_filename = "surface_given_output.nc"
130 self.grid = shallow_grid()
131 self.geometry = PISM.Geometry(self.grid)
132
133 self.T = 272.15
134 self.M = 1001.0
135
136 output = PISM.util.prepare_output(self.filename)
137 ice_surface_temp(self.grid, self.T).write(output)
138 climatic_mass_balance(self.grid, self.M).write(output)
139 output.close()
140
141 def test_surface_given(self):
142 "Model 'given'"
143 atmosphere = PISM.AtmosphereUniform(self.grid)
144
145 config.set_string("surface.given.file", self.filename)
146
147 model = PISM.SurfaceGiven(self.grid, atmosphere)
148
149 model.init(self.geometry)
150
151 model.update(self.geometry, 0, 1)
152
153 check_model(model, self.T, 0.0, self.M, accumulation=self.M)
154
155 write_state(model, self.output_filename)
156 probe_interface(model)
157
158 def tearDown(self):
159 os.remove(self.filename)
160 os.remove(self.output_filename)
161
162 class DeltaT(TestCase):
163 def setUp(self):
164 self.filename = "surface_delta_T_input.nc"
165 self.output_filename = "surface_delta_T_output.nc"
166 self.grid = shallow_grid()
167 self.model = surface_simple(self.grid)
168 self.dT = -5.0
169 self.geometry = PISM.Geometry(self.grid)
170
171 create_scalar_forcing(self.filename, "delta_T", "Kelvin", [self.dT], [0])
172
173 def test_surface_delta_t(self):
174 "Modifier 'delta_T'"
175
176 modifier = PISM.SurfaceDeltaT(self.grid, self.model)
177
178 config.set_string("surface.delta_T.file", self.filename)
179
180 modifier.init(self.geometry)
181 modifier.update(self.geometry, 0, 1)
182
183 check_modifier(self.model, modifier, T=self.dT)
184
185 write_state(modifier, self.output_filename)
186 probe_interface(modifier)
187
188 def tearDown(self):
189 os.remove(self.filename)
190 os.remove(self.output_filename)
191
192 class ElevationChange(TestCase):
193 def setUp(self):
194 self.filename = "surface_reference_surface.nc"
195 self.output_filename = "surface_lapse_rates_output.nc"
196 self.grid = shallow_grid()
197 self.dTdz = 1.0 # 1 Kelvin per km
198 self.dSMBdz = 2.0 # m year-1 per km
199 self.dz = 1500.0 # m
200
201 self.geometry = PISM.Geometry(self.grid)
202
203 # save current surface elevation to use it as a "reference" surface elevation
204 self.geometry.ice_surface_elevation.dump(self.filename)
205
206 config.set_string("surface.elevation_change.file", self.filename)
207 config.set_number("surface.elevation_change.temperature_lapse_rate", self.dTdz)
208 config.set_number("surface.elevation_change.smb.lapse_rate", self.dSMBdz)
209
210 self.dT = self.dz * convert(-self.dTdz, "Kelvin / km", "Kelvin / m")
211
212 def test_elevation_change_shift(self):
213 "Modifier elevation_change: shift"
214
215 config.set_string("surface.elevation_change.smb.method", "shift")
216
217 model = surface_simple(self.grid)
218 modifier = PISM.SurfaceElevationChange(self.grid, model)
219
220 modifier.init(self.geometry)
221
222 # change surface elevation
223 self.geometry.ice_surface_elevation.shift(self.dz)
224
225 # check changes in outputs
226 modifier.update(self.geometry, 0, 1)
227
228 ice_density = config.get_number("constants.ice.density")
229
230 dSMB = self.dz * ice_density * convert(-self.dSMBdz,
231 "kg m-2 year-1 / km",
232 "kg m-2 s-1 / m")
233 # shifting SMB by dSMB reduced accumulation to zero
234 dA = 0.0 - sample(model.accumulation())
235 dM = dA - dSMB
236 dR = dM
237
238 check_modifier(model, modifier, T=self.dT, SMB=dSMB,
239 accumulation=dA, melt=dM, runoff=dR)
240
241 write_state(modifier, self.output_filename)
242 probe_interface(modifier)
243
244 def test_elevation_change_scale(self):
245 "Modifier elevation_change: scale"
246
247 config.set_string("surface.elevation_change.smb.method", "scale")
248 config.set_number("surface.elevation_change.smb.exp_factor", 1.0)
249
250 model = surface_simple(self.grid)
251 modifier = PISM.SurfaceElevationChange(self.grid, model)
252
253 modifier.init(self.geometry)
254
255 # change surface elevation
256 self.geometry.ice_surface_elevation.shift(self.dz)
257
258 # check changes in outputs
259 modifier.update(self.geometry, 0, 1)
260
261 SMB = sample(model.mass_flux())
262 C = config.get_number("surface.elevation_change.smb.exp_factor")
263
264 dSMB = np.exp(C * self.dT) * SMB - SMB
265 # SMB stays positive, so all SMB corresponds to accumulation
266 dA = dSMB
267 dM = dA - dSMB
268 dR = dM
269
270 check_modifier(model, modifier, T=self.dT, SMB=dSMB,
271 accumulation=dA, melt=dM, runoff=dR)
272
273 write_state(modifier, self.output_filename)
274 probe_interface(modifier)
275
276 def tearDown(self):
277 os.remove(self.filename)
278 os.remove(self.output_filename)
279
280 class Elevation(TestCase):
281 def setUp(self):
282 self.grid = shallow_grid()
283 self.geometry = PISM.Geometry(self.grid)
284 self.output_filename = "surface_elevation_output.nc"
285
286 # change geometry just to make this a bit more interesting
287 self.geometry.ice_thickness.set(1000.0)
288 self.geometry.ensure_consistency(0.0)
289
290 def elevation_1_test(self):
291 "Model 'elevation', test 1"
292 model = PISM.SurfaceElevation(self.grid, PISM.AtmosphereUniform(self.grid))
293
294 model.init(self.geometry)
295
296 model.update(self.geometry, 0, 1)
297
298 T = 268.15
299 omega = 0.0
300 SMB = -8.651032746943449e-05
301 accumulation = 0.0
302 melt = -SMB
303 runoff = melt
304
305 check_model(model, T, omega, SMB,
306 accumulation=accumulation, melt=melt, runoff=runoff)
307
308 probe_interface(model)
309
310 def elevation_2_test(self):
311 "Model 'elevation', test 2"
312 T_min = -5.0
313 T_max = 0.0
314 z_min = 1000.0
315 z_ela = 1100.0
316 z_max = 1500.0
317 M_min = -1.0
318 M_max = 5.0
319
320 self.geometry.ice_thickness.set(0.5 * (z_min + z_max))
321 self.geometry.ensure_consistency(0.0)
322
323 options.setValue("-ice_surface_temp", "{},{},{},{}".format(T_min, T_max, z_min, z_max))
324 options.setValue("-climatic_mass_balance",
325 "{},{},{},{},{}".format(M_min, M_max, z_min, z_ela, z_max))
326
327 T = PISM.util.convert(0.5 * (T_min + T_max), "Celsius", "Kelvin")
328 SMB = PISM.util.convert(1.87504, "m/year", "m/s") * config.get_number("constants.ice.density")
329
330 model = PISM.SurfaceElevation(self.grid, PISM.AtmosphereUniform(self.grid))
331
332 model.init(self.geometry)
333
334 model.update(self.geometry, 0, 1)
335
336 check_model(model, T=T, SMB=SMB, omega=0, mass=0, thickness=0, accumulation=SMB)
337
338 class TemperatureIndex1(TestCase):
339 def setUp(self):
340 self.grid = shallow_grid()
341 self.geometry = PISM.Geometry(self.grid)
342 self.atmosphere = PISM.AtmosphereUniform(self.grid)
343 self.output_filename = "surface_pdd_output.nc"
344
345 def pdd_test(self):
346 "Model 'pdd', test 1"
347 config.set_string("surface.pdd.method", "expectation_integral")
348
349 model = PISM.SurfaceTemperatureIndex(self.grid, self.atmosphere)
350
351 model.init(self.geometry)
352
353 model.update(self.geometry, 0, 1)
354
355 T = config.get_number("atmosphere.uniform.temperature", "Kelvin")
356 omega = 0.0
357
358 accumulation = config.get_number("atmosphere.uniform.precipitation", "kg m-2 second-1")
359 melt = accumulation
360 runoff = melt * (1.0 - config.get_number("surface.pdd.refreeze"))
361 SMB = accumulation - runoff
362
363 check_model(model, T, omega, SMB, accumulation=accumulation, melt=melt,
364 runoff=runoff)
365
366 write_state(model, self.output_filename)
367 probe_interface(model)
368
369 def tearDown(self):
370 os.remove(self.output_filename)
371
372 class TemperatureIndex2(TestCase):
373 def setUp(self):
374 self.air_temp = config.get_number("atmosphere.uniform.temperature")
375 self.precip = config.get_number("atmosphere.uniform.precipitation")
376
377 self.grid = shallow_grid()
378
379 self.geometry = PISM.Geometry(self.grid)
380 # make sure that there's ice to melt
381 self.geometry.ice_thickness.set(1000.0)
382
383 T_above_zero = 1
384 dt_days = 5
385 self.T = 273.15 + T_above_zero
386 self.dt = dt_days * 86400
387
388 ice_density = config.get_number("constants.ice.density")
389 beta_ice = config.get_number("surface.pdd.factor_ice")
390 refreeze_fraction = config.get_number("surface.pdd.refreeze")
391 PDD = dt_days * T_above_zero
392 ice_melted = PDD * beta_ice
393 refreeze = ice_melted * refreeze_fraction
394 self.SMB = -(ice_melted - refreeze) * ice_density / self.dt
395
396 config.set_number("atmosphere.uniform.temperature", self.T)
397 # disable daily variability so that we can compute the number of PDDs exactly
398 config.set_number("surface.pdd.std_dev", 0.0)
399 # no precipitation
400 config.set_number("atmosphere.uniform.precipitation", 0)
401
402 def tearDown(self):
403 config.set_number("atmosphere.uniform.temperature", self.air_temp)
404 config.set_number("atmosphere.uniform.precipitation", self.precip)
405
406 def test_surface_pdd(self):
407 "Model 'pdd'"
408 model = PISM.SurfaceTemperatureIndex(self.grid, PISM.AtmosphereUniform(self.grid))
409
410 model.init(self.geometry)
411
412 model.update(self.geometry, 0, self.dt)
413
414 check_model(model, T=self.T, SMB=self.SMB, omega=0.0, mass=0.0, thickness=0.0,
415 melt=40, runoff=16)
416
417 class PIK(TestCase):
418 def setUp(self):
419 self.filename = "surface_pik_input.nc"
420 self.output_filename = "surface_pik_output.nc"
421 self.grid = shallow_grid()
422 self.geometry = PISM.Geometry(self.grid)
423
424 self.M = 1001.0
425 self.T = 233.13
426
427 self.geometry.latitude.set(-80.0)
428 self.geometry.ice_thickness.set(2000.0)
429 self.geometry.ensure_consistency(0.0)
430
431 climatic_mass_balance(self.grid, self.M).dump(self.filename)
432
433 config.set_string("input.file", self.filename)
434
435 def surface_pik_test(self):
436 "Model 'pik'"
437 model = PISM.SurfacePIK(self.grid, PISM.AtmosphereUniform(self.grid))
438
439 model.init(self.geometry)
440
441 model.update(self.geometry, 0, 1)
442
443 check_model(model, self.T, 0.0, self.M, accumulation=self.M)
444
445 write_state(model, self.output_filename)
446 probe_interface(model)
447
448 def tearDown(self):
449 os.remove(self.filename)
450 os.remove(self.output_filename)
451 config.set_string("input.file", "")
452
453 class Simple(TestCase):
454 def setUp(self):
455 self.grid = shallow_grid()
456 self.output_filename = "surface_simple_output.nc"
457 self.atmosphere = PISM.AtmosphereUniform(self.grid)
458 self.geometry = PISM.Geometry(self.grid)
459
460 def simple_test(self):
461 "Model 'simple'"
462 atmosphere = self.atmosphere
463
464 model = PISM.SurfaceSimple(self.grid, atmosphere)
465
466 model.init(self.geometry)
467
468 model.update(self.geometry, 0, 1)
469
470 T = sample(atmosphere.mean_annual_temp())
471 M = sample(atmosphere.mean_precipitation())
472
473 check_model(model, T, 0.0, M, accumulation=M)
474
475 write_state(model, self.output_filename)
476 probe_interface(model)
477
478 def tearDown(self):
479 os.remove(self.output_filename)
480
481 class Anomaly(TestCase):
482 def setUp(self):
483 self.filename = "surface_anomaly_input.nc"
484 self.output_filename = "surface_anomaly_output.nc"
485 self.grid = shallow_grid()
486 self.geometry = PISM.Geometry(self.grid)
487 self.model = surface_simple(self.grid)
488 self.dSMB = -(config.get_number("atmosphere.uniform.precipitation", "kg m-2 s-1") + 5.0)
489 self.dT = 2.0
490
491 PISM.util.prepare_output(self.filename)
492
493 delta_SMB = PISM.IceModelVec2S(self.grid, "climatic_mass_balance_anomaly",
494 PISM.WITHOUT_GHOSTS)
495 delta_SMB.set_attrs("climate_forcing",
496 "2D surface mass flux anomaly", "kg m-2 s-1", "kg m-2 s-1", "", 0)
497 delta_SMB.set(self.dSMB)
498
499 delta_SMB.write(self.filename)
500
501 delta_T = PISM.IceModelVec2S(self.grid, "ice_surface_temp_anomaly",
502 PISM.WITHOUT_GHOSTS)
503 delta_T.set_attrs("climate_forcing",
504 "2D surface temperature anomaly", "Kelvin", "Kelvin", "", 0)
505 delta_T.set(self.dT)
506
507 delta_T.write(self.filename)
508
509 def anomaly_test(self):
510 "Modifier 'anomaly'"
511
512 config.set_string("surface.anomaly.file", self.filename)
513
514 modifier = PISM.SurfaceAnomaly(self.grid, self.model)
515
516 modifier.init(self.geometry)
517 modifier.update(self.geometry, 0, 1)
518
519 # once anomaly is applied the SMB is negative, so the new accumulation is zero
520 dA = 0.0 - sample(self.model.accumulation())
521 dM = dA - self.dSMB
522 dR = dM
523
524 check_modifier(self.model, modifier, T=self.dT, SMB=self.dSMB,
525 accumulation=dA, melt=dM, runoff=dR)
526
527 write_state(modifier, self.output_filename)
528 probe_interface(modifier)
529
530 def tearDown(self):
531 os.remove(self.filename)
532 os.remove(self.output_filename)
533
534 class Cache(TestCase):
535 def setUp(self):
536 self.filename = "surface_dT.nc"
537 self.output_filename = "surface_cache_output.nc"
538 self.grid = shallow_grid()
539 self.geometry = PISM.Geometry(self.grid)
540
541 self.simple = surface_simple(self.grid)
542 self.delta_T = PISM.SurfaceDeltaT(self.grid, self.simple)
543
544 time_bounds = np.array([0, 1, 1, 2, 2, 3, 3, 4]) * seconds_per_year
545 create_scalar_forcing(self.filename, "delta_T", "Kelvin", [1, 2, 3, 4],
546 times=None, time_bounds=time_bounds)
547
548 config.set_string("surface.delta_T.file", self.filename)
549
550 config.set_number("surface.cache.update_interval", 2.0)
551
552 def test_surface_cache(self):
553 "Modifier 'cache'"
554
555 modifier = PISM.SurfaceCache(self.grid, self.delta_T)
556
557 modifier.init(self.geometry)
558
559 dt = seconds_per_year
560
561 N = 4
562 ts = np.arange(float(N)) * dt
563 diff = []
564 for t in ts:
565 modifier.update(self.geometry, t, dt)
566
567 original = sample(self.simple.temperature())
568 cached = sample(modifier.temperature())
569
570 diff.append(cached - original)
571
572 np.testing.assert_almost_equal(diff, [1, 1, 3, 3])
573
574 write_state(modifier, self.output_filename)
575 probe_interface(modifier)
576
577 def tearDown(self):
578 os.remove(self.filename)
579 os.remove(self.output_filename)
580
581 class ForceThickness(TestCase):
582 def setUp(self):
583 self.grid = shallow_grid()
584 self.geometry = PISM.Geometry(self.grid)
585 self.model = surface_simple(self.grid)
586 self.filename = "surface_force_to_thickness_input.nc"
587 self.output_filename = "surface_force_to_thickness_output.nc"
588
589 self.H = 1000.0
590 self.dH = 1000.0
591
592 self.geometry.ice_thickness.set(self.H)
593
594 # save ice thickness to a file to use as the target thickness
595 PISM.util.prepare_output(self.filename)
596 self.geometry.ice_thickness.write(self.filename)
597
598 ftt_mask = PISM.IceModelVec2S(self.grid, "ftt_mask", PISM.WITHOUT_GHOSTS)
599 ftt_mask.set(1.0)
600 ftt_mask.write(self.filename)
601
602 alpha = 10.0
603 ice_density = config.get_number("constants.ice.density")
604 self.dSMB = -ice_density * alpha * self.dH
605
606 config.set_string("surface.force_to_thickness_file", self.filename)
607 config.set_number("surface.force_to_thickness.alpha", convert(alpha, "1/s", "1/year"))
608
609 def forcing_test(self):
610 "Modifier ForceThickness"
611 modifier = PISM.SurfaceForceThickness(self.grid, self.model)
612
613 modifier.init(self.geometry)
614
615 self.geometry.ice_thickness.set(self.H + self.dH)
616
617 modifier.update(self.geometry, 0, 1)
618
619 dA = 0.0 - sample(self.model.accumulation())
620 dM = dA - self.dSMB
621 dR = dM
622
623 check_modifier(self.model, modifier, SMB=self.dSMB,
624 accumulation=dA, melt=dM, runoff=dR)
625
626 write_state(modifier, self.output_filename)
627 probe_interface(modifier)
628
629 def tearDown(self):
630 os.remove(self.filename)
631 os.remove(self.output_filename)
632
633 class EISMINTII(TestCase):
634 def setUp(self):
635 self.grid = shallow_grid()
636 self.geometry = PISM.Geometry(self.grid)
637 self.output_filename = "surface_eismint_output.nc"
638
639 def eismintii_test(self):
640 "Model EISMINTII: define and write model state; get diagnostics"
641
642 for experiment in "ABCDEFGHIJKL":
643 model = PISM.SurfaceEISMINTII(self.grid, ord(experiment))
644
645 model.init(self.geometry)
646
647 model.update(self.geometry, 0, 1)
648
649 write_state(model, self.output_filename)
650 probe_interface(model)
651
652 os.remove(self.output_filename)
653
654 def tearDown(self):
655 try:
656 os.remove(self.output_filename)
657 except:
658 pass
659
660 class Initialization(TestCase):
661 def setUp(self):
662 self.grid = shallow_grid()
663 self.geometry = PISM.Geometry(self.grid)
664 self.output_filename = "surface_init_output.nc"
665 self.model = surface_simple(self.grid)
666
667 def initialization_test(self):
668 "Modifier InitializationHelper"
669
670 modifier = PISM.SurfaceInitialization(self.grid, self.model)
671
672 modifier.init(self.geometry)
673
674 modifier.update(self.geometry, 0, 1)
675
676 write_state(modifier, self.output_filename)
677 probe_interface(modifier)
678
679 def tearDown(self):
680 os.remove(self.output_filename)
681
682 class Factory(TestCase):
683 def setUp(self):
684 self.grid = shallow_grid()
685 self.geometry = PISM.Geometry(self.grid)
686
687 def factory_test(self):
688 "Surface model factory"
689 atmosphere = PISM.AtmosphereUniform(self.grid)
690
691 factory = PISM.SurfaceFactory(self.grid, atmosphere)
692
693 simple = factory.create("simple")
694
695 model = factory.create("simple,cache")
696
697 try:
698 factory.create("invalid_model")
699 return False
700 except RuntimeError:
701 pass
702
703 try:
704 factory.create("simple,invalid_modifier")
705 return False
706 except RuntimeError:
707 pass
708
709 def tearDown(self):
710 pass
711
712
713 class ISMIP6(TestCase):
714 def prepare_reference_data(self, grid, filename):
715
716 usurf = PISM.model.createIceSurfaceVec(grid)
717 SMB_ref = PISM.IceModelVec2S(grid, "climatic_mass_balance", PISM.WITHOUT_GHOSTS)
718 T_ref = PISM.IceModelVec2S(grid, "ice_surface_temp", PISM.WITHOUT_GHOSTS)
719
720 usurf.metadata(0).set_string("units", "m")
721
722 SMB_ref.set_attrs("climate_forcing", "reference SMB", "kg m-2 s-1", "kg m-2 s-1",
723 "land_ice_surface_specific_mass_balance_flux", 0)
724
725 T_ref.metadata(0).set_string("units", "K")
726
727 out = PISM.util.prepare_output(filename, append_time=True)
728
729 T_ref.set(260.0)
730 SMB_ref.set(0.0)
731 usurf.set(0.0)
732
733 # write time-independent fields
734 for v in [usurf, SMB_ref, T_ref]:
735 v.set_time_independent(True)
736 v.write(out)
737
738 out.close()
739
740 def prepare_climate_forcing(self, grid, filename):
741
742 aSMB = PISM.IceModelVec2S(grid, "climatic_mass_balance_anomaly", PISM.WITHOUT_GHOSTS)
743 aSMB.set_attrs("climate_forcing", "SMB anomaly", "kg m-2 s-1", "kg m-2 s-1", "", 0)
744
745 dSMBdz = PISM.IceModelVec2S(grid, "climatic_mass_balance_gradient", PISM.WITHOUT_GHOSTS)
746 dSMBdz.set_attrs("climate_forcing", "SMB gradient", "kg m-2 s-1 m-1", "kg m-2 s-1 m-1", "", 0)
747
748 aT = PISM.IceModelVec2S(grid, "ice_surface_temp_anomaly", PISM.WITHOUT_GHOSTS)
749 aT.set_attrs("climate_forcing", "temperature anomaly", "Kelvin", "Kelvin", "", 0)
750
751 dTdz = PISM.IceModelVec2S(grid, "ice_surface_temp_gradient", PISM.WITHOUT_GHOSTS)
752 dTdz.set_attrs("climate_forcing", "surface temperature gradient", "K m-1", "K m-1", "", 0)
753
754 out = PISM.util.prepare_output(filename, append_time=False)
755
756 bounds = PISM.TimeBoundsMetadata("time_bounds", "time", self.ctx.unit_system)
757
758 SMB_anomaly = 1.0
759 T_anomaly = 1.0
760 SMB_gradient = 1.0
761 T_gradient = 1.0
762
763 # monthly steps
764 dt = (365 * 86400) / 12.0
765
766 for j in range(12):
767
768 t = self.ctx.time.current() + j * dt
769
770 PISM.append_time(out, self.ctx.config, t)
771
772 PISM.write_time_bounds(out, bounds, j, [t, t + dt])
773
774 aSMB.set(t * SMB_anomaly)
775
776 dSMBdz.set(t * SMB_gradient)
777
778 aT.set(t * T_anomaly)
779
780 dTdz.set(t * T_gradient)
781
782 for v in [aSMB, dSMBdz, aT, dTdz]:
783 v.write(out)
784
785 out.redef()
786 out.write_attribute("time", "bounds", "time_bounds")
787 out.write_attribute("time", "units", "seconds since 2000-1-1")
788
789 out.close()
790
791 def setUp(self):
792
793 self.ctx = PISM.Context()
794
795 self.grid = shallow_grid()
796 self.geometry = PISM.Geometry(self.grid)
797
798 self.geometry.ice_surface_elevation.set(100.0)
799
800 self.forcing_file = "surface_ismip6_forcing.nc"
801 self.reference_file = "surface_ismip6_reference.nc"
802
803 self.prepare_reference_data(self.grid, self.reference_file)
804 self.prepare_climate_forcing(self.grid, self.forcing_file)
805
806 self.ctx.config.set_string("surface.ismip6.file", self.forcing_file)
807 self.ctx.config.set_string("surface.ismip6.reference_file", self.reference_file)
808
809 def ismip6_test(self):
810 "Surface model ISMIP6"
811
812 atmosphere = PISM.AtmosphereUniform(self.grid)
813
814 model = PISM.SurfaceISMIP6(self.grid, atmosphere)
815
816 model.init(self.geometry)
817
818 t = self.ctx.time.current()
819 dt = model.max_timestep(t).value()
820
821 model.update(self.geometry, t, dt)
822
823 def tearDown(self):
824 os.remove(self.reference_file)
825 os.remove(self.forcing_file)
826
827 if __name__ == "__main__":
828
829 t = ISMIP6()
830
831 t.setUp()
832 t.ismip6_test()
833 t.tearDown()