subroutine spline( x, y, n, i_bc1, i_bc2, param1, param2, y2 )

c  given arrays x(1:n) and y(1:n) containing a tabulated function, i.e., 
c  yi = f(xi), with x1 < x2 < . . . < xn, this routine returns an 
c  array y2(1:n) of length n which contains the second derivatives of 
c  the interpolating function at the tabulated points xi. i_bc1 and i_bc2
c  are boundary condition flags, used as follows:
c    i_bc1 = 0  natural spline with y''(1) = 0
c    i_bc1 = 1  second derivative specified via y''(1) = param1*y''(2)
c    i_bc1 = 2  first derivative specified via y'(1) = param1
c    i_bc1 = 3  periodic
c
c  parameter nmax is the largest anticipated value of n.

      integer n, nmax
      real param1, param2, p1, p2, x(n), y(n), y2(n)
      parameter( nmax=1000 )
      integer i, k
      real a(nmax), b(nmax), c(nmax)

c   *** check inputs

      if( n .gt. nmax ) then
         print *, 'ERROR: n is greater than nmax'
         stop
      end if

      if( (i_bc1 .eq. 3 .and. i_bc2 .ne. 3) .or.
     &    (i_bc2 .eq. 3 .and. i_bc2 .ne. 3)      ) then
         print *, 'ERROR: both i_bc1 and i_bc2 must be set to 3 ',
     &            'for periodic bcs'
         stop
      end if

      p1 = param1
      p2 = param2

      if( i_bc1 .eq. 0 ) p1 = 0.0
      if( i_bc2 .eq. 0 ) p2 = 0.0

c   *** compute and store the tridiagonal matrix elements and rhs
c   *** do this at the endpoints using periodicity even for non-periodic
c   *** cases.  The endpoint values will be rewritten below for non-periodic
c   *** cases.

      c3 = 1.0/3.0
      c6 = 1.0/6.0

      do i=1, n
         ip1 = i + 1
         im1 = i - 1
         if( i .eq. 1 ) then
            im1 = n-1
            delta_m = x(n) - x(n-1)
            delta_0 = x(ip1) - x(i)
         else if( i .eq. n ) then
            ip1 = 2
            delta_m = x(i) - x(im1)
            delta_0 = x(2) - x(1)
         else
            delta_m = x(i) - x(im1)
            delta_0 = x(ip1) - x(i)
         end if
         a(i) = c6*delta_m
         c(i) = c6*delta_0
         b(i) = 2.0*( a(i) + c(i) )
         y2(i) = ( y(ip1) - y(i) )/delta_0 - ( y(i) - y(im1) )/delta_m
      end do

c   *** set boundary conditions

      select case( i_bc1 )
      case( 0, 1 )
         a( 1) = 0.0
         b( 1) = 1.0
         c( 1) = -p1
         y2(1) = 0.0
      case( 2 )
         delta_0 = x(2) - x(1)
         a( 1) = 0.0
         b( 1) = c3*delta_0
         c( 1) = c6*delta_0
         y2(1) = ( y(2) - y(1) )/delta_0 - p1
      case( 3 )
         beta = a(1)
      end select

      select case( i_bc2 )
      case( 0, 1 )
         a( n) = -p2
         b( n) = 1.0
         c( n) = 0.0
         y2(n) = 0.0
      case( 2 )
         delta_0 = x(n) - x(n-1)
         a( n) = c6*delta_0
         b( n) = c3*delta_0
         c( n) = 0.0
         y2(n) = p2 - ( y(n) - y(n-1) )/delta_0
      case( 3 )
         alpha = c(n-1)
      end select
         
c      do i=1, n
c         write(4,40) i, a(i), b(i), c(i), y2(i)
c40       format(i5,1p4e12.4)
c      end do

      if( i_bc1 .eq. 3 ) then
        call tridag_p( a, b, c, alpha, beta, y2, n-1 ) 
        y2(n) = y2(1)
      else
         call tridag( a, b, c, y2, n )
      end if

      return
      end


      subroutine splint(xa,ya,y2a,n,x,y,dy,d2y)

c  given the arrays xa(1:n) and ya(1:n) of length n, which tabulate a 
c  function (with the xai ’s in order), and given the array y2a(1:n), 
c  which is the output from spline above, and given a value of x, this 
c  routine returns a cubic-spline interpolated value y.

      integer n
      real x,xa(n),y2a(n),ya(n)
      integer k,khi,klo
      real a,b,h

      if( x .lt. xa(1) ) then
c         print *, 'Warning x < xa(1) in splint x, xa(1) = ', x, xa(1)
         klo = 1
         khi = 2
         goto 2
      end if

      klo=1
      khi=n
1     if (khi-klo.gt.1) then
         k=(khi+klo)/2
         if(xa(k).gt.x)then
            khi=k
         else
            klo=k
         end if
         goto 1
      endif

2     h=xa(khi)-xa(klo)
      if (h.eq.0.) then
         print *, 'bad x input in splint'
         stop
      end if
      fac1 = h/6.0
      fac2 = h*fac1
      fac3 = 1.0/h

      a =  (xa(khi)-x)*fac3
      b =  (x-xa(klo))*fac3
      c =  (a**3-a)*fac2
      d =  (b**3-b)*fac2
      e = -fac3
      f =  fac3
      g = -(3.0*a**2-1.0)*fac1
      h =  (3.0*b**2-1.0)*fac1

      y   = a*ya(klo) + b*ya(khi) + c*y2a(klo) + d*y2a(khi)
      dy  = e*ya(klo) + f*ya(khi) + g*y2a(klo) + h*y2a(khi)
      d2y =                         a*y2a(klo) + b*y2a(khi)

c      write(88,88) klo, a, b, g, h, y2a(klo), y2a(khi)
c88    format(i5,1p,6e12.4)

      return
      end


      subroutine tridag( a, b, c, r, n )

c     Tri-diagonal matrix solver from Numerical Recipies.

      real a(n), b(n), c(n), r(n)
      real w(n), u(n)

      b_inv = 1.0/b(1)

      u(1) = r(1)*b_inv
      do i=2, n
         w(i) = c(i-1)*b_inv
         b_inv = 1.0/( b(i) - a(i)*w(i) )
         u(i) = ( r(i) - a(i)*u(i-1) )*b_inv
      end do

      r(n) = u(n)
      do i=n-1, 1, -1
         r(i) = u(i) - w(i+1)*r(i+1)
      end do
      
      return
      end


      subroutine tridag_p( a, b, c, alpha, beta, r, n )

c     Periodic tri-diagonal matrix solver from Numerical Recipies.

      real a(n), b(n), c(n), r(n)
      real bb(n), u(n)

      gamma = -b(1)
      gamma_inv = 1.0/gamma

      bb = b
      bb(1) = b(1) - gamma
      bb(n) = b(n) - alpha*beta*gamma_inv

      call tridag( a, bb, c, r, n )

      u = 0.0
      u(1) = gamma
      u(n) = alpha

      call tridag( a, bb, c, u, n )

      fact = (       r(1) + beta*r(n)*gamma_inv ) /
     &       ( 1.0 + u(1) + beta*u(n)*gamma_inv )

      do i=1, n
         r(i) = r(i) - fact*u(i)
      end do

      return
      end

      subroutine spline_coef( x_dat, y_dat, dydx, N )

      real x_dat(N), y_dat(N), dydx(N)

      do i=1, N-1
         dydx(i) = ( y_dat(i+1) - y_dat(i) )/( x_dat(i+1) - x_dat(i) )
      end do

      return
      end

      subroutine lin_spline( x_dat, y_dat, dydx, N, x, y )

      real x_dat(N), y_dat(N), dydx(N) 

      do i=1, N-1
         if( x .ge. x_dat(i) .and. x .lt. x_dat(i+1) ) goto 5
      end do

      if( abs(x-x_dat(1)) .lt. 1.0e-12 ) then
         i=1
         goto 5
      end if
      if( abs(x-x_dat(N)) .lt. 1.0e-12 ) then
         i=N-1
         goto 5
      end if

      print *, 'ERROR: bad value x to lin_spline'
      stop

5     m = i

      y = y_dat(i) + dydx(i)*( x - x_dat(i) )

      return
      end

      function smooth ( x, y, dy, npoint, s, v, a )
c  from  * a practical guide to splines *  by c. de boor    
c  calls  setupq, chol1d
c
c  constructs the cubic smoothing spline  f  to given data  (x(i),y(i)),
c  i=1,...,npoint, which has as small a second derivative as possible
c  while
c  s(f) = sum( ((y(i)-f(x(i)))/dy(i))**2 , i=1,...,npoint ) .le. s .
c
c******  i n p u t  ******
c  x(1),...,x(npoint)   data abscissae,  a s s u m e d  to be strictly
c        increasing .
c  y(1),...,y(npoint)     corresponding data ordinates .
c  dy(1),...,dy(npoint)     estimate of uncertainty in data,  a s s u m-
c        e d  to be positive .
c  npoint.....number of data points,  a s s u m e d  .gt. 1
c  s.....upper bound on the discrete weighted mean square distance of
c        the approximation  f  from the data .
c
c******  w o r k  a r r a y s  *****
c  v.....of size (npoint,7)
c  a.....of size (npoint,4)
c
c*****  o u t p u t  *****
c  a(.,1).....contains the sequence of smoothed ordinates .
c  a(i,j) = f^(j-1)(x(i)), j=2,3,4, i=1,...,npoint-1 ,  i.e., the
c        first three derivatives of the smoothing spline  f  at the
c        left end of each of the data intervals .
c     w a r n i n g . . .   a  would have to be transposed before it
c        could be used in  ppvalu .
c
c******  m e t h o d  ******
c     The matrices  Q-transp*d  and  Q-transp*D**2*Q  are constructed in
c   s e t u p q  from  x  and  dy , as is the vector  qty = Q-transp*y .
c  Then, for given  p , the vector  u  is determined in  c h o l 1 d  as
c  the solution of the linear system
c               (6(1-p)Q-transp*D**2*Q + p*Q)u  = qty  .
c  From  u , the smoothing spline  f  (for this choice of smoothing par-
c  ameter  p ) is obtained in the sense that
c                        f(x(.))  =  y - 6(1-p)D**2*Q*u        and
c                      f''(x(.))  =  6*p*u                      .
c     The smoothing parameter  p  is found (if possible) so that
c                sf(p)  =  s ,
c  with  sf(p) = s(f) , where  f  is the smoothing spline as it depends
c  on  p .  if  s = 0, then p = 1 . if  sf(0) .le. s , then p = 0 .
c  Otherwise, the secant method is used to locate an appropriate  p  in
c  the open interval  (0,1) . However, straightforward application of
c  the secant method, as done in the original version of this program,
c  can be very slow and is influenced by the units in which  x  and  y 
c  are measured, as C. Reinsch has pointed out. Instead, on recommend-
c  ation from C. Reinsch, the secant method is applied to the function
c           g:q |--> 1/sqrt{sfq(q)} - 1/sqrt{s} ,
c  with  sfq(q) := sf(q/(1+q)), since  1/sqrt{sfq}  is monotone increasing
c  and close to linear for larger  q . One starts at  q = 0  with a
c  Newton step, i.e., 
c                q_0 = 0,  q_1 = -g(0)/g'(0)
c  with  g'(0) = -(1/2) sfq(0)^{-3/2} dsfq, where dsfq = -12*u-transp*r*u ,
c  and  u  as obtained for  p = 0 . Iteration terminates as soon as 
c   abs(sf - s) .le. .01*s .
c
c     logical test 
c     parameter (test = .true.)
c     integer itercnt

      integer npoint,   i,npm1
      real a(npoint,4),dy(npoint),s,v(npoint,7),x(npoint),y(npoint)
     *     ,change,ooss,oosf,p,prevsf,prevq,q,sfq,sixp,six1mp,utru
      call setupq(x,dy,y,npoint,v,a(1,4))
      if (s .gt. 0.)                    go to 20
   10 p = 1.
      six1mp = 0.
      call chol1d(p,v,a(1,4),npoint,1,a(1,3),a(1,1))
      sfq = 0.
                                        go to 60
   20 p = 0.
      six1mp = 6.
      call chol1d(p,v,a(1,4),npoint,1,a(1,3),a(1,1))
      sfq = 0.
      do 21 i=1,npoint
   21    sfq = sfq + (a(i,1)*dy(i))**2
      sfq = sfq*36.
      if (sfq .le. s)                   go to 60
      utru = 0.
      do 25 i=2,npoint
   25    utru = utru + v(i-1,4)*(a(i-1,3)*(a(i-1,3)+a(i,3))+a(i,3)**2)
      ooss = 1./sqrt(s)
      oosf = 1./sqrt(sfq)
      q = -(oosf-ooss)*sfq/(6.*utru*oosf)
c  secant iteration for the determination of p starts here.
c     itercnt = 0
      prevq = 0.
      prevsf = oosf
   30 call chol1d(q/(1.+q),v,a(1,4),npoint,1,a(1,3),a(1,1))
      sfq = 0.
      do 35 i=1,npoint
   35    sfq = sfq + (a(i,1)*dy(i))**2
      sfq = sfq*36./(1.+q)**2
      if (abs(sfq-s) .le. .01*s)        go to 59
      oosf = 1./sqrt(sfq)
      change = (q-prevq)/(oosf-prevsf)*(oosf-ooss)
      prevq = q
      q = q - change
      prevsf = oosf
c     itercnt = itercnt + 1
                                        go to 30
   59 p = q/(1.+q)
      six1mp = 6./(1.+q)
correct value of p has been found.
compute pol.coefficients from  Q*u (in a(.,1)).
   60 smooth = sfq
c     if (test) then
c        print *, 'number of iterations = ', itercnt
c     end if
      do 61 i=1,npoint
   61    a(i,1) = y(i) - six1mp*dy(i)**2*a(i,1)
      sixp = 6.*p
      do 62 i=1,npoint
   62    a(i,3) = a(i,3)*sixp
      npm1 = npoint - 1
      do 63 i=1,npm1
         a(i,4) = (a(i+1,3)-a(i,3))/v(i,4)
   63    a(i,2) = (a(i+1,1)-a(i,1))/v(i,4)
     *                      - (a(i,3)+a(i,4)/3.*v(i,4))/2.*v(i,4)
                                        return
      end

      subroutine setupq ( x, dx, y, npoint, v, qty )
c  from  * a practical guide to splines *  by c. de boor    
c  to be called in  s m o o t h
c  put  delx = x(.+1) - x(.)  into  v(.,4),
c  put  the three bands of  q-transp*d  into  v(.,1-3), and
c  put the three bands of  (d*q)-transp*(d*q)  at and above the diagonal
c     into  v(.,5-7) .
c     here,  q is  the tridiagonal matrix of order (npoint-2,npoint)
c  with general row  1/delx(i) , -1/delx(i) - 1/delx(i+1) , 1/delx(i+1)
c  and   d  is the diagonal matrix  with general row  dx(i) .

      integer npoint,   i,npm1
      real dx(npoint),qty(npoint),v(npoint,7),x(npoint),y(npoint),
     *                                                       diff,prev
      npm1 = npoint - 1
      v(1,4) = x(2) - x(1)
      do 11 i=2,npm1
         v(i,4) = x(i+1) - x(i)
         v(i,1) = dx(i-1)/v(i-1,4)
         v(i,2) = - dx(i)/v(i,4) - dx(i)/v(i-1,4)
   11    v(i,3) = dx(i+1)/v(i,4)
       v(npoint,1) = 0.
      do 12 i=2,npm1
   12    v(i,5) = v(i,1)**2 + v(i,2)**2 + v(i,3)**2
      if (npm1 .lt. 3)                  go to 14
      do 13 i=3,npm1
   13    v(i-1,6) = v(i-1,2)*v(i,1) + v(i-1,3)*v(i,2)
   14 v(npm1,6) = 0.
      if (npm1 .lt. 4)                  go to 16
      do 15 i=4,npm1
   15    v(i-2,7) = v(i-2,3)*v(i,1)
   16 v(npm1-1,7) = 0.
      v(npm1,7) = 0.
construct  q-transp. * y  in  qty.
      prev = (y(2) - y(1))/v(1,4)
      do 21 i=2,npm1
         diff = (y(i+1)-y(i))/v(i,4)
         qty(i) = diff - prev
   21    prev = diff
                                        return
      end

      subroutine chol1d ( p, v, qty, npoint, ncol, u, qu )
c  from  * a practical guide to splines *  by c. de boor    
c  from  * a practical guide to splines *  by c. de boor    
c to be called in  s m o o t h
constructs the upper three diags. in v(i,j), i=2,,npoint-1, j=1,3, of
c  the matrix  6*(1-p)*q-transp.*(d**2)*q + p*r, then computes its
c  l*l-transp. decomposition and stores it also in v, then applies
c  forward and backsubstitution to the right side q-transp.*y in  qty
c  to obtain the solution in  u .

      integer ncol,npoint,   i,npm1,npm2
      real p,qty(npoint),qu(npoint),u(npoint),v(npoint,7), prev,ratio
     *    ,six1mp,twop
      npm1 = npoint - 1
c     construct 6*(1-p)*q-transp.*(d**2)*q  +  p*r
      six1mp = 6.*(1.-p)
      twop = 2.*p
      do 2 i=2,npm1
         v(i,1) = six1mp*v(i,5) + twop*(v(i-1,4)+v(i,4))
         v(i,2) = six1mp*v(i,6) + p*v(i,4)
    2    v(i,3) = six1mp*v(i,7)
      npm2 = npoint - 2
      if (npm2 .ge. 2)                  go to 10
      u(1) = 0.
      u(2) = qty(2)/v(2,1)
      u(3) = 0.
                                        go to 41
c  factorization
   10 do 20 i=2,npm2
         ratio = v(i,2)/v(i,1)
         v(i+1,1) = v(i+1,1) - ratio*v(i,2)
         v(i+1,2) = v(i+1,2) - ratio*v(i,3)
         v(i,2) = ratio
         ratio = v(i,3)/v(i,1)
         v(i+2,1) = v(i+2,1) - ratio*v(i,3)
   20    v(i,3) = ratio
c
c  forward substitution
      u(1) = 0.
      v(1,3) = 0.
      u(2) = qty(2)
      do 30 i=2,npm2
   30    u(i+1) = qty(i+1) - v(i,2)*u(i) - v(i-1,3)*u(i-1)
c  back substitution
      u(npoint) = 0.
      u(npm1) = u(npm1)/v(npm1,1)
      i = npm2    
   40    u(i) = u(i)/v(i,1)-u(i+1)*v(i,2)-u(i+2)*v(i,3)
         i = i - 1
         if (i .gt. 1)                  go to 40
c  construct q*u
   41 prev = 0.
      do 50 i=2,npoint
         qu(i) = (u(i) - u(i-1))/v(i-1,4)
         qu(i-1) = qu(i) - prev
   50    prev = qu(i)
      qu(npoint) = -qu(npoint)
                                        return
      end

      subroutine my_sort( a, n, i_order )

      real a(n)
      integer i_order(n)

      do i=1, n
         i_order(i) = i
c         print *, i, a(i), i_order(i)
      end do

      do j=2, n
         a1 = a(j)
         j1 = i_order(j)
         do i=j-1, 1, -1
            if( a(i) .le. a1 ) goto 10
            a(i+1) = a(i)
            i_order(i+1) = i_order(i)
         end do
         i = 0
10       a(i+1) = a1
         i_order(i+1) = j1
c         print *, ' ' 
c         do i=1, n
c            print *, i, a(i), i_order(i)
c         end do
      end do

      return
      end

      subroutine auto_smooth(z, u, work, nz, i_method, dev_est, target)

c   *** if i_method=1:
c   *** Smooth the data vector u so that abs(dev(d2udz2)/avg(d2udz2))<=target.
c   *** if i_method=0:
c   *** Smooth the data vector u using a simple call to the smooth function
c   *** with the measurement deviation set to dev_est.
c   *** dev_est is an estimate of the standard deviation of the input data.
c   *** The smoothed data are returned in work(1:nz).  The first and second
c   *** derivatives are in work(nz+1:2*nz) and work(2*nz+1:3*nz) respecively.
c   *** Evidently the smooth function assumes d2udz2=0 at the endpoints.  This
c   *** can cause problems for data that actually have significant second
c   *** derivative values at the boundaries.  

      real z(nz), u(nz), us(nz), work(nz*12)
      real sum1(3)

      pi = acos(-1.0)
      wave_z = 2*pi/float(nz-1)

      ic =  2*nz + 1
      id =  4*nz + 1
      iw =  5*nz + 1

      s_fac = nz + sqrt(2.0*nz)

c   *** Make sure we have at least 3 points

      if( nz .lt. 3 ) then
         print *, 'Error in auto_smooth nz must be at least 3'
         stop
      end if

c   *** Use simple smoothing if called for, or if nz <= nz_min

      nz_min = 6
      if( i_method .eq. 0 .or. nz .le. nz_min ) goto 90

c   *** Perform the adptive smoothing

      dev_meas = 1.0e-10

      do L=1, 20

         do k=id, id+nz
            work(k) = dev_meas
         end do

         sfp = smooth( z, u, work(id), nz, s_fac, work(iw), work(1) )

c      *** Compute the mean and variance of the second derivatives

         sum1 = 0.0
         do k=1, nz-1
            d2udz2 = work(2*nz+k)
            sum1(1) = sum1(1) + d2udz2
            sum1(2) = sum1(2) + d2udz2**2
         end do

         scale = 1.0/float(nz-1)
         avg     = sum1(1)*scale
c         var     = sum1(2)*scale - avg**2
         var     = sum1(2)*scale
         dev = sqrt(var)

c      *** Compute the contribution to the standard deviation in
c      *** the upper 1/4 of the spectrum

c      print *, 'nz, k_start, k_stop = ', 
c     & nz, nint(float(3*(nz-1))/8.0), (nz-1)/2

         do m=nint(float(3*(nz-1))/8.0), (nz-1)/2
            sum2 = 0.0
            sum3 = 0.0
            do k=1, nz-1
               arg = float(k-1)*float(m)*wave_z
               c1 = scale*cos(arg)
               s1 = scale*sin(arg)
               d2udz2 = work(2*nz+k)
               sum2 = sum2 + d2udz2*c1
               sum3 = sum3 + d2udz2*s1
            end do
            fac = 2.0
            if( m .eq. (nz-1)/2 .and. mod(m,2) .eq. 0.0 ) fac = 1.0
            sum1(3) = sum1(3) + fac*( sum2**2 + sum3**2 )
         end do

         dev_hw = sqrt( sum1(3) )

c      *** The quality is the fraction of high wavenumber deviation to
c      *** the total deviation

         quality = abs(dev_hw/dev)
         f = quality - target

c         print *, 'L, dev, dev_hw = ', L, dev, dev_hw
c         print *, 'dev_hw/dev, dev_meas = ', dev_hw/dev, dev_meas

         if( f .lt. 0.1*target ) return

         if( quality .lt. 5*target .and. L .gt. 5 ) then
            dfdx = ( f - f_old )/( dev_meas - dev_old )
            dev_old = dev_meas
            dev_meas = dev_meas - f/dfdx
         else
            dev_old = dev_meas
            dev_meas = dev_meas + 0.2*dev_est
         end if

         f_old = f

      end do

c   *** If flow gets here the adaptive smoothing iteration failed
c   *** In this case we issue a warning and then use the simple smoothing
c   *** approach

      print *, 'WARNING: Iteration did not converge in auto_smooth.'
      print *, 'Using simple method instead.'

90    do k=id, id+nz
         work(k) = dev_est
      end do

      sfp = smooth( z, u, work(id), nz, s_fac, work(iw), work(1) )

      if( nz .le. nz_min ) then
         print *, 'Warning: not enough points for auto smoothing ',
     &            'using simple method'
      end if

      return
      end