tslr_show.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
---
tslr_show.py (3943B)
---
1 #!/usr/bin/env python3
2
3 # Copyright (C) 2011-2012, 2014, 2016, 2017 The PISM Authors
4
5 # script to generate figure: results from SeaRISE experiments
6 # usage: if UAFX_G_D3_C?_??.nc are result NetCDF files then do
7 # $ slr_show.py -m UAFX
8
9 # try different netCDF modules
10 try:
11 from netCDF4 import Dataset as CDF
12 except:
13 print("netCDF4 is not installed!")
14 sys.exit(1)
15
16 from numpy import zeros
17 import pylab as plt
18 from optparse import OptionParser
19
20 parser = OptionParser()
21 parser.usage = "usage: %prog [options]"
22 parser.description = "A script for PISM output files to show time series plots using pylab."
23 parser.add_option("-a", dest="t_a", type="int",
24 help="start year, in years since 2004, default = 0", default=0)
25 parser.add_option("-e", dest="t_e", type="int",
26 help="end year, in years since 2004, default = 500", default=500)
27 parser.add_option("-m", "--model", dest="model",
28 help="choose experiment, default UAF1", default="UAF1")
29
30
31 (options, args) = parser.parse_args()
32 model = options.model
33 t_a = options.t_a
34 t_e = options.t_e
35
36 # first name in this list is CONTROL
37 NCNAMES = [model + "_G_D3_C1_E0.nc", model + "_G_D3_C2_E0.nc", model + "_G_D3_C3_E0.nc", model + "_G_D3_C4_E0.nc", model + "_G_D3_C1_S1.nc", model +
38 "_G_D3_C1_S2.nc", model + "_G_D3_C1_S3.nc", model + "_G_D3_C1_M1.nc", model + "_G_D3_C1_M2.nc", model + "_G_D3_C1_M3.nc", model + "_G_D3_C1_T1.nc"]
39
40 # labels
41 labels = ["AR4 A1B", "AR4 A1B 1.5x", "AR4 A1B 2x", "2x basal sliding", "2.5x basal sliding",
42 "3x basal sliding", "2 m/a bmr", "20 m/a bmr", "200 m/a bmr", "AR4 A1B + 2x sliding"]
43 # line colors
44 colors = ['#984EA3', # violet
45 '#984EA3', # violet
46 '#984EA3', # violet
47 '#FF7F00', # orange
48 '#FF7F00', # orange
49 '#FF7F00', # orange
50 '#377EB8', # light blue
51 '#377EB8', # light blue
52 '#377EB8', # light blue
53 '#4DAF4A'] # green
54
55 dashes = ['-', '--', '-.', '-', '--', '-.', '-', '--', '-.', '-']
56
57 print("control run name is " + NCNAMES[0])
58
59 n = len(NCNAMES)
60 nc0 = CDF(NCNAMES[0], 'r')
61 try:
62 t_units = nc0.variables['tseries'].units
63 t = nc0.variables['tseries'][t_a:t_e]
64 except:
65 t_units = nc0.variables['time'].units
66 t = nc0.variables['time'][t_a:t_e]
67 nc0.close()
68
69 # convert to years if t is in seconds
70 if (t_units.split()[0] == ('seconds' or 's')):
71 t /= 3.15569259747e7
72
73 ice_volume_glacierized = zeros((len(t), n))
74 ivolshift = zeros((len(t), n - 1))
75
76 for j in range(n):
77 nc = CDF(NCNAMES[j], 'r')
78 ice_volume_glacierized[:, j] = nc.variables['ice_volume_glacierized'][t_a:t_e]
79 nc.close()
80
81 for j in range(n - 1):
82 ivolshift[:, j] = ice_volume_glacierized[:, j + 1] - ice_volume_glacierized[:, 0]
83
84 # "2,850,000 km3 of ice were to melt, global sea levels would rise 7.2 m"
85 scale = 7.2 / 2.850e6
86
87
88 # screen plot with high contrast
89 fig = plt.figure()
90 ax = fig.add_subplot(111, axisbg='0.15')
91 for j in range(n - 1):
92 ax.plot(t, -(ivolshift[:, j] / 1.0e9) * scale, dashes[j], color=colors[j], linewidth=3)
93 ax.set_xlabel('years from 2004')
94 ax.set_ylabel('sea level rise relative to control (m)')
95 ax.legend(labels, loc='upper left')
96 ax.grid(True, color='w')
97
98 plt.show()
99
100
101 # line colors
102 colors = ['#984EA3', # violet
103 '#984EA3', # violet
104 '#984EA3', # violet
105 '#FF7F00', # orange
106 '#FF7F00', # orange
107 '#FF7F00', # orange
108 '#084594', # dark blue
109 '#084594', # dark blue
110 '#084594', # dark blue
111 '#4DAF4A'] # green
112
113 # print plot with white background
114 fig = plt.figure()
115 ax = fig.add_subplot(111)
116 for j in range(n - 1):
117 ax.plot(t, -(ivolshift[:, j] / 1.0e9) * scale, dashes[j], color=colors[j], linewidth=2)
118 ax.set_xlabel('years from 2004')
119 ax.set_ylabel('sea level rise relative to control (m)')
120 ax.legend(labels, loc='upper left')
121 ax.grid(True)
122 plt.savefig(model + '_slr.pdf')