module refindex

  use vtype_module
  use ni3_module

  implicit NONE

  private

  public :: RefIndexInit
  public :: RefIndexFinalize
  public :: RefIndexCalc

  integer, save        :: NRITbl
  real(8), allocatable :: w_RITblWl  (:)
  real(8), allocatable :: w_RITblReal(:)
  real(8), allocatable :: w_RITblImag(:)

contains

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

  subroutine RefIndexInit( &
    & NcFN &
    & )

    character(*), intent(in) :: NcFN


    ! Local variables
    !
    character(stdstr) :: Mode
    integer           :: NcID


    Mode = "read"
    call ni3_open( NcFN, Mode, NcID )
    call ni3_inq_dimlen( NcID, "wl", NRITbl )
    allocate( w_RITblWl  (NRITbl) )
    allocate( w_RITblReal(NRITbl) )
    allocate( w_RITblImag(NRITbl) )
    call ni3_get_var( NcID, "wl", w_RITblWl )
    call ni3_get_var( NcID, "refindexreal", w_RITblReal )
    call ni3_get_var( NcID, "refindeximag", w_RITblImag )
    call ni3_close( NcID )


  end subroutine RefIndexInit

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

  subroutine RefIndexFinalize


    deallocate( w_RITblWl   )
    deallocate( w_RITblReal )
    deallocate( w_RITblImag )


  end subroutine RefIndexFinalize

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

  subroutine RefIndexCalc( &
    & Wl,                  & ! (in)
    & RIReal, RIImag       & ! (out)
    & )

    real(8), intent(in ) :: Wl
    real(8), intent(out) :: RIReal
    real(8), intent(out) :: RIImag


    integer :: l


    if ( Wl < w_RITblWl(1) ) then
      write( 6, * ) 'Specified wavelength is out of range.'
      write( 6, * ) 'Wavelength has to be larger than ', &
        & w_RITblWl(1)
      stop
    end if
    if ( Wl > w_RITblWl(NRITbl) ) then
      write( 6, * ) 'Specified wavelength is out of range.'
      write( 6, * ) 'Wavelength has to be less than ', &
        & w_RITblWl(NRITbl)
      stop
    end if

    do l = 2, NRITbl
      if ( Wl <= w_RITblWl(l) ) exit
    end do

    ! see http://www.atmos.washington.edu/ice_optical_constants/
    ! for selection of interpolation method
!!$    RIReal =                                    &
!!$      &   ( w_RITblReal(l) - w_RITblReal(l-1) ) &
!!$      & / ( w_RITblWl  (l) - w_RITblWl  (l-1) ) &
!!$      & * ( Wl             - w_RITblWl  (l-1) ) &
!!$      & + w_RITblReal(l-1)
!!$    RIImag=                                     &
!!$      &   ( w_RITblImag(l) - w_RITblImag(l-1) ) &
!!$      & / ( w_RITblWl  (l) - w_RITblWl  (l-1) ) &
!!$      & * ( Wl             - w_RITblWl  (l-1) ) &
!!$      & + w_RITblImag(l-1)
    RIReal =                                       &
      &      ( w_RITblReal(l) - w_RITblReal(l-1) ) &
      & / log( w_RITblWl  (l) / w_RITblWl  (l-1) ) &
      & * log( Wl             / w_RITblWl  (l-1) ) &
      & + w_RITblReal(l-1)
    RIImag =                                       &
      &   log( w_RITblImag(l) / w_RITblImag(l-1) ) &
      & / log( w_RITblWl  (l) / w_RITblWl  (l-1) ) &
      & * log( Wl             / w_RITblWl  (l-1) ) &
      & + log( w_RITblImag(l-1) )
    RIImag = exp( RIImag )


  end subroutine RefIndexCalc

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

end module refindex
