module opt_prop_kddevelop

  use vtype_module

  implicit none

  private

  public :: OptPropKDDevelopInit
  public :: OptPropKDDevelopFinalize
  public :: OptPropKDDevelopGetNWN
  public :: OptPropKDDevelopGetNPress
  public :: OptPropKDDevelopGetNMol
  public :: OptPropKDDevelopGetNBand
  public :: OptPropKDDevelopGetDelWN
  public :: OptPropKDDevelopGetBandNum
  public :: OptPropKDDevelopGetBandWNBnds
  public :: OptPropKDDevelopGetBandBinIndexBnds
!!$  public :: OptPropKDDevelopGetTotalStrFlux
  public :: OptPropKDDevelopGetAbsCoefProf
  public :: OptPropKDDevelopGetRayScatCoef
  public :: OptPropKDDevelopGetStrBandAvePF
  public :: OptPropKDDevelopGetPFIntedRatio
  public :: OptPropKDDevelopGetStrPFIntedRatio


  ! Common
  character(128)       , save :: PressAxisName
  character(128)       , save :: MolAxisName
  character(128)       , save :: PressForPFAxisName
  !
  integer              , save :: NPressSave
  integer              , save :: NMolSave
  integer              , save :: NPressForPFSave
  !
  real(DP), allocatable, save :: p_PressSave(:)
  integer , allocatable, save :: m_MolNumSave(:)

  ! LBL
  integer              , save :: NWaveNumSave

!!$  real(DP), allocatable, save :: b_BandAveStrPFSave(:)
  real(DP), allocatable, save :: b_RatioBandIntFluxtoTotFluxSave(:)

  real(DP), allocatable, save :: aa_BandWaveNumBndsSave     (:,:)

  ! KDDevelop
  integer , parameter         :: NBndSave = 2
  integer              , save :: NBinSave
  integer              , save :: NBandSave
  integer , allocatable, save :: a_BinNumSave(:)
!!$  integer , allocatable, save :: aa_BandWaveNumIndexBndsSave(:,:)
  integer , allocatable, save :: aa_BandBinIndexBndsSave    (:,:)
  !
  character(128), save :: WaveNumAxisName
  character(128), save :: AbsCoefVarName
  character(128), save :: AtmPFRatioVarName
  character(128), save :: StrPFRatioVarName
  character(128), save :: StrSpeVarName
  !
  real(DP), allocatable, save :: a_DelWaveNum(:)
  !
  ! KDDevelop only
  real(DP), allocatable, save :: ra_AtmPFRatio      (:,:)
  real(DP), allocatable, save :: zma_AbsCoef     (:,:,:)

  real(DP), allocatable, save :: b_BandAveRayScatCoefNonRadAct(:)
  real(DP), allocatable, save :: bm_BandAveRayScatCoef(:,:)

  integer , allocatable, save :: ma_IDAbsCoefType(:,:)

  real(DP), allocatable, save :: a_StrPFRatio(:)

  character(128), save :: ModuleName = 'opt_prop_kddevelop'

  logical, save :: FlagInited = .false.

contains

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

  subroutine OptPropKDDevelopInit( &
    & OptPropNcFN, RayScatCoefNcFN, StrSpeNcFN &
    & )

    use ni3_module


    character(*), intent(in) :: OptPropNcFN
    character(*), intent(in) :: RayScatCoefNcFN
    character(*), intent(in) :: StrSpeNcFN


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

    real(DP), allocatable :: ra_StrPFRatio(:,:)
!!$    real(DP), allocatable :: pb_BandAveStrPF(:,:)
    real(DP), allocatable :: rb_RatioBandIntFluxtoTotFlux(:,:)

    integer :: iWaveNum
    integer :: iBand

    integer              :: TmpNum
    integer, allocatable :: m_MolNumTmp(:)

    integer :: m


!!$    if ( .not. FlagInited ) then
!!$      write( 6, * ) trim( ModuleName ), " is not initialized."
!!$      stop
!!$    end if


    NPressSave      = -1
    NMolSave        = -1
    NWaveNumSave    = -1
    NBinSave        = -1
    NBandSave       = -1

    NPressForPFSave = -1


    PressAxisName      = "Press"
    MolAxisName        = "MolNum"
    PressForPFAxisName = "PressForPF"

    Mode = "read"
    call ni3_open( OptPropNcFN, Mode, NcID )
    call ni3_inq_dimlen( NcID, PressAxisName, NPressSave )
    call ni3_inq_dimlen( NcID, MolAxisName  , NMolSave   )
    allocate( p_PressSave (1:NPressSave) )
    allocate( m_MolNumSave(1:NMolSave  ) )
    call ni3_get_var( NcID, PressAxisName, p_PressSave  )
    call ni3_get_var( NcID, MolAxisName  , m_MolNumSave )

    WaveNumAxisName = 'BinNum'
    !
    call ni3_inq_dimlen( NcID, WaveNumAxisName, NBinSave  )
    call ni3_inq_dimlen( NcID, "BandNum"      , NBandSave )
    allocate( a_BinNumSave               (NBinSave) )
    allocate( aa_BandWaveNumBndsSave     (NBndSave,NBandSave) )
    allocate( aa_BandBinIndexBndsSave    (NBndSave,NBandSave) )

    call ni3_get_var( NcID, WaveNumAxisName   , a_BinNumSave            )
    call ni3_get_var( NcID, 'BandWaveNumBnds' , aa_BandWaveNumBndsSave  )
    call ni3_get_var( NcID, 'BandBinIndexBnds', aa_BandBinIndexBndsSave )
    !
    AbsCoefVarName    = 'BinAbsCoef'
    AtmPFRatioVarName = 'BinAtmPFRatio'
    StrPFRatioVarName = 'BinStrPFRatio'

    allocate( a_DelWaveNum(NBinSave) )
    call ni3_get_var( NcID, "BinWeight", a_DelWaveNum )

    allocate( zma_AbsCoef (NPressSave,NMolSave,NBinSave) )
    call ni3_get_var( NcID, AbsCoefVarName, zma_AbsCoef )

    allocate( ma_IDAbsCoefType(NMolSave,NBinSave) )
    call ni3_get_var( NcID, "BinIDAbsCoefType", ma_IDAbsCoefType )

    call ni3_inq_dimlen( NcID, PressForPFAxisName, NPressForPFSave )
    allocate( ra_AtmPFRatio(NPressForPFSave,NBinSave) )
    call ni3_get_var( NcID, AtmPFRatioVarName, ra_AtmPFRatio )
    allocate( ra_StrPFRatio(NPressForPFSave,NBinSave) )
    call ni3_get_var( NcID, StrPFRatioVarName, ra_StrPFRatio )
    allocate( a_StrPFRatio(NBinSave) )
    do iWaveNum = 1, NBinSave
      a_StrPFRatio(iWaveNum) = ra_StrPFRatio(1,iWaveNum)
    end do
    deallocate( ra_StrPFRatio )

    call ni3_close( NcID )


    if ( len_trim( RayScatCoefNcFN ) > 0 ) then
      Mode = "read"
      call ni3_open( RayScatCoefNcFN, Mode, NcID )
      !
      call ni3_inq_dimlen( NcID, "BandNum", TmpNum )
      if ( TmpNum /= NBandSave ) then
        write( 6, * ) 'Array size for stellar spectrum is inappropriate: '
        write( 6, * ) '  ', trim( StrSpeNcFN  ), ': ', TmpNum
        write( 6, * ) '  ', trim( OptPropNcFN ), ': ', NBandSave
        stop
      end if
      call ni3_inq_dimlen( NcID, MolAxisName, TmpNum )
      if ( TmpNum /= NMolSave ) then
        write( 6, * ) 'MolNum size for Rayleigh scattering coefficient is inappropriate:'
        write( 6, * ) '  ', trim( RayScatCoefNcFN ), ': ', TmpNum
        write( 6, * ) '  ', trim( OptPropNcFN     ), ': ', NMolSave
        stop
      end if
      !
      allocate( m_MolNumTmp(1:NMolSave) )
      call ni3_get_var( NcID, MolAxisName, m_MolNumTmp )
      do m = 1, NMolSave
        if ( m_MolNumTmp(m) /= m_MolNumSave(m) ) then
          write( 6, * ) 'Molecular number ', m, ' is different between files for gas and particles.'
          write( 6, * ) '  ', trim( RayScatCoefNcFN ), ': ', m_MolNumTmp
          write( 6, * ) '  ', trim( OptPropNcFN     ), ': ', m_MolNumSave
          stop
        end if
      end do
      deallocate( m_MolNumTmp )
      !
      allocate( b_BandAveRayScatCoefNonRadAct(NBandSave) )
      allocate( bm_BandAveRayScatCoef(NBandSave,NMolSave) )
      call ni3_get_var( NcID, "BandAveRayScatCoefNonRadAct", b_BandAveRayScatCoefNonRadAct )
      call ni3_get_var( NcID, "BandAveRayScatCoef"         , bm_BandAveRayScatCoef )
      !
      call ni3_close( NcID )
    else
      allocate( b_BandAveRayScatCoefNonRadAct(NBandSave) )
      allocate( bm_BandAveRayScatCoef        (NBandSave,NMolSave) )
      b_BandAveRayScatCoefNonRadAct = 0.0d0
      bm_BandAveRayScatCoef         = 0.0d0
    end if


!!$    allocate( b_BandAveStrPFSave( NBandSave ) )
!!$    allocate( pb_BandAveStrPF   ( NPressSave, NBandSave ) )
    allocate( b_RatioBandIntFluxtoTotFluxSave( NBandSave ) )
    allocate( rb_RatioBandIntFluxtoTotFlux   ( NPressForPFSave, NBandSave ) )

!!$    StrSpeVarName   = 'BandAveStrPF'
    StrSpeVarName   = 'RatioBandIntFluxtoTotFlux'
    !
    Mode = "read"
    call ni3_open( StrSpeNcFN, Mode, NcID )
    call ni3_inq_dimlen( NcID, "BandNum", TmpNum )
    if ( TmpNum /= NBandSave ) then
      write( 6, * ) 'Array size for stellar spectrum is inappropriate: '
      write( 6, * ) '  ', trim( StrSpeNcFN  ), ': ', TmpNum
      write( 6, * ) '  ', trim( OptPropNcFN ), ': ', NBandSave
      stop
    end if
!!$    call ni3_get_var( NcID, "BandAveStrPF", pb_BandAveStrPF )
    call ni3_get_var( NcID, "RatioBandIntFluxtoTotFlux", rb_RatioBandIntFluxtoTotFlux )
    call ni3_close( NcID )
    !
    do iBand = 1, NBandSave
!!$      b_BandAveStrPFSave(iBand) = pb_BandAveStrPF(1,iBand)
      b_RatioBandIntFluxtoTotFluxSave(iBand) = rb_RatioBandIntFluxtoTotFlux(1,iBand)
    enddo
    !
!!$    deallocate( pb_BandAveStrPF )
    deallocate( rb_RatioBandIntFluxtoTotFlux )


    FlagInited = .true.


  end subroutine OptPropKDDevelopInit

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

  subroutine OptPropKDDevelopFinalize


    ! local variables
    !


    if ( .not. FlagInited ) then
      write( 6, * ) trim( ModuleName ), " is not initialized."
      stop
    end if


    NPressSave      = -1
    NMolSave        = -1
    NWaveNumSave    = -1
    NBinSave        = -1
    NBandSave       = -1

    NPressForPFSave = -1


    deallocate( p_PressSave  )
    deallocate( m_MolNumSave )

    NBinSave    = 0
    deallocate( a_BinNumSave                )
    deallocate( aa_BandWaveNumBndsSave      )
!!$      deallocate( aa_BandWaveNumIndexBndsSave )
    deallocate( aa_BandBinIndexBndsSave     )

    deallocate( a_DelWaveNum )

    deallocate( zma_AbsCoef )

    deallocate( b_BandAveRayScatCoefNonRadAct )
    deallocate( bm_BandAveRayScatCoef         )

    deallocate( ma_IDAbsCoefType )

    deallocate( ra_AtmPFRatio )
    deallocate( a_StrPFRatio )

!!$    deallocate( b_BandAveStrPFSave )
    deallocate( b_RatioBandIntFluxtoTotFluxSave )


  end subroutine OptPropKDDevelopFinalize

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

  function OptPropKDDevelopGetNWN() result( NWN )


    integer :: NWN

    ! local variables
    !


    if ( .not. FlagInited ) then
      write( 6, * ) trim( ModuleName ), " is not initialized."
      stop
    end if


    NWN = NBinSave


  end function OptPropKDDevelopGetNWN

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

  function OptPropKDDevelopGetNPress() result( NPress )


    integer :: NPress

    ! local variables
    !


    if ( .not. FlagInited ) then
      write( 6, * ) trim( ModuleName ), " is not initialized."
      stop
    end if


    NPress = NPressSave


  end function OptPropKDDevelopGetNPress

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

  function OptPropKDDevelopGetNMol() result( NMol )


    integer :: NMol

    ! local variables
    !


    if ( .not. FlagInited ) then
      write( 6, * ) trim( ModuleName ), " is not initialized."
      stop
    end if


    NMol = NMolSave


  end function OptPropKDDevelopGetNMol

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

  function OptPropKDDevelopGetNBand() result( NBand )


    integer :: NBand

    ! local variables
    !


    if ( .not. FlagInited ) then
      write( 6, * ) trim( ModuleName ), " is not initialized."
      stop
    end if


    NBand = NBandSave


  end function OptPropKDDevelopGetNBand

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

  function OptPropKDDevelopGetDelWN( iWaveNum ) result( DelWaveNum )


    integer, intent(in) :: iWaveNum

    real(DP) :: DelWaveNum

    ! local variables
    !


    if ( .not. FlagInited ) then
      write( 6, * ) trim( ModuleName ), " is not initialized."
      stop
    end if


    DelWaveNum = a_DelWaveNum(iWaveNum)


  end function OptPropKDDevelopGetDelWN

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

  function OptPropKDDevelopGetBandNum( iWaveNum ) result( BandNum )


    integer, intent(in) :: iWaveNum

    integer :: BandNum


    ! local variables
    !
    integer :: n


    if ( .not. FlagInited ) then
      write( 6, * ) trim( ModuleName ), " is not initialized."
      stop
    end if


    do n = 1, NBandSave
      if ( ( aa_BandBinIndexBndsSave(1,n) <= iWaveNum ) .and. &
        &  ( iWaveNum <= aa_BandBinIndexBndsSave(2,n) ) ) &
        exit
    end do

    BandNum = n


  end function OptPropKDDevelopGetBandNum

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

  subroutine OptPropKDDevelopGetBandWNBnds( &
    & NBand, &
    & aa_BandWNBnds &
    & )


    integer , intent(in ) :: NBand
    real(DP), intent(out) :: aa_BandWNBnds(2,NBand)

    ! local variables
    !


    if ( .not. FlagInited ) then
      write( 6, * ) trim( ModuleName ), " is not initialized."
      stop
    end if


    ! check
    if ( NBand /= NBandSave ) then
      write( 6, * ) 'Number of bands is inappropriate: ', NBand, NBandSave
      stop
    end if
    aa_BandWNBnds = aa_BandWaveNumBndsSave


  end subroutine OptPropKDDevelopGetBandWNBnds

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

  subroutine OptPropKDDevelopGetBandBinIndexBnds( &
    & NBand, &
    & aa_BandBinIndexBnds &
    & )

    integer, intent(in ) :: NBand
    integer, intent(out) :: aa_BandBinIndexBnds(2,NBand)


    ! local variables
    !


    if ( .not. FlagInited ) then
      write( 6, * ) trim( ModuleName ), " is not initialized."
      stop
    end if


    ! check
    if ( NBand /= NBandSave ) then
      write( 6, * ) 'Number of bands is appropriate: ', NBand, NBandSave
      stop
    end if
    aa_BandBinIndexBnds = aa_BandBinIndexBndsSave


  end subroutine OptPropKDDevelopGetBandBinIndexBnds

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

!!$  subroutine OptPropKDDevelopGetTotalStrFlux( &
!!$    & TotalStrFlux &
!!$    & )
!!$
!!$    real(DP), intent(out) :: TotalStrFlux
!!$
!!$
!!$    ! local variables
!!$    !
!!$    real(DP) :: DelWaveNum
!!$
!!$    integer  :: l
!!$
!!$
!!$    TotalStrFlux = 0.0d0
!!$    do l = 1, NBandSave
!!$      DelWaveNum = aa_BandWaveNumBndsSave(2,l) - aa_BandWaveNumBndsSave(1,l)
!!$      TotalStrFlux = TotalStrFlux &
!!$        & + b_BandAveStrPFSave( l ) * DelWaveNum
!!$    end do
!!$
!!$
!!$  end subroutine OptPropKDDevelopGetTotalStrFlux

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

  subroutine OptPropKDDevelopGetAbsCoefProf( &
    & iWaveNum, kmax, NMol, m_MolNum, pm_VMR, &
    & pm_AbsCoef &
    & )


    integer     , intent(in ) :: iWaveNum
    integer     , intent(in ) :: kmax
    integer     , intent(in ) :: NMol
    integer     , intent(in ) :: m_MolNum  (1:NMol)
    real(DP)    , intent(in ) :: pm_VMR    (1:kmax,1:NMol)
    real(DP)    , intent(out) :: pm_AbsCoef(1:kmax,1:NMol)


    ! local variables
    !
!!$    integer        :: k
    integer        :: m


    if ( .not. FlagInited ) then
      write( 6, * ) trim( ModuleName ), " is not initialized."
      stop
    end if


    ! check
    if ( kmax /= NPressSave ) then
      write( 6, * ) 'In ', trim( ModuleName )
      write( 6, * ) '  Number of layers does not match with number of pressure levels in optical properties: ', &
        & kmax, NPressSave
      stop
    end if
!!$      do k = 1, kmax
!!$        if ( abs( p_Press(k) - p_PressSave(k) ) > p_PressSave(k) * 1.0d-6 ) then
!!$          write( 6, * ) 'In ', trim( ModuleName )
!!$          write( 6, * ) '  Pressure is inappropriate: ', p_PressSave(k), p_Press(k)
!!$          stop
!!$        end if
!!$      end do


    if ( NMol /= NMolSave ) then
      write( 6, * ) 'NMol in arguments does not match with NMol in absorption coefficient file.'
      write( 6, * ) NMol, NMolSave
      stop
    end if
    do m = 1, NMol
      if ( m_MolNum(m) /= m_MolNumSave(m) ) then
        write( 6, * ) 'Molecular number is inappropriate: ', &
          & m, m_MolNum(m), m_MolNumSave(m)
        stop
      end if
    end do


    pm_AbsCoef = zma_AbsCoef(:,:,iWaveNum)

    do m = 1, NMolSave
      select case ( ma_IDAbsCoefType(m,iWaveNum) )
      case ( -1 )
        pm_AbsCoef(:,m) = 0.0d0
      case ( 0 )
        pm_AbsCoef(:,m) = pm_AbsCoef(:,m) * pm_VMR(:,m)
!!$          case ( 1 )
!!$            pm_AbsCoef(:,m) = pm_AbsCoef(:,m)
      end select
    end do


  end subroutine OptPropKDDevelopGetAbsCoefProf

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

  subroutine OptPropKDDevelopGetRayScatCoef( &
    & iWaveNum, &
    & NMol, m_MolNum, &
    & RayScatCoefNonRadAct, m_RayScatCoef &
    & )


    integer     , intent(in ) :: iWaveNum
    integer     , intent(in ) :: NMol
    integer     , intent(in ) :: m_MolNum(NMol)
    real(DP)    , intent(out) :: RayScatCoefNonRadAct
    real(DP)    , intent(out) :: m_RayScatCoef(NMol)


    ! local variables
    !
    integer :: iBand

    integer :: m


    if ( .not. FlagInited ) then
      write( 6, * ) trim( ModuleName ), " is not initialized."
      stop
    end if


    if ( NMol /= NMolSave ) then
      write( 6, * ) 'NMol in arguments does not match with NMol in absorption coefficient file.'
      write( 6, * ) NMol, NMolSave
      stop
    end if
    do m = 1, NMol
      if ( m_MolNum(m) /= m_MolNumSave(m) ) then
        write( 6, * ) 'Molecular number is inappropriate: ', &
          & m, m_MolNum(m), m_MolNumSave(m)
        stop
      end if
    end do


    iBand = OptPropKDDevelopGetBandNum( iWaveNum )
    RayScatCoefNonRadAct = b_BandAveRayScatCoefNonRadAct( iBand )
    do m = 1, NMol
      m_RayScatCoef(m) = bm_BandAveRayScatCoef( iBand, m )
    end do


  end subroutine OptPropKDDevelopGetRayScatCoef

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

  subroutine OptPropKDDevelopGetStrBandAvePF( &
    & iBand, StrFluxTOA, &
    & StrBandAvePF &
    & )

    integer , intent(in ) :: iBand
    real(DP), intent(in ) :: StrFluxTOA
    real(DP), intent(out) :: StrBandAvePF


    ! local variables
    !


!!$    StrBandAvePF = b_BandAveStrPFSave( iBand )
    StrBandAvePF = b_RatioBandIntFluxtoTotFluxSave( iBand )
    StrBandAvePF = StrBandAvePF * StrFluxTOA
    StrBandAvePF = StrBandAvePF &
      & / ( aa_BandWaveNumBndsSave(2,iBand) - aa_BandWaveNumBndsSave(1,iBand) )


  end subroutine OptPropKDDevelopGetStrBandAvePF

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

  subroutine OptPropKDDevelopGetPFIntedRatio( &
    & iWaveNum, kmax, &
    & r_AtmPFRatio, SurfPFRatio &
    & )


    integer     , intent(in ) :: iWaveNum
    integer     , intent(in ) :: kmax
    real(DP)    , intent(out) :: r_AtmPFRatio(1:kmax)
    real(DP)    , intent(out) :: SurfPFRatio


    ! local variables
    !

!!$    real(DP), allocatable :: za_PFRatio(:,:)


    if ( .not. FlagInited ) then
      write( 6, * ) trim( ModuleName ), " is not initialized."
      stop
    end if


    ! check
    if ( kmax /= NPressForPFSave ) then
      write( 6, * ) 'In ', trim( ModuleName )
      write( 6, * ) '  Number of layers does not match with number of pressure levels in optical properties: ', &
        & kmax, NPressSave
      stop
    end if


    r_AtmPFRatio = ra_AtmPFRatio(:,iWaveNum)


!!$    allocate( za_PFRatio( NPressSave, NWaveNumSave ) )
!!$    call HistoryGet( OptPropNcFNSave, PFRatioVarName, za_PFRatio, &
!!$      & quiet = .true. )
!!$    z_PFRatio = za_PFRatio(:,iWaveNum)
!!$    deallocate( za_PFRatio )


    SurfPFRatio = r_AtmPFRatio(1)


  end subroutine OptPropKDDevelopGetPFIntedRatio

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

  subroutine OptPropKDDevelopGetStrPFIntedRatio( &
    & iWaveNum, &
    & StrPFRatio &
    & )


    integer     , intent(in ) :: iWaveNum
    real(DP)    , intent(out) :: StrPFRatio


    ! local variables
    !

!!$    real(DP), allocatable :: za_PFRatio(:,:)


    if ( .not. FlagInited ) then
      write( 6, * ) trim( ModuleName ), " is not initialized."
      stop
    end if


    StrPFRatio = a_StrPFRatio(iWaveNum)


  end subroutine OptPropKDDevelopGetStrPFIntedRatio

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

end module opt_prop_kddevelop
