!> @brief
!> splintt calculates an interpolated value from function
!> using cubic spline interpolation.
!> Use spline first for derivative calculation.

!> @param[in] n | length of xa, ya, y2a
!> @param[in] xa | locations for ya=f(xa)
!> @param[in] ya | values for ya=f(xa)
!> @param[in] y2a | array of derivatives calcualted by spline
!> @param[in] x | location for interpolated value
!> @param[inout] tension | tension between 0-1.0
!> @param[out] y | interpolated value valid at x
SUBROUTINE splintt(xa,ya,y2a,n,x,tension,y)
  implicit none
  ! tension added to the numerical recipies routine splint
  ! tension should be between 0. and 1.0
  !
  INTEGER:: n
  REAL, INTENT(IN), DIMENSION(n) :: xa, ya, y2a
  REAL, INTENT(INOUT) :: tension
  REAL, INTENT(IN) :: x
  REAL, INTENT(OUT):: y
  INTEGER :: k,khi,klo
  REAL:: a,b,h

  if (tension < 0.0) tension=1.0
  if (tension > 1.0) tension=1.0

  klo=1
  khi=n
  do while (khi-klo.gt.1)
     k=(khi+klo)/2
     if(xa(k).gt.x)then
        khi=k
     else
        klo=k
     endif
  enddo
  h=xa(khi)-xa(klo)
  ! if (h.eq.0.) pause 'bad xa input in splint'
  a=(xa(khi)-x)/h
  b=(x-xa(klo))/h
  y=a*ya(klo)+b*ya(khi)+ tension*((a**3-a)*y2a(klo)+  &
       (b**3-b)*y2a(khi))*(h**2)/6.
  return
END SUBROUTINE splintt

!> @brief
!> This calculates the cubic spline second derivative of y for fast
!> repeated use in splintt.

!> @param[in] x | array of locations for function
!> @param[in] y | array of function y = f(x)
!> @param[in] n | length of x, y, y2
!> @param[in] yp1 | first derivative of y at start of array
!> @param[in] ypn | first derivative of y at end of array
!> @param[out] y2 | cubic spline derivative
SUBROUTINE spline(x,y,n,yp1,ypn,y2)
  implicit none
  INTEGER :: i,k
  INTEGER:: n
  REAL, INTENT(IN):: yp1,ypn
  REAL, INTENT(IN), Dimension(n) :: x,y
  REAL, INTENT(INOUT),Dimension(n):: y2
  REAL:: p,qn,sig,un
  REAL,dimension(n)::u

  if (yp1.gt..99e30) then
     y2(1)=0.
     u(1)=0.
  else
     y2(1)=-0.5
     u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
  endif
  do  i=2,n-1
     sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
     p=sig*y2(i-1)+2.
     y2(i)=(sig-1.)/p
     u(i)=(6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))&
          /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
  end do
  if (ypn.gt..99e30) then
     qn=0.
     un=0.
  else
     qn=0.5
     un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
  endif
  y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
  do k=n-1,1,-1
     y2(k)=y2(k)*y2(k+1)+u(k)
  end do
  return
END SUBROUTINE spline


subroutine linte(f,dtl,idtl,nt)
  !
  !     This routine linearly interpolates f to fill in
  !     missing forecasts.
  !
  !     If idtl=1, then f is extrapolated
  !     rather than interpolated if the distance to land array
  !     (dtl) is negative for the point at t+12, but positive
  !     at t and t-12. This option is normally only used for
  !     interpolation of intensities.
  !
  implicit none
  integer, intent(in) :: idtl, nt
  real, dimension(0:nt), intent(inout) :: f, dtl
  ! local variables
  integer k
  real fp, fm, fk, dp, dm, dk

  if (nt .lt. 2) return

  if (idtl .eq. 1) then
      do k=2,nt-1
        fp = f(k+1)
        fm = f(k-1)
        fk = f(k)

        if (fk .le. 0.0) then
          dp = dtl(k+1)
          dm = dtl(k-1)
          dk = dtl(k)

          if (dk .gt. 0.0 .and. dp .lt. 0.0) then
              if (fm .gt. 0.0) f(k) = fm
          else
              if (fm .gt. 0.0 .and. fp .gt. 0.0) then
                f(k) = 0.5*(fp+fm)
              endif
          endif
        end if
      enddo
  else
      do k=2,nt-1
        fp = f(k+1)
        fm = f(k-1)
        fk = f(k)

        if ((fk .le. 0.0 .and. fp .gt. 0.0 &
                        .and. fm .gt. 0.0) .or. &
            (fk .ge. 0.0 .and. fp .lt. 0.0 &
                        .and. fm .lt. 0.0)) then
            f(k) = 0.5*(fp+fm)
        endif
      enddo
  endif

  return
end