module sortac_module

  use dc_types

  implicit none

  private

  public :: sortac

  character(128), save :: AbsCoefVarName = 'rmw_TAbsorpK'

contains

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

  subroutine SortAC( &
    & FlagSort, &
    & InNcFN, OutNcFN, &
    & FlagChkConstH2OMR, FlagBinning, &
    & FlagKeyMolACAverage, FlagNonKeyMolACAverage, &
    & FlagPFRatioOne, &
    & FlagSortPFRatioAverage, &
    & FlagCorrelation, &
    & FlagTwoMolAbsOverlap, &
    & NBand, a_BandWNBnds, a_BandKeyMolNum, a_ArgBandNBin &
    & )

    use gtool_history                                 ! モジュール指定
    use sort_module, only : sort_heap
    use rad_utils, only : InqDimLength, InqPFIntGQNum

    use planck_func_wrapper

    use constants0, only: PI
    use planetconst, only: Grav

    logical     , intent(in ) :: FlagSort
    character(*), intent(in ) :: InNcFN
    character(*), intent(in ) :: OutNcFN

    logical, intent(in ) :: FlagChkConstH2OMR
    logical, intent(in ) :: FlagBinning
    logical, intent(in ) :: FlagKeyMolACAverage
    logical, intent(in ) :: FlagNonKeyMolACAverage
    logical, intent(in ) :: FlagPFRatioOne
    logical, intent(in ) :: FlagSortPFRatioAverage
    logical, intent(in ) :: FlagCorrelation
    logical, intent(in ) :: FlagTwoMolAbsOverlap
    integer, intent(in ) :: NBand
    real(8), intent(in ) :: a_BandWNBnds(0:NBand)
                                    ! Wavenumber bounds of bands
    integer, intent(in ) :: a_BandKeyMolNum(1:NBand)
    integer, intent(in ) :: a_ArgBandNBin  (1:NBand)
                           ! Number of bins in bands


    ! local variables
    !
    integer :: IndexH2OVap
    integer :: IndexKeyMol

    integer                :: a_BandNBin  (1:NBand)

    integer                :: NPress
    real(8), allocatable   :: r_Press(:)

    integer                :: NTemp
    real(8), allocatable   :: r_Temp(:)

    integer                :: NWaveNum
    real(8), allocatable   :: w_WaveNum(:)

    integer                :: NVMRH2O

    integer                :: NMol
    integer, allocatable   :: m_MolNum(:)             ! Molecular number
    real(8), allocatable   :: m_MolWt(:)              ! Molecular weight

    real(8), allocatable   :: m_MixRatio(:)           ! volume mixing ratio
    real(8), allocatable   :: rm_MixRatio(:,:)        ! volume mixing ratio
    real(8), allocatable   :: rmw_TAbsorpK(:,:,:)     ! Absorption coefficient
    real(8), allocatable   :: wm_TAbsorpK(:,:)
    real(8), allocatable   :: rw_PFRatio(:,:)       ! Planck func. ratio
    real(8), allocatable   :: w_PFRatio(:)
    real(8), allocatable   :: w_AbsorpK(:)
    integer, allocatable   :: w_SortedIndex   (:)
!!$    integer, allocatable   :: w_SortedIndexTmp(:)
    integer, allocatable   :: wr_SortedIndexSave(:,:)

    real(8), allocatable   :: rm_MMR(:,:)


!!$  real(8) :: PressMin
!!$  real(8) :: PressMax
!!$  integer :: NPressPer10
!!$
!!$  real(8) :: TempMin
!!$  real(8) :: TempMax
!!$  real(8) :: TempInc

!!$  real(8) :: Press
!!$  real(8) :: Temp

    integer              :: a_BndNum(2)

    integer, allocatable :: a_BandNum(:)
    ! Number of bands
    real(8), allocatable :: aa_BandWNBnds(:,:)
                           ! Wavenumber bounds of bands
    integer, allocatable :: aa_BandWNIndexBnds(:,:)
                           ! Wavenumber index bounds of bands
    integer, allocatable :: aa_BandBinIndexBnds(:,:)
                           ! Index bounds of bins in bands

    integer              :: NBin
                           ! Number of bins
    integer, allocatable :: a_BinNum    (:)
                           ! Bin number
    real(8), allocatable :: a_BinPointWN(:)
                           ! Bin point wavenumber
    real(8), allocatable :: a_BinWeight (:)
                           ! Bin weight
    integer, allocatable :: aa_BinWNIndexBnds(:,:)
                           ! Wavenumber index bounds of bins
    real(8), allocatable :: rma_BinTAbsorpK(:,:,:)     ! Absorption coefficient
    real(8), allocatable :: ra_BinPFRatio (:,:)        ! Planck func. ratio

    real(8), allocatable :: a_BinWNBnds(:)

!!$    real(8), allocatable :: r_DelMass(:)
!!$    real(8) :: DelAbsMass
    real(8) :: Sum

!!$    integer :: kp
!!$    integer :: kf

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

    real(DP), allocatable :: ra_BandAvePF(:,:)

    real(DP) :: MeanMolWt

    integer :: j
    integer :: k
    integer :: l
    integer :: m
    integer :: n
    integer :: iTemp
    integer :: iVMRH2O

    integer :: ls
    integer :: le

    integer :: ks
    integer :: ke



    a_BndNum(1) = 1
    a_BndNum(2) = 2

    allocate( aa_BandWNBnds     (2,NBand) )
    allocate( aa_BandWNIndexBnds(2,NBand) )

    allocate( a_BandNum(NBand) )
    do n = 1, NBand
      a_BandNum(n) = n
    end do

    do n = 1, NBand
      aa_BandWNBnds(1,n) = a_BandWNBnds(n-1)
      aa_BandWNBnds(2,n) = a_BandWNBnds(n  )
    end do



    NVMRH2O = 1
    NTemp  = 1

    do iVMRH2O = 1, NVMRH2O

      do iTemp = 1, NTemp

!!$    Temp = TempMin + TempInc * ( n - 1 )
!!$
!!$    r_Temp = Temp
!!$
!!$    write( ncfn, '(a,i4.4,a)' ) 'out/Profile_T', int( Temp ), '.nc'

        NWaveNum = InqDimLength( InNcFN, "w_WaveNum" )
        NMol     = InqDimLength( InNcFN, "m_MolNum" )
        NPress   = InqDimLength( InNcFN, "r_Press" )

        allocate( w_WaveNum(NWaveNum) )
        allocate( m_MolNum(NMol) )
        allocate( m_MolWt (NMol) )
        allocate( r_Press(NPress) )
        allocate( r_Temp (NPress) )

        allocate( rm_MixRatio (NPress, NMol) )
        allocate( rmw_TAbsorpK(NPress, NMol, NWaveNum) )
        allocate( rw_PFRatio  (NPress, NWaveNum) )

        call HistoryGet( InNcFN, "r_Press" , r_Press         , quiet = .true. )
        call HistoryGet( InNcFN, "m_MolNum", m_MolNum        , quiet = .true. )
        call HistoryGet( InNcFN, "m_MolWt" , m_MolWt         , quiet = .true. )
        call HistoryGet( InNcFN, "w_WaveNum", w_WaveNum      , quiet = .true. )
        call HistoryGet( InNcFN, "r_Temp"      , r_Temp      , quiet = .true. )
        call HistoryGet( InNcFN, "rm_MixRatio" , rm_MixRatio , quiet = .true. )
        call HistoryGet( InNcFN, "rmw_TAbsorpK", rmw_TAbsorpK, quiet = .true. )


        allocate( rm_MMR(NPress,NMol) )
        MeanMolWt = m_MolWt(1)
        do m = 1, NMol
          rm_MMR(:,m) = rm_MixRatio(:,m) * m_MolWt(m) / MeanMolWt
        end do


        allocate( m_MixRatio  (NMol) )
        k = 1
        do m = 1, NMol
          m_MixRatio(m) = rm_MixRatio(k,m)
        end do
        if ( FlagChkConstH2OMR ) then
          do m = 1, NMol
            do k = 1+1, NPress
              if ( m_MixRatio(m) /= rm_MixRatio(k,m) ) then
                write( 6, * ) 'ERROR: Mixing ratio is not constant over pressure.'
                write( 6, * ) m, k, rm_MixRatio(k,m), m_MixRatio(m)
                stop
              end if
            end do
          end do
        end if


        ! Prepare wavenumber integrated band Planck function
        allocate( ra_BandAvePF(NPress,NBand) )
        do n = 1, NBand
          WNs = a_BandWNBnds(n-1)
          WNe = a_BandWNBnds(n  )
!!$          GQNum = int( ( WNe - WNs ) / 100.0d0 ) + 1
          GQNum = InqPFIntGQNum( WNe - WNs )
          call Integ_PF_GQ_Array1D( &
            & WNs, WNe, GQNum, &
            & 1, NPress, &
            & r_Temp, &
            & ra_BandAvePF(:,n) &
            & )
          ra_BandAvePF(:,n) = PI * ra_BandAvePF(:,n) / ( WNe - WNs )
        end do


        ! Adjustment of absorption coefficient
        ! MEMO:
        !   w * B * exp(-tau) = B * w * exp( -tau ) = B * exp( -tau* )
        !   w * exp(-tau) = exp(-tau*)
        !   -tau* = log(w) - tau
        !   tau* = tau - log(w)
!!$        allocate( r_DelMass(1:kmax) )
!!$        do k = 1, kmax
!!$          if ( k == 1 ) then
!!$            kp = 1
!!$            kf = k+1
!!$          else
!!$            kp = k-1
!!$            kf = k
!!$          end if
!!$          r_DelMass(k) = ( r_Press(kp) - r_Press(kf) ) / Grav
!!$        end do
!!$        do l = 1, NWaveNum
!!$          do m = 1, NMol
!!$            do k = 1, kmax
!!$              rmw_PFRatio(k,m,l) = &
!!$                & PI * integral_planck(w_WaveNum(l)-0.5d0, w_WaveNum(l)+0.5d0, r_Temp(k)) &
!!$                & / GetPFInted( k, w_WaveNum(l) )
!!$              DelAbsMass = r_DelMass(k) * rm_MixRatio(k,m)
!!$              rmw_TAbsorpK(k,m,l) = rmw_TAbsorpK(k,m,l) &
!!$                & - log( PFRatio ) / DelAbsMass
!!$            end do
!!$          end do
!!$        end do
!!$        deallocate( r_DelMass )


        ! Look for bound index
        loop_search_bound_index_1 : do l = 1, NWaveNum
          if ( w_WaveNum(l) >= a_BandWNBnds(0) ) exit loop_search_bound_index_1
        end do loop_search_bound_index_1
        aa_BandWNIndexBnds(1,1) = l
        !
        do n = 1, NBand
          loop_search_bound_index : do l = aa_BandWNIndexBnds(1,n), NWaveNum
            if ( w_WaveNum(l) > a_BandWNBnds(n) ) exit loop_search_bound_index
          end do loop_search_bound_index
          aa_BandWNIndexBnds(2,n) = l-1
          if ( n < NBand ) then
            aa_BandWNIndexBnds(1,n+1) = aa_BandWNIndexBnds(2,n) + 1
          end if
        end do


        ! check band
        do n = 1, NBand
          if ( aa_BandWNIndexBnds(2,n) - aa_BandWNIndexBnds(1,n) + 1 <= 0 ) then
            write( 6, * ) 'ERROR: Band ', n, ' has no absorption coefficient.'
            stop
          end if
        end do

        if ( FlagPFRatioOne ) then

          do l = 1, NWaveNum
            rw_PFRatio(:,l) = 1.0d0
          end do

        else

          do l = 1, aa_BandWNIndexBnds(1,1)-1
            rw_PFRatio(:,l) = 0.0d0
          end do
          do l = aa_BandWNIndexBnds(1,1), aa_BandWNIndexBnds(2,NBand)
            DelWaveNum = w_WaveNum(2) - w_WaveNum(1)
!!$          WNs   = w_WaveNum(l)-1.0d0/2.0d0
!!$          WNe   = w_WaveNum(l)+1.0d0/2.0d0
            WNs   = w_WaveNum(l)-DelWaveNum/2.0d0
            WNe   = w_WaveNum(l)+DelWaveNum/2.0d0
!!$          GQNum = 1
!!$          GQNum = int( DelWaveNum / 10.0d0 ) + 1
            GQNum = InqPFIntGQNum( DelWaveNum )
            call Integ_PF_GQ_Array1D( &
              & WNs, WNe, GQNum, &
              & 1, NPress, &
              & r_Temp, &
              & rw_PFRatio(:,l) &
              & )

            rw_PFRatio(:,l) = PI * rw_PFRatio(:,l)
            rw_PFRatio(:,l) = rw_PFRatio(:,l) / ( WNe - WNs )
            n = GetBandNum( NBand, a_BandWNBnds, w_WaveNum(l) )
            rw_PFRatio(:,l) = rw_PFRatio(:,l) / ra_BandAvePF(:,n)

          end do
          do l = aa_BandWNIndexBnds(2,NBand)+1, NWaveNum
            rw_PFRatio(:,l) = 0.0d0
          end do

        end if



        ! modify absorption coefficient
        !
        if ( FlagTwoMolAbsOverlap ) then
          IndexH2OVap = -1
          do m = 1, NMol
            if ( m_MolNum(m) == 1 ) then
              IndexH2OVap = m
            end if
          end do
          do n = 1, NBand
            if ( a_BandKeyMolNum(n) /= 1 ) then

              IndexKeyMol = -1
              do m = 1, NMol
                if ( m_MolNum(m) == a_BandKeyMolNum(n) ) then
                  IndexKeyMol = m
                end if
              end do
              m = IndexKeyMol
              do l = aa_BandWNIndexBnds(1,n), aa_BandWNIndexBnds(2,n)
                rmw_TAbsorpK(:,m,l) = &
                  &   (   rmw_TAbsorpK(:,IndexH2OVap,l) &
                  &         * rm_MMR(:,IndexH2OVap)     &
                  &     + rmw_TAbsorpK(:,m          ,l) &
                  &         * rm_MMR(:,m          ) )   &
                  & / rm_MMR(:,m)
!!$                rmw_TAbsorpK(:,m,l) = &
!!$                  &   rmw_TAbsorpK(:,IndexH2OVap,l)            &
!!$                  &     * rm_MMR(:,IndexH2OVap) / rm_MMR(:,m)  &
!!$                  & + rmw_TAbsorpK(:,m          ,l)
                  !
                rmw_TAbsorpK(:,IndexH2OVap,l) = 0.0d0
              end do

            end if
          end do
        end if





        ! Index for a key specie
        allocate( wm_TAbsorpK(NWaveNum, NMol) )
        allocate( w_PFRatio  (NWaveNum) )
        allocate( w_AbsorpK( NWaveNum ) )
        allocate( w_SortedIndex   ( NWaveNum ) )
!!$        allocate( w_SortedIndexTmp( NWaveNum ) )
        allocate( wr_SortedIndexSave( NWaveNum, NPress ) )


        if ( FlagSort ) then

          if ( FlagCorrelation ) then
            ks = 1
            ke = NPress
          else
            ks = 1
            ke = 1
          end if

          do k = ks, ke

            ! Work variable
            do m = 1, NMol
              do l = 1, NWaveNum
                wm_TAbsorpK(l,m) = rmw_TAbsorpK(k,m,l)
              end do
            end do

            ! Initialization of index
            !   Index range which is sorted below will be overwritted.
            do l = 1, NWaveNum
              w_SortedIndex( l ) = l
            end do

            ! loop for band
            do n = 1, NBand

              ! sort absorption coefficient in a band
              if ( (aa_BandWNIndexBnds(2,n)-aa_BandWNIndexBnds(1,n)+1) > 0 ) then

                ! pick up a key specie in a band
                !   initialize not to miss an unset variable
                IndexKeyMol = -1
                do m = 1, NMol
                  if ( m_MolNum(m) == a_BandKeyMolNum(n) ) then
                    IndexKeyMol = m
                  end if
                end do

                m = IndexKeyMol

                call sort_heap( &
                  & (aa_BandWNIndexBnds(2,n)-aa_BandWNIndexBnds(1,n)+1), &
                  & wm_TAbsorpK(aa_BandWNIndexBnds(1,n):aa_BandWNIndexBnds(2,n),m), &
                  & w_SortedIndex(aa_BandWNIndexBnds(1,n):aa_BandWNIndexBnds(2,n)) &
                  & )
                w_SortedIndex(aa_BandWNIndexBnds(1,n):aa_BandWNIndexBnds(2,n)) = &
                  &   w_SortedIndex(aa_BandWNIndexBnds(1,n):aa_BandWNIndexBnds(2,n)) &
                  & + aa_BandWNIndexBnds(1,n) - 1


                ! Order is reversed.
!!$                do l = aa_BandWNIndexBnds(1,n), aa_BandWNIndexBnds(2,n)
!!$                  w_SortedIndexTmp(l) = w_SortedIndex(l)
!!$                end do
!!$                do l = aa_BandWNIndexBnds(1,n), aa_BandWNIndexBnds(2,n)
!!$                  w_SortedIndex(l) = w_SortedIndexTmp(aa_BandWNIndexBnds(2,n)-(l-aa_BandWNIndexBnds(1,n)))
!!$                end do


              end if

            end do

            do l = 1, NWaveNum
              wr_SortedIndexSave(l,k) = w_SortedIndex(l)
            end do
          end do

          if ( FlagCorrelation ) then
          else
            do k = 1, NPress
              if ( k == ks ) cycle
              do l = 1, NWaveNum
                wr_SortedIndexSave(l,k) = wr_SortedIndexSave(l,ks)
              end do
            end do
          end if

        else

          do k = 1, NPress
            do l = 1, NWaveNum
              wr_SortedIndexSave(l,k) = l
            end do
          end do

        end if



        if ( FlagSort ) then

          do k = 1, NPress

            ! TESTTESTTESTTESTTEST
!!$            if ( k == 1 ) then
!!$              do n = 1, NBand
!!$                m = 1+1
!!$                do l = aa_BandWNIndexBnds(1,n), aa_BandWNIndexBnds(2,n)
!!$                  write( 80, * ) l, rmw_TAbsorpK(k,m,l), rw_PFRatio(k,l )
!!$                end do
!!$              end do
!!$            end if
            ! TESTTESTTESTTESTTEST



            ! sort ac based on order of key specie
            do m = 1, NMol
              do l = 1, NWaveNum
                wm_TAbsorpK(l,m) = rmw_TAbsorpK(k,m,wr_SortedIndexSave(l,k))
              end do
            end do
            ! sort PFRatio based on order of key specie
            do l = 1, NWaveNum
              w_PFRatio(l) = rw_PFRatio(k,wr_SortedIndexSave(l,k))
            end do

            do m = 1, NMol
              do l = 1, NWaveNum
                rmw_TAbsorpK(k,m,l) = wm_TAbsorpK(l,m)
              end do
            end do
            do l = 1, NWaveNum
              rw_PFRatio(k,l) = w_PFRatio(l)
            end do

            ! TESTTESTTESTTESTTEST
!!$            if ( k == 1 ) then
!!$              do n = 1, NBand
!!$                m = 1+1
!!$                do l = aa_BandWNIndexBnds(1,n), aa_BandWNIndexBnds(2,n)
!!$                  write( 81, * ) l, rmw_TAbsorpK(k,m,l), rw_PFRatio(k,l )
!!$                end do
!!$              end do
!!$            end if
            ! TESTTESTTESTTESTTEST

          end do

        end if


        deallocate( wr_SortedIndexSave )
        deallocate( w_SortedIndex    )
!!$        deallocate( w_SortedIndexTmp )
        deallocate( w_AbsorpK )


        if ( .not. FlagBinning ) then

          do n = 1, NBand
            a_BandNBin(n) = aa_BandWNIndexBnds(2,n) - aa_BandWNIndexBnds(1,n) + 1
          end do

          NBin = 0
          do n = 1, NBand
            NBin = NBin + a_BandNBin(n)
          end do

          ! this if block exists for a debug purpose
!!$          NBin = NWaveNum

          allocate( a_BinNum    (NBin) )
          allocate( aa_BandBinIndexBnds(2,NBand) )
          allocate( a_BinWeight (NBin) )
          allocate( a_BinPointWN(NBin) )
          allocate( aa_BinWNIndexBnds(2,NBin) )


          do j = 1, NBin
            a_BinNum(j) = j
          end do
          le = 0
          do n = 1, NBand
            ls = le + 1
            le = ls + a_BandNBin(n) - 1
            aa_BandBinIndexBnds(1,n) = ls
            aa_BandBinIndexBnds(2,n) = le
          end do

          a_BinPointWN(1:NBin) = w_WaveNum(aa_BandWNIndexBnds(1,1):aa_BandWNIndexBnds(2,NBand))
          a_BinWeight (1:NBin) = a_BinPointWN(2) - a_BinPointWN(1)

!!$          a_BinWeight  = 0.0d0
!!$          do n = 1, NBand
!!$            if ( aa_BandWNIndexBnds(2,n) - aa_BandWNIndexBnds(1,n) >= 0 ) then
!!$!              do l = aa_BandWNIndexBnds(1,n), aa_BandWNIndexBnds(2,n)
!!$              do l = aa_BandBinIndexBnds(1,n), aa_BandBinIndexBnds(2,n)
!!$                a_BinWeight(l) = a_BinPointWN(2) - a_BinPointWN(1)
!!$              end do
!!$            else
!!$!              do l = aa_BandWNIndexBnds(1,n), aa_BandWNIndexBnds(2,n)
!!$              do l = aa_BandBinIndexBnds(1,n), aa_BandBinIndexBnds(2,n)
!!$                a_BinWeight(l) = 1.0d100
!!$              end do
!!$            end if
!!$          end do


          ! Calculate Planck function ratio
          !
          if ( FlagSortPFRatioAverage ) then

            allocate( a_BinWNBnds(0:NBin) )

            ! Actually, bin is not bounded by wavenumber, but is bounded by 
            ! cumulative probability.
            ! But, here, wavenumber bound is considered for clarity.
            !
            a_BinWNBnds(0) = aa_BandWNBnds(1,1)
            do j = 1, NBin
              a_BinWNBnds(j) = a_BinWNBnds(j-1) + a_BinWeight(j)
            end do

            ! Look for bin bound index
            search_bin_bound_index_1 : do l = 1+1, NWaveNum
              if ( w_WaveNum(l) >= a_BinWNBnds(0) ) exit search_bin_bound_index_1
            end do search_bin_bound_index_1
            aa_BinWNIndexBnds(1,1) = l
            !
            do j = 1, NBin
              search_bin_bound_index : do l = aa_BinWNIndexBnds(1,j), NWaveNum
                if ( w_WaveNum(l) > a_BinWNBnds(j) ) exit search_bin_bound_index
              end do search_bin_bound_index
              aa_BinWNIndexBnds(2,j) = l-1
              if ( j < NBin ) then
                aa_BinWNIndexBnds(1,j+1) = aa_BinWNIndexBnds(2,j) + 1
              end if
            end do




            do k = 1, NPress
              do n = 1, NBand

                do j = aa_BandBinIndexBnds(1,n), aa_BandBinIndexBnds(2,n)
                  ls = aa_BinWNIndexBnds(1,j)
                  le = aa_BinWNIndexBnds(2,j)
                  Sum = 0.0d0
                  do l = ls, le
                    Sum = Sum + rw_PFRatio(k,l)
                  end do
                  if ( ls <= le ) then
                    Sum = Sum / ( le - ls + 1 )
                  else
                    Sum = 1.0d0
                  end if
                  do l = ls, le
                    rw_PFRatio(k,l) = Sum
                  end do
                end do

              end do
            end do

            deallocate( a_BinWNBnds )

          end if


          allocate( rma_BinTAbsorpK(NPress, NMol, NBin) )
          allocate( ra_BinPFRatio  (NPress, NBin) )

!!$          rma_BinTAbsorpK = rmw_TAbsorpK
!!$          ra_BinPFRatio   = rw_PFRatio
          rma_BinTAbsorpK(:,:,1:NBin) = rmw_TAbsorpK(:,:,aa_BandWNIndexBnds(1,1):aa_BandWNIndexBnds(2,NBand))
          ra_BinPFRatio  (:,1:NBin)   = rw_PFRatio  (:,aa_BandWNIndexBnds(1,1):aa_BandWNIndexBnds(2,NBand))

        else

          a_BandNBin = a_ArgBandNBin

          NBin = 0
          do n = 1, NBand
            NBin = NBin + a_BandNBin(n)
          end do

          allocate( a_BinNum    (NBin) )
          allocate( aa_BandBinIndexBnds(2,NBand) )
          allocate( a_BinWeight (NBin) )
          allocate( a_BinPointWN(NBin) )
          allocate( aa_BinWNIndexBnds(2,NBin) )

          call SetupBin( &
            & NWaveNum, w_WaveNum,  &
            & NBand, a_BandNBin, aa_BandWNBnds, &
            & NBin, a_BinNum, aa_BandBinIndexBnds, a_BinPointWN, a_BinWeight, &
            & aa_BinWNIndexBnds &
            & )


          allocate( rma_BinTAbsorpK(NPress, NMol, NBin) )
          allocate( ra_BinPFRatio  (NPress, NBin) )


          do k = 1, NPress
            do n = 1, NBand

              ! Calculate absorption coefficient of a key molecule
              !
              !   pick up a key specie in a band
              !
              !     initialize not to miss an unset variable
              IndexKeyMol = -1
              do m = 1, NMol
                if ( m_MolNum(m) == a_BandKeyMolNum(n) ) then
                  IndexKeyMol = m
                end if
              end do
              !
              if ( FlagKeyMolACAverage ) then

                m = IndexKeyMol
                !
                do j = aa_BandBinIndexBnds(1,n), aa_BandBinIndexBnds(2,n)
                  ls = aa_BinWNIndexBnds(1,j)
                  le = aa_BinWNIndexBnds(2,j)
                  if ( ls <= le ) then
                    Sum = 0.0d0
                    do l = ls, le
                      Sum = Sum + rmw_TAbsorpK(k,m,l)
                    end do
                    Sum = Sum / ( le - ls + 1 )
                  else
                    Sum = 0.0d0
                  end if
                  rma_BinTAbsorpK(k,m,j) = Sum
                end do

              else

                m = IndexKeyMol
                !


                ! sort absorption coefficient in a bin
                if ( (aa_BandBinIndexBnds(2,n)-aa_BandBinIndexBnds(1,n)+1) > 0 ) then
                  do l = 1, aa_BandBinIndexBnds(1,n)-1
                    wm_TAbsorpK(l,m) = 0.0d0
                  end do
                  do l = aa_BandBinIndexBnds(1,n), aa_BandBinIndexBnds(2,n)
                    wm_TAbsorpK(l,m) = rmw_TAbsorpK(k,m,l)
                  end do
                  do l = aa_BandBinIndexBnds(2,n)+1, NWaveNum
                    wm_TAbsorpK(l,m) = 0.0d0
                  end do
                  call sort_heap( &
                    & (aa_BandBinIndexBnds(2,n)-aa_BandBinIndexBnds(1,n)+1), &
                    & wm_TAbsorpK(aa_BandBinIndexBnds(1,n):aa_BandBinIndexBnds(2,n),m) &
                    & )
                  do l = aa_BandBinIndexBnds(1,n), aa_BandBinIndexBnds(2,n)
                    rmw_TAbsorpK(k,m,l) = wm_TAbsorpK(l,m)
                  end do
                  !
                  ls = 2
                  do j = aa_BandBinIndexBnds(1,n), aa_BandBinIndexBnds(2,n)
                    search_gausspoint : do l = ls, NWaveNum-1
                      if ( a_BinPointWN(j) < w_WaveNum(l) ) then
                        exit search_gausspoint
                      end if
                    end do search_gausspoint
                    rma_BinTAbsorpK(k,m,j) = &
                      &   ( rmw_TAbsorpK(k,m,l) - rmw_TAbsorpK(k,m,l-1) ) &
                      & / ( w_WaveNum(l)        - w_WaveNum(l-1)        ) &
                      & * ( a_BinPointWN(j)     - w_WaveNum(l-1)        ) &
                      & + rmw_TAbsorpK(k,m,l-1)
                    ls = l
                  end do
                end if

              end if

              ! Calculate absorption coefficient of non-key molecules by average
              !
              do m = 1, NMol
                if ( m == IndexKeyMol ) cycle

                if ( FlagNonKeyMolACAverage ) then

                  do j = aa_BandBinIndexBnds(1,n), aa_BandBinIndexBnds(2,n)
                    ls = aa_BinWNIndexBnds(1,j)
                    le = aa_BinWNIndexBnds(2,j)
                    if ( ls <= le ) then
                      Sum = 0.0d0
                      do l = ls, le
                        Sum = Sum + rmw_TAbsorpK(k,m,l)
                      end do
                      Sum = Sum / ( le - ls + 1 )
                    else
                      Sum = 0.0d0
                    end if
                    rma_BinTAbsorpK(k,m,j) = Sum
                  end do

                else

                  ! sort absorption coefficient in a bin
                  if ( (aa_BandBinIndexBnds(2,n)-aa_BandBinIndexBnds(1,n)+1) > 0 ) then
                    do l = 1, aa_BandBinIndexBnds(1,n)-1
                      wm_TAbsorpK(l,m) = 0.0d0
                    end do
                    do l = aa_BandBinIndexBnds(1,n), aa_BandBinIndexBnds(2,n)
                      wm_TAbsorpK(l,m) = rmw_TAbsorpK(k,m,l)
                    end do
                    do l = aa_BandBinIndexBnds(2,n)+1, NWaveNum
                      wm_TAbsorpK(l,m) = 0.0d0
                    end do
                    call sort_heap( &
                      & (aa_BandBinIndexBnds(2,n)-aa_BandBinIndexBnds(1,n)+1), &
                      & wm_TAbsorpK(aa_BandBinIndexBnds(1,n):aa_BandBinIndexBnds(2,n),m) &
                      & )
                    do l = aa_BandBinIndexBnds(1,n), aa_BandBinIndexBnds(2,n)
                      rmw_TAbsorpK(k,m,l) = wm_TAbsorpK(l,m)
                    end do
                    !
                    ls = 2
                    do j = aa_BandBinIndexBnds(1,n), aa_BandBinIndexBnds(2,n)
                      search_gausspoint_nonkey : do l = ls, NWaveNum-1
                        if ( a_BinPointWN(j) < w_WaveNum(l) ) then
                          exit search_gausspoint_nonkey
                        end if
                      end do search_gausspoint_nonkey
                      rma_BinTAbsorpK(k,m,j) = &
                        &   ( rmw_TAbsorpK(k,m,l) - rmw_TAbsorpK(k,m,l-1) ) &
                        & / ( w_WaveNum(l)        - w_WaveNum(l-1)        ) &
                        & * ( a_BinPointWN(j)     - w_WaveNum(l-1)        ) &
                        & + rmw_TAbsorpK(k,m,l-1)
                      ls = l
                    end do
                  end if

                end if

              end do



              ! Calculate Planck function ratio
              !
              if ( FlagSortPFRatioAverage ) then

                do j = aa_BandBinIndexBnds(1,n), aa_BandBinIndexBnds(2,n)
                  ls = aa_BinWNIndexBnds(1,j)
                  le = aa_BinWNIndexBnds(2,j)
                  Sum = 0.0d0
                  do l = ls, le
                    Sum = Sum + rw_PFRatio(k,l)
                  end do
                  if ( ls <= le ) then
                    Sum = Sum / ( le - ls + 1 )
                  else
                    Sum = 1.0d0
                  end if
                  ra_BinPFRatio(k,j) = Sum
                end do

              else

                ! sort absorption coefficient in a bin
                if ( (aa_BandBinIndexBnds(2,n)-aa_BandBinIndexBnds(1,n)+1) > 0 ) then
                  call sort_heap( &
                    & (aa_BandBinIndexBnds(2,n)-aa_BandBinIndexBnds(1,n)+1), &
                    & w_PFRatio(aa_BandBinIndexBnds(1,n):aa_BandBinIndexBnds(2,n)) &
                    & )
                  !
                  ls = 2
                  do j = aa_BandBinIndexBnds(1,n), aa_BandBinIndexBnds(2,n)
                    search_gausspoint_pf : do l = ls, NWaveNum-1
                      if ( a_BinPointWN(j) < w_WaveNum(l) ) then
                        exit search_gausspoint_pf
                      end if
                    end do search_gausspoint_pf
                    ra_BinPFRatio(k,j) = &
                      &   ( rw_PFRatio(k,l) - rw_PFRatio(k,l-1) ) &
                      & / ( w_WaveNum(l)    - w_WaveNum(l-1)    ) &
                      & * ( a_BinPointWN(j) - w_WaveNum(l-1)    ) &
                      & + rw_PFRatio(k,l-1)
                    ls = l
                  end do
                end if

              end if




            end do
          end do

        end if

        deallocate( wm_TAbsorpK )
        deallocate( w_PFRatio   )

        call HistoryCreate( &                             ! ヒストリー作成
          & quiet = .true., &
          & file=OutNcFn, title='Absorption coefficient', &
          & source='Program for sorting',   &
          & institution='GFD_Dennou Club dcmodel project',       &
          & dims=(/'r_Press   ', 'm_MolNum  ', 'w_WaveNum ', 'BandNum   ', 'BndNum    ', 'BinNum    ', 'BinWaveNum'/), &
          & dimsizes=(/NPress, NMol, NWaveNum, NBand, 2, NBin, NBin/), &
          & xtypes=(/'double ', 'double ', 'double ', 'integer', 'integer', 'integer', 'double '/), &
!!$      & xtypes=(/'double ', 'double ', 'double ', 'integer', 'integer', 'double '/), &
          & longnames=(/'Pressure        ', 'molecular number', &
          &             'wavenumber      ', 'band number     ', &
          &             'bound number    ', 'bin number      ', &
          &             'wavenumber      '/), &
          & units=(/'Pa ', '1  ', 'm-1', '1  ', '1  ', '1  ', 'm-1'/) )

        call HistoryPut('r_Press'  , r_Press  )                ! 次元変数出力
        call HistoryPut('m_MolNum' , m_MolNum )                ! 次元変数出力
        if ( .not. FlagBinning ) then
          call HistoryPut('w_WaveNum', w_WaveNum)                ! 次元変数出力
        end if
        call HistoryPut('BandNum'  , a_BandNum)                ! 次元変数出力
        call HistoryPut('BndNum'   , a_BndNum )                ! 次元変数出力
        call HistoryPut('BinNum'   , a_BinNum )                ! 次元変数出力
!!$      call HistoryPut('BinNum'   , a_BinPointWN )            ! 次元変数出力

!!$        call HistoryAddVariable( &                          ! 変数定義
!!$          & varname='BinPointWN', &
!!$          & dims=(/'BinNum'/), &
!!$          & longname='bin point', &
!!$          & units='1', xtype='double')
        call HistoryPut('BinWaveNum', a_BinPointWN)      ! 変数出力

        call HistoryAddVariable( &                        ! 変数定義
          & varname='r_Temp', &
          & dims=(/'r_Press'/), &
          & longname='temperature', &
          & units='K', xtype='double')
        call HistoryPut('r_Temp', r_Temp)          ! 変数出力

        call HistoryAddVariable( &                        ! 変数定義
          & varname='m_MolWt', &
          & dims=(/'m_MolNum'/), &
          & longname='molecular weight', &
          & units='kg mol-1', xtype='double')
        call HistoryPut('m_MolWt', m_MolWt)          ! 変数出力

        call HistoryAddVariable( &                        ! 変数定義
          & varname='rm_MixRatio', &
          & dims=(/'r_Press ', 'm_MolNum'/), &
          & longname='volume mixing ratio', &
          & units='1', xtype='double')
        call HistoryPut('rm_MixRatio', rm_MixRatio)          ! 変数出力

        call HistoryAddVariable( &                        ! 変数定義
          & varname='m_MixRatio', &
          & dims=(/'m_MolNum'/), &
          & longname='volume mixing ratio', &
          & units='1', xtype='double')
        call HistoryPut('m_MixRatio', m_MixRatio)          ! 変数出力

        if ( .not. FlagBinning ) then
          call HistoryAddVariable( &                        ! 変数定義
            & varname=AbsCoefVarName, &
            & dims=(/'r_Press  ', 'm_MolNum ', 'w_WaveNum'/), &
            & longname='line & continuum absorption coefficient', &
            & units='m2 kg-1', xtype='double')
          call HistoryPut(AbsCoefVarName, rmw_TAbsorpK)          ! 変数出力
!!$          call HistoryPut(AbsCoefVarName, rma_BinTAbsorpK)          ! 変数出力
          call HistoryAddVariable( &                        ! 変数定義
            & varname='rw_PFRatio', &
            & dims=(/'r_Press  ', 'w_WaveNum'/), &
            & longname='Planck function ratio', &
            & units='1', xtype='double')
          call HistoryPut('rw_PFRatio', rw_PFRatio)          ! 変数出力
        end if

        call HistoryAddVariable( &                          ! 変数定義
          & varname='BandWaveNumBnds', &
          & dims=(/'BndNum ', 'BandNum'/), &
          & longname='band wavenumber bound', &
          & units='m-1', xtype='double')
        call HistoryPut('BandWaveNumBnds', aa_BandWNBnds)          ! 変数出力

        if ( .not. FlagBinning ) then
          call HistoryAddVariable( &                          ! 変数定義
            & varname='BandWaveNumIndexBnds', &
            & dims=(/'BndNum ', 'BandNum'/), &
            & longname='band index bound', &
            & units='1', xtype='integer')
          call HistoryPut('BandWaveNumIndexBnds', aa_BandWNIndexBnds)     ! 変数出力
        end if

        call HistoryAddVariable( &                          ! 変数定義
          & varname='BandNBin', &
          & dims=(/'BandNum'/), &
          & longname='number of bins in a band', &
          & units='1', xtype='integer')
        call HistoryPut('BandNBin', a_BandNBin) ! 変数出力

        call HistoryAddVariable( &                          ! 変数定義
          & varname='BandBinIndexBnds', &
          & dims=(/'BndNum ', 'BandNum'/), &
          & longname='index bound of bins in a band', &
          & units='1', xtype='integer')
        call HistoryPut('BandBinIndexBnds', aa_BandBinIndexBnds) ! 変数出力

        call HistoryAddVariable( &                          ! 変数定義
          & varname='BandKeyMolNum', &
          & dims=(/'BandNum'/), &
          & longname='band key molecular number', &
          & units='1', xtype='integer')
        call HistoryPut('BandKeyMolNum', a_BandKeyMolNum) ! 変数出力

        call HistoryAddVariable( &                          ! 変数定義
          & varname='BinWeight', &
          & dims=(/'BinNum'/), &
          & longname='bin weight', &
          & units='1', xtype='double')
        call HistoryPut('BinWeight', a_BinWeight)          ! 変数出力

        call HistoryAddVariable( &                        ! 変数定義
          & varname='rma_TAbsorpK', &
!!$          & dims=(/'r_Press ', 'm_MolNum', 'BinNum  '/), &
          & dims=(/'r_Press   ', 'm_MolNum  ', 'BinWaveNum'/), &
          & longname='line & continuum absorption coefficient', &
          & units='m2 kg-1', xtype='double')
        call HistoryAddAttr( 'rma_TAbsorpK', 'missing_value', 0.0d0 )
        call HistoryPut('rma_TAbsorpK', rma_BinTAbsorpK)          ! 変数出力

        call HistoryAddVariable( &                        ! 変数定義
          & varname='ra_PFRatio', &
!!$          & dims=(/'r_Press ', 'BinNum  '/), &
          & dims=(/'r_Press   ', 'BinWaveNum'/), &
          & longname='Planck function ratio', &
          & units='1', xtype='double')
        call HistoryPut('ra_PFRatio', ra_BinPFRatio)          ! 変数出力

        call HistoryClose( quiet = .true. )


        deallocate( ra_BandAvePF )
        deallocate( ra_BinPFRatio   )
        deallocate( rma_BinTAbsorpK )
        deallocate( aa_BinWNIndexBnds )
        deallocate( a_BinPointWN    )
        deallocate( a_BinWeight     )
        deallocate( aa_BandBinIndexBnds )
        deallocate( a_BinNum        )

        deallocate( rm_MMR )

        deallocate( w_WaveNum )
        deallocate( m_MolNum )
        deallocate( m_MolWt )
        deallocate( r_Press )
        deallocate( r_Temp  )

        deallocate( rmw_TAbsorpK )
        deallocate( rw_PFRatio   )

      end do
    end do


  end subroutine SortAC

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

  function GetBandNum( NBand, a_WNBnds, WaveNum ) result( BandNum )

    integer , intent(in ) :: NBand
    real(DP), intent(in ) :: a_WNBnds(0:NBand)
    real(DP), intent(in ) :: WaveNum

    integer :: BandNum

    integer :: n


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

    search_Band : do n = 1, NBand
       if ( WaveNum <= a_WNBnds(n) ) exit search_Band
    end do search_Band

    if ( n >= NBand + 1 ) then
!!$      write( 6, * ) 'Given wavenumber is out of given bands: ', &
!!$        & WaveNum, a_WNBnds(NBand)
!!$      stop
      n = -999
    end if

    BandNum = n

  end function GetBandNum

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

  subroutine SetupBin( &
    & NWaveNum, w_WaveNum,  &
    & NBand, a_BandNBin, aa_BandWNBnds, &
    & NBin, a_BinNum, aa_BandBinIndexBnds, a_BinPointWN, a_BinWeight, &
    & aa_BinWNIndexBnds &
    & )

    use gauss_quad, only : GauLeg

    integer, intent(in ) :: NWaveNum
    real(8), intent(in ) :: w_WaveNum(NWaveNum)
    integer, intent(in ) :: NBand
    integer, intent(in ) :: a_BandNBin(NBand)
                           ! Number of bins in bands
    real(8), intent(in ) :: aa_BandWNBnds(2,NBand)
                           ! Wavenumber bounds of bands
    integer, intent(in ) :: NBin
                           ! Number of bins
    integer, intent(out) :: a_BinNum    (NBin)
                           ! Bin number
    integer, intent(out) :: aa_BandBinIndexBnds(2,NBand)
                           ! Index bounds of bins in bands
    real(8), intent(out) :: a_BinPointWN(NBin)
                           ! Bin point wavenumber
    real(8), intent(out) :: a_BinWeight (NBin)
                           ! Bin weight
    integer, intent(out) :: aa_BinWNIndexBnds(2,NBin)
                           ! Wavenumber index bounds of bins


    ! Local variables
    !
    real(8) :: a_BinWNBnds      (0:NBin)
                           ! Wavenumber bounds of bins

    real(8) :: WNs
    real(8) :: WNe
    integer :: NSubBin

    integer :: j
    integer :: l
    integer :: n

    integer :: ls
    integer :: le


    do j = 1, NBin
      a_BinNum(j) = j
    end do
    le = 0
    do n = 1, NBand
      ls = le + 1
      le = ls + a_BandNBin(n) - 1
      aa_BandBinIndexBnds(1,n) = ls
      aa_BandBinIndexBnds(2,n) = le
    end do
    !
    a_BinWeight = 0.0d0
    do n = 1, NBand

      ! 1 region
      !
!!$      ls = aa_BandBinIndexBnds(1,n)
!!$      le = aa_BandBinIndexBnds(2,n)
!!$      call GauLeg( aa_BandWNBnds(1,n), aa_BandWNBnds(2,n), &
!!$        & a_BandNBin(n), a_BinPointWN(ls:le), a_BinWeight(ls:le) )

      ! 2 regions
      !   1st: from 0 to 0.95
      !   2nd: from 0.95 to 0.95
      !
      ls = aa_BandBinIndexBnds(1,n)
      le = ls + ( aa_BandBinIndexBnds(2,n) - aa_BandBinIndexBnds(1,n) ) / 2
      WNs = aa_BandWNBnds(1,n)
      WNe = WNs + ( aa_BandWNBnds(2,n) - aa_BandWNBnds(1,n) ) * 0.95d0
      NSubBin = a_BandNBin(n) / 2
      call GauLeg( WNs, WNe, &
        & NSubBin, a_BinPointWN(ls:le), a_BinWeight(ls:le) )
      ls = le + 1
      le = aa_BandBinIndexBnds(2,n)
      WNs = WNe
      WNe = aa_BandWNBnds(2,n)
      NSubBin = a_BandNBin(n) - NSubBin
      call GauLeg( WNs, WNe, &
        & NSubBin, a_BinPointWN(ls:le), a_BinWeight(ls:le) )

    end do

    ! Actually, bin is not bounded by wavenumber, but is bounded by 
    ! cumulative probability.
    ! But, here, wavenumber bound is considered for clarity.
    !
    a_BinWNBnds(0) = aa_BandWNBnds(1,1)
    do j = 1, NBin
      a_BinWNBnds(j) = a_BinWNBnds(j-1) + a_BinWeight(j)
    end do

    ! Look for bin bound index
    search_bin_bound_index_1 : do l = 1+1, NWaveNum
      if ( w_WaveNum(l) >= a_BinWNBnds(0) ) exit search_bin_bound_index_1
    end do search_bin_bound_index_1
    aa_BinWNIndexBnds(1,1) = l
    !
    do j = 1, NBin
      search_bin_bound_index : do l = aa_BinWNIndexBnds(1,j), NWaveNum
        if ( w_WaveNum(l) > a_BinWNBnds(j) ) exit search_bin_bound_index
      end do search_bin_bound_index
      aa_BinWNIndexBnds(2,j) = l-1
      if ( j < NBin ) then
        aa_BinWNIndexBnds(1,j+1) = aa_BinWNIndexBnds(2,j) + 1
      end if
    end do


  end subroutine SetupBin

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

end module sortac_module
