Hello Dennis
i used these functions before but i got wrong result and when i asked the wrf
help about the problem, they told me this functions was developed for wrf model
and there is a large change that this will not work for other model output. i
work in MM5 model.
Thanks
Ali
Hi group
i'm trying to calculate the CAPE from skew T digram , all my variables 4D
(time, level,lat,lon) i attachment the figure, it work just in (eps) and i
can't open it as (x11) or (pdf),but i got the figure without the red line of the
parcel lapse rate like in example no. 2 in skew examples because i want to
calculate the CAPE, if anyone can help me please.
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/skewt_func.ncl"
begin
colmin = 40
colmax = 47
rowmin = 33
rowmax = 40
class = 6
;1999
diri = "/srv/science/maths1/z3236814/data02/swap/co2/fc/"
f = addfile(diri+"evb_minclass_200noev_2095.nc","r") ; open data netcdf file
topo = f->terrain
w = f->w(class-1,:,:,:,:)
u = f->u(class-1,:,:,:,:)
v = f->v(class-1,:,:,:,:)
t = f->t(class-1,:,:,:,:)
q = f->q(class-1,:,:,:,:)
pp = f->pp(class-1,:,:,:,:)
sigma = f->sigma_level
ps = f->pstarcrs
time = f->time
lat = f->latitcrs
lon = f->longicrs
ptop = f->ptop
dsize = dimsizes(q)
pres = new(dsize,float)
;to calculate the pressure for each grid point
pres = pp + ptop +conform(pp,ps,(/2,3/))*conform(pp,sigma,1)
rel = relhum(t,q,pres)
;to calculates the dewpoint temperature
dewpoint= dewtemp_trh(t,rel)
dewpoint = dewpoint -273 ;returned depoint temperature in (C)
; to define the dewpoint
dewpoint!0 = "time"
dewpoint!1 = "sigma_level"
dewpoint!2 = "i_cross"
dewpoint!3 = "j_cross"
pres = pres*0.01 ; to converte it to (hPa)
;calculate virtual temperature
vir_temp = new(dsize,float)
vir_temp = t*(1.+q*0.61)
; to define the pres
pres!0 = "time"
pres!1 = "sigma_level"
pres!2 = "i_cross"
pres!3 = "j_cross"
; to define the virtual temperature
vir_temp!0 = "time"
vir_temp!1 = "sigma_level"
vir_temp!2 = "i_cross"
vir_temp!3 = "j_cross"
t = t-273 ; to converte it to (C)
gh = hydro(pres(time|:,i_cross|:,j_cross|:,sigma_level|::-1), \
vir_temp(time|:,i_cross|:,j_cross|:,sigma_level|::-1), \
conform(vir_temp(:,0,:,:),topo,(/1,2/)))
; to define the geopotential height
gh!0 = "time"
gh!1 = "i_cross"
gh!2 = "j_cross"
gh!3 = "sigma_level"
;to calculates the meteorological wind direction
r2d = 45.0/atan(1.0) ;to conversion from radians to degrees because the return
value in radians
wdir = atan2(u,v)* r2d + 180
wdir!0 = "time"
wdir!1 = "sigma_level"
wdir!2 = "i_cross"
wdir!3 = "j_cross"
wks = gsn_open_wks ("eps", "skewclass62095")
pres1=
dim_avg(dim_avg(pres(time|:,sigma_level|:,i_cross|rowmin:rowmax,j_cross|colmin:colmax)))
z =
dim_avg(dim_avg(gh(time|:,sigma_level|:,i_cross|rowmin:rowmax,j_cross|colmin:colmax)))
wind =
dim_avg(dim_avg(v(time|:,sigma_level|:,i_cross|rowmin:rowmax,j_cross|colmin:colmax)))
windd=
dim_avg(dim_avg(wdir(time|:,sigma_level|:,i_cross|rowmin:rowmax,j_cross|colmin:colmax)))
temp =
dim_avg(dim_avg(t(time|:,sigma_level|:,i_cross|rowmin:rowmax,j_cross|colmin:colmax)))
dew1 =
dim_avg(dim_avg(dewpoint(time|:,sigma_level|:,i_cross|rowmin:rowmax,j_cross|colmin:colmax)))
skewtOpts = True
skewtOpts@DrawFahrenheit = False
skewtOpts@DrawWind = True
skewtOpts@DrawHeightScale = True
dataOpts = True ; options describing data and ploting
dataOpts@linePatternDewPt = 1
dataOpts@colDewPt = "red"
;add legend
legend = create "Legend" legendClass wks
"vpXF" : 0.15
"vpYF" : 0.88
"vpWidthF" : 0.24 ; width
"vpHeightF" : 0.09 ; height
"lgPerimOn" : True ; perimeter
"lgItemCount" : 2 ; how many
"lgLabelStrings" : (/"Temperature","Dew point Temperature"/)
"lgLineThicknessF" : 3.0
"lgBoxMinorExtentF" : 0.2
end create
skewt_bkgd = skewT_BackGround (wks, skewtOpts)
skewt_data=skewT_PlotData(wks,skewt_bkgd,pres1(15,:),temp(15,:),dew1(15,:),z(15,:),wind(15,:),windd(15,:),dataOpts)
draw (skewt_data)
draw (skewt_bkgd)
draw(legend)
frame(wks)
end
----- End forwarded message -----
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
This archive was generated by hypermail 2.1.8 : Thu Feb 25 2010 - 09:38:32 MST