c-----------------------------------------------------------------------
c
c  Name:
c      SUBROUTINE INTERP
c
c  Purpose:
c      Linearly interpolate an input array.
c
c  Usage:
c      CALL INTERP( NOLD, XOLD, YOLD, NNEW, XNEW, YNEW )
c 
c  Input:
c      NOLD    Number of elements in the input arrays
c      XOLD    Abscissa values for the input array (must be monotonic)
c      YOLD    Values of the input array
c      NNEW    Number of elements in the output arrays
c      XNEW    Abscissa values for the output array
c
c  Output:
c      YNEW    Linearly interpolated values of the output array
c
c  Revised:
c      Liam.Gumley@ssec.wisc.edu
c      $Id: interp.f,v 1.3 1999/11/11 20:04:45 gumley Exp $
c
c  Notes:
c      Values outside the range of the input abscissae are obtained
c      via extrapolation.
c
c-----------------------------------------------------------------------

      subroutine interp( nold, xold, yold, nnew, xnew, ynew )

      implicit none
      
      integer nold, nnew
      real xold( nold ), yold( nold ), xnew( nnew ), ynew (nnew )

      real slope, intercept
      integer lo, hi, j, init
      
      lo = 1
      hi = 2
      init = 1
      
      do j = 1, nnew
        
20      continue

c ...   check for monotonicity
      
        if( xold( lo ) .ge. xold( hi ) ) then
          write(*,*) 'Fatal Error in subroutine INTERP'
          write(*,*) 'Input abscissa array XOLD is not monotonic'
          write(*,*) 'Offending XOLD values are', xold( lo ), xold( hi )
          stop
        endif

c ...   check if output point falls between current input points
                
        if( xnew( j ) .gt. xold( hi ) ) then
          if( hi .lt. nold ) then
            lo = lo + 1
            hi = hi + 1
            init = 1
            goto 20
          endif
        endif

c ...   compute slope and intercept only when necessary

        if( init .eq. 1 ) then
          slope = ( yold( hi ) - yold( lo ) ) /
     &      ( xold( hi ) - xold( lo ) )
          intercept = yold( lo ) - slope * xold( lo )
          init = 0
        endif

c ...   compute output value

        ynew( j ) = slope * xnew( j ) + intercept

      end do
        
      end