[502] | 1 | from common import * |
---|
| 2 | |
---|
| 3 | def plot_T850(lon,lat,T850): # T850 at final time |
---|
| 4 | print 'Reading data ...' |
---|
| 5 | lon, lat, T850 = lon[:], lat[:], T850[-1, :, :] |
---|
| 6 | print '... done.' |
---|
| 7 | plt.figure(figsize=(12,6)) |
---|
| 8 | plt.contourf(lon,lat,T850) |
---|
| 9 | plt.colorbar() |
---|
| 10 | plt.title('T850') |
---|
| 11 | axis_longitude() |
---|
| 12 | axis_latitude() |
---|
| 13 | plt.savefig('T850.png') |
---|
| 14 | |
---|
[512] | 15 | def plot_dT(nlon,nlat,llm, lon,T,p,Phi): # perturbation temp, final time |
---|
[502] | 16 | # vertical slice at final time |
---|
| 17 | print 'Reading data ...' |
---|
[512] | 18 | T,p, Phi = T[-1,:,nlat/2,:], p[-1,:,nlat/2,:], Phi[-1,:,nlat/2,:] |
---|
[502] | 19 | print '... done.' |
---|
| 20 | Cpd, kappa, g = 1004.5, 0.2857143, 9.80616 |
---|
| 21 | N, Teq, peq = 0.01, 300., 1e5 |
---|
| 22 | N2, g2 = N*N, g*g |
---|
| 23 | G = g2/(N2*Cpd) |
---|
| 24 | |
---|
| 25 | lon2, z = np.zeros((llm,nlon)), np.zeros((llm,nlon)) |
---|
| 26 | for lev in range(llm): |
---|
| 27 | z[lev,:] = (.5/g)*(Phi[lev,:]+Phi[lev+1,:]) # average from interfaces to full levels |
---|
| 28 | lon2[lev,:] = lon[:] |
---|
| 29 | |
---|
[512] | 30 | theta = T*((peq/p)**kappa) |
---|
| 31 | Thetab = Teq*np.exp(N2*z/g) |
---|
[502] | 32 | Tb = G + (Teq-G)*np.exp(N2*z/g) # background temperature |
---|
| 33 | plt.figure(figsize=(12,6)) |
---|
[512] | 34 | plt.contourf(lon2,z,theta-Thetab, levels=np.arange(-0.12,0.12,0.02) ) |
---|
[502] | 35 | plt.colorbar() |
---|
[512] | 36 | plt.title('$\\Theta\'$') |
---|
[502] | 37 | axis_longitude() |
---|
| 38 | plt.ylabel('z (m)') |
---|
| 39 | plt.yticks(np.arange(0, 10001, 1000)) |
---|
| 40 | plt.savefig('dT.png') |
---|
| 41 | |
---|
| 42 | gridfile = 'netcdf/output_dcmip2016_regular.nc' |
---|
| 43 | nc = cdf.Dataset(gridfile, "r") |
---|
| 44 | llm, nlon, nlat, ntime = getdims(nc, 'lev','lon','lat','time_counter') |
---|
[512] | 45 | lon, lat, T850, T, Phi, p = getvars(nc, 'lon','lat','T850', 'T', 'PHI','P') |
---|
| 46 | plot_dT(nlon,nlat,llm, lon,T,p,Phi) |
---|
[502] | 47 | plot_T850(lon,lat,T850) |
---|