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