module rad_flux

  use dc_types
  use gtool_history

  implicit none

  private

  public :: RadFluxCalcFlux
  public :: RadFluxInit
  public :: RadFluxFinalize
  public :: RadFluxOutputInit
  public :: RadFluxOutputFinalize


  integer, save      :: IDMethod
  integer, parameter :: IDMethodKD         = 0
  integer, parameter :: IDMethodLBL        = 1
  integer, parameter :: IDMethodLBLConstPF = 2

  real(8), save :: Grav

  integer              , save :: NBand
  real(DP), allocatable, save :: aa_BandWNBnds      (:,:)
  integer , allocatable, save :: aa_BandBinIndexBnds(:,:)


  logical         , save :: FlagOutputInited = .false.
  type(GT_HISTORY), save :: GTHistOut


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

  logical, save :: FlagInited       = .false.


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

contains

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

  subroutine RadFluxCalcFlux( &
    & kmax, &
    & r_Press, &
    & SurfTemp, r_Temp, r_MMRH2OVap, r_MMRCO2, r_MMRO3, &
    & WNIntegStart, WNIntegEnd, &
    & r_RadUwFlux, r_RadDwFlux &
    & )

    use constants0, only : PI

    use rad_utils, only : MkRangeString, InqPFIntGQNum

    use rad_rte_two_stream_app_wrapper, only : RadRTETwoStreamAppLWWrapper

    use opt_prop, only : &
      & OptPropGetNPress, &
      & OptPropGetNMol, &
      & OptPropGetNWN, &
      & OptPropGetNBand, &
      & OptPropGetDelWN, &
      & OptPropGetWN, &
      & OptPropGetBandWNBnds, &
      & OptPropGetBandBinIndexBnds, &
      & OptPropGetAbsCoefProf, &
      & OptPropGetPFInted, &
      & OptPropGetPFIntedRatio

    use planck_func_wrapper

    integer , intent(in ) :: kmax
    real(DP), intent(in ) :: r_Press     (0:kmax)
    real(DP), intent(in ) :: SurfTemp
    real(DP), intent(in ) :: r_Temp      (0:kmax)
    real(DP), intent(in ) :: r_MMRH2OVap (0:kmax)
    real(DP), intent(in ) :: r_MMRCO2    (0:kmax)
    real(DP), intent(in ) :: r_MMRO3     (0:kmax)
    real(DP), intent(in ) :: WNIntegStart
    real(DP), intent(in ) :: WNIntegEnd
    real(DP), intent(out) :: r_RadUwFlux (0:kmax)
    real(DP), intent(out) :: r_RadDwFlux (0:kmax)


    ! local variables
    !
    real(DP), allocatable :: a_BandSurfAvePF   (:)
    real(DP), allocatable :: ra_BandAvePF      (:,:)
    real(DP), allocatable :: a_BandSurfAveDPFDT(:)

    integer :: NPress
    integer :: NMol
    integer :: NWaveNum

    integer :: iBand
    integer :: iBandStart
    integer :: iBandEnd

    integer :: iWaveNum
    integer :: iWaveNumStart
    integer :: iWaveNumEnd
    real(8) :: WaveNum
    real(8) :: DelWaveNum

    real(8), allocatable :: rm_AbsCoef(:,:)

    real(8) :: z_DelAtmMass   (1:kmax)
    real(8) :: r_MMR          (0:kmax)
    real(8) :: DelOptDep

    real(DP) :: z_SSA          (1:kmax)
    real(DP) :: z_AF           (1:kmax)
    real(DP) :: r_OptDep       (0:kmax)
    real(DP) :: SurfAlbedo
    real(DP) :: r_PFInted      (0:kmax)
    real(DP) :: SurfPFInted
    real(DP) :: SurfDPFDTInted
    real(DP) :: r_RadUwFluxEach    (0:kmax)
    real(DP) :: r_RadDwFluxEach    (0:kmax)
    real(DP) :: ra_DelRadUwFluxEach(0:kmax,0:1)
    real(DP) :: ra_DelRadDwFluxEach(0:kmax,0:1)

    real(DP) :: r_PFIntedRatio (0:kmax)
    real(DP) :: SurfPFIntedRatio

    real(DP) :: WNs
    real(DP) :: WNe
    integer  :: GQNum

    character(128) :: String

    integer :: k
    integer :: m
    integer :: n


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


    select case ( IDMethod )
    case ( IDMethodLBLConstPF, IDMethodKD )

      allocate( a_BandSurfAvePF   (NBand) )
      allocate( ra_BandAvePF      (0:kmax,NBand) )
      allocate( a_BandSurfAveDPFDT(NBand) )
!!$      allocate( a_SolarFluxTOA(NBand) )

      do n = 1, NBand
        WNs = aa_BandWNBnds(1,n)
        WNe = aa_BandWNBnds(2,n)
        GQNum = InqPFIntGQNum( WNe - WNs )
        call Integ_PF_GQ_Array0D( &
          & WNs, WNe, GQNum, &
          & SurfTemp, &
          & a_BandSurfAvePF(n) &
          & )
        a_BandSurfAvePF(n) = PI * a_BandSurfAvePF(n) / ( WNe - WNs )
        call Integ_PF_GQ_Array1D( &
          & WNs, WNe, GQNum, &
          & 0, kmax, &
          & r_Temp, &
          & ra_BandAvePF(:,n) &
          & )
        ra_BandAvePF(:,n) = PI * ra_BandAvePF(:,n) / ( WNe - WNs )
        call Integ_DPFDT_GQ_Array0D( &
          & WNs, WNe, GQNum,                 & ! (in )
          & SurfTemp,                & ! (in )
          & a_BandSurfAveDPFDT(n) &
          & )
        a_BandSurfAveDPFDT(n) = PI * a_BandSurfAveDPFDT(n) / ( WNe - WNs )
      end do

    end select



    do k = 1, kmax
      z_DelAtmMass(k) = ( r_Press(k-1) - r_Press(k) ) / Grav
    end do


    NPress     = OptPropGetNPress()
    NMol       = OptPropGetNMol()
    NWaveNum   = OptPropGetNWN()
    ! check
    if ( kmax /= NPress-1 ) then
      stop 'Number of layers does not match with number of pressure levels in optical properties.'
    end if

    allocate( rm_AbsCoef(0:kmax, 0:NMol-1) )


    r_RadUwFlux = 0.0d0
    r_RadDwFlux = 0.0d0


    select case ( IDMethod )
    case ( IDMethodKD )
      if ( ( WNIntegStart < 0.0d0 ) .or. ( WNIntegEnd < 0.0d0 ) ) then
        iWaveNumStart = 1
        iWaveNumEnd   = NWaveNum
      else
        do iBand = 1, NBand
          if ( aa_BandWNBnds(1,iBand) >= WNIntegStart ) exit
        end do
        iBandStart = iBand
        iBandStart = min( iBandStart, NBand )
        do iBand = iBandStart, NBand
          if ( aa_BandWNBnds(2,iBand) >= WNIntegEnd ) exit
        end do
        iBandEnd = iBand
        iWaveNumStart = aa_BandBinIndexBnds(1,iBandStart)
        iWaveNumEnd   = aa_BandBinIndexBnds(2,iBandEnd  )
      end if
    case ( IDMethodLBL, IDMethodLBLConstPF )
      if ( ( WNIntegStart < 0.0d0 ) .or. ( WNIntegEnd < 0.0d0 ) ) then
        iWaveNumStart = 1
        iWaveNumEnd   = NWaveNum
      else
        loop_search_wnstart: do iWaveNum = 1, NWaveNum
          if ( OptPropGetWN( iWaveNum ) >= WNIntegStart ) exit
        end do loop_search_wnstart
        iWaveNumStart = iWaveNum
        iWaveNumStart = min( iWaveNumStart, NWaveNum )
        loop_search_wnend: do iWaveNum = iWaveNumStart, NWaveNum
          if ( OptPropGetWN( iWaveNum ) > WNIntegEnd ) exit
        end do loop_search_wnend
        iWaveNumEnd = iWaveNum - 1
      end if
    case default
      write( 6, * ) 'In ', trim( ModuleName )
      write( 6, * ) '  Unexpected case default.'
      stop
    end select


!!$    do iWaveNum = 1, NWaveNum
    do iWaveNum = iWaveNumStart, iWaveNumEnd

      DelWaveNum = OptPropGetDelWN( iWaveNum )

      z_SSA = 0.0d0
      z_AF  = 0.0d0
      !
      call OptPropGetAbsCoefProf( &
        & iWaveNum, kmax+1, NMol, &
        & rm_AbsCoef &
        & )


      r_OptDep(kmax) = 0.0d0
      do k = kmax, 1, -1
        DelOptDep = 0.0d0
        do m = 0, NMol-1
          select case ( m )
          case ( 1 )
            r_MMR = r_MMRH2OVap
          case ( 2 )
            r_MMR = r_MMRCO2
          case ( 3 )
            r_MMR = r_MMRO3
          case default
            r_MMR = 0.0d0
          end select
          DelOptDep = DelOptDep &
            & + (   rm_AbsCoef(k-1,m) * r_MMR(k-1)           &
            &     + rm_AbsCoef(k  ,m) * r_MMR(k  ) ) / 2.0d0 &
            & * z_DelAtmMass(k)
        end do
        r_OptDep(k-1) = r_OptDep(k) + DelOptDep
      end do
      !
      SurfAlbedo = 0.0d0


      select case ( IDMethod )
      case ( IDMethodLBLConstPF, IDMethodKD )

        select case ( IDMethod )
        case ( IDMethodLBLConstPF )
          WaveNum = OptPropGetWN( iWaveNum )
          search_Band : do n = 1, NBand
            if ( WaveNum <= aa_BandWNBnds(2,n) ) exit search_Band
          end do search_Band
        case ( IDMethodKD )
          do n = 1, NBand
            if ( ( aa_BandBinIndexBnds(1,n) <= iWaveNum ) .and. &
              &  ( iWaveNum <= aa_BandBinIndexBnds(2,n) ) ) &
              exit
          end do
        case default
          write( 6, * ) 'In ', trim( ModuleName )
          write( 6, * ) '  Unexpected case default'
          stop
        end select
        if ( n >= NBand + 1 ) then
          write( 6, * ) 'In ', trim( ModuleName )
          write( 6, * ) '  Unable to find appropriate band number: ', n
          stop
        end if

        r_PFInted      = ra_BandAvePF(:,n)
        SurfPFInted    = a_BandSurfAvePF(n)
        SurfDPFDTInted = a_BandSurfAveDPFDT(n)
        !
        call OptPropGetPFIntedRatio( &
          & iWaveNum, kmax+1, &
          & r_PFIntedRatio, SurfPFIntedRatio &
          & )
        !
        r_PFInted      = r_PFInted * r_PFIntedRatio
        SurfPFInted    = SurfPFInted * SurfPFIntedRatio
        SurfDPFDTInted = SurfDPFDTInted * SurfPFIntedRatio

        r_PFInted      = r_PFInted      * DelWaveNum
        SurfPFInted    = SurfPFInted    * DelWaveNum
        SurfDPFDTInted = SurfDPFDTInted * DelWaveNum

      case default
        call OptPropGetPFInted( &
          & iWaveNum, &
          & kmax, &
          & r_Temp, SurfTemp, &
          & r_PFInted, SurfPFInted, SurfDPFDTInted &
          & )
      end select

      call RadRTETwoStreamAppLWWrapper(      &
        & z_SSA, z_AF,                             & ! (in)
        & r_OptDep,                                & ! (in)
        & SurfAlbedo,                              & ! (in)
        & r_PFInted, SurfPFInted, SurfDPFDTInted,  & ! (in)
        & r_RadUwFluxEach, r_RadDwFluxEach,        & ! (out)
        & ra_DelRadUwFluxEach, ra_DelRadDwFluxEach & ! (out)
        & )

      if ( FlagOutputInited ) then
        String = MkRangeString( "w_WaveNum", iWaveNum )
        call HistoryPut('rw_RadFluxUp',r_RadUwFluxEach/DelWaveNum,range=String,history=GTHistOut)
        call HistoryPut('rw_RadFluxDn',r_RadDwFluxEach/DelWaveNum,range=String,history=GTHistOut)
      end if

      r_RadUwFlux = r_RadUwFlux + r_RadUwFluxEach
      r_RadDwFlux = r_RadDwFlux + r_RadDwFluxEach

    end do

    deallocate( rm_AbsCoef )


    select case ( IDMethod )
    case ( IDMethodLBLConstPF, IDMethodKD )
      deallocate( a_BandSurfAvePF    )
      deallocate( ra_BandAvePF       )
      deallocate( a_BandSurfAveDPFDT )
    end select


  end subroutine RadFluxCalcFlux

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

  subroutine RadFluxInit( &
    & ArgGrav, &
    & OptPropNcFN, &
    & Method, &
    & a_ArgBandWNBnds &
    & )

    use opt_prop, only : &
      & OptPropInit, &
      & OptPropGetNBand, &
      & OptPropGetBandWNBnds, &
      & OptPropGetBandBinIndexBnds

    use rad_rte_two_stream_app_wrapper, only : RadRTETwoStreamAppLWWrapperInit

    real(8)     , intent(in )           :: ArgGrav
    character(*), intent(in )           :: OptPropNcFN
    character(*), intent(in ), optional :: Method
    real(DP)    , intent(in ), optional :: a_ArgBandWNBnds(:)


    ! local variables
    !
    integer :: n


    if ( FlagInited ) return


    if ( present( Method ) ) then
      select case ( Method )
      case ( 'KD' )
        IDMethod = IDMethodKD
      case ( 'LBL' )
        IDMethod = IDMethodLBL
      case ( 'LBLConstPF' )
        IDMethod = IDMethodLBLConstPF
      case default
        write( 6, * ) "In module, ", trim( ModuleName )
        write( 6, * ) "  Unexpected method: ", trim( Method )
        stop
      end select
    else
      IDMethod = IDMethodKD
    end if


    ! Check arguments
    !
    select case ( IDMethod )
    case ( IDMethodLBLConstPF )
      if ( present( a_ArgBandWNBnds ) ) then
      else
        write( 6, * ) "In module, ", trim( ModuleName )
        write( 6, * ) 'a_BandWNBnds has to be present when IDMethod is IDMethodLBLConstPF.'
        stop
      end if
    end select


    select case ( IDMethod )
    case ( IDMethodKD )
      call OptPropInit( &
        & 'KD', &
        & OptPropNcFN &
        & )
    case ( IDMethodLBL )
      call OptPropInit( &
        & 'LBL', &
        & OptPropNcFN &
        & )
    case ( IDMethodLBLConstPF )
      call OptPropInit( &
        & 'LBLConstPF', &
        & OptPropNcFN &
        & )
    end select


    call RadRTETwoStreamAppLWWrapperInit


    Grav = ArgGrav


    select case ( IDMethod )
    case ( IDMethodKD )

      NBand = OptPropGetNBand()
      allocate( aa_BandWNBnds      (2,NBand) )
      allocate( aa_BandBinIndexBnds(2,NBand) )
      aa_BandWNBnds       = OptPropGetBandWNBnds()
      aa_BandBinIndexBnds = OptPropGetBandBinIndexBnds()

    case ( IDMethodLBLConstPF )

      NBand = size( a_ArgBandWNBnds ) - 1
      allocate( aa_BandWNBnds(2,NBand) )

      do n = 1, NBand
        aa_BandWNBnds(1,n) = a_ArgBandWNBnds(n  )
        aa_BandWNBnds(2,n) = a_ArgBandWNBnds(n+1)
      end do

    end select



    FlagInited = .true.


  end subroutine RadFluxInit

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

  subroutine RadFluxFinalize

    use opt_prop, only : OptPropFinalize


    ! local variables
    !


    call OptPropFinalize

    if ( FlagOutputInited ) then
      call RadFluxOutputFinalize
    end if


    select case ( IDMethod )
    case ( IDMethodKD )
      deallocate( aa_BandWNBnds       )
      deallocate( aa_BandBinIndexBnds )
    case ( IDMethodLBLConstPF )
      deallocate( aa_BandWNBnds )
    end select


    FlagInited = .false.

  end subroutine RadFluxFinalize

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

  subroutine RadFluxOutputInit( &
    & OptPropNcFN, &
    & FluxOutNcFn, &
    & kmax, r_Press &
    & )

    use rad_utils, only : InqDimLength

    character(*), intent(in) :: OptPropNcFN
    character(*), intent(in) :: FluxOutNcFn
    integer     , intent(in) :: kmax
    real(8)     , intent(in) :: r_Press    (0:kmax)


    ! local variables
    !
    integer :: LBLNWN
    real(DP), allocatable :: w_WaveNum(:)


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


    LBLNWN = InqDimLength( OptPropNcFN, 'w_WaveNum' )
    allocate( w_WaveNum( 1:LBLNWN ) )
    call HistoryGet( OptPropNcFN, 'w_WaveNum', w_WaveNum, quiet = .true. )


    ! 出力ファイルの設定
    ! setting of output file
    !
    call HistoryCreate(                                    &
      & quiet = .true., &
      & file=FluxOutNcFn, &
      & title='Flux', &
      & source='flux.f90', &
      & institution='DCRTM',       &
      & dims=(/'r_Press  ', 'w_WaveNum'/), dimsizes=(/kmax+1,LBLNWN/), &
      & longnames=(/'pressure  ', 'wavenumber'/),       &
      & units=(/'Pa ', 'm-1'/), &
      & history=GTHistOut &
      & )

    call HistoryPut('r_Press  ',r_Press, history=GTHistOut)
    call HistoryPut('w_WaveNum',w_WaveNum, history=GTHistOut)

    call HistoryAddVariable( &
      & varname='rw_RadFluxUp', dims=(/'r_Press  ', 'w_WaveNum'/), &
      & longname='upward flux(planet)', units='W m-2 m', xtype='double', &
      & history=GTHistOut &
      & )

    call HistoryAddVariable( &
      & varname='rw_RadFluxDn', dims=(/'r_Press  ', 'w_WaveNum'/), &
      & longname='downward flux(planet)', units='W m-2 m', xtype='double', &
      & history=GTHistOut &
      & )


    deallocate( w_WaveNum )


    FlagOutputInited = .true.

  end subroutine RadFluxOutputInit

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

  subroutine RadFluxOutputFinalize


    ! local variables
    !


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


    if ( FlagOutputInited ) then
      call HistoryClose( history=GTHistOut, quiet = .true. )
      FlagOutputInited = .false.
    end if


  end subroutine RadFluxOutputFinalize

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

end module rad_flux
