vertical profile of WRF output in NCL

From: Joseph Brodie <jbrodie_at_nyahnyahspammersnyahnyah>
Date: Tue Aug 21 2012 - 08:33:31 MDT

Hi all,

 

I'm trying to get a lower atmosphere vertical profile from some WRF output
data I have (I only want the first 200 meters of the atmosphere) and for
some reason, when I use the wrf_user_intrp3d function, and then plot the
data, it removes the lowest 100 m of the atmosphere from the data! I know
there is data there, because I have model layers starting as low as 8m, and
if I plot horizontal data, I have no problem getting the data from the
layers below 100 m. Any guidance would be appreciated! The code is below.

 

Cheers,

Joe Brodie

 

-------- NCL CODE

; WRF windfarm output processesing - VERTICAL CROSS SECTION

; J. F. Brodie

; University of Delaware

; 8 August 2012

 

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"

load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

 

begin

  a = addfile("wrfout_farm1_newlays.nc","r")

; a = addfile("wrfout_d01_0001-01-03_00:00:00.nc","r")

 

  type = "x11" ; X11 visual output

; type = "pdf" ; PDF output

 

  plane = (/ 17.0, 17.0/) ; x,y coordinate for vertical cross

  angle = 90.0 ; 90 degrees = W-E cross section

  FirstTime = True

  

  wks = gsn_open_wks(type,"WRF_AWEA_farm1")

  gsn_define_colormap(wks,"precip3_16lev")

  gsn_reverse_colormap(wks)

 

; plot options

  res = True

  res@MainTitle = "WRF WINDFARM - AWEA"

 

  pltres = True

 

  height = 200. ; maximum height I would like displayed, meters

 

; let's get all the times

  times = wrf_user_getvar(a,"times",-1)

  ntimes = dimsizes(times)

 

  do it = 0,ntimes-1,1

 

    print("Working on time " + times(it) )

    res@TimeLabel = times(it)

 

    u = wrf_user_getvar(a,"ua",it)

    v = wrf_user_getvar(a,"va",it) ; U and V velocities averaged to mass pts

    z = wrf_user_getvar(a, "z",it) ; grid point height

 

    if ( FirstTime ) then

      zz = wrf_user_intrp3d(z,z,"v",plane,angle,False)

      b = ind(zz(:,0) .gt. height )

      zmax_pos = b(0) - 1

      if ( abs(zz(zmax_pos,0)-height) .lt. abs(zz(zmax_pos+1,0)-height) )
then

        zspan = b(0) - 1

      else

        zspan = b(0)

      end if

      delete(zz)

      delete(b)

      FirstTime = False

    end if

 

    u_plane = wrf_user_intrp3d(u,z,"v",plane,angle,False)

    v_plane = wrf_user_intrp3d(v,z,"v",plane,angle,False)

    u_plane@units = "m/s"

    v_plane@units = "m/s"

 

    spd = (u_plane*u_plane + v_plane*v_plane)^(0.5) ; convert to speed

    spd@description = "Wind Speed"

    spd@units = "m/s"

 

; contour of speed

    opts = res

    opts@cnFillOn = True

    opts@trYMinF = 0

    opts@trYMaxF = zspan

; contour = wrf_contour(a,wks,u_plane(0:zmax_pos,:),opts)

    contour = wrf_contour(a,wks,spd(0:zmax_pos,:),opts)

    delete(opts)

 

; overlay plots together

    plot = wrf_overlays(a,wks,(/contour/),pltres)

 

    delete(u)

    delete(u_plane)

    delete(v)

   delete(v_plane)

    delete(z)

    delete(spd)

 

  end do

 

end

 

--
Joseph F. Brodie
Graduate Research Assistant, Center for Carbon-free Power Integration
School of Marine Science and Policy
College of Earth, Ocean, and Environment
University of Delaware
007A Robinson Hall
Newark, DE 19716 USA
Phone: +1 302-831-0694
Fax: +1 302-831-6838
Email: jbrodie@udel.edu
 

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Tue Aug 21 08:33:56 2012

This archive was generated by hypermail 2.1.8 : Thu Aug 23 2012 - 16:16:15 MDT