Re: Vertical 1D 90m Windspeed interpolation to ASCII file

Hi Justin,
First thing you need to do is create an array that will house the 90m
wind speed for each ensemble member. You can do this by putting an if
statement in, say right before your wind speed calculation at line 50.
Then fill in the array:

v3 = (v1 + v2)/2
if (i.eq.1) then
    finarr = new((/10/),typeof(v3))
    finarr!0 = "ensemble_member"
    finarr&ensemble_member = ispan(1,10,1)
    finarr@long_name = "90m wind speed"
    finarr@units = v3@units
end if
spd = sqrt(u3*u3+v3*v3)
finarr(i) = (/ spd(5) /) ; I'm not sure which index refers to 90m, I
                            ; just put (5)

In the above I assumed that the wind speed you want resides in spd(5),
so you'll want to change that to whatever index refers to the 90m wind
speed.

After your do while loop is complete, write out the .nc file:

    i = i+1
end do
b = addfile("windspeed_90m.nc","c")
b@source = systemfunc("pwd")+get_script_name
b@title = "Title"
b->wspd = finarr

Hope that helps...
Adam

On 10/04/2010 01:07 PM, Justin Traiteur wrote:
> 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

-- 
__________________________________________________
Adam Phillips 
asphilli@ucar.edu
National Center for Atmospheric Research   tel: (303) 497-1726
Climate and Global Dynamics Division         fax: (303) 497-1333
P.O. Box 3000				
Boulder, CO 80307-3000    http://www.cgd.ucar.edu/cas/asphilli
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk