function zinterp(field,zgrid,zlev)
*----------------------------------------------------------------------
* 
* Bob Hart (rhart@fsu.edu) /  FSU Meteorology
* 3/4/1999
*
* GrADS function to interpolate within a 3-D grid to a specified
* height level.  Can also be used on non-pressure level data, such
* as sigma or eta-coordinate output where height is a function
* of time and grid level.
* 
* Advantages:  Easy to use, no UDFs.  Disadvantages:  Can take 3-10 secs.
*
* Arguments:
*    field = name of 3-D grid to interpolate
*
*    zgrid = name of 3-D grid holding height values at each gridpoint
*            Units of measurement are arbitrary.
*
*    zlev  = height level at which to interpolate (having same units as zgrid)
*
* Function returns:  defined grid interp holding interpolated values
*
* NOTE:  Areas having zlev below bottom level or above upper level 
*        in output will be undefined.
*
* NOTE:  No distinction in the function is made between height above
*        sea level, height above ground surface, or geopotential. The
*        function will give output regardless of which is sent.  
*        It is up to the user to be aware which height variable is
*        being passed to the function and treat the output accordingly.
*
* Example function calls:
*
*      "d "zinterp(temp,height,5000)
*
* Would display a temperature field interpolated to 5000.
*      
*      "define t1000="zinterp(temp,height,1000)
*
* Would define a new variable, t1000, as a temp field at 1000.
*
*      "d p1000="zinterp(lev,height,1000)
*
* Would display a field of the pressure at a height of 1000.
*
* PROBLEMS:  Send email to Bob Hart (rhart@fsu.edu)
* 
*-----------------------------------------------------------------------
*-------------------- BEGINNING OF FUNCTION ----------------------------
*-----------------------------------------------------------------------
*
* Get initial dimensions of dataset so that exit dimensions will be
* same

"q dims"
rec=sublin(result,4)
ztype=subwrd(rec,3)
if (ztype = "fixed") 
   zmin=subwrd(rec,9)
   zmax=zmin
else
   zmin=subwrd(rec,11)
   zmax=subwrd(rec,13)
endif

* Get full z-dimensions of dataset.

"q file"
rec=sublin(result,5)
zsize=subwrd(rec,9)

* Determine spatially varying bounding height levels for height surface
* zabove = height-value at level above ; zbelow = height value at level
* below for each gridpoint

"set z 2 "zsize
"define zabove=0.5*maskout("zgrid","zgrid"-"zlev")+0.5*maskout("zgrid","zlev"-"zgrid"(z-1))"
"set z 1 "zsize-1
"define zbelow=0.5*maskout("zgrid","zgrid"(z+1)-"zlev")+0.5*maskout("zgrid","zlev"-"zgrid")"

* Isolate field values at bounding height levels
* fabove = requested field value above height surface
* fbelow = requested field value below height surface

"set z 2 "zsize
"define fabove=zabove*0+"field
"set z 1 "zsize-1
"define fbelow=zbelow*0+"field

* Turn this 3-D grid of values (mostly undefined) into a 2-D height layer
* mean is used here below for simplicity, since mean ignores undefined
* values.

"set z 1"
"define zabove=mean(zabove,z=2,z="zsize")"
"define fabove=mean(fabove,z=2,z="zsize")"
"define zbelow=mean(zbelow,z=1,z="zsize-1")"
"define fbelow=mean(fbelow,z=1,z="zsize-1")"

* Finally, interpolate linearly in height and create surface.

"set z "zmin " " zmax

"define slope=(fabove-fbelow)/(zabove-zbelow)"
"define b=fbelow-slope*zbelow"
"define interp=slope*"zlev"+b"

* variable interp now holds height field and its named it returned
* for use by the user.

say "Done.  Newly defined variable interp has "zlev" "field"-field."

Return(interp)