program main

  ! モジュール引用 ; USE statements
  !

  use dc_types

  ! netcdf データ入出力
  ! netcdf data input or output
  !
  use gtool_history

  use gridset, only : GridSetInit

  use rad_utils, only : InqDimLength
  use interpolate, only : Interpolate1D


  ! 物理定数設定1
  ! physical constants setteing 1
  !
  use constants0, only: PI, &
    &                   StB

  use sortac_module, only : sortac

!!$  use rad_ckd_utils, only : GetDimLen

  use rad_flux, only : RadFluxCalcFlux, RadFluxInit, RadFluxFinalize, RadFluxOutputInit, RadFluxOutputFinalize

  use planetconst, only : Grav


  ! 宣言文 ; Declaration statements
  !
  implicit none


  character(128) :: ctlfn = 'sortac.conf'

  integer :: FileUnit
  logical :: FlagEx
  logical :: FlagOpened
  integer :: iostat

  logical :: FlagChkConstH2OMR
  logical :: FlagKeyMolACAverage
  logical :: FlagNonKeyMolACAverage
  logical :: FlagSortPFRatioAverage

  integer, parameter :: MaxNmlArySize = 1024
  real(8)            :: WNBnds(MaxNmlArySize)
  integer            :: BandNBin(MaxNmlArySize)
  integer            :: KeyMolNum(MaxNmlArySize)

  integer :: NBand
  real(8), allocatable :: a_BandWNBnds   (:)
  integer, allocatable :: a_BandKeyMolNum(:)
  integer, allocatable :: a_BandNBin     (:)

  integer              :: NBandCur = 1
  real(8), allocatable :: a_BandCurWNBnds   (:)
  integer, allocatable :: a_BandCurKeyMolNum(:)
  integer, allocatable :: a_BandCurNBin     (:)

  integer :: n


  ! 設定変数, 他
  ! setteing parameter, values
  !
  character(128) :: InOrgOptPropNcFN
                          ! 
                          ! input file name for original optical properties
  logical        :: FlagOutSpeFlux
  character(128) :: OutSpeFluxNcFN
                          ! 
                          ! 
  character(128) :: OutFluxNcFN
                          ! 出力ファイル名
                          ! output file name

  character(128) :: InProfNcFn

!!$  real(8) :: SolTemp
                          ! 中心星表面温度 [K]
                          ! surface temperature of central star [K]
!!$  real(8) :: SolConst
                          ! 太陽定数 [W m-2]
                          ! solar constant [W m-2]
!!$  real(8) :: InAngle
                          ! sec (太陽光入射角) [1]
                          ! sec (angle of incidence) [1]
!!$  real(8) :: SurfAlbedo
                          ! 地表面アルベド [1]
                          ! surface albedo [1]

  real(8), allocatable :: r_Press(:)
  real(8), allocatable :: r_Temp (:)
  real(8), allocatable :: r_MMRH2OVap(:)
  real(8), allocatable :: r_MMRCO2   (:)
  real(8), allocatable :: r_MMRO3    (:)

  real(8), allocatable :: z_Press(:)
  real(8)              :: SurfTemp
  real(8), allocatable :: z_Temp (:)
  real(8), allocatable :: zm_MMR(:,:)
  real(8), allocatable :: z_MMRH2OVap(:)
  real(8), allocatable :: z_MMRCO2   (:)
  real(8), allocatable :: z_MMRO3    (:)

  real(8), allocatable :: r_RadUwFluxTot (:)
                          ! upward flux
  real(8), allocatable :: r_RadDwFluxTot (:)
                          ! downward flux

  real(8), allocatable :: r_RadUwFluxBandLBL(:)
                          ! upward flux
  real(8), allocatable :: r_RadDwFluxBandLBL(:)
                          ! downward flux

  real(8), allocatable :: r_RadUwFlux    (:)
                          ! upward flux
  real(8), allocatable :: r_RadDwFlux    (:)
                          ! downward flux

!!$  real(8) :: RatioSolarConst

!!$  real(8), pointer :: a_Axis(:) => null()

  real(8) :: RadUwFluxErrMax
  real(8) :: RadDwFluxErrMax
  real(8) :: RadUwFluxRatioMax
  real(8) :: RadDwFluxRatioMax
  real(8) :: RadNetFluxErrMax

  character(128) :: InNcFN
  character(128) :: OutNcFN

  integer :: kmax
  integer :: NMol

  logical :: FlagSort
  logical :: FlagBinning
  logical :: FlagCorrelation
  logical :: FlagTwoMolAbsOverlap

  character(128) :: Method

  integer :: iLoop
  integer :: iBand

  integer :: k


  ! namelist file の設定
  ! setting of namelist file
  !
!!$  namelist /flux_nml/ &
!!$    & SurfAlbedo,      & 
!!$    & SolTemp,         &
!!$    & SolConst,        & 
!!$    & InAngle


  ! namelist file の読み込み
  ! read namelist file


  namelist /flag_nml/ &
!!$    & FlagChkConstH2OMR, FlagKeyMolACAverage, FlagNonKeyMolACAverage
    & FlagChkConstH2OMR, FlagKeyMolACAverage, FlagNonKeyMolACAverage, &
    & FlagSortPFRatioAverage
  namelist /file_nml/ &
    & InProfNcFn, &
!!$    & InNcFN, OutNcFN
    & InOrgOptPropNcFN, &
    & FlagOutSpeFlux, &
    & OutSpeFluxNcFN, &
    & OutFluxNcFN,    &
    & OutNcFN
  namelist /band_nml/ &
    & NBand, WNBnds, BandNBin, KeyMolNum


  InProfNcFn = ''

  FlagChkConstH2OMR      = .true.
  FlagKeyMolACAverage    = .false.
!!$  FlagNonKeyMolACAverage = .true.
  FlagNonKeyMolACAverage = .false.

  FlagSortPFRatioAverage     = .true.

  InNcFN  = ''
  OutNcFN = ''

  NBand = 1
  WNBnds    = -999.0d0
  WNBnds(1) = 0.0d0
  WNBnds(2) = 1.0d100
  KeyMolNum    = -1
  KeyMolNum(1) = 1
  BandNBin     = -1
  BandNBin(1)  = 10


  ! Read a namelist file

  FileUnit = GetFileUnit( ctlfn )

  open( FileUnit, file=ctlfn, status='unknown' )
  !
  rewind( FileUnit )
  read( FileUnit, nml = file_nml, iostat = iostat )
  if ( iostat > 0 ) then
    stop 'Namelist in file_nml is something wrong.'
  end if
  write( 6, nml = file_nml )
  !
  open( FileUnit, file=ctlfn, status='unknown' )
  rewind( FileUnit )
  read( FileUnit, nml = flag_nml, iostat = iostat )
  if ( iostat > 0 ) then
    stop 'Namelist is something wrong.'
  end if
  write( 6, nml = flag_nml )
  !
  rewind( FileUnit )
  read( FileUnit, nml = band_nml, iostat = iostat )
  if ( iostat > 0 ) then
    stop 'Namelist is something wrong.'
  end if
  write( 6, nml = band_nml )
  !
  close( FileUnit )


  kmax = InqDimLength( InProfNcFN, 'z_Temp' )
  NMol = InqDimLength( InProfNcFN, 'm_MolNum' )

  allocate( r_Press    (0:kmax) )
  !
  allocate( z_Press    (1:kmax) )
  allocate( z_Temp     (1:kmax) )
  allocate( z_MMRH2OVap(1:kmax) )
  allocate( z_MMRCO2   (1:kmax) )
  allocate( z_MMRO3    (1:kmax) )

  call HistoryGet( InProfNcFN, 'r_Press'    , r_Press    , quiet = .true. )
  !
  call HistoryGet( InProfNcFN, 'z_Press'    , z_Press    , quiet = .true. )
  call HistoryGet( InProfNcFN, 'SurfTemp'   , SurfTemp   , quiet = .true. )
  call HistoryGet( InProfNcFN, 'z_Temp'     , z_Temp     , quiet = .true. )
  !
  allocate( zm_MMR(1:kmax,NMol) )
  call HistoryGet( InProfNcFN, 'zm_MassMixRatio', zm_MMR, quiet = .true. )
  z_MMRH2OVap = zm_MMR(:,1)
  z_MMRCO2    = zm_MMR(:,2)
  z_MMRO3     = zm_MMR(:,3)
  deallocate( zm_MMR )


  allocate( r_Temp     (0:kmax) )
  allocate( r_MMRH2OVap(0:kmax) )
  allocate( r_MMRCO2   (0:kmax) )
  allocate( r_MMRO3    (0:kmax) )

  call Interpolate1D( &
    & kmax, z_Press, z_Temp,         &
    & kmax+1, r_Press,                  &
    & r_Temp                         &
    & )
  call Interpolate1D( &
    & kmax, z_Press, z_MMRH2OVap,         &
    & kmax+1, r_Press,                  &
    & r_MMRH2OVap                         &
    & )
  call Interpolate1D( &
    & kmax, z_Press, z_MMRCO2,         &
    & kmax+1, r_Press,                  &
    & r_MMRCO2                         &
    & )
  call Interpolate1D( &
    & kmax, z_Press, z_MMRO3,         &
    & kmax+1, r_Press,                  &
    & r_MMRO3                         &
    & )


!!$  write( 6, * ) r_Press
!!$  write( 6, * ) z_Temp
!!$  write( 6, * ) z_MMRH2OVap
!!$  write( 6, * ) z_MMRCO2
!!$  write( 6, * ) z_MMRO3


  allocate( r_RadUwFlux(0:kmax) )
  allocate( r_RadDwFlux(0:kmax) )


  call GridSetInit( &
    & kmax &
    & )

  Method = 'LBL'
  call RadFluxInit( &
    & Grav, &
    & InOrgOptPropNcFN, &
    & Method &
    & )
  if ( FlagOutSpeFlux ) then
    call RadFluxOutputInit( &
      & InOrgOptPropNcFN, &
      & OutSpeFluxNcFN, &
      & kmax, r_Press &
      & )
  end if
  call RadFluxCalcFlux( &
    & kmax, &
    & r_Press, &
    & SurfTemp, r_Temp, r_MMRH2OVap, r_MMRCO2, r_MMRO3, &
    & -1.0d0, -1.0d0, &
    & r_RadUwFlux, r_RadDwFlux &
    & )
  if ( FlagOutSpeFlux ) then
    call RadFluxOutputFinalize
  end if
  call RadFluxFinalize

  call OutputFlux( &
    & OutFluxNcFN,    &
    & kmax, r_Press, &
    & r_RadUwFlux, r_RadDwFlux &
    & )

  allocate( r_RadUwFluxTot (0:kmax) )
  allocate( r_RadDwFluxTot (0:kmax) )

  r_RadUwFluxTot = r_RadUwFlux
  r_RadDwFluxTot = r_RadDwFlux







  InNcFN = InOrgOptPropNcFN



  !TESTTESTTESTTESTTESTTESTTESTTESTTESTTEST
!!$  NBand = 3000.0d2 / 100.0d2
!!$  do n = 1, NBand+1
!!$    WNBnds(n) = 100.0d2 * (n-1)
!!$  end do
!!$  BandNBin = 10
!!$  KeyMolNum = 1
  !TESTTESTTESTTESTTESTTESTTESTTESTTESTTEST



  if ( NBand+1 > MaxNmlArySize ) then
    stop 'Array size is too large.'
  end if
  do n = 1, NBand
    if ( WNBnds(n) > WNBnds(n+1) ) then
      stop 'Band bounds does not increase monotonicaly.'
    end if
  end do


!!$  RatioSolarConst = SolConst/(StB*SolTemp*SolTemp*SolTemp*SolTemp)


  allocate( a_BandWNBnds   (0:NBand) )
  allocate( a_BandKeyMolNum(1:NBand) )
  allocate( a_BandNBin     (1:NBand) )

  a_BandWNBnds   (0:NBand) = WNBnds   (1:NBand+1)
  a_BandKeyMolNum(1:NBand) = KeyMolNum(1:NBand  )
  a_BandNBin     (1:NBand) = BandNBin (1:NBand  )



  allocate( r_RadUwFluxBandLBL(0:kmax) )
  allocate( r_RadDwFluxBandLBL(0:kmax) )


  do iBand = 0, NBand
!!$  do iBand = 0, 0
!!$  do iBand = 1, 1
!!$  do iBand = 4, 4

    if ( iBand == 0 ) then
      NBandCur = NBand
      allocate( a_BandCurWNBnds   (0:NBandCur) )
      allocate( a_BandCurKeyMolNum(1:NBandCur) )
      allocate( a_BandCurNBin     (1:NBandCur) )

      a_BandCurWNBnds   (0:NBandCur) = a_BandWNBnds   (0:NBand)
      a_BandCurKeyMolNum(1:NBandCur) = a_BandKeyMolNum(1:NBand)
      a_BandCurNBin     (1:NBandCur) = a_BandNBin     (1:NBand)
    else
      NBandCur = 1
      allocate( a_BandCurWNBnds   (0:NBandCur) )
      allocate( a_BandCurKeyMolNum(1:NBandCur) )
      allocate( a_BandCurNBin     (1:NBandCur) )

      a_BandCurWNBnds   (0:NBandCur) = a_BandWNBnds   (iBand-1:iBand)
      a_BandCurKeyMolNum(1:NBandCur) = a_BandKeyMolNum(iBand  :iBand)
      a_BandCurNBin     (1:NBandCur) = a_BandNBin     (iBand  :iBand)
    end if



    call RadFluxInit( &
      & Grav, &
      & InOrgOptPropNcFN, &
      & "LBL" &
      & )
    call RadFluxCalcFlux( &
      & kmax, &
      & r_Press, &
      & SurfTemp, r_Temp, r_MMRH2OVap, r_MMRCO2, r_MMRO3, &
      & a_BandCurWNBnds(0), a_BandCurWNBnds(NBandCur), &
      & r_RadUwFlux, r_RadDwFlux &
      & )
    call RadFluxFinalize

    r_RadUwFluxBandLBL = r_RadUwFlux
    r_RadDwFluxBandLBL = r_RadDwFlux



!!$    do iLoop = 0, 0
!!$    do iLoop = 0, 1
!!$    do iLoop = 0, 3
!!$    do iLoop = 0, 4
    do iLoop = 0, 5
!!$    do iLoop = 1, 1
!    do iLoop = 2, 2
!!$    do iLoop = 3, 3
!!$    do iLoop = 4, 4
      write( 6, * ) 'iBand = ', iBand
      write( 6, * ) 'iLoop = ', iLoop


      select case ( iLoop )
      case ( 0 )
        FlagSort             = .false.
        Method               = 'LBL'
        FlagCorrelation      = .false.
        FlagTwoMolAbsOverlap = .false.
        FlagBinning          = .false.
      case ( 1 )
        FlagSort             = .false.
        Method               = 'LBLConstPF'
        FlagCorrelation      = .false.
        FlagTwoMolAbsOverlap = .false.
        FlagBinning          = .false.
      case ( 2 )
        FlagSort             = .true.
        Method               = 'LBLConstPF'
        FlagCorrelation      = .false.
        FlagTwoMolAbsOverlap = .false.
        FlagBinning          = .false.
      case ( 3 )
        FlagSort             = .true.
        Method               = 'LBLConstPF'
        FlagCorrelation      = .false.
        FlagTwoMolAbsOverlap = .true.
        FlagBinning          = .false.
      case ( 4 )
        FlagSort             = .true.
        Method               = 'KD'
        FlagCorrelation      = .false.
        FlagTwoMolAbsOverlap = .true.
        FlagBinning          = .false.
      case ( 5 )
        FlagSort             = .true.
        Method               = 'KD'
        FlagCorrelation      = .true. ! test
        FlagTwoMolAbsOverlap = .true.
        FlagBinning          = .false.
      case ( 6 )
        FlagSort             = .true.
        Method               = 'KD'
        FlagCorrelation      = .false.
        FlagTwoMolAbsOverlap = .true.
        FlagBinning          = .true.
      case default
        stop 'Unexpected value of iLoop'
      end select

      write( 6, * ) "FlagSort             : ", FlagSort
      write( 6, * ) "Method               : ", trim( Method )
      write( 6, * ) "FlagCorrelation      : ", FlagCorrelation
      write( 6, * ) "FlagTwoMolAbsOverlap : ", FlagTwoMolAbsOverlap
      write( 6, * ) "FlagBinning          : ", FlagBinning

      InNcFN = InOrgOptPropNcFN
      call sortac( &
        & FlagSort, &
        & InNcFN, OutNcFN, &
        & FlagChkConstH2OMR, FlagBinning, &
        & FlagKeyMolACAverage, FlagNonKeyMolACAverage, &
        & FlagSortPFRatioAverage, &
!        & NBand, a_BandWNBnds, a_BandKeyMolNum, a_BandNBin &
        & FlagCorrelation, &
        & FlagTwoMolAbsOverlap, &
        & NBandCur, a_BandCurWNBnds, a_BandCurKeyMolNum, a_BandCurNBin &
        & )

      InNcFN = OutNcFN
      !
      call RadFluxInit( &
        & Grav, &
        & InNcFN, &
        & Method, &
        & a_BandCurWNBnds &
        & )
      call RadFluxCalcFlux( &
        & kmax, &
        & r_Press, &
        & SurfTemp, r_Temp, r_MMRH2OVap, r_MMRCO2, r_MMRO3, &
        & a_BandCurWNBnds(0), a_BandCurWNBnds(NBandCur), &
        & r_RadUwFlux, r_RadDwFlux &
        & )
      call RadFluxFinalize

      RadUwFluxErrMax = 0.0d0
      RadDwFluxErrMax = 0.0d0
      do k = kmax, 0, -1
        RadUwFluxErrMax = &
          & max( RadUwFluxErrMax, &
          &   abs( ( r_RadUwFlux(k) - r_RadUwFluxBandLBL(k) ) / r_RadUwFluxBandLBL(k) ) )
        RadDwFluxErrMax = &
          & max( RadDwFluxErrMax, &
          &   abs( ( r_RadDwFlux(k) - r_RadDwFluxBandLBL(k) ) / ( r_RadDwFluxBandLBL(k) + 1.0d-10 ) ) )
      end do

      RadNetFluxErrMax = 0.0d0
      do k = kmax, 0, -1
        RadNetFluxErrMax = &
          & max( RadNetFluxErrMax, &
          &   abs( ( ( r_RadUwFlux(k) - r_RadDwFlux(k) ) &
          &        - ( r_RadUwFluxBandLBL(k) - r_RadDwFluxBandLBL(k) ) ) &
          &        / ( r_RadUwFluxBandLBL(k) - r_RadDwFluxBandLBL(k) ) ) )
      end do

      RadUwFluxRatioMax = 0.0d0
      RadDwFluxRatioMax = 0.0d0
      do k = kmax, 0, -1
        RadUwFluxRatioMax = &
          & max( RadUwFluxRatioMax, &
          &   abs( r_RadUwFlux(k) / r_RadUwFluxTot(k) ) )
        RadDwFluxRatioMax = &
          & max( RadDwFluxRatioMax, &
          &   abs( r_RadDwFlux(k) / ( r_RadDwFluxTot(k) + 1.0d-10 ) ) )
      end do

      write( 6, * ) iBand, iLoop, a_BandCurWnBnds(0), '||', &
        & RadNetFluxErrMax, '|', RadDwFluxErrMax, RadDwFluxRatioMax, '|', RadUwFluxErrMax, RadUwFluxRatioMax
      write( 60, * ) iBand, iLoop, int(a_BandCurWnBnds(0)), int(a_BandCurWnBnds(NBandCur)), '||', &
        & RadNetFluxErrMax, '|', RadDwFluxErrMax, RadDwFluxRatioMax, '|', RadUwFluxErrMax, RadUwFluxRatioMax
      call flush( 60 )


    end do


    deallocate( a_BandCurWNBnds    )
    deallocate( a_BandCurKeyMolNum )
    deallocate( a_BandCurNBin      )

  end do


  deallocate( r_RadUwFlux )
  deallocate( r_RadDwFlux )

  deallocate( r_Press     )
  deallocate( z_Press     )
  deallocate( z_Temp      )
  deallocate( r_Temp      )
  deallocate( z_MMRH2OVap )
  deallocate( z_MMRCO2    )
  deallocate( z_MMRO3     )


  deallocate( a_BandWNBnds    )
  deallocate( a_BandKeyMolNum )
  deallocate( a_BandNBin      )

  deallocate( r_RadUwFluxBandLBL )
  deallocate( r_RadDwFluxBandLBL )

  deallocate( r_RadUwFluxTot )
  deallocate( r_RadDwFluxTot )

contains

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

  function GetFileUnit( InFileName ) result( FileUnit )

    character(*), intent(in) :: InFileName

    integer :: FileUnit


    ! Read a namelist file

    inquire( file = InFileName, opened = FlagOpened, exist = FlagEx )
    if ( .not. FlagEx ) then
      stop 'Configure file does not exist.'
    end if
    if ( FlagOpened ) then
      stop 'Configure file is opened.'
    end if
    search_fileunit : do FileUnit = 50, 99
      inquire( unit = FileUnit, opened = FlagOpened, iostat = iostat )
      if ( .not. FlagOpened ) exit search_fileunit
    end do search_fileunit

  end function GetFileUnit

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

  subroutine OutputFlux( &
    & OutNcFN, &
    & kmax, r_Press, &
    & r_RadUwFlux, r_RadDwFlux &
    & )

    character(*), intent(in ) :: OutNcFN
    integer     , intent(in ) :: kmax
    real(DP)    , intent(in ) :: r_Press(0:kmax)
    real(DP)    , intent(in ) :: r_RadUwFlux(0:kmax)
    real(DP)    , intent(in ) :: r_RadDwFlux(0:kmax)


    ! Local variables
    !
    type(GT_HISTORY), save :: GTHistOut


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

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

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

    call HistoryPut('r_RadFluxUp',r_RadUwFlux,history=GTHistOut)
    call HistoryPut('r_RadFluxDn',r_RadDwFlux,history=GTHistOut)


    call HistoryClose( history=GTHistOut, quiet = .true. )

  end subroutine OutputFlux

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

end program main
