module hunting

contains

  ! for spees reasond we decided nt to use optional i parameter but
  ! rather have two suboutines as they are rather simple

  function huntpoli(t, tt, nt, data, nv, i) result(val)

    use typedef, only: &
         int32, real64

    implicit none

    real(kind=real64), intent(in) :: &
         t
    real(kind=real64), dimension(nt), intent(in) :: &
         tt
    real(kind=real64), dimension(nv, nt), intent(in) :: &
         data

    integer(kind=int32), intent(in) :: &
         nt, nv
    integer(kind=int32), intent(inout) :: &
         i
    real(kind=real64), dimension(nv) :: &
         val

    real(kind=real64) :: &
         f

    call hunt(t, tt, nt, i, f)
    val(:) = data(:,i+1) * f + data(:,i) * (1.d0 - f)

  end function huntpoli


  function huntpol(t, tt, nt, data, nv) result(val)

    use typedef, only: &
         int32, real64

    implicit none

    real(kind=real64), intent(in) :: &
         t
    real(kind=real64), dimension(nt), intent(in) :: &
         tt
    real(kind=real64), dimension(nv, nt), intent(in) :: &
         data

    integer(kind=int32), intent(in) :: &
         nt, nv
    real(kind=real64), dimension(nv) :: &
         val

    integer(kind=int32), save :: &
         i = -1
    real(kind=real64) :: &
         f

    call hunt(t, tt, nt, i, f)
    val(:) = data(:,i+1) * f + data(:,i) * (1.d0 - f)

  end function huntpol


  subroutine hunt(t, tt, nt, i, f)

    use typedef, only: &
         int32, real64

    implicit none

    real(kind=real64), intent(in) :: &
         t
    real(kind=real64), dimension(nt), intent(in) :: &
         tt
    integer(kind=int32), intent(in) :: &
         nt

    integer(kind=int32), intent(inout) :: &
         i
    real(kind=real64), intent(out) :: &
         f

    integer(kind=int32) :: &
         i0, i1

    if ((i >= 1).and.(i < nt)) then
       if ((tt(i) < t).and.(tt(i+1) > t)) then
          goto 1000
       endif
    endif

    if (t <= tt(1)) then
       i = 1
       f = 0.d0
       return
    endif
    if (t >= tt(nt)) then
       i = nt - 1
       f = 1.d0
       return
    endif
    if (nt == 2) then
       i = 1
       goto 1000
    endif
    i0 = 1
    i1 = nt
    do while (i0 < i1 - 1)
       i = (i0 + i1) / 2
       if (tt(i) < t) then
          i0 = i
       else
          i1 = i
       endif
    enddo
    if (tt(i) > t) then
       i = i - 1
    endif

1000 continue
    f = (t - tt(i)) / (tt(i+1) - tt(i))

  end subroutine hunt

end module hunting