subroutine splineslab (z,y,nh,km,dydz1,dydzkm,wk,d2ydz2)
c
c=======================================================================
c                                                                    ===
c  This routine computes the second derivative used by the slab      ===
c  data, vertical cubic spline.                                      ===
c                                                                    ===
c  On Input:                                                         ===
c                                                                    ===
c     DYDZ1    first derivative at the surface (real array)          ===
c     DYDZKM   first derivative at the bottom (real array)           ===
c     KM       number of vertical levels (integer)                   ===
c     NH       number of horizontal grid points (integer)            ===
c     WK       Work space (real array)                               ===
c     Y        volume field to be interpolated (real array)          ===
c     Z        depths of the volume grid points (real array)         ===
c                                                                    ===
c  On Output:                                                        ===
c                                                                    ===
c     D2YDZ2   second derivative in the volume (real array)          ===
c                                                                    ===
c  Calls:  none                                                      ===
c                                                                    ===
c  Notes on the boundary conditions (special values):                ===
c                                                                    ===
c     if DYDZ1(i) or DYDZKM(i) is greater than 9.9e29 (spval), then  ===
c        the routine is signalled to set the second derivative to    ===
c        zero at that point.  ("natural" spline)                     ===
c                                                                    ===
c     if DYDZ1(i) or DYDZKM(i) is less than -9.9e29 (spval2), then   ===
c        the routine is signalled to compute the first derivative at ===
c        that point from an interpolating cubic polynomial.          ===
c        ("fitted" spline)  This requires at least 4 levels in the   ===
c        vertical.                                                   ===
c                                                                    ===
c=======================================================================
c
#include <cdefs.h>
#include <pconst.h>
c
c-----------------------------------------------------------------------
c  Define local data.
c-----------------------------------------------------------------------
c
      integer i,k,km,nh
      FLOAT
     *      d21,d31,d41,d32,d42,d43,p,qn,sig,spval,spval2,un,yp
      FLOAT
     *      dydz1(nh),dydzkm(nh),d2ydz2(nh,km),wk(nh,km),y(nh,km),
     *      z(nh,km)
      parameter (spval=c9p9e29,spval2=-c9p9e29)
c
c=======================================================================
c  Begin executable code.
c=======================================================================
c
c-----------------------------------------------------------------------
c  Set surface boundary condition.
c-----------------------------------------------------------------------
c
      do 10 i=1,nh
        if(dydz1(i).gt.spval) then
          d2ydz2(i,1)=c0
          wk(i,1)=c0
        else if (dydz1(i).ge.spval2) then
          d2ydz2(i,1)=mp5
          wk(i,1)=(c3/(z(i,2)-z(i,1)))*((y(i,2)-y(i,1))
     *            /(z(i,2)-z(i,1))-dydz1(i))
        else
          d21=z(i,2)-z(i,1)
          d31=z(i,3)-z(i,1)
          d41=z(i,4)-z(i,1)
          d32=z(i,3)-z(i,2)
          d42=z(i,4)-z(i,2)
          d43=z(i,4)-z(i,3)
          yp=(d31*d41)/(d21*d32*d42)*y(i,2)
     *       +(d21*d31)/(d41*d42*d43)*y(i,4)
     *       -(d21*d41)/(d31*d32*d43)*y(i,3)
     *                 -(c1/d21+c1/d31+c1/d41)*y(i,1)
          d2ydz2(i,1)=mp5
          wk(i,1)=(c3/d21)*((y(i,2)-y(i,1))/d21-yp)
        endif
  10  continue
c
c-----------------------------------------------------------------------
c  Decompose for tridiagonal solver for second derivatives.
c-----------------------------------------------------------------------
c
      do 20 k=2,km-1
      do 20 i=1,nh
        sig=(z(i,k)-z(i,k-1))/(z(i,k+1)-z(i,k-1))
        p=sig*d2ydz2(i,k-1)+c2
        d2ydz2(i,k)=(sig-c1)/p
        wk(i,k)=(c6*((y(i,k+1)-y(i,k))/(z(i,k+1)-z(i,k))-
     *               (y(i,k)-y(i,k-1))/(z(i,k)-z(i,k-1)))
     *              /(z(i,k+1)-z(i,k-1))-sig*wk(i,k-1))/p
  20  continue
c
c-----------------------------------------------------------------------
c  Set bottom boundary condition.  Find second derivative at the
c  bottom.
c-----------------------------------------------------------------------
c
      do 30 i=1,nh
        if(dydzkm(i).gt.spval) then
          qn=c0
          un=c0
        else if (dydzkm(i).ge.spval2) then
          qn=p5
          un=(c3/(z(i,km)-z(i,km-1)))*(dydzkm(i)-
     *       (y(i,km)-y(i,km-1))/(z(i,km)-z(i,km-1)))
        else
          d21=z(i,km)-z(i,km-1)
          d31=z(i,km)-z(i,km-2)
          d41=z(i,km)-z(i,km-3)
          d32=z(i,km-1)-z(i,km-2)
          d42=z(i,km-1)-z(i,km-3)
          d43=z(i,km-2)-z(i,km-3)
          yp=(c1/d21+c1/d31+c1/d41)*y(i,km)
     *       +(d21*d41)/(d31*d32*d43)*y(i,km-2)
     *       -(d31*d41)/(d21*d32*d42)*y(i,km-1)
     *       -(d21*d31)/(d41*d42*d43)*y(i,km-3)
          qn=p5
          un=(c3/d21)*(yp-(y(i,km)-y(i,km-1))/d21)
        endif
        d2ydz2(i,km)=(un-qn*wk(i,km-1))/(qn*d2ydz2(i,km-1)+c1)
  30  continue
c
c-----------------------------------------------------------------------
c  Back substitute to find the remaining second derivatives.
c-----------------------------------------------------------------------
c
      do 40 k=km-1,1,-1
      do 40 i=1,nh
        d2ydz2(i,k)=d2ydz2(i,k)*d2ydz2(i,k+1)+wk(i,k)
  40  continue
      return
      end