c-----------------------------------------------------------------------
c
c  Name:
c      SUBROUTINE HMWSRF
c
c  Purpose:
c      Read a spectral response data file and linearly interpolate
c      the spectral response to 0.1 inverse centimeter resolution.
c
c  Usage:
c      CALL HMWSRF( FILE, N, W1, W2, SRF )
c 
c  Input:
c      FILE    Name of spectral response data file
c              The input spectral response data file should be ASCII
c              format with two columns. The first column should be
c              wavelength in microns, and the second column should
c              be spectral response (normalized to 1.0 or un-normalized
c              are both acceptable).
c              Wavelength values in the range [0.4-20.0] are accepted.
c              Wavelength values outside this range will cause an exit.      
c              Negative response values in the range [-0.01,0.0] will be
c              set to zero. Negative values less than -0.01 will cause
c              an exit.
c
c  Output:
c      N       Number of points in interpolated output array
c      W1      Lower wavenumber limit of interpolated output array
c              (inverse centimeters)
c      W2      Upper wavenumber limit of interpolated output array
c              (inverse centimeters)
c      SRF     Linearly interpolated spectral response at
c              0.1 inverse centimeter resolution, normalized so that
c              maximum respose is 1.0
c
c  Revised:
c      Liam.Gumley@ssec.wisc.edu
c      $Id: hmwsrf.f,v 1.3 1999/11/11 20:04:44 gumley Exp $
c
c  Notes:
c      (1) Calls the linear interpolation subroutine INTERP.
c
c-----------------------------------------------------------------------

      subroutine hmwsrf( file, n, w1, w2, srf )
      
      implicit none

c ... arguments
      
      character*(*) file
      integer n
      real w1, w2, srf( * )

c ... local variables
      
      real xold( 10000 ), yold( 10000 ), xnew( 10000 ), temp( 10000 ),
     &     big
      integer nold, ios, i, j, lun, getlun
      external getlun

c ... open input file
      
      lun = getlun()
      if( lun .le. 0 ) then
        write(*,*) 'Fatal Error in subroutine HMWSRF'
        write(*,*) 'Could not get a free logical unit number'
        stop
      endif
      open( unit = lun, file = file, status = 'old', iostat = ios )
      if( ios .ne. 0 ) then
        write(*,*) 'Fatal Error in subroutine HMWSRF'
        write(*,*) 'Could not open input file'
        write(*,*) 'Input filename is ', file
        stop
      endif

c ... read wavelength and spectral response values and check ranges
      
      nold = 0
20    continue
      read( unit = lun, fmt = *, iostat = ios )
     &  xold( nold + 1 ), yold( nold + 1 )
      if( ios .ne. 0 ) goto 40
      
      nold = nold + 1
      
      if( xold( nold ) .le. 0.4 .or. xold( nold ) .gt. 20.0 ) then
        write(*,*) 'Fatal Error in subroutine HMWSRF'
        write(*,*) 'Input wavelength is outside the range [0.4,20.0]'
        write(*,*) 'Offending value is at element ', nold,
     &             ', value is ', xold( nold )  
        write(*,*) 'Input filename is ', file
        stop
      endif
      
      if( yold( nold ) .lt. -0.01 ) then
        write(*,*) 'Fatal Error in subroutine HMWSRF'
        write(*,*) 'Input response is less than -0.01'
        write(*,*) 'Offending value is at element ', nold,
     &             ', value is ', yold( nold )  
        write(*,*) 'Input filename is ', file
        stop
      endif
      yold( nold ) = max( yold( nold ), 0.0 )
      
      goto 20

40    continue

      close( unit = lun )
      
c ... reverse order of spectral response data

      j = 1
      do i = nold, 1, -1
        temp( j ) = xold( i )
        j = j + 1
      end do
      do i = 1, nold
        xold( i ) = temp( i )
      end do
      
      j = 1
      do i = nold, 1, -1
        temp( j ) = yold( i )
        j = j + 1
      end do
      do i = 1, nold
        yold( i ) = temp( i )
      end do

c ... convert wavelength (microns) to wavenumber

      do i = 1, nold
        xold( i ) = 1.0e4 / xold( i )
      end do
      
c ... construct new wavenumber array at 0.1 inverse centimeter intervals

      w1 = xold( 1 ) - mod( xold( 1 ), 0.1 ) - 0.05
      w2 = xold( nold ) - mod( xold( nold ), 0.1 ) + 0.15
      n = int( ( w2 - w1 ) / 0.1 ) + 1
      do i = 1, n
        xnew( i ) = ( w2 - w1 ) * real( i - 1 ) / real( n - 1 ) + w1
      end do
      
c ... interpolate to new wavenumber array spacing

      call interp( nold, xold, yold, n, xnew, srf )
      srf( 1 ) = 0.0
      srf( n ) = 0.0

c ... normalize to maximum value of 1.0

      big = -1.0e8
      do i = 1, n
        big = max( big, srf( i ) )
      end do
      do i = 1, n
        srf( i ) = srf( i ) / big
      end do
      
      end