module flux

  use dc_types
  use gtool_history

  implicit none

  private

  public :: FluxCalcFlux
  public :: FluxInit
  public :: FluxFinalize
  public :: FluxOutputInit
  public :: FluxOutputFinalize


  real(8), save :: Grav

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

  type(GT_HISTORY), save :: GTHistOut


  logical, save :: FlagInited       = .false.
  logical, save :: FlagOutputInited = .false.


contains

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

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

    use utils, only : MkRangeString

    use constants0, only : PI

    use planck_func_wrapper

    use rad_rte_two_stream_app_wrapper, only : RadRTETwoStreamAppLWWrapper

    use opt_prop, only : &
      & OptPropLBLGetNPress, &
      & OptPropLBLGetNMol, &
      & OptPropLBLGetNWN, &
      & OptPropLBLGetDelWN, &
      & OptPropLBLGetWN, &
      & OptPropLBLGetAbsCoefProf

    use band, only : &
      & GetBandNum, &
      & r_GetPFInted, &
      & GetSurfPFInted, &
      & GetSurfDPFDTInted


    logical , intent(in ) :: FlagConstPF
    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
    !
    integer :: GQNum = 5

    real(DP) :: z_MMRH2OVap (1:kmax)
    real(DP) :: z_MMRCO2    (1:kmax)
    real(DP) :: z_MMRO3     (1:kmax)

    integer :: NPress
    integer :: NMol
    integer :: NWaveNum
    integer :: iWaveNum
    integer :: iWaveNumStart
    integer :: iWaveNumEnd
    real(8) :: WaveNum
    real(8) :: DelWaveNum

    real(8) :: WNs
    real(8) :: WNe

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

    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)

    character(128) :: String

    integer :: k
    integer :: kp
    integer :: kf
    integer :: m


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



    NPress     = OptPropLBLGetNPress()
    NMol       = OptPropLBLGetNMol()
    NWaveNum   = OptPropLBLGetNWN()
    DelWaveNum = OptPropLBLGetDelWN()
    ! 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:NPress-1, 0:NMol-1) )


    if ( ( WNIntegStart < 0.0d0 ) .or. ( WNIntegEnd < 0.0d0 ) ) then
      iWaveNumStart = 1
      iWaveNumEnd   = NWaveNum
    else
      loop_search_wnstart: do iWaveNum = 1, NWaveNum
        if ( WNIntegStart < OptPropLBLGetWN( iWaveNum ) ) exit
      end do loop_search_wnstart
      iWaveNumStart = iWaveNum
      iWaveNumStart = min( iWaveNumStart, NWaveNum )
      loop_search_wnend: do iWaveNum = iWaveNumStart, NWaveNum
        if ( WNIntegEnd < OptPropLBLGetWN( iWaveNum ) ) exit
      end do loop_search_wnend
      iWaveNumEnd = iWaveNum - 1
    end if

    r_RadUwFlux = 0.0d0
    r_RadDwFlux = 0.0d0

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

      z_SSA = 0.0d0
      z_AF  = 0.0d0
      !
      call OptPropLBLGetAbsCoefProf( &
        & iWaveNum, NPress, 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) * r_Press(k-1) &
            &     + rm_AbsCoef(k  ,m) * r_MMR(k  ) * r_Press(k  ) ) / 2.0d0 &
            & / Grav
        end do
        r_OptDep(k-1) = r_OptDep(k) + DelOptDep
      end do
      !
      SurfAlbedo = 0.0d0
      !

      WaveNum = OptPropLBLGetWN( iWaveNum )
      if ( FlagConstPF ) then
        r_PFInted = r_GetPFInted( GetBandNum( WaveNum ) )
        SurfPFInted = GetSurfPFInted( GetBandNum( WaveNum ) )
        SurfDPFDTInted = GetSurfDPFDTInted( GetBandNum( WaveNum ) )
      else
        WNs = WaveNum - DelWaveNum / 2.0d0
        WNe = WaveNum + DelWaveNum / 2.0d0
        GQNum = 5
        !
        call Integ_PF_GQ_Array1D( &
          & WNs, WNe, GQNum, &
          & 0, kmax, &
          & r_Temp, &
          & r_PFInted &
          & )
        r_PFInted = PI * r_PFInted
        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 if


!!$      write( 6, * ) r_OptDep(0), zm_AbsCoef(1,1+1), z_DelH2OVapMass(1)
      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 )


  end subroutine FluxCalcFlux

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

  subroutine FluxInit( &
    & ArgGrav, &
    & OptPropNcFN &
    & )

    use opt_prop, only : OptPropInit

    use rad_rte_two_stream_app_wrapper, only : RadRTETwoStreamAppLWWrapperInit

    real(8)     , intent(in) :: ArgGrav
    character(*), intent(in) :: OptPropNcFN


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


    if ( FlagInited ) return


    Grav = ArgGrav


    call OptPropInit( &
      & OptPropNcFN &
      & )

    call RadRTETwoStreamAppLWWrapperInit


    FlagInited = .true.

  end subroutine FluxInit

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

  subroutine FluxFinalize

    use opt_prop, only : OptPropFinalize


    ! local variables
    !


    call OptPropFinalize

    if ( FlagOutputInited ) then
      call FluxOutputFinalize
    end if



    FlagInited = .false.

  end subroutine FluxFinalize

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

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

    use 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 FluxOutputInit

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

  subroutine FluxOutputFinalize


    ! 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 FluxOutputFinalize

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

end module flux
