simple cubic spline module

    ! =======================================================================
    ! Module
    ! =======================================================================
    module spline_module
      ! [specification part]
      implicit none
      integer, parameter, private :: dp = kind(1.0d0)

      private

      public :: spline,splint,ratint,polint

      ! [execution part]

    contains
      ! [your functions]
    ! =======================================================================
    ! subroutine
    ! =======================================================================
      subroutine spline(x,y,yp1,ypn,y2)
        implicit none
        real (dp), dimension (:), intent (in) :: x, y
        real (dp), intent (in) :: yp1, ypn
        real (dp), dimension (:), intent (inout) :: y2
        integer :: n
        real (dp), dimension (size(x)) :: a, b, c, r

        n = assert_eq3(size(x),size(y),size(y2),'spline')
        c(1:n-1) = x(2:n) - x(1:n-1)
        r(1:n-1) = 6.0_dp*((y(2:n)-y(1:n-1))/c(1:n-1))
        r(2:n-1) = r(2:n-1) - r(1:n-2)
        a(2:n-1) = c(1:n-2)
        b(2:n-1) = 2.0_dp*(c(2:n-1)+a(2:n-1))
        b(1) = 1.0_dp
        b(n) = 1.0_dp
        if (yp1>0.99d30) then
          r(1) = 0.0_dp
          c(1) = 0.0_dp
        else
          r(1) = (3.0_dp/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
          c(1) = 0.5_dp
        end if
        if (ypn>0.99d30) then
          r(n) = 0.0_dp
          a(n) = 0.0_dp
        else
          r(n) = (-3.0_dp/(x(n)-x(n-1)))*((y(n)-y(n-1))/(x(n)-x(n-1))-ypn)
          a(n) = 0.5_dp
        end if
        call tridag_par(a(2:n),b(1:n),c(1:n-1),r(1:n),y2(1:n))
      end subroutine spline

    ! =======================================================================
    ! subroutine
    ! =======================================================================
      function splint(xa,ya,y2a,x)

        implicit none
        real (dp), dimension (:), intent (in) :: xa, ya, y2a
        real (dp), intent (in) :: x
        real (dp) :: splint
        integer :: khi, klo, n
        real (dp) :: a, b, h

        n = assert_eq3(size(xa),size(ya),size(y2a),'splint')
        klo = max(min(locate(xa,x),n-1),1)
        khi = klo + 1
        h = xa(khi) - xa(klo)
        if (h==0.0_dp) call nrerror('bad xa input in splint') ! MOD
        a = (xa(khi)-x)/h
        b = (x-xa(klo))/h
        splint = a*ya(klo) + b*ya(khi) + ((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi)) *(h**2)/6.0_dp
      end function splint



    ! =======================================================================
    ! subroutine
    ! =======================================================================
      subroutine polint(xa,ya,x,y,dy)
        implicit none
        real (dp), dimension (:), intent (in) :: xa, ya
        real (dp), intent (in) :: x
        real (dp), intent (out) :: y, dy
        integer :: m, n, ns
        real (dp), dimension (size(xa)) :: c, d, den, ho

        n = assert_eq2(size(xa),size(ya),'polint')
        c = ya
        d = ya
        ho = xa - x
        ns = iminloc(abs(xa-x))
        y = ya(ns)
        ns = ns - 1
        do m = 1, n - 1
          den(1:n-m) = ho(1:n-m) - ho(1+m:n)
          if (any(den(1:n-m)==0.0_dp)) call nrerror('polint: calculation failure') ! MOD
          den(1:n-m) = (c(2:n-m+1)-d(1:n-m))/den(1:n-m)
          d(1:n-m) = ho(1+m:n)*den(1:n-m)
          c(1:n-m) = ho(1:n-m)*den(1:n-m)
          if (2*ns<n-m) then
            dy = c(ns+1)
          else
            dy = d(ns)
            ns = ns - 1
          end if
          y = y + dy
        end do
      end subroutine polint

    ! =======================================================================
    ! subroutine
    ! =======================================================================
      subroutine ratint(xa,ya,x,y,dy)
        implicit none
        real (dp), dimension (:), intent (in) :: xa, ya
        real (dp), intent (in) :: x
        real (dp), intent (out) :: y, dy
        integer :: m, n, ns
        real (dp), dimension (size(xa)) :: c, d, dd, h, t
        real (dp), parameter :: MTINY = 1.0d-25

        n = assert_eq2(size(xa),size(ya),'ratint')
        h = xa - x
        ns = iminloc(abs(h))
        y = ya(ns)
        if (x==xa(ns)) then        ! TODO
          dy = 0.0_dp
          return
        end if
        c = ya
        d = ya + MTINY
        ns = ns - 1
        do m = 1, n - 1
          t(1:n-m) = (xa(1:n-m)-x)*d(1:n-m)/h(1+m:n)
          dd(1:n-m) = t(1:n-m) - c(2:n-m+1)
          if (any(dd(1:n-m)==0.0_dp)) call nrerror('failure in ratint') ! MOD
          dd(1:n-m) = (c(2:n-m+1)-d(1:n-m))/dd(1:n-m)
          d(1:n-m) = c(2:n-m+1)*dd(1:n-m)
          c(1:n-m) = t(1:n-m)*dd(1:n-m)
          if (2*ns<n-m) then
            dy = c(ns+1)
          else
            dy = d(ns)
            ns = ns - 1
          end if
          y = y + dy
        end do
      end subroutine ratint


    ! =======================================================================
    ! subroutine
    ! =======================================================================
      subroutine tridag_ser(a,b,c,r,u)
        implicit none
        real (dp), dimension (:), intent (in) :: a, b, c, r
        real (dp), dimension (:), intent (inout) :: u
        real (dp), dimension (size(b)) :: gam
        integer :: n, j
        real (dp) :: bet

        n = assert_eqn((/size(a)+1,size(b),size(c)+1,size(r),size(u)/), 'tridag_ser')
        bet = b(1)
        if (bet==0.0_dp) call nrerror('tridag_ser: error at code stage 1')
        ! MOD
        u(1) = r(1)/bet
        do j = 2, n
          gam(j) = c(j-1)/bet
          bet = b(j) - a(j-1)*gam(j)
          if (bet==0.0_dp) call nrerror('tridag_ser: error at code stage 2')
          ! MOD
          u(j) = (r(j)-a(j-1)*u(j-1))/bet
        end do
        do j = n - 1, 1, -1
          u(j) = u(j) - gam(j+1)*u(j+1)
        end do
      end subroutine tridag_ser

    ! =======================================================================
    ! subroutine
    ! =======================================================================
      recursive subroutine tridag_par(a,b,c,r,u)
        implicit none
        real (dp), dimension (:), intent (in) :: a, b, c, r
        real (dp), dimension (:), intent (inout) :: u
        integer, parameter :: NPAR_TRIDAG = 4
        integer :: n, n2, nm, nx
        real (dp), dimension (size(b)/2) :: y, q, piva
        real (dp), dimension (size(b)/2-1) :: x, z
        real (dp), dimension (size(a)/2) :: pivc

        n = assert_eqn((/size(a)+1,size(b),size(c)+1,size(r),size(u)/), 'tridag_par')
        if (n<NPAR_TRIDAG) then
          call tridag_ser(a,b,c,r,u)
        else
          if (maxval(abs(b(1:n)))==0.0_dp) call nrerror('tridag_par: possible singular matrix') ! MOD
          n2 = size(y)
          nm = size(pivc)
          nx = size(x)
          piva = a(1:n-1:2)/b(1:n-1:2)
          pivc = c(2:n-1:2)/b(3:n:2)
          y(1:nm) = b(2:n-1:2) - piva(1:nm)*c(1:n-2:2) - pivc*a(2:n-1:2)
          q(1:nm) = r(2:n-1:2) - piva(1:nm)*r(1:n-2:2) - pivc*r(3:n:2)
          if (nm<n2) then
            y(n2) = b(n) - piva(n2)*c(n-1)
            q(n2) = r(n) - piva(n2)*r(n-1)
          end if
          x = -piva(2:n2)*a(2:n-2:2)
          z = -pivc(1:nx)*c(3:n-1:2)
          call tridag_par(x,y,z,q,u(2:n:2))
          u(1) = (r(1)-c(1)*u(2))/b(1)
          u(3:n-1:2) = (r(3:n-1:2)-a(2:n-2:2)*u(2:n-2:2)-c(3:n-1:2)*u(4:n:2))/b(3:n-1:2)
          if (nm==n2) u(n) = (r(n)-a(n-1)*u(n-1))/b(n)
        end if
      end subroutine tridag_par

    ! =======================================================================
    ! function
    ! =======================================================================
      function locate(xx,x)
        implicit none
        real (dp), dimension (:), intent (in) :: xx
        real (dp), intent (in) :: x
        integer :: locate
        integer :: n, jl, jm, ju
        logical :: ascnd

        n = size(xx)
        ascnd = (xx(n)>=xx(1))
        jl = 0
        ju = n + 1
        do
          if (ju-jl<=1) exit
          jm = (ju+jl)/2
          if (ascnd .eqv. (x>=xx(jm))) then
            jl = jm
          else
            ju = jm
          end if
        end do
        if (x==xx(1)) then         ! TODO
          locate = 1
        else if (x==xx(n)) then    ! TODO
          locate = n - 1
        else
          locate = jl
        end if
      end function locate

    ! =======================================================================
    ! function
    ! =======================================================================
      function assert_eq2(n1,n2,string)
        implicit none
        character (len=*), intent (in) :: string
        integer, intent (in) :: n1, n2
        integer :: assert_eq2

        if (n1==n2) then
          assert_eq2 = n1
        else
          write (*,*) 'nrerror: an assert_eq failed with this tag:', string
          stop 'program terminated by assert_eq2'
        end if
      end function assert_eq2

    ! =======================================================================
    ! function
    ! =======================================================================
      function assert_eq3(n1,n2,n3,string)
        implicit none
        character (len=*), intent (in) :: string
        integer, intent (in) :: n1, n2, n3
        integer :: assert_eq3

        if (n1==n2 .and. n2==n3) then
          assert_eq3 = n1
        else
          write (*,*) 'nrerror: an assert_eq failed with this tag:', string
          stop 'program terminated by assert_eq3'
        end if
      end function assert_eq3

    ! =======================================================================
    ! function
    ! =======================================================================
      function assert_eqn(nn,string)
        implicit none
        character (len=*), intent (in) :: string
        integer, dimension (:), intent (in) :: nn
        integer :: assert_eqn

        if (all(nn(2:)==nn(1))) then
          assert_eqn = nn(1)
        else
          write (*,*) 'nrerror: an assert_eq failed with this tag:', string
          stop 'program terminated by assert_eqn'
        end if
      end function assert_eqn

    ! =======================================================================
    ! subroutine
    ! =======================================================================
      subroutine nrerror(string)
        character (len=*), intent (in) :: string

        write (*,*) 'nrerror: ', string
        stop 'program terminated by nrerror'
      end subroutine nrerror

    ! =======================================================================
    ! function
    ! =======================================================================
      function iminloc(arr)
        implicit none
        real (dp), dimension (:), intent (in) :: arr
        integer, dimension (1) :: imin
        integer :: iminloc

        imin = minloc(arr(:))
        iminloc = imin(1)
      end function iminloc

    end module spline_module