source: trunk/UTIL/PYTHON/2dplot.py @ 445

Last change on this file since 445 was 310, checked in by aslmd, 13 years ago

MESOSCALE:
-- Code in storm scenario corresponds to 'OMEGA' reference case [Julien Faure]
-- Easier settings for dust lifting without the need to recompile [see below]
-- A few modifications to plot and dust-devil detection PYTHON routines

29/09/11 == AS

--> To easily explore sensitivity to lifting thresholds: in dustlift.F, ustar_seuil=sqrt(stress/rho)

and alpha_lift[dust_mass] can be prescribed through an external stress.def parameter file.
--- alpha_lift[dust_number] is computed from alpha_lift[dust_mass] as in initracer.F
--- ustar_seuil is more user-friendly than stress because direct comparison with ustar from model

--> For the moment this is MESOSCALE only, but potentially useful to everyone

  • Property svn:executable set to *
File size: 1.7 KB
Line 
1#! /usr/bin/env python
2
3from netCDF4 import Dataset
4import matplotlib.pyplot as mpl
5import numpy as np
6from myplot import reducefield,getfield,getcoorddef,calculate_bounds,bounds,fmtvar,ptitle,makeplotres
7from matplotlib.pyplot import contourf,colorbar,show,xlabel,ylabel
8from matplotlib.cm import get_cmap
9
10name = "wrfout_d01_9999-09-09_09:00:00_z"
11itime = 12
12#itime = 1
13ndiv = 10
14zey = 0
15var = "W"
16vmin = -1.
17vmax = 1.
18title = "Vertical velocity"
19var = "Um"
20vmin = -2.
21vmax = 18.
22title = "Horizontal velocity"
23#var = "tk"
24#vmin = 150.
25#vmax = 170.
26#title = "Atmospheric temperature"
27
28name = "wrfout_d01_2024-03-22_01:00:00_z"
29itime = 12
30ndiv = 10
31zey = 120
32var = "VMR_ICE"
33title = "Volume mixing ratio of water ice [ppm]"
34vmin = -0.5
35vmax = 400.
36
37nc = Dataset(name)
38
39what_I_plot, error = reducefield( getfield(nc,var), d4=itime, d2=zey )
40
41
42y = nc.variables["vert"][:]
43
44horinp = len(what_I_plot[0,:])
45x = np.linspace(0.,horinp*500.,horinp) / 1000.
46xlabel("Horizontal coordinate (km)")
47
48horinp = len(what_I_plot[0,:])
49x = np.linspace(0.,horinp,horinp) 
50xlabel("x grid points")
51
52
53zevmin, zevmax = calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax)
54#if colorb in ["def","nobar"]:   palette = get_cmap(name=defcolorb(fvar))
55#else:                           palette = get_cmap(name=colorb)
56palette = get_cmap(name="jet")
57what_I_plot = bounds(what_I_plot,zevmin,zevmax)
58zelevels = np.linspace(zevmin,zevmax)
59contourf( x, y, what_I_plot, zelevels, cmap = palette )
60colorbar(fraction=0.05,pad=0.03,format=fmtvar(var),\
61                       ticks=np.linspace(zevmin,zevmax,ndiv+1),\
62                       extend='neither',spacing='proportional')
63ptitle(title)
64ylabel("Altitude (m)")
65makeplotres(var+str(itime),res=200.,disp=False)
66#show()
Note: See TracBrowser for help on using the repository browser.