Question involving vertical interpolation

I put together a ncl code displaying WRF SCM output. This code displays height (y-axis) vs potential temperature (x-axis) and height (y-axis) vs wind speed (x-axis) for the center of the grid. I know that WRF uses eta coordinates for the vertical. My question is whether the way I derived height from the model and plotted it against the two variables aformentioned is correct or do I need to perform a wrf_interp function instead. Below is a copy of the code.

; SCM: Single Column Model
; WRF SCM: time-eta cross section.
;********************************************************
; %> ncl_filedump wrfout_SCM.nc | less
; netcdf wrfout_CSM {
; dimensions:
; Time = UNLIMITED ; // (60 currently)
; DateStrLen = 19 ;
; west_east = 2 ;
; south_north = 2 ;
; bottom_top = 59 ;
; bottom_top_stag = 60 ;
; soil_layers_stag = 4 ;
; west_east_stag = 3 ;
; force_layers = 8 ;
; south_north_stag = 3 ;
; [SNIP]
; 2x2 west-east/south_north values
; 3x3 west_east_stag/south_north_stag values
; 60 time steps
; 59 or 60 vertical levels
;********************************************************
; Variables are on different spatial grids.
;
; User *must* look!
;
; float U(Time, bottom_top, south_north, west_east_stag)
; float V(Time, bottom_top, south_north_stag, west_east)
; float W(Time, bottom_top_stag, south_north, west_east)
; float PH(Time, bottom_top_stag, south_north, west_east)
;********************************************************
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

  diri = "../"
  fili = "wrfout.nc"
  f = addfile (diri+fili, "r") ;Load file
  
  pltType = "x11" ;Use x11 for plotting
  wks = gsn_open_wks(pltType,"plt_Surface1") ;Open workstation

  res = True
  res@MainTitle = "SCM WRF"
  res@Footer = False
  pltres = True
  mpres = True
  
;********************************************************
  ;res@TimeLabel = times(it)

  times = wrf_user_list_times(f)
  it = 20
  res@TimeLabel = times(it)

;********************************************************
  plot = new (2,"graphic")
  theta = wrf_user_getvar(f,"theta",it)
  thetan1 = theta((/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40/),0,0)
  thetan2 = theta((/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40/),0,1)
  thetan3 = theta((/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40/),1,0)
  thetan4 = theta((/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40/),1,1)
  
  thetan5 = (thetan1 + thetan2)/2
  thetan6 = (thetan3 + thetan4)/2
  thetan7 = (thetan5 + thetan6)/2

  z = wrf_user_getvar(f,"z",it)
  printVarSummary(z)
  zn = z((/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40/),0,0)
  printVarSummary(zn)
  printVarSummary(theta)
  u1 = f->U(it,(/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40/),0,1)
  u2 = f->U(it,(/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40/),1,1)

  v1 = f->V(it,(/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40/),1,0)
  v2 = f->V(it,(/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40/),1,1)
  u3 = (u1 + u2)/2
  v3 = (v1 + v2)/2
  spd = sqrt(u3*u3+v3*v3)
  printVarSummary(spd)

  res = True ; plot mods desired
  res@tiMainString = "WRF SCM THETA" ; title
  res@tiYAxisString = "Height (m)" ; y axis title
  res@tiXAxisString = "Potential Temperature (K)" ; x axis title

  res@xyLineColors = (/"red"/) ; line colors
  res@xyLineThicknesses = 3.0 ; line thicknesses
  res@TimeLabel = times(it)
  

  opts = res
  opts@cnLineColor = "Red"
  plot(0) = gsn_xy(wks,thetan7,zn,res)
  delete(opts)
  delete(res)
  

  res = True ; plot mods desired
  res@tiMainString = "WRF SCM Wind Speed" ; title
  res@tiYAxisString = "Height (m)" ; y axis title
  res@tiXAxisString = "Speed (m/s)" ; x axis title

  res@xyLineColors = (/"blue"/) ; line colors
  res@xyLineThicknesses = 3.0 ; line thicknesses
  res@TimeLabel = times(it)

  plot(1) = gsn_xy(wks,spd,zn,res)
  gsn_panel(wks,plot,(/1,2/),res)

end
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed May 19 15:04:42 2010