Hello everyone,
Below I have created a NCL code to plot several ensemble members of
the WRF Single Column Model on the same plot. Now I need to just get
the 90m windspeed from each ensemble member and write it to a file.
Any suggestions would be appreciated.
Thanks,
Justin
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/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl" ;
WRF_Times2Udunits_c
begin
i = 1
do while (i .le. 10)
date = "_08021300"
diri = "./"
ens = i
fili = "wrfout_ens"
nc = ".nc"
f = addfile (diri+fili+ens+date+nc,
"r") ; Load file
pltType = "png" ;
Use x11 for plotting
wks = gsn_open_wks(pltType,"Stable_1300") ;
Open workstation
res = True
res@MainTitle = "SCM WRF"
pltres = True
mpres = True
times = wrf_user_list_times(f)
it = 1
res@TimeLabel = times(it)
eta_i = ispan(1,14,1)
plot = new (1,"graphic") ;
Plot 2 maps
theta = wrf_user_getvar(f,"theta",it)
thetan1 = theta(eta_i,0,0)
thetan2 = theta(eta_i,0,1)
thetan3 = theta(eta_i,1,0)
thetan4 = theta(eta_i,1,1)
thetan5 = (thetan1 + thetan2)/2 ;
Calculate average Theta between U points
thetan6 = (thetan3 + thetan4)/2 ;
Calculate average Theta between V poitns
thetan7 = (thetan5 + thetan6)/2
z = wrf_user_getvar(f,"z",it) ;
Model height
zn = z(eta_i,0,0) ;
Only first 40 levs
u1 = f->U(it,eta_i,0,1)
u2 = f->U(it,eta_i,1,1)
v1 = f->V(it,eta_i,1,0)
v2 = f->V(it,eta_i,1,1)
u3 = (u1 + u2)/2 ;
Calculate average between U grid points
v3 = (v1 + v2)/2 ;
Calculate average between V grid points
spd = sqrt(u3*u3+v3*v3) ;
Calculate wind speed
res = True ;
plot mods desired
; res@tiMainString = "WRF SCM WIND SPEED ~C~ ~Z75~ 1-Hour
Forecast Valid 14Z ~C~ ~Z75~
Unstable" ; title
res@tiYAxisString = "Height (m)" ; y
axis title
res@tiXAxisString = "Wind Magnitude (m/s)" ; x axis
title
res@trXMaxF = 6
res@trXMinF = 0
res@gsMarkerSizeF = 20
res@gsMarkerThicknessF= 2
res@gsnDraw = False
res@gsnFrame = False
colors =
(/"black
","red
","blue
","green
","brown
","orange","blueviolet","yellow","tomato","aquamarine","slategray2"/)
res@xyLineThicknesses = 3.0 ;
line thicknesses
res@TimeLabel = times(it)
opts = res
opts@cnLineColor = "Red"
FramePlot = False
gsnDraw = False
gsnFrame = False
if(i .eq. 1) then
spd1 = spd
z1 = zn
end if
if(i .eq. 2) then
spd2 = spd
z2 = zn
end if
if(i .eq. 3) then
spd3 = spd
z3 = zn
end if
if(i .eq. 4) then
spd4 = spd
z4 = zn
end if
if(i .eq. 5) then
spd5 = spd
z5 = zn
end if
if(i .eq. 6) then
spd6 = spd
z6 = zn
end if
if(i .eq. 7) then
spd7 = spd
z7 = zn
end if
if(i .eq. 8) then
spd8 = spd
z8 = zn
end if
if(i .eq. 9) then
spd9 = spd
z9 = zn
end if
if(i .eq. 10) then
spd10 = spd
z10 = zn
end if
if(i .eq. 11) then
spd11 = spd
z11 = zn
end if
i = i+1
end do
plot1 = gsn_csm_xy(wks,spd1,z1,res)
res@xyLineColor = colors(1)
plot2 = gsn_csm_xy(wks,spd2,z2,res)
res@xyLineColor = colors(2)
plot3 = gsn_csm_xy(wks,spd3,z3,res)
res@xyLineColor = colors(3)
plot4 = gsn_csm_xy(wks,spd4,z4,res)
res@xyLineColor = colors(4)
plot5 = gsn_csm_xy(wks,spd5,z5,res)
res@xyLineColor = colors(5)
plot6 = gsn_csm_xy(wks,spd6,z6,res)
res@xyLineColor = colors(6)
res@gsMarkerIndex = 4
plot7 = gsn_add_polymarker(wks,plot1,3.0,90,res)
plot12 = gsn_add_polymarker(wks,plot1,2.9,90,res)
plot8 = gsn_csm_xy(wks,spd7,z7,res)
res@xyLineColor = colors(7)
plot9 = gsn_csm_xy(wks,spd8,z8,res)
res@xyLineColor = colors(8)
plot10 = gsn_csm_xy(wks,spd9,z9,res)
res@xyLineColor = colors(9)
plot11 = gsn_csm_xy(wks,spd10,z10,res)
res@xyLineColor = colors(10)
overlay(plot1,plot2)
overlay(plot1,plot3)
overlay(plot1,plot4)
overlay(plot1,plot5)
overlay(plot1,plot6)
overlay(plot1,plot8)
overlay(plot1,plot9)
overlay(plot1,plot10)
overlay(plot1,plot11)
overlay(plot1,plot12)
number = (/1., 2., 3., 4., 5., 6., 7., 8., 9., 10./)
lgres = True
lgres@lgLineColors = colors
lgres@lgItemType = "Lines"
lgres@lgLabelFontHeightF = .08
lgres@vpWidthF = 0.18
lgres@vpHeightF = 0.16
lgres@perimThicknessF = 2.0
lgres@lgMonoDashIndex= True
lgres@lgDashIndex = 0
labels = " Ensemble "+number
legend = gsn_create_legend (wks,10,labels,lgres)
amres = True
amres@amJust = "TopRight"
amres@amParallelPosF = -0.1
amres@amOrthogonalPosF = -0.4
annoid = gsn_add_annotation(plot1,legend,amres)
draw(plot1)
frame(wks)
end
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Oct 4 13:06:03 2010
This archive was generated by hypermail 2.1.8 : Wed Oct 06 2010 - 09:53:35 MDT