module particle_ec_calc

  use vtype_module

  implicit none

  private

  public :: ParticleECCalc

  logical, parameter :: FlagChkRangeAndStop = .true.
!!$  logical, parameter :: FlagChkRangeAndStop = .false.

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

contains

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

  subroutine ParticleECCalc( &
    & FlagReadQExtSSAAFFromFile, &
    & DelWaveNum, WaveNumS, WaveNumE, &
    & NMolNum, m_MolNum, &
    & NParticle, a_ParticleName, a_ParticleDataFN, &
    & a_ParticleNRadius, a_ParticleRadiusS, a_ParticleRadiusE, &
    & OutputFileName &
    & )

    use rayscat, only : CalcRayScatCoef

    use refindex, only : &
      & RefIndexInit, &
      & RefIndexCalc, &
      & RefIndexFinalize

    use mie_theory, only : CALLBH3

    use particle_ec_io_module, only : &
      & particle_ec_io_output_init, &
      & particle_ec_io_output_end, &
      & particle_ec_io_output


    logical     , intent(in ) :: FlagReadQExtSSAAFFromFile
    real(dp)    , intent(in ) :: DelWaveNum
    real(dp)    , intent(in ) :: WaveNumS
    real(dp)    , intent(in ) :: WaveNumE
    integer     , intent(in ) :: NMolNum
    integer     , intent(in ) :: m_MolNum(NMolNum)
    integer     , intent(in ) :: NParticle
    character(*), intent(in ) :: a_ParticleName   (NParticle)
    character(*), intent(in ) :: a_ParticleDataFN (NParticle)
    integer     , intent(in ) :: a_ParticleNRadius(NParticle)
    real(dp)    , intent(in ) :: a_ParticleRadiusS(NParticle)
    real(dp)    , intent(in ) :: a_ParticleRadiusE(NParticle)
    character(*), intent(in ) :: OutputFileName

!!$    allocate( z_Press ( NLev    ) )
!!$    allocate( m_MolNum( NMolNum ) )
!!$    allocate( z_Temp  ( NLev    ) )
!!$    allocate( zm_VMR  ( NLev, NMolNum ) )


    ! Local variables
    !
    real(dp)    , parameter   :: MissingValue = -1.0d100

    integer(i4b)              :: NWaveNum
    real(dp)    , allocatable :: w_WaveNum   (:)
    real(dp)    , allocatable :: w_Wavelength(:)

    real(dp)                  :: DelWaveNumMod


    real(dp)    , allocatable :: w_RayScatCoefNonRadAct(:)
    real(dp)    , allocatable :: wm_RayScatCoef(:,:)
    real(dp)    , allocatable :: wra_QExt (:,:,:)
    real(dp)    , allocatable :: wra_SSA  (:,:,:)
    real(dp)    , allocatable :: wra_AF   (:,:,:)

    integer     , parameter   :: NRadiusDivide = 1
    integer                   :: iRadius
    real(DP)                  :: RadiusCenter
!!$    real(DP)                  :: RadiusWidth
    real(DP)                  :: QExt
    real(DP)                  :: SSA
    real(DP)                  :: SSCoA
    real(DP)                  :: AF

!!$    real(DP)                  :: a_ParticleRadiusStd(NParticle)
!!$    real(DP)                  :: Weight

!!$    real(dp)    , parameter   :: Rmed = 1.00045D0
    real(dp)    , parameter   :: Rmed = 1.0d0

    real(DP)                  :: RIReal
    real(DP)                  :: RIImag

    real(DP)                  :: Radius
    real(DP)                  :: RadiusS
    real(DP)                  :: RadiusE

    integer                   :: NRadius
    integer                   :: NRadiusMax
    real(DP)    , allocatable :: ra_ParticleRadius(:,:)

    integer                   :: iWaveNum
    integer                   :: l
    integer                   :: m
    integer                   :: n


    NWaveNum   = int( ( WaveNumE - WaveNumS ) / DelWaveNum ) + 1
    DelWaveNumMod = ( WaveNumE - WaveNumS ) / dble( NWaveNum )

    allocate( w_WaveNum   ( NWaveNum ) )
    allocate( w_Wavelength( NWaveNum ) )

    allocate( w_RayScatCoefNonRadAct( NWaveNum ) )
    allocate( wm_RayScatCoef( NWaveNum, NMolNum ) )


    NRadiusMax = 0
    do n = 1, NParticle
      NRadius = a_ParticleNRadius(n)
      NRadiusMax = max( NRadiusMax, NRadius )
    end do
    allocate( ra_ParticleRadius( NRadiusMax, NParticle ) )
    do n = 1, NParticle
      NRadius = a_ParticleNRadius(n)
      RadiusS = a_ParticleRadiusS(n)
      RadiusE = a_ParticleRadiusE(n)
      do l = 1, NRadius
        if ( NRadius == 1 ) then
          Radius = RadiusS
        else
          Radius = &
            &   ( RadiusE - RadiusS ) &
            & / ( NRadius - 1       ) &
            & * ( l       - 1       ) &
            & + RadiusS
        end if
        ra_ParticleRadius(l,n) = Radius
      end do
      do l = NRadius+1, NRadiusMax
        ra_ParticleRadius(l,n) = MissingValue
      end do
    end do

    allocate( wra_QExt  ( NWaveNum, NRadiusMax, NParticle ) )
    allocate( wra_SSA   ( NWaveNum, NRadiusMax, NParticle ) )
    allocate( wra_AF    ( NWaveNum, NRadiusMax, NParticle ) )


    call particle_ec_io_output_init( &
      & OutputFileName, NWaveNum, &
      & NMolNum, m_MolNum, &
      & NParticle, a_ParticleName, a_ParticleNRadius, &
      & NRadiusMax, ra_ParticleRadius, &
      & MissingValue, DelWaveNum &
      & )


    !
    ! loop for wavenumber
    !

    do iWaveNum = 1, NWaveNum
      w_WaveNum( iWaveNum ) = &
        & WaveNumS + DelWaveNumMod * dble( iWaveNum - 1 )
    end do
    w_Wavelength = 1.0d0 / w_WaveNum

    m = 0
    do iWaveNum = 1, NWaveNum
      call CalcRayScatCoef( &
        & m, w_Wavelength(iWaveNum), &
        & w_RayScatCoefNonRadAct(iWaveNum) &
        & )
    end do
    do m = 1, NMolNum
      do iWaveNum = 1, NWaveNum
        call CalcRayScatCoef( &
          & m_MolNum(m), w_Wavelength(iWaveNum), &
          & wm_RayScatCoef(iWaveNum,m) &
          & )
      end do
    end do

    wra_QExt = MissingValue
    wra_SSA  = MissingValue
    wra_AF   = MissingValue
    !
    if ( FlagReadQExtSSAAFFromFile ) then

      do n = 1, NParticle
        do l = 1, a_ParticleNRadius(n)
          if ( l /= 1 ) then
            write( 6, * ) 'ParticleNRadius has to be 1, when FlagReadOptPropFromFile is .true..'
            stop
          end if
        end do
      end do

      do n = 1, NParticle
        do l = 1, a_ParticleNRadius(n)
          call ParticleECCalcSetOptProp( &
            & a_ParticleDataFN(n), &
            & n, &
            & NWaveNum, w_Wavelength, &
            & wra_QExt(:,l,n), wra_SSA(:,l,n), wra_AF(:,l,n) &
            & )
        end do
      end do

    else

      do n = 1, NParticle
        call RefIndexInit( &
          & a_ParticleDataFN(n), &
          & FlagChkRangeAndStop &
          & )
        do l = 1, a_ParticleNRadius(n)

          RadiusCenter = ra_ParticleRadius(l,n)

!!$          RadiusWidth = 2.0d0 * a_ParticleRadiusStd(n)

          do iWaveNum = 1, NWaveNum
            call RefIndexCalc( &
              & w_Wavelength(iWaveNum),                  & ! (in)
              & RIReal, RIImag       & ! (out)
              & )

            wra_QExt(iWaveNum,l,n) = 0.0d0
            wra_SSA (iWaveNum,l,n) = 0.0d0
            wra_AF  (iWaveNum,l,n) = 0.0d0
            do iRadius = 1, NRadiusDivide
!!$              if ( NRadiusDivide == 1 ) then
                Radius = RadiusCenter
!!$              else
!!$                Radius = &
!!$                  &   RadiusCenter - RadiusWidth &
!!$                  & + 2.0d0 * RadiusWidth / dble( NRadiusDivide-1 ) * ( iRadius-1 )
!!$              end if
              call CALLBH3(                &
                & Rmed, RIreal, RIimag, Radius, w_Wavelength(iWaveNum), &
                & QExt, SSA, SSCoA, AF &
                & )
!!$              Weight = 1.0d0 / Radius * exp( - ( log( Radius / RadiusMode )**2 / 2.0d0  )
              wra_QExt(iWaveNum,l,n) = wra_QExt(iWaveNum,l,n) + QExt
              wra_SSA (iWaveNum,l,n) = wra_SSA (iWaveNum,l,n) + SSA
              wra_AF  (iWaveNum,l,n) = wra_AF  (iWaveNum,l,n) + AF
            end do
            wra_QExt(iWaveNum,l,n) = wra_QExt(iWaveNum,l,n) / NRadiusDivide
            wra_SSA (iWaveNum,l,n) = wra_SSA (iWaveNum,l,n) / NRadiusDivide
            wra_AF  (iWaveNum,l,n) = wra_AF  (iWaveNum,l,n) / NRadiusDivide
!!$          write( 6, * ) w_Wavelength(iWaveNum), RIReal, RIImag, wra_QExt(iWaveNum,l,n)
          end do
        end do
        call RefIndexFinalize
      end do

    end if


    call particle_ec_io_output( &
      & NWaveNum, w_WaveNum, w_RayScatCoefNonRadAct, NMolNum, wm_RayScatCoef, &
      & NRadiusMax, NParticle, a_ParticleName, wra_QExt, wra_SSA, wra_AF &
      & )

    call particle_ec_io_output_end


    deallocate( ra_ParticleRadius )

    deallocate( wm_RayScatCoef )
    deallocate( wra_QExt  )
    deallocate( wra_SSA   )
    deallocate( wra_AF    )


  end subroutine ParticleECCalc

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

  subroutine ParticleECCalcSetOptProp( &
    & InNcFN, &
    & iParticle, &
    & NWavelen, w_Wavelen, &
    & w_QExt, w_SSA, w_AF &
    & )

    use ni3_module

    character(*), intent(in ) :: InNcFN
    integer     , intent(in ) :: iParticle
    integer     , intent(in ) :: NWavelen
    real(dp)    , intent(in ) :: w_Wavelen( NWavelen )
    real(dp)    , intent(out) :: w_QExt   ( NWavelen )
    real(dp)    , intent(out) :: w_SSA    ( NWavelen )
    real(dp)    , intent(out) :: w_AF     ( NWavelen )


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

    integer               :: NPtcl
    integer               :: NData
    real(dp), allocatable :: w_InWavelen(:)
    real(dp), allocatable :: wa_InQExt  (:,:)
    real(dp), allocatable :: wa_InSSA   (:,:)
    real(dp), allocatable :: wa_InAF    (:,:)

    real(dp), allocatable :: w_Tmp       (:)

    integer :: l
    integer :: m
    integer :: n


    Mode = "read"
    call ni3_open( InNcFN, Mode, NcID )

    Name = "wl"
    call ni3_inq_dimlen( NcID, Name, NData )

    Name = "PtclNum"
    call ni3_inq_dimlen( NcID, Name, NPtcl )

    ! check
    if ( ( iParticle < 1 ) .or. ( NPtcl < iParticle ) ) then
      write( 6, * ) 'Number of particles in ', trim( InNcFN ), ' is inappropriate: ', iParticle, NPtcl
      stop
    end if

    allocate( w_InWaveLen( NData )        )
    allocate( wa_InQExt  ( NData, NPtcl ) )
    allocate( wa_InSSA   ( NData, NPtcl ) )
    allocate( wa_InAF    ( NData, NPtcl ) )

    Name = "wl"
    call ni3_get_var( NcID, Name, w_InWaveLen )
    Name = "QExt"
    call ni3_get_var( NcID, Name, wa_InQExt   )
    Name = "SSA"
    call ni3_get_var( NcID, Name, wa_InSSA    )
    Name = "AF"
    call ni3_get_var( NcID, Name, wa_InAF     )

    call ni3_close( NcID )


    if ( w_InWavelen(1) > w_InWavelen(2) ) then
      allocate( w_Tmp( NData ) )
      !
      w_Tmp = w_InWavelen(:)
      do l = 1, NData
        w_InWavelen(l) = w_Tmp(NData-l+1)
      end do
      !
      do n = 1, NPtcl
        w_Tmp = wa_InQExt(:,n)
        do l = 1, NData
          wa_InQExt(l,n) = w_Tmp(NData-l+1)
        end do
        w_Tmp = wa_InSSA(:,n)
        do l = 1, NData
          wa_InSSA(l,n) = w_Tmp(NData-l+1)
        end do
        w_Tmp = wa_InAF(:,n)
        do l = 1, NData
          wa_InAF(l,n) = w_Tmp(NData-l+1)
        end do
      end do
      !
      deallocate( w_Tmp )
    end if

    do l = 1, NWavelen
      if ( w_Wavelen(l) <= w_InWavelen(1) ) then
        m = 1
        w_QExt(l) = wa_InQExt(m,iParticle)
        w_SSA (l) = wa_InSSA (m,iParticle)
        w_AF  (l) = wa_InAF  (m,iParticle)
      else if ( w_InWavelen(NData) <= w_Wavelen(l) ) then
        m = NData
        w_QExt(l) = wa_InQExt(m,iParticle)
        w_SSA (l) = wa_InSSA (m,iParticle)
        w_AF  (l) = wa_InAF  (m,iParticle)
      else
        do m = 1+1, NData
          if ( w_Wavelen(l) <= w_InWavelen(m) ) exit
        end do
        w_QExt(l) = &
          &   ( wa_InQExt  (m,iParticle) - wa_InQExt  (m-1,iParticle) ) &
          & / ( w_InWavelen(m)           - w_InWavelen(m-1)           ) &
          & * ( w_Wavelen  (l)           - w_InWavelen(m-1)           ) &
          & + wa_InQExt(m-1,iParticle)
        w_SSA (l) = &
          &   ( wa_InSSA   (m,iParticle) - wa_InSSA   (m-1,iParticle) ) &
          & / ( w_InWavelen(m)           - w_InWavelen(m-1)           ) &
          & * ( w_Wavelen  (l)           - w_InWavelen(m-1)           ) &
          & + wa_InSSA (m-1,iParticle)
        w_AF  (l) = &
          &   ( wa_InAF    (m,iParticle) - wa_InAF    (m-1,iParticle) ) &
          & / ( w_InWavelen(m)           - w_InWavelen(m-1)           ) &
          & * ( w_Wavelen  (l)           - w_InWavelen(m-1)           ) &
          & + wa_InAF  (m-1,iParticle)
      end if
    end do


    deallocate( w_InWavelen  )
    deallocate( wa_InQExt    )
    deallocate( wa_InSSA     )
    deallocate( wa_InAF      )


  end subroutine ParticleECCalcSetOptProp

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

end module particle_ec_calc
