Vertical 1D 90m Windspeed interpolation to ASCII file

From: Justin Traiteur <jtraite2_at_nyahnyahspammersnyahnyah>
Date: Mon Oct 04 2010 - 13:07:46 MDT

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