tb15a.f90

Go to the documentation of this file.

subroutine tb15a(n ,n_bnd, x, f, d, w, lp)
!------------------------------------------------------------------
! HSL routine for cubic spline with periodic boundary conditions
! first point must be the same as last : f(1) = f(n)
!    n : number of points
!    x : coordinate (input)
!    f : the function values to be splined (input)
!    d : the derivatives at the points (output)
!    w : workspace (dimension 3n)
!   lp : unit number for output
!------------------------------------------------------------------

  use itm_types

  implicit none

  integer(itm_i4) :: lp, n, n_bnd
  real(r8) :: d(n_bnd + 2), f(n_bnd + 2), x(n_bnd + 2)
  real(r8) :: w(3 * n_bnd + 6)
  real(r8) :: a3n1, f1, f2, h1, h2, p
  integer(itm_i4) :: i, k

!.  write(lp, *) f(1), f(n)
  if (n < 4) then
    write (lp, '(a39)')  'return from tb15a because n too small'
    w(1) = 1._r8
    return
  end if
  do i = 2, n
    if (x(i) <= x(i - 1)) then
      write (lp, '(a29, i3, a13)') 'return from tb15a because ', i, &
       ' out of order'
      w(1) = 2._r8
      return
    end if
  end do
  if (f(1) /= f(n)) then
    write (lp, '(a40)')  'return from tb15ad because f(1) /= f(n)'
    w(1) = 3._r8
    return
  end if
  do i = 2, n
    h1 = 1._r8 / (x(i) - x(i - 1))
    f1 = f(i - 1)
    if (i == n) then
      h2 = 1._r8 / (x(2) - x(1))
      f2 = f(2)
    else
      h2 = 1._r8 / (x(i + 1) - x(i))
      f2 = f(i + 1)
    end if
    w(3 * i - 2) = h1
    w(3 * i - 1) = 2._r8 * (h1 + h2)
    w(3 * i) = h2
    d(i) = 3._r8 * (f2 * h2 * h2 + f(i) * (h1 * h1 - h2 * h2) - f1 * h1 &
     * h1)
  end do
  k = 5
  a3n1 = w(3 * n - 1)
  do i = 2, n - 2
    p = w(k + 2) / w(k)
    w(k + 3) = w(k + 3) - p * w(k + 1)
    d(i + 1) = d(i + 1) - p * d(i)
    w(k + 2) = -p * w(k - 1)
    p = w(k - 1) / w(k)
    a3n1 = -p * w(k - 1) + a3n1
    d(n) = d(n) - p * d(i)
    k = k + 3
  end do
  p = (w(k + 2) + w(k - 1)) / w(k)
  a3n1 = a3n1 - p * (w(k + 1) + w(k - 1))
  d(n) = (d(n) - p * d(n - 1)) / a3n1
  do i = 3, n
    k = n + 2 - i
    d(k) = (d(k) - w(3 * k) * d(k + 1) - w(3 * k - 2) * d(n)) / w(3 * k - 1)
  end do
  d(1) = d(n)
  w(1) = 0._r8

  return

end subroutine tb15a