Problems defining an array to be plotted in gsn_csm_xy

Hi,
I have what is probably an incredibly basic question, but I'm stuck.
I am trying to: read in a grid, interpolate it to height coordinates,
sum values over all gridpoints per each individual pressure level, and
then write those sum values vs. height on an xy plot.
The problem is that I can't seem to correctly define the
sum-values-at-each-height-level array: attempts have resulted in index
errors, only one value getting written (at the top height level,
presumably because I'm overwriting it every time through the loop), and
too many other problems to list.
Again, I think this is basic array construction, but I can't figure it
out for this purpose.
Any help would be incredibly appreciated!
Thanks in advance,
Kelly Mahoney

************************************************************************************************************************************************

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/csm/shea_util.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

begin
setvalues NhlGetWorkspaceObjectId()
         "wsMaximumSize" : 100000000000
end setvalues
setfileoption("nc", "SuppressClose", False)

imgdir="./"
diri="./"

fils = systemfunc("ls ./wrfout_d02_2040*")+".nc"
nfil = dimsizes(fils)

fils_past = systemfunc("ls ./wrfout_d02_1970*")+".nc"

   do nf=0,nfil-1
   f=addfile(fils(nf),"r")
   times = wrf_user_list_times(f) ; get time(s) in the file
   ntimes = dimsizes(times)

   f_past=addfile(fils_past(nf),"r")
   times_past = wrf_user_list_times(f_past) ; get time(s) in the file

lat=f->XLAT(0,:,:)
lon=f->XLONG(0,:,:)
hgt=f->HGT(0,:,:)

;KM read in 3D graupel:
;****** PAST ********
qgraup_past = wrf_user_getvar(f_past,"QGRAUP",-1)
p_past = wrf_user_getvar(f_past, "pressure",-1) ; pressure is our
vertical coordinate
z_past = wrf_user_getvar(f_past, "z",-1) ; grid point height

; The specific height levels that we want the data interpolated to.
; And interpolate to these levels
   height_levels = (/0.,200.,400.,600.,800.,1000., \ ; height
levels - in meters
                     1200.,1400.,1600.,1800.,2000.,\
                     2200.,2400.,2600.,2800.,3000.,\
             3200.,3400.,3600.,3800.,4000.,\
             4200.,4400.,4600.,4800.,5000.,\
             5200.,5400.,5600.,5800.,6000.,\
             6200.,6400.,6600.,6800.,7000.,\
             7200.,7400.,7600.,7800.,8000.,\
             9200.,9400.,9600.,9800.,10000.,\
             10200.,10400.,10600.,10800.,11000.,\
             11200.,11400.,11600.,11800.,12000.,\
             12200.,12400.,12600.,12800.,13000.,\
             13200.,13400.,13600.,13800.,14000.,\
             14200.,14400.,14600.,14800.,15000./)

   nlevels = dimsizes(height_levels) ; number of height levels

;****** INTERPOLATE TO HEIGHT: PAST ********
qgraup_plane_past = int2p_n_Wrap(z_past,qgraup_past,height_levels,1,1)
qgraup_plane_past@_FillValue = -999.
qgraup_plane_past@lat = lat
qgraup_plane_past@lon = lon

;****** SUM/AVG/MEDIAN ETC AMT OF GRAUPEL PER HEIGHT LEVEL ********

do level = 0,nlevels-1 ; LOOP OVER LEVELS

    height = height_levels(level)

sum_graup_past = sum(qgraup_plane_past(0,level,:,:))

printVarSummary(sum_graup_past)

;km sum_graup_past_ht_array = new((/sum_graup_past,height(level)/,float))
;km sum_graup_past_ht_array = (/(/sum_graup_past/),(/height/)/)
;km this is just a sample of many things I've tried...
sum_graup_past_ht_array = sum_graup_past(level)
;KM ????? How to define an array of sum_graup_past *and* height so that
it can be plotted on an xy plot ?????

print("At "+height+"m sum_graup_past from all gridpoints = "+sum_graup_past)

end do ; file loop

end do ; height loop

print(sum_graup_past_ht_array)
print(height_levels)

;***** TRY TO PLOT HEIGHT(YAXIS) VS. SUM_GRAUP_PAST (XAXIS): ************

x = sum_graup_past_ht_array
y = height_levels

wks = gsn_open_wks("ps","xy_test")

res = True
res@trYMinF = 0.
res@trYMaxF = 15000.

res@trXMinF = 0.
res@trXMaxF = 30.

plot = gsn_csm_xy(wks,x,y,res)

end

exit

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Fri Feb 18 16:42:04 2011