Idealized WRF simulation - Integrate with NCL

I am trying to use example 3 in the Vibeta function (
http://www.ncl.ucar.edu/Document/Functions/Built-in/vibeta.shtml) to
integrate temperature but am confused about what some of the necessary
inputs are for this command even after visiting (
http://www.ncl.ucar.edu/Document/Functions/Built-in/pres_hybrid_ccm.shtml).

I'd like to do this for a column and then for a plane and I'm interested in
temperature.

I have a 525 x 250 grid with 120 levels irregularly spaced.

I'd like to integrate the temperature for a column/point (100,100, nlevles)
(XYZ)
nlevels should extend from the 101500 Pa (surface p) and max out at 5000 Pa
(arbitrary upper boundary)

Then I'd like to do the same for (100, 100:180, nlevels) (XYZ)

Below is my attempt at the column integration, I am getting an error at
psfc saying that it at least needs to be two-dimensional but when I do that
I still get the same error. Surface pressure in WRF is defined as an input
and I don't think that it's explicitly stated as a variable (that I can
find). In the example they grab the variable from the input file.

Any help would be appreciated, thanks!

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
;Column integration of T
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

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

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

  times = wrf_user_list_times(a) ; get times in the file
  ntimes = dimsizes(times) ; number of times in the file

; time loop
; do it = 2,ntimes-1,2 ;every other time step
; do it = 1,ntimes-1 ;every time step

; time = it
     time = 6 ; t=6 hours, steady state reached for the simulation

  p_in = wrf_user_getvar(a,"p",time)
  T_in = wrf_user_getvar(a,"tk",time)

  hyam = p_in(:,100,100) ; read from a file the mid-layer coef
  hybm = T_in(:,100,100) ; read from a file
  psfc = (101500., 100, 100) ; surface pressure [Pa]
  p0 = 101500. ; since ps is in Pa or [ f->P0]

  x = T_in(:,100,100) ; (time,lev,lat,lon)
  pm = x ; simple way to create a variable with
                 ; appropriate metadata
  pm = pres_hybrid_ccm (psfc, p0, hyam, hybm)
; pm(ntim,klvl,nlat,mlon) top-to-bottom
  pm@long_name = "pressure"
  pm@units = "Pa"

  pbot = 101500. ; arbitrary lower level
  ptop = 5000. ; upper level
                  ; create array with appropriate metadata
  vint = x(time|:,lat|:,lon|:,lev|0) ; create a variable with metadata
                                      ; reverse pressure order via ::-1
syntax
                                      ; pm and x will span bottom-to-top
  vint = vibeta (pm(time|:,lat|:,lon|:,lev|::-1), \
                  x(time|:,lat|:,lon|:,lev|::-1),linlog,psfc,pbot,ptop)

  vint = vint/(pbot-ptop) ; normalize [ original units ]
; returns vint(time,lat,lon)

end

Received on Thu Jun 05 07:05:15 2014