tg02a.f90

Go to the documentation of this file.

subroutine tg02a(ix, n, n_bnd, u, s, d, x, v)
!------------------------------------------------------------------
! HSL subroutine to calculate splined values
!   n   : number of points
!   ix  : negative 0 -> no initial guess for where xi is
!         positive -> guess for index close to value x
!   u   : the coordinates of the spline points
!   s   : the function values of the spline points
!   d   : the derivatives on the spline points
!   x   : the coordinate where the output is wanted
!   v(1-4) : value and derivatives of the spline interpolation
!------------------------------------------------------------------

  use itm_types

  implicit none

  real(r8), intent(in) :: x
  integer(itm_i4), intent(in) :: ix, n, n_bnd
  real(r8), intent(in) :: d(n_bnd + 2), s(n_bnd + 2), u(n_bnd + 2)
  real(r8), intent(out) :: v(4)

  real(r8), parameter :: eps = 1e-33_r8
  real(r8) :: a, b, c, c3, d0, d1
  real(r8) :: gama, h, hr, hrr, phi, s0, s1, t, theta
  integer(itm_i4) :: i, j, k
  integer(itm_i4) :: iflg

  k = 0
  iflg = 0
  if (x < u(1)) go to 990
  if (x > u(n)) go to 991
  if (ix < 0 .or. iflg == 0) go to 12
  if (x > u(j + 1)) go to 1
  if (x >= u(j)) go to 18
  go to 2

    1 j = j + 1
   11 if (x > u(j + 1)) go to 1
      go to 7
   12 j = abs(x - u(1)) / (u(n) - u(1)) * (n - 1) + 1
      j = min(j, n - 1)
      iflg = 1
      if (x >= u(j)) go to 11
    2 j = j - 1
      if (x < u(j)) go to 2
    7 k = j
      h = u(j + 1) - u(j)
      hr = 1._r8 / h
      hrr = 2._r8 * hr * hr
      s0 = s(j)
      s1 = s(j + 1)
      d0 = d(j)
      d1 = d(j + 1)
      a = s1 - s0
      b = a - h * d1
      a = a - h * d0
      c = a + b
      c3 = c * 3._r8
   18 theta = (x - u(j)) * hr
      phi = 1._r8 - theta
      t = theta * phi
      gama = theta * b - phi * a
      v(1) = theta * s1 + phi * s0 + t * gama
      v(2) = theta * d1 + phi * d0 + t * c3 * hr
      v(3) = (c * (phi - theta) - gama) * hrr
      v(4) = -c3 * hrr * hr
      return
  990 if (x <= u(1) - eps * max(abs(u(1)), abs(u(n)))) go to 99
      j = 1
      go to 7
  991 if (x >= u(n) + eps * max(abs(u(1)), abs(u(n)))) go to 995
      j = n - 1
      go to 7
  995 k = n
   99 iflg = 0
  do i = 1, 4
    v(i) = 0._r8
  end do

  return

end subroutine tg02a