! $Id$
module lmdz_interpolation

  ! From Press et al., 1996, version 2.10a
  ! B3 Interpolation and Extrapolation

  IMPLICIT NONE; PRIVATE
  PUBLIC locate, hunt

CONTAINS

  pure FUNCTION locate(xx, x)

    REAL, DIMENSION(:), INTENT(IN) :: xx
    REAL, INTENT(IN) :: x
    INTEGER  locate

    ! Given an array xx(1:N), and given a value x, returns a value j,
    ! between 0 and N, such that x is between xx(j) and xx(j + 1). xx
    ! must be monotonic, either increasing or decreasing. j = 0 or j =
    ! N is returned to indicate that x is out of range. This
    ! procedure should not be called with a zero-sized array argument.
    ! See notes.

    INTEGER  n, jl, jm, ju
    LOGICAL  ascnd

    !----------------------------

    n = size(xx)
    ascnd = (xx(n) >= xx(1))
    ! (True if ascending order of table, false otherwise.)
    ! Initialize lower and upper limits:
    jl = 0
    ju = n + 1
    do while (ju - jl > 1)
      jm = (ju + jl) / 2 ! Compute a midpoint,
      IF (ascnd .eqv. (x >= xx(jm))) THEN
        jl = jm ! and replace either the lower limit
      else
        ju = jm ! or the upper limit, as appropriate.
      end if
    END DO
    ! {ju == jl + 1}

    ! {(ascnd .AND. xx(jl) <= x < xx(jl+1))
    !  .neqv. 
    !  (.NOT. ascnd .AND. xx(jl+1) <= x < xx(jl))}

    ! Then set the output, being careful with the endpoints:
    IF (x == xx(1)) THEN
      locate = 1
    ELSE IF (x == xx(n)) THEN
      locate = n - 1
    else
      locate = jl
    end if

  END FUNCTION locate

  !***************************

  pure SUBROUTINE hunt(xx, x, jlo)

    ! Given an array xx(1:N ), and given a value x, returns a value
    ! jlo such that x is between xx(jlo) and xx(jlo+1). xx must be
    ! monotonic, either increasing or decreasing. jlo = 0 or jlo = N is
    ! returned to indicate that x is out of range. jlo on input is taken as
    ! the initial guess for jlo on output.
    ! Modified so that it uses the information "jlo = 0" on input.

    INTEGER, INTENT(INOUT) :: jlo
    REAL, INTENT(IN) :: x
    REAL, DIMENSION(:), INTENT(IN) :: xx
    INTEGER  n, inc, jhi, jm
    LOGICAL  ascnd, hunt_up

    !-----------------------------------------------------

    n = size(xx)
    ascnd = (xx(n) >= xx(1))
    ! (True if ascending order of table, false otherwise.)
    IF (jlo < 0 .OR. jlo > n) THEN
      ! Input guess not useful. Go immediately to bisection.
      jlo = 0
      jhi = n + 1
    else
      inc = 1 ! Set the hunting increment.
      IF (jlo == 0) THEN
        hunt_up = .TRUE.
      else
        hunt_up = x >= xx(jlo) .eqv. ascnd
      end if
      IF (hunt_up) then ! Hunt up:
        DO
          jhi = jlo + inc
          IF (jhi > n) then ! Done hunting, since off end of table.
            jhi = n + 1
            exit
          else
            IF (x < xx(jhi) .eqv. ascnd) exit
            jlo = jhi ! Not done hunting,
            inc = inc + inc ! so double the increment
          end if
        END DO ! and try again.
      else ! Hunt down:
        jhi = jlo
        DO
          jlo = jhi - inc
          IF (jlo < 1) then ! Done hunting, since off end of table.
            jlo = 0
            exit
          else
            IF (x >= xx(jlo) .eqv. ascnd) exit
            jhi = jlo ! Not done hunting,
            inc = inc + inc ! so double the increment
          end if
        END DO ! and try again.
      end if
    end if ! Done hunting, value bracketed.

    do ! Hunt is done, so begin the final bisection phase:
      IF (jhi - jlo <= 1) THEN
        IF (x == xx(n)) jlo = n - 1
        IF (x == xx(1)) jlo = 1
        exit
      else
        jm = (jhi + jlo) / 2
        IF (x >= xx(jm) .eqv. ascnd) THEN
          jlo = jm
        else
          jhi = jm
        end if
      end if
    END DO

  END SUBROUTINE hunt

END MODULE lmdz_interpolation