module opt_prop_lbl

  use vtype_module

  implicit none

  private

  public :: OptPropLBLInit
  public :: OptPropLBLFinalize
  public :: OptPropLBLGetNWN
  public :: OptPropLBLGetNPress
  public :: OptPropLBLGetNMol
  public :: OptPropLBLGetDelWN
  public :: OptPropLBLGetWN
  public :: OptPropLBLGetTotalStrFlux
  public :: OptPropLBLGetAbsCoefProf
  public :: OptPropLBLGetRayScatCoef
  public :: OptPropLBLGetPtclParam
  public :: OptPropLBLGetPFInted
  public :: OptPropLBLGetStrPFInted
  public :: OptPropLBLGetPFIntedRatio
  public :: OptPropLBLGetStrPFIntedRatio

  public :: OptPropLBLGetInputedOptDepProf


  ! Common
  integer              , save :: NPressSave
  integer              , save :: NMolSave
  !
  real(DP), allocatable, save :: p_PressSave(:)
  integer , allocatable, save :: m_MolNumSave(:)

  ! LBL
  integer              , save :: NWaveNumSave
  real(DP), allocatable, save :: w_WaveNumSave(:)
  real(DP), allocatable, save :: w_StrSpeSave(:)
  real(DP)             , save :: StrSpeIntedSave

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

  real(DP), allocatable, save :: w_RayScatCoefNonRadActSave(:)
  real(DP), allocatable, save :: wm_RayScatCoefSave(:,:)

  integer(i4b)                  , save :: NPtclSave
  character(stdstr), allocatable, save :: a_PtclNameSave(:)
  integer(i4b)                  , save :: PtclNRadiusMaxSave
  integer          , allocatable, save :: a_PtclNRadiusSave(:)
  real(dp)         , allocatable, save :: ra_PtclRadiusSave(:,:)
  real(dp)         , allocatable, save :: wra_QExtSave(:,:,:)
  real(dp)         , allocatable, save :: wra_SSASave (:,:,:)
  real(dp)         , allocatable, save :: wra_AFSave  (:,:,:)


  ! Variables for input optical depth
  logical , parameter         :: FlagInputDelOptDep = .false.
  real(DP), allocatable, save :: za_DelOptDep(:,:)



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

  logical, save :: FlagInited = .false.

contains

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

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

    use ni3_module


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

    ! local variables
    !
    character(stdstr) :: PressAxisName
    character(stdstr) :: MolAxisName
    character(stdstr) :: WaveNumAxisName
    character(stdstr) :: AbsCoefVarName
    character(stdstr) :: AtmPFRatioVarName
    character(stdstr) :: StrSpeVarName
    character(stdstr) :: StrPFRatioVarName
    character(stdstr) :: Mode
    integer           :: NcID

    integer              :: TmpNum
    real(DP)             :: DelWaveNum
    integer, allocatable :: m_MolNumTmp(:)

    integer  :: l
    integer  :: m


    ! Variables for input optical depth
    integer           :: NTmp
    character(extstr) :: CTmp



    NPressSave   = -1
    NMolSave     = -1
    NWaveNumSave = -1


    PressAxisName = "Press"
    MolAxisName   = "MolNum"

    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 = 'WaveNum'
    AbsCoefVarName  = 'AbsCoef'
    AtmPFRatioVarName = 'AtmPFRatio'
    StrPFRatioVarName = 'StrPFRatio'
    !
    call ni3_inq_dimlen( NcID, WaveNumAxisName, NWaveNumSave )
    allocate( w_WaveNumSave( 1:NWaveNumSave ) )
    call ni3_get_var( NcID, WaveNumAxisName, w_WaveNumSave )

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

    call ni3_close( NcID )


    allocate( w_RayScatCoefNonRadActSave(NWaveNumSave) )
    allocate( wm_RayScatCoefSave(NWaveNumSave,NMolSave) )
    !
    if ( len_trim( RayScatCoefNcFN ) >= 1 ) then
      Mode = "read"
      call ni3_open( RayScatCoefNcFN, Mode, NcID )
      !
      call ni3_inq_dimlen( NcID, WaveNumAxisName, TmpNum )
      if ( TmpNum /= NWaveNumSave ) then
        write( 6, * ) 'Array size for Rayleigh scattering coefficient is inappropriate:'
        write( 6, * ) '  ', trim( RayScatCoefNcFN ), ': ', TmpNum
        write( 6, * ) '  ', trim( OptPropNcFN     ), ': ', NWaveNumSave
        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 )
      !
      call ni3_get_var( NcID, "RayScatCoefNonRadAct", w_RayScatCoefNonRadActSave )
      call ni3_get_var( NcID, "RayScatCoef"         , wm_RayScatCoefSave )
      !
      if ( ni3_chk_var( NcID, "PtclNum" ) ) then
        call ni3_inq_dimlen( NcID, "PtclNum"      , NPtclSave          )
        call ni3_inq_dimlen( NcID, "PtclRadiusNum", PtclNRadiusMaxSave )
        allocate( a_PtclNameSave   (NPtclSave) )
        allocate( a_PtclNRadiusSave(NPtclSave) )
        allocate( ra_PtclRadiusSave(PtclNRadiusMaxSave,NPtclSave) )
        allocate( wra_QExtSave(NWaveNumSave,PtclNRadiusMaxSave,NPtclSave) )
        allocate( wra_SSASave (NWaveNumSave,PtclNRadiusMaxSave,NPtclSave) )
        allocate( wra_AFSave  (NWaveNumSave,PtclNRadiusMaxSave,NPtclSave) )
        call ni3_get_var( NcID, "PtclName"   , a_PtclNameSave    )
        call ni3_get_var( NcID, "PtclNRadius", a_PtclNRadiusSave )
        call ni3_get_var( NcID, "PtclRadius" , ra_PtclRadiusSave )
        call ni3_get_var( NcID, "QExt"           , wra_QExtSave      )
        call ni3_get_var( NcID, "SSA"            , wra_SSASave       )
        call ni3_get_var( NcID, "AF"             , wra_AFSave        )
      else
        NPtclSave          = 0
        PtclNRadiusMaxSave = 0
        allocate( a_PtclNameSave   (NPtclSave) )
        allocate( a_PtclNRadiusSave(NPtclSave) )
        allocate( ra_PtclRadiusSave(PtclNRadiusMaxSave,NPtclSave) )
        allocate( wra_QExtSave(NWaveNumSave,PtclNRadiusMaxSave,NPtclSave) )
        allocate( wra_SSASave (NWaveNumSave,PtclNRadiusMaxSave,NPtclSave) )
        allocate( wra_AFSave  (NWaveNumSave,PtclNRadiusMaxSave,NPtclSave) )
      end if
      !
      call ni3_close( NcID )
    else
      w_RayScatCoefNonRadActSave = 0.0d0
      wm_RayScatCoefSave         = 0.0d0

      NPtclSave          = 0
      PtclNRadiusMaxSave = 0
      allocate( a_PtclNameSave   (NPtclSave) )
      allocate( a_PtclNRadiusSave(NPtclSave) )
      allocate( ra_PtclRadiusSave(PtclNRadiusMaxSave,NPtclSave) )
      allocate( wra_QExtSave(NWaveNumSave,PtclNRadiusMaxSave,NPtclSave) )
      allocate( wra_SSASave (NWaveNumSave,PtclNRadiusMaxSave,NPtclSave) )
      allocate( wra_AFSave  (NWaveNumSave,PtclNRadiusMaxSave,NPtclSave) )
    end if


    allocate( w_StrSpeSave( 1:NWaveNumSave ) )
    !
    StrSpeVarName   = 'InFlux'
    !
    Mode = "read"
    call ni3_open( StrSpeNcFN, Mode, NcID )
    call ni3_inq_dimlen( NcID, WaveNumAxisName, TmpNum )
    if ( TmpNum /= NWaveNumSave ) then
      write( 6, * ) 'Array size for stellar spectrum is inappropriate: '
      write( 6, * ) '  ', trim( StrSpeNcFN  ), ': ', TmpNum
      write( 6, * ) '  ', trim( OptPropNcFN ), ': ', NWaveNumSave
      stop
    end if
    call ni3_get_var( NcID, StrSpeVarName, w_StrSpeSave )
    call ni3_close( NcID )


    DelWaveNum = w_WaveNumSave(2) - w_WaveNumSave(1)
    StrSpeIntedSave = 0.0d0
    do l = 1, NWaveNumSave
      StrSpeIntedSave = StrSpeIntedSave + w_StrSpeSave(l) * DelWaveNum
    end do


    ! Lines below are for input optical depth
    if ( FlagInputDelOptDep ) then
      CTmp = ""
      Mode = "read"
      call ni3_open( CTmp, Mode, NcID )
      call ni3_inq_dimlen( NcID, "level", NTmp )
      if ( NTmp /= NPressSave-1 ) then
        write( 6, * ) 'Dimension length for level in ', trim( CTmp ), ' in inappropriate: ', NTmp, NPressSave-1
        stop
      end if
      call ni3_inq_dimlen( NcID, "wn", NTmp )
      if ( NTmp /= NWaveNumSave ) then
        write( 6, * ) 'Dimension length for wn in ', trim( CTmp ), ' in inappropriate: ', NTmp, NWaveNumSave
        stop
      end if
      allocate( za_DelOptDep(NPressSave-1,NWaveNumSave) )
      call ni3_get_var( NcID, "value", za_DelOptDep )
      call ni3_close( NcID )
    end if




    FlagInited = .true.

  end subroutine OptPropLBLInit

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

  subroutine OptPropLBLFinalize


    ! local variables
    !


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


    NPressSave   = -1
    NMolSave     = -1
    NWaveNumSave = -1


    deallocate( p_PressSave  )
    deallocate( m_MolNumSave )

    NWaveNumSave    = 0
    deallocate( w_WaveNumSave )

    deallocate( zma_AbsCoef       )

    deallocate( w_RayScatCoefNonRadActSave )
    deallocate( wm_RayScatCoefSave         )
    deallocate( a_PtclNameSave    )
    deallocate( a_PtclNRadiusSave )
    deallocate( ra_PtclRadiusSave )
    deallocate( wra_QExtSave      )
    deallocate( wra_SSASave       )
    deallocate( wra_AFSave        )

    deallocate( w_StrSpeSave )


    ! Lines below are for input optical depth
    if ( FlagInputDelOptDep ) then
      deallocate( za_DelOptDep )
    end if


  end subroutine OptPropLBLFinalize

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

  function OptPropLBLGetNWN() result( NWN )


    integer :: NWN

    ! local variables
    !


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


    NWN = NWaveNumSave


  end function OptPropLBLGetNWN

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

  function OptPropLBLGetNPress() result( NPress )


    integer :: NPress

    ! local variables
    !


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


    NPress = NPressSave


  end function OptPropLBLGetNPress

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

  function OptPropLBLGetNMol() result( NMol )


    integer :: NMol

    ! local variables
    !


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


    NMol = NMolSave


  end function OptPropLBLGetNMol

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

  function OptPropLBLGetDelWN() result( DelWaveNum )


    real(DP) :: DelWaveNum

    ! local variables
    !


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


    DelWaveNum = w_WaveNumSave( 2 ) - w_WaveNumSave( 1 )


  end function OptPropLBLGetDelWN

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

  function OptPropLBLGetWN( iWaveNum ) result( WaveNum )


    integer, intent(in) :: iWaveNum

    real(DP) :: WaveNum

    ! local variables
    !


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


    WaveNum = w_WaveNumSave( iWaveNum )


  end function OptPropLBLGetWN

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

  subroutine OptPropLBLGetTotalStrFlux( &
    & TotalStrFlux &
    & )


    real(DP), intent(out) :: TotalStrFlux


    ! local variables
    !
    real(DP) :: DelWaveNum

    integer  :: l


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


    DelWaveNum = w_WaveNumSave(2) - w_WaveNumSave(1)

    TotalStrFlux = 0.0d0
    do l = 1, NWaveNumSave
      TotalStrFlux = TotalStrFlux &
        & + w_StrSpeSave( l ) * DelWaveNum
    end do


  end subroutine OptPropLBLGetTotalStrFlux

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

  subroutine OptPropLBLGetAbsCoefProf( &
    & 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 :: 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)

    pm_AbsCoef = pm_AbsCoef * pm_VMR


  end subroutine OptPropLBLGetAbsCoefProf

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

  subroutine OptPropLBLGetInputedOptDepProf( &
    & iWaveNum, kmax, &
    & p_DelOptDep &
    & )


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


    ! local variables
    !


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


    ! Lines below are for input optical depth
    if ( FlagInputDelOptDep ) then

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

      p_DelOptDep = za_DelOptDep(:,iWaveNum)

    else

      write( 6, * ) 'FlagInputDelOptDep is ', FlagInputDelOptDep
      stop

    end if


  end subroutine OptPropLBLGetInputedOptDepProf

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

  subroutine OptPropLBLGetRayScatCoef( &
    & 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 :: 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


    RayScatCoefNonRadAct  = w_RayScatCoefNonRadActSave(iWaveNum)
    do m = 1, NMol
      m_RayScatCoef(m) = wm_RayScatCoefSave(iWaveNum,m)
    end do


  end subroutine OptPropLBLGetRayScatCoef

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

  subroutine OptPropLBLGetPtclParam( &
    & iWaveNum, &
    & kmax, NPtcl, a_PtclName, za_PtclRadius, &
    & za_QExt, za_SSA, za_AF &
    & )


    integer     , intent(in ) :: iWaveNum
    integer     , intent(in ) :: kmax
    integer     , intent(in ) :: NPtcl
    character(*), intent(in ) :: a_PtclName(NPtcl)
    real(DP)    , intent(in ) :: za_PtclRadius(kmax,NPtcl)
    real(DP)    , intent(out) :: za_QExt(kmax,NPtcl)
    real(DP)    , intent(out) :: za_SSA (kmax,NPtcl)
    real(DP)    , intent(out) :: za_AF  (kmax,NPtcl)


    ! local variables
    !
    integer  :: NRadius
    real(DP) :: Radius
    real(DP) :: a_X(2)
    real(DP) :: a_Val(2)

    integer  :: iPtcl
    integer  :: i
    integer  :: k
    integer  :: l
    integer  :: n


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

    i = iWaveNum

    do n = 1, NPtcl

      do iPtcl = 1, NPtclSave
        if ( trim( a_PtclName(n) ) == trim( a_PtclNameSave(iPtcl) ) ) exit
      end do
      if ( iPtcl > NPtclSave ) then
        write( 6, * ) trim( a_PtclName(n) ) // ' is not found.'
        stop
      end if

      NRadius = a_PtclNRadiusSave(iPtcl)

      if ( NRadius == 1 ) then
        l = 1
        do k = 1, kmax
          za_QExt(k,n) = wra_QExtSave(i,l,iPtcl)
          za_SSA (k,n) = wra_SSASave (i,l,iPtcl)
          za_AF  (k,n) = wra_AFSave  (i,l,iPtcl)
        end do
      else
        do k = 1, kmax
          Radius = za_PtclRadius(k,n)
          if ( Radius <= ra_PtclRadiusSave(1,n) ) then
            l = 1
            za_QExt(k,n) = wra_QExtSave(i,l,iPtcl)
            za_SSA (k,n) = wra_SSASave (i,l,iPtcl)
            za_AF  (k,n) = wra_AFSave  (i,l,iPtcl)
          else if ( ra_PtclRadiusSave(NRadius,n) <= Radius ) then
            l = NRadius
            za_QExt(k,n) = wra_QExtSave(i,l,iPtcl)
            za_SSA (k,n) = wra_SSASave (i,l,iPtcl)
            za_AF  (k,n) = wra_AFSave  (i,l,iPtcl)
          else
            do l = 1+1, NRadius
              if ( Radius <= ra_PtclRadiusSave(l,n) ) exit
            end do
            a_X(1)   = ra_PtclRadiusSave(l-1,iPtcl)
            a_X(2)   = ra_PtclRadiusSave(l  ,iPtcl)
            !
            a_Val(1) = wra_QExtSave(i,l-1,iPtcl)
            a_Val(2) = wra_QExtSave(i,l  ,iPtcl)
            za_QExt(k,n) = &
              &   ( a_Val(2) - a_Val(1) ) &
              & / ( a_X  (2) - a_X  (1) ) &
              & * ( Radius   - a_X  (1) ) &
              & + a_Val(1)
            !
            a_Val(1) = wra_SSASave (i,l-1,iPtcl)
            a_Val(2) = wra_SSASave (i,l  ,iPtcl)
            za_SSA (k,n) = &
              &   ( a_Val(2) - a_Val(1) ) &
              & / ( a_X  (2) - a_X  (1) ) &
              & * ( Radius   - a_X  (1) ) &
              & + a_Val(1)
            !
            a_Val(1) = wra_AFSave  (i,l-1,iPtcl)
            a_Val(2) = wra_AFSave  (i,l  ,iPtcl)
            za_AF  (k,n) = &
              &   ( a_Val(2) - a_Val(1) ) &
              & / ( a_X  (2) - a_X  (1) ) &
              & * ( Radius   - a_X  (1) ) &
              & + a_Val(1)
          end if
        end do
      end if
    end do


  end subroutine OptPropLBLGetPtclParam

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

  subroutine OptPropLBLGetPFInted( &
    & iWaveNum, &
    & kmax, &
    & r_Temp, SurfTemp, &
    & r_AtmPFInted, SurfPFInted, SurfDPFDTInted &
    & )

    use constants0, only : PI

    use planck_func_wrapper, only : &
      & Integ_PF_GQ_Array1D   , &
      & Integ_PF_GQ_Array0D   , &
      & Integ_DPFDT_GQ_Array0D

    use rad_kd_utils, only : InqPFIntGQNum


    integer , intent(in ) :: iWaveNum
    integer , intent(in ) :: kmax
    real(DP), intent(in ) :: r_Temp   (1:kmax)
    real(DP), intent(in ) :: SurfTemp
    real(DP), intent(out) :: r_AtmPFInted(1:kmax)
    real(DP), intent(out) :: SurfPFInted
    real(DP), intent(out) :: SurfDPFDTInted


    ! local variables
    !
    real(DP) :: WaveNum
    real(DP) :: DelWaveNum
    real(DP) :: WNs
    real(DP) :: WNe
    integer  :: GQNum


    WaveNum    = OptPropLBLGetWN( iWaveNum )
    DelWaveNum = OptPropLBLGetDelWN()
    WNs = WaveNum - DelWaveNum / 2.0d0
    WNe = WaveNum + DelWaveNum / 2.0d0

!!$    GQNum = 1
!!$    GQNum = int( DelWaveNum / 10.0d0 ) + 1
    GQNum = InqPFIntGQNum( DelWaveNum )

    call Integ_PF_GQ_Array1D( &
      & WNs, WNe, GQNum, &
      & 1, kmax, &
      & r_Temp, &
      & r_AtmPFInted &
      & )
    r_AtmPFInted = PI * r_AtmPFInted
    call Integ_PF_GQ_Array0D( &
      & WNs, WNe, GQNum, &
      & SurfTemp, &
      & SurfPFInted &
      & )
    SurfPFInted = PI * SurfPFInted
    call Integ_DPFDT_GQ_Array0D( &
      & WNs, WNe, GQNum,                 & ! (in )
      & SurfTemp,                & ! (in )
      & SurfDPFDTInted                  & ! (out)
      & )
    SurfDPFDTInted = PI * SurfDPFDTInted


  end subroutine OptPropLBLGetPFInted

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

  subroutine OptPropLBLGetStrPFInted( &
    & iWaveNum, StrFluxTOA, &
    & StrPFInted &
    & )

    use constants0, only : PI

    use rad_kd_utils, only : InqPFIntGQNum

    integer , intent(in ) :: iWaveNum
    real(DP), intent(in ) :: StrFluxTOA
    real(DP), intent(out) :: StrPFInted


    ! local variables
    !
    real(DP) :: DelWaveNum


    DelWaveNum = w_WaveNumSave(2) - w_WaveNumSave(1)
    StrPFInted = w_StrSpeSave( iWaveNum ) * DelWaveNum
    StrPFInted = StrPFInted * StrFluxTOA / StrSpeIntedSave


  end subroutine OptPropLBLGetStrPFInted

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

  subroutine OptPropLBLGetPFIntedRatio( &
    & kmax, &
    & r_AtmPFRatio, SurfPFRatio &
    & )

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


    ! local variables
    !


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


    r_AtmPFRatio = 1.0d0
    SurfPFRatio  = 1.0d0


  end subroutine OptPropLBLGetPFIntedRatio

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

  subroutine OptPropLBLGetStrPFIntedRatio( &
    & StrPFRatio &
    & )


    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 = 1.0d0


  end subroutine OptPropLBLGetStrPFIntedRatio

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

end module opt_prop_lbl

