1 | import numpy as np |
---|
2 | import netCDF4 as cdf |
---|
3 | # select non-interactive backend, cf http://stackoverflow.com/questions/4931376/generating-matplotlib-graphs-without-a-running-x-server |
---|
4 | import matplotlib |
---|
5 | matplotlib.use('Agg') |
---|
6 | import matplotlib.pyplot as plt |
---|
7 | |
---|
8 | def getdims(nc, *names): return [len(nc.dimensions[name]) for name in names] |
---|
9 | def getvars(nc, *names): return [nc.variables[name] for name in names] |
---|
10 | |
---|
11 | def axis_longitude(): |
---|
12 | plt.xlim((0,360)) |
---|
13 | plt.xlabel('longitude (degrees)') |
---|
14 | plt.xticks(np.arange(0,361,30)) |
---|
15 | |
---|
16 | def axis_latitude(): |
---|
17 | plt.ylim((-90,90)) |
---|
18 | plt.ylabel('latitude (degrees)') |
---|
19 | plt.yticks(np.arange(-90,91,30)) |
---|
20 | |
---|
21 | def slice_lon(nlon,llm,lon,Phi): |
---|
22 | lon2, z = np.zeros((llm,nlon)), np.zeros((llm,nlon)) |
---|
23 | for lev in range(llm): |
---|
24 | # average from interfaces to full levels |
---|
25 | z[lev,:] = (.5/g)*(Phi[lev,:]+Phi[lev+1,:]) |
---|
26 | lon2[lev,:] = lon[:] |
---|
27 | return lon2, z |
---|
28 | |
---|
29 | #--------------------------- DCMIP21 --------------------------- |
---|
30 | |
---|
31 | def post_DCMIP21(): |
---|
32 | def plot_var(nlon,nlat,llm, lon,Phi,var,varname): |
---|
33 | # vertical slice at final time |
---|
34 | print 'Reading data ...' |
---|
35 | var, Phi = var[-1,:,nlat/2,:], Phi[-1,:,nlat/2,:] |
---|
36 | print '... done.' |
---|
37 | lon2, z = slice_lon(nlon,llm,lon,Phi) |
---|
38 | plt.figure(figsize=(12,6)) |
---|
39 | plt.contourf(lon2,z,var) |
---|
40 | plt.colorbar() |
---|
41 | plt.title('%s at final time'%varname) |
---|
42 | axis_longitude() |
---|
43 | plt.ylabel('z (m)') |
---|
44 | # plt.yticks(np.arange(0, 10001, 1000)) |
---|
45 | plt.savefig('%s.png'%varname) |
---|
46 | lon, lat, Omega,T,u,Phi = getvars(nc, 'lon','lat','OMEGA', 'T', 'U', 'PHI') |
---|
47 | plot_var(nlon,nlat,llm, lon,Phi,Omega,'Omega') |
---|
48 | plot_var(nlon,nlat,llm, lon,Phi,u,'u') |
---|
49 | plot_var(nlon,nlat,llm, lon,Phi,T,'T') |
---|
50 | |
---|
51 | #--------------------------- DCMIP31 --------------------------- |
---|
52 | |
---|
53 | def post_DCMIP31(): |
---|
54 | |
---|
55 | def plot_T850(lon,lat,T850): # T850 at final time |
---|
56 | print 'Reading data ...' |
---|
57 | lon, lat, T850 = lon[:], lat[:], T850[-1, :, :] |
---|
58 | print '... done.' |
---|
59 | plt.figure(figsize=(12,6)) |
---|
60 | plt.contourf(lon,lat,T850) |
---|
61 | plt.colorbar() |
---|
62 | plt.title('T850') |
---|
63 | axis_longitude() |
---|
64 | axis_latitude() |
---|
65 | plt.savefig('T850.png') |
---|
66 | |
---|
67 | def plot_dT(nlon,nlat,llm, lon,T,p,Phi): # perturbation temp, final time |
---|
68 | # vertical slice at final time |
---|
69 | print 'Reading data ...' |
---|
70 | T,p, Phi = T[-1,:,nlat/2,:], p[-1,:,nlat/2,:], Phi[-1,:,nlat/2,:] |
---|
71 | print '... done.' |
---|
72 | N, Teq, peq = 0.01, 300., 1e5 |
---|
73 | N2, g2 = N*N, g*g |
---|
74 | G = g2/(N2*Cpd) |
---|
75 | |
---|
76 | lon2, z = slice_lon(nlon,llm,lon,Phi) |
---|
77 | theta = T*((peq/p)**kappa) |
---|
78 | Thetab = Teq*np.exp(N2*z/g) |
---|
79 | plt.figure(figsize=(12,6)) |
---|
80 | plt.contourf(lon2,z,theta-Thetab, levels=np.arange(-0.12,0.12,0.02) ) |
---|
81 | plt.colorbar() |
---|
82 | plt.title('$\\Theta\'$') |
---|
83 | axis_longitude() |
---|
84 | plt.ylabel('z (m)') |
---|
85 | plt.yticks(np.arange(0, 10001, 1000)) |
---|
86 | plt.savefig('dT.png') |
---|
87 | |
---|
88 | lon, lat, T850, T, Phi, p = getvars(nc, 'lon','lat','T850', 'T', 'PHI','P') |
---|
89 | plot_dT(nlon,nlat,llm, lon,T,p,Phi) |
---|
90 | plot_T850(lon,lat,T850) |
---|
91 | |
---|
92 | #--------------------------- DCMIP41 --------------------------- |
---|
93 | |
---|
94 | def post_DCMIP41(): |
---|
95 | def plot_T850(day, lon,lat,T850): # T850 at final time |
---|
96 | print 'Reading data ...' |
---|
97 | lon, lat, T850 = lon[:], lat[:], T850[day-1, :, :] |
---|
98 | print '... done.' |
---|
99 | plt.figure(figsize=(12,5)) |
---|
100 | plt.contourf(lon,lat,T850) |
---|
101 | plt.colorbar() |
---|
102 | plt.title('T850 at day %d'%day) |
---|
103 | axis_longitude() |
---|
104 | axis_latitude() |
---|
105 | plt.savefig('T850_day%02d.png'%day) |
---|
106 | |
---|
107 | lon, lat, T850, T, Phi = getvars(nc, 'lon','lat','T850', 'T', 'PHI') |
---|
108 | plot_T850(7,lon,lat,T850) |
---|
109 | plot_T850(10,lon,lat,T850) |
---|
110 | plot_T850(30,lon,lat,T850) |
---|
111 | |
---|
112 | #------------------------ Held & Suarez ---------------------- |
---|
113 | |
---|
114 | post_held_suarez = post_DCMIP41 |
---|
115 | |
---|
116 | #--------------------------- MAIN ---------------------------- |
---|
117 | |
---|
118 | gridfile = 'netcdf/output_dcmip2016_regular.nc' |
---|
119 | nc = cdf.Dataset(gridfile, "r") |
---|
120 | llm, nlon, nlat, ntime = getdims(nc, 'lev','lon','lat','time_counter') |
---|
121 | Cpd, kappa, g = 1004.5, 0.2857143, 9.80616 |
---|
122 | |
---|
123 | # Now call a routine post_XXX() |
---|