subroutine splintslab (zra,yra,d2ydz2,nh,km,indx,klow,khi,z,y,xtp)
c
c=======================================================================
c                                                                    ===
c  This routines performs the vertical cubic spline interpolation    ===
c  with  slab data.  This  routines  assumes  that the  calls are    ===
c  increasing in depth.  It also performs constant extrapolations.   ===
c                                                                    ===
c  On Input:                                                         ===
c                                                                    ===
c     D2YDZ2   second derivatives for the cubic interpolation        ===
c              (real array)                                          ===
c     INDX     horizontal location for interpolation (integer)       ===
c     KHI      data level deeper than the desired depth (integer)    ===
c     KLOW     data level shallower than the desired depth (integer) ===
c     KM       number of vertical levels (integer)                   ===
c     NH       number of horizontal grid points (integer)            ===
c     YRA      field to be interpolated (real array)                 ===
c     Z        depth at which to interpolate (real)                  ===
c     ZRA      depths of the tabulated field (real array)            ===
c                                                                    ===
c  On Output:                                                        ===
c                                                                    ===
c     XTP      flag to indicate if extrapolation was used (logical): ===
c              XTP = .true.   =>  extrapolation was used             ===
c              XTP = .false.  =>  extrapolation was NOT used         ===
c     KHI      data level deeper than the desired depth (integer)    ===
c     KLOW     data level shallower than the desired depth (integer) ===
c     Y        interpolated value (real)                             ===
c                                                                    ===
c  Calls:  none                                                      ===
c                                                                    ===
c=======================================================================
c
c-----------------------------------------------------------------------
c  Define global data.
c-----------------------------------------------------------------------
c
#include <cdefs.h>
#include <pconst.h>
c
c-----------------------------------------------------------------------
c  Define local data.
c-----------------------------------------------------------------------
c
      logical xtp
      integer indx,khi,klow,km,nh
      FLOAT
     *      a,b,d2ydz2(nh,km),h,y,yra(nh,km),z,zra(nh,km)
c
c=======================================================================
c  Begin executable code.
c=======================================================================
c
c-----------------------------------------------------------------------
c  Search for bracketting levels, starting from previous ones.
c-----------------------------------------------------------------------
c
      do 10 while ((z.gt.zra(indx,khi)).and.(khi.lt.km))
        khi=khi+1
        klow=klow+1
  10  continue
c
c-----------------------------------------------------------------------
c  Check to see if interpolation or extrapolation is appropriate.
c-----------------------------------------------------------------------
c
      if((z.ge.zra(indx,1)).and.(z.le.zra(indx,km))) then
c
c-----------------------------------------------------------------------
c  Z is an interior point, interpolate.
c-----------------------------------------------------------------------
c
        h=zra(indx,khi)-zra(indx,klow)
        a=(zra(indx,khi)-z)/h
        b=(z-zra(indx,klow))/h
        if (abs(a).le.abs(b)) then
          y=yra(indx,khi)+a*(yra(indx,klow)-yra(indx,khi))
        else
          y=yra(indx,klow)+b*(yra(indx,khi)-yra(indx,klow))
        endif
        y=y+
     *    ((a**3-a)*d2ydz2(indx,klow)+(b**3-b)*d2ydz2(indx,khi))*h*h*r6
        xtp=.false.
      elseif (z.lt.zra(indx,1)) then
c
c-----------------------------------------------------------------------
c  Z is too shallow, use shallowest tabulated value.
c-----------------------------------------------------------------------
c
        y=yra(indx,1)
        xtp=.true.
      else
c
c-----------------------------------------------------------------------
c  Z is too deep, use deepest tabulated value.
c-----------------------------------------------------------------------
c
        y=yra(indx,km)
        xtp=.true.
      endif
      return
      end