begin

  a = addfile("wrfout_d04_2007-09-24_05:00:00.nc","r")

  
  wks = gsn_open_wks("pdf","vertical_test")
  pltres                             = True
  mpres                              = True

;***********************************************
; plot control
;***********************************************

  res          = True
  res@gsnDraw  = False
  res@gsnFrame = False
  wres         = res

  wres@FieldTitle = "Wind"
  wres@vcRefMagnitudeF         = 3                 ; make vectors larger
  wres@vcRefLengthF            = 0.01               ; ref vec length
  wres@vcGlyphStyle            = "LineArrow"        ; select wind barbs
  wres@vcMinDistanceF          = 0.025              ; thin out windbarbs
  wres@vcLineArrowThicknessF   = 1

  res@cnFillOn             = True               ; turns on the color
  res@cnLinesOn            = False              ; turn off contour lines
  res@cnLevelSelectionMode = "ManualLevels"     ; set manual contour levels
  res@cnMinLevelValF       = 0.                 ; set min contour level
  res@cnMaxLevelValF       = 16.                ; set max contour level
  res@cnLevelSpacingF      = 1.                 ; set contour spacing
  
  xlat  = wrf_user_getvar(a, "XLAT",0) 
  xlon  = wrf_user_getvar(a, "XLONG",0) 
  ter   = wrf_user_getvar(a, "HGT",0) 

;**********************************************
; for a single time
;**********************************************

  time    = 3
  z       = wrf_user_getvar(a,"z",time)      ; height in mass points
  tc      = wrf_user_getvar(a,"tc",time)     ; temperature in degrees c        
  w       = wrf_user_getvar(a,"wa",time)     
  u       = wrf_user_getvar(a,"ua",time)

  printVarSummary(z)
  printVarSummary(tc) 
 
  plane   = (/  500,530 , 600,530  /)

;--- Y-axis labels 
  zmax = 6000. ; We only want to see the first 6 km 
  zz = wrf_user_intrp3d(z,z,"v",plane,90.0,True) 
  dims = dimsizes(zz) 
  do imax = 0,dims(0)-1 
  if ( .not.ismissing(zz(imax,0)) .and. zz(imax,0) .lt. zmax ) then 
    zmax_pos = imax 
  end if 
  end do 
  zspan = zmax_pos 
  zmin = z(0,0,0) 
  zmax = zz(zmax_pos,0)
  print(zmax) 
  zmin=zmin/1000. 
  zmax=zmax/1000. 
  nzlabs = floattoint(zmax + 1) 

;--- X-axis labels
  X_plane = wrf_user_intrp2d(xlon,plane,90.0,True) 
  X_desc = "longitude" 

  dimsX = dimsizes(X_plane) 
  xmin = X_plane(0) 
  xmax = X_plane(dimsX(0)-1) 
  xspan = dimsX(0)-1
  nxlabs = floattoint( (xmax-xmin)/2 + 1) 

  print(dimsX)
  print(xmin)
  print(xmax)
  print(xspan)

  res@tmYLMode = "Explicit" 
  res@tmYLValues = fspan(0,zspan,nzlabs)                     ; Create the correct tick marks 
  res@tmYLLabels = sprintf("%.1f",fspan(zmin,zmax,nzlabs))   ; Create labels 

  res@tmXBMode = "Explicit" 
  res@tmXBValues = fspan(0,xspan,nxlabs)                     ; Create the correct tick marks 
  res@tmXBLabels = sprintf("%.1f",fspan(xmin,xmax,nxlabs))   ; Create labels 

  w_plane  = wrf_user_intrp3d(w,z,"v",plane,90.0,True)
  u_plane  = wrf_user_intrp3d(u,z,"v",plane,90.0,True) 
  tc_plane = wrf_user_intrp3d(tc,z,"v",plane,90.0,True) 

  w_plane2 = w_plane(0:10,:)
  u_plane2 = u_plane(0:10,:)
  tc_plane2 = tc_plane(0:10,:)  

  uw_plane_vector  = wrf_vector(a,wks,u_plane2,w_plane2,wres)
  tc_plane_contour = wrf_contour(a,wks,tc_plane2,res) 

  plot = wrf_overlays(a,wks,(/tc_plane_contour,uw_plane_vector/),pltres)

end