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, Cp


  ! 宣言文 ; Declaration statements
  !
  implicit none


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

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

  logical :: FlagChkConstH2OMR
  logical :: FlagKeyMolACAverage
  logical :: FlagNonKeyMolACAverage
  logical :: FlagPFRatioOne
  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) :: OutOptPropNcFN
  character(128) :: OutFluxNcFN

  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_RadUwFluxBandLBL(:)
                          ! upward flux
  real(8), allocatable :: r_RadDwFluxBandLBL(:)
                          ! downward flux

  real(8), allocatable :: z_DTempDt      (:)
                          ! tendency
  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 :: FlagBinningDemonstration
  logical :: FlagCorrelation
  logical :: FlagTwoMolAbsOverlap

  integer              :: NGaussReg
  real(8), allocatable :: a_GaussRegBnds(:)
  real(8)              :: GaussRegBnds(MaxNmlArySize)


  character(128) :: Method

  integer :: iLoop
  integer :: iBand

  integer :: iBands
  integer :: iBande
  integer :: iModes
  integer :: iModee

  integer :: k


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


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


  namelist /band_mode_nml/ &
    & iBands, iBande, iModes, iModee
  namelist /flag_nml/ &
!!$    & FlagChkConstH2OMR, FlagKeyMolACAverage, FlagNonKeyMolACAverage
    & FlagChkConstH2OMR
  namelist /file_nml/ &
    & InProfNcFn, &
!!$    & InNcFN, OutNcFN
    & InOrgOptPropNcFN, &
    & FlagOutSpeFlux, &
    & OutSpeFluxNcFN, &
    & OutFluxNcFN,    &
    & OutOptPropNcFN
  namelist /band_nml/ &
    & NBand, WNBnds, BandNBin, KeyMolNum

  namelist /gauss_region_nml/ &
    & NGaussReg, GaussRegBnds


  InProfNcFn = ''

  FlagChkConstH2OMR      = .true.

  InNcFN  = ''
  OutOptPropNcFN = ''

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

  GaussRegBnds = 0.0d0


  ! 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 = band_mode_nml, iostat = iostat )
  write( 6, * ) iostat
  if ( iostat /= 0 ) then
    stop 'Namelist is something wrong.'
  end if
  write( 6, nml = band_mode_nml )
  !
  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 )
  !
  rewind( FileUnit )
  read( FileUnit, nml = gauss_region_nml, iostat = iostat )
  if ( iostat /= 0 ) then
    stop 'Namelist is something wrong.'
  end if
  write( 6, nml = gauss_region_nml )
  !
  close( FileUnit )


  allocate( a_GaussRegBnds(0:NGaussReg) )
  a_GaussRegBnds = GaussRegBnds(1:NGaussReg+1)
  ! check
  if ( NGaussReg <= 0 ) then
    write( 6, * ) "NGaussReg is not appropriate: ", NGaussReg
    stop
  end if
  if ( ( a_GaussRegBnds(0) /= 0.0d0 ) .or. ( a_GaussRegBnds(NGaussReg) /= 1.0d0 ) ) then
    write( 6, * ) "GaussRegBnds is not appropriate: ", a_GaussRegBnds
    stop
  end if


  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( z_DTempDt  (1:kmax) )
  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 CalcTendency( &
    & Grav, Cp, &
    & kmax, r_Press, &
    & r_RadUwFlux, r_RadDwFlux, &
    & z_DTempDt &
    & )
  call OutputFluxTendency( &
    & OutFluxNcFN,    &
    & kmax, r_Press, z_Press, &
    & r_RadUwFlux, r_RadDwFlux, z_DTempDt &
    & )






  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 = max( iBands, 0 ), min( iBande, NBand+1 )
!!$  do iBand = 1, NBand+1
!!$  do iBand = NBand+1, NBand+1
!!$  do iBand = 0, 0
!!$  do iBand = 1, 1
!!$  do iBand = 4, 4
!!$  do iBand = 10, 10

!!$    if ( iBand == 0 ) then
    if ( ( iBand == 0 ) .or. ( iBand == NBand+1 ) ) 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 = iModes, iModee
!!$    do iLoop = 0, 0
!!$    do iLoop = 0, 1
!!$    do iLoop = 0, 3
!!$    do iLoop = 0, 4
!!$    do iLoop = 0, 5
!!$    do iLoop = 0, 6
!!$    do iLoop = 0, 7
!!$    do iLoop = 0, 9
!!$    do iLoop = 1, 1
!!$    do iLoop = 2, 2
!!$    do iLoop = 3, 3
!!$    do iLoop = 4, 4
!!$    do iLoop = 5, 5
!!$    do iLoop = 6, 6
!!$    do iLoop = 7, 7
!!$    do iLoop = 7, 9
!!$    do iLoop = 8, 9
      write( 6, * ) 'iBand = ', iBand
      write( 6, * ) 'iLoop = ', iLoop


      select case ( iLoop )
      case ( 0 )
        FlagSort                 = .false.
        Method                   = 'LBL'
        FlagCorrelation          = .false.
        FlagTwoMolAbsOverlap     = .false.
        FlagPFRatioOne           = .false.
        FlagSortPFRatioAverage   = .false.
        FlagBinning              = .false.
        FlagBinningDemonstration = .false.
        FlagKeyMolACAverage      = .false.
        FlagNonKeyMolACAverage   = .true.
      case ( 1 )
        FlagSort               = .false.
        Method                 = 'LBLConstPF'
        FlagCorrelation        = .false.
        FlagTwoMolAbsOverlap   = .false.
        FlagPFRatioOne         = .false.
        FlagSortPFRatioAverage = .false.
        FlagBinning            = .false.
        FlagBinningDemonstration = .false.
        FlagKeyMolACAverage      = .false.
        FlagNonKeyMolACAverage   = .true.
      case ( 2 )
        FlagSort               = .true.
        Method                 = 'LBLConstPF'
        FlagCorrelation        = .false.
        FlagTwoMolAbsOverlap   = .false.
        FlagPFRatioOne         = .false.
        FlagSortPFRatioAverage = .false.
        FlagBinning            = .false.
        FlagBinningDemonstration = .false.
        FlagKeyMolACAverage      = .false.
        FlagNonKeyMolACAverage   = .true.
      case ( 3 )
        FlagSort               = .true.
        Method                 = 'LBLConstPF'
        FlagCorrelation        = .false.
        FlagTwoMolAbsOverlap   = .false.
        FlagPFRatioOne         = .true.
        FlagSortPFRatioAverage = .false.
        FlagBinning            = .false.
        FlagBinningDemonstration = .false.
        FlagKeyMolACAverage      = .false.
        FlagNonKeyMolACAverage   = .true.
      case ( 4 )
        FlagSort               = .true.
        Method                 = 'LBLConstPF'
        FlagCorrelation        = .true.
        FlagTwoMolAbsOverlap   = .true.
        FlagPFRatioOne         = .false.
        FlagSortPFRatioAverage = .false.
        FlagBinning            = .false.
        FlagBinningDemonstration = .false.
        FlagKeyMolACAverage      = .false.
        FlagNonKeyMolACAverage   = .true.
      case ( 5 )
        FlagSort               = .true.
        Method                 = 'LBLConstPF'
        FlagCorrelation        = .true.
        FlagTwoMolAbsOverlap   = .true.
        FlagPFRatioOne         = .true.
        FlagSortPFRatioAverage = .false.
        FlagBinning            = .false.
        FlagBinningDemonstration = .false.
        FlagKeyMolACAverage      = .false.
        FlagNonKeyMolACAverage   = .true.
      case ( 6 ) ! no effect (meaningless) now
        FlagSort                 = .true.
        Method                   = 'LBLConstPF'
        FlagCorrelation          = .false.
        FlagTwoMolAbsOverlap     = .true.
        FlagPFRatioOne           = .false.
        FlagSortPFRatioAverage   = .true.
        FlagBinning              = .false.
        FlagBinningDemonstration = .false.
        FlagKeyMolACAverage      = .false.
        FlagNonKeyMolACAverage   = .true.
      case ( 7 ) ! the "best" with correlation before binning (2016/12/14)
        FlagSort                 = .true.
        Method                   = 'KDDevelop'
        FlagCorrelation          = .true.
        FlagTwoMolAbsOverlap     = .true.
        FlagPFRatioOne           = .false.
        FlagSortPFRatioAverage   = .false.
        FlagBinning              = .false.
        FlagBinningDemonstration = .false.
        FlagKeyMolACAverage      = .false.
        FlagNonKeyMolACAverage   = .true.
      case ( 8 )
        FlagSort                 = .true.
        Method                   = 'KDDevelop'
        FlagCorrelation          = .true.
        FlagTwoMolAbsOverlap     = .true.
        FlagPFRatioOne           = .false.
        FlagSortPFRatioAverage   = .true.
        FlagBinning              = .true.
        FlagBinningDemonstration = .true.
        FlagKeyMolACAverage      = .false.
        FlagNonKeyMolACAverage   = .true.
      case ( 9 ) ! the "best" with correlation and binning (2016/12/16)
        FlagSort                 = .true.
        Method                   = 'KDDevelop'
        FlagCorrelation          = .true.
        FlagTwoMolAbsOverlap     = .true.
        FlagPFRatioOne           = .false.
        FlagSortPFRatioAverage   = .true.
        FlagBinning              = .true.
        FlagBinningDemonstration = .false.
        FlagKeyMolACAverage      = .false.
        FlagNonKeyMolACAverage   = .true.
      case ( 10 )
        FlagSort               = .true.
        Method                 = 'KDDevelop'
        FlagCorrelation        = .true.
        FlagTwoMolAbsOverlap   = .true.
        FlagPFRatioOne         = .false.
        FlagSortPFRatioAverage = .false.
        FlagBinning            = .false.
        FlagBinningDemonstration = .false.
        FlagKeyMolACAverage      = .false.
        FlagNonKeyMolACAverage   = .true.
      case ( 11 )
        FlagSort               = .true.
        Method                 = 'KDDevelop'
        FlagCorrelation        = .true.
        FlagTwoMolAbsOverlap   = .true.
        FlagPFRatioOne         = .false.
        FlagSortPFRatioAverage = .false.
        FlagBinning            = .true.
        FlagBinningDemonstration = .false.
        FlagKeyMolACAverage      = .false.
        FlagNonKeyMolACAverage   = .true.
      case default
        stop 'Unexpected value of iLoop'
      end select

      write( 6, * ) "KeyMolNum                : ", a_BandCurKeyMolNum(1:NBandCur)
      write( 6, * ) "FlagSort                 : ", FlagSort
      write( 6, * ) "Method                   : ", trim( Method )
      write( 6, * ) "FlagCorrelation          : ", FlagCorrelation
      write( 6, * ) "FlagTwoMolAbsOverlap     : ", FlagTwoMolAbsOverlap
      write( 6, * ) "FlagPFRatioOne           : ", FlagPFRatioOne
      write( 6, * ) "FlagSortPFRatioAverage   : ", FlagSortPFRatioAverage
      write( 6, * ) "FlagBinning              : ", FlagBinning
      write( 6, * ) "FlagBinningDemonstration : ", FlagBinningDemonstration
      write( 6, * ) "FlagKeyMolACAverage      : ", FlagKeyMolACAverage
      write( 6, * ) "FlagNonKeyMolACAverage   : ", FlagNonKeyMolACAverage

      InNcFN = InOrgOptPropNcFN
      write( OutNcFN, '(a,a,i2.2,a,i2.2,a)' ) &
        & OutOptPropNcFN(1:len_trim(OutOptPropNcFN)-3), '-', iband, '-', iloop, '.nc'
      call sortac( &
        & FlagSort, &
        & InNcFN, OutNcFN, &
        & FlagChkConstH2OMR, FlagBinning, FlagBinningDemonstration, &
        & FlagKeyMolACAverage, FlagNonKeyMolACAverage, &
        & FlagPFRatioOne, &
        & FlagSortPFRatioAverage, &
!        & NBand, a_BandWNBnds, a_BandKeyMolNum, a_BandNBin &
        & FlagCorrelation, &
        & FlagTwoMolAbsOverlap, &
        & NBandCur, a_BandCurWNBnds, a_BandCurKeyMolNum, a_BandCurNBin, &
        & NGaussReg, a_GaussRegBnds &
        & )

      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


      write( OutNcFN, '(a,a,i2.2,a,i2.2,a)' ) &
        & OutFluxNcFN(1:len_trim(OutFluxNcFN)-3), '-', iband, '-', iloop, '.nc'
      call CalcTendency( &
        & Grav, Cp, &
        & kmax, r_Press, &
        & r_RadUwFlux, r_RadDwFlux, &
        & z_DTempDt &
        & )
      call OutputFluxTendency( &
        & OutNcFN,    &
        & kmax, r_Press, z_Press, &
        & r_RadUwFlux, r_RadDwFlux, z_DTempDt &
        & )


      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_RadUwFluxBandLBL(k) ) )
        RadDwFluxRatioMax = &
          & max( RadDwFluxRatioMax, &
          &   abs( r_RadDwFlux(k) / ( r_RadDwFluxBandLBL(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( z_DTempDt   )
  deallocate( r_RadUwFlux )
  deallocate( r_RadDwFlux )

  deallocate( a_GaussRegBnds )

  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 )


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 CalcTendency( &
    & Grav, Cp, &
    & kmax, r_Press, &
    & r_RadUwFlux, r_RadDwFlux, &
    & z_DTempDt &
    & )

    real(DP)    , intent(in ) :: Grav
    real(DP)    , intent(in ) :: Cp
    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)
    real(DP)    , intent(out) :: z_DTempDt  (1:kmax)


    ! Local variables
    !
    real(DP) :: r_RadNetFlux(0:kmax)
    integer  :: k


    r_RadNetFlux = r_RadUwFlux - r_RadDwFlux

    do k = 1, kmax
      z_DTempDt(k) = &
        &   ( r_RadNetFlux(k-1) - r_RadNetFlux(k) ) &
        & / ( r_Press     (k-1) - r_Press     (k) ) &
        & * Grav / Cp
    end do


  end subroutine CalcTendency

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

  subroutine OutputFluxTendency( &
    & OutNcFN, &
    & kmax, r_Press, z_Press, &
    & r_RadUwFlux, r_RadDwFlux, z_DTempDt &
    & )

    character(*), intent(in ) :: OutNcFN
    integer     , intent(in ) :: kmax
    real(DP)    , intent(in ) :: r_Press    (0:kmax)
    real(DP)    , intent(in ) :: z_Press    (1:kmax)
    real(DP)    , intent(in ) :: r_RadUwFlux(0:kmax)
    real(DP)    , intent(in ) :: r_RadDwFlux(0:kmax)
    real(DP)    , intent(in ) :: z_DTempDt  (1: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', 'z_Press'/), dimsizes=(/kmax+1, kmax/), &
      & longnames=(/'pressure', 'pressure'/),       &
      & units=(/'Pa', 'Pa'/), &
      & history=GTHistOut &
      & )

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

    call HistoryAddVariable( &
      & varname='r_RadFluxUp', dims=(/'r_Press'/), &
      & longname='upward flux(planet)', units='W m-2', xtype='double', &
      & history=GTHistOut &
      & )
    call HistoryPut('r_RadFluxUp',r_RadUwFlux,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_RadFluxDn',r_RadDwFlux,history=GTHistOut)
    call HistoryAddVariable( &
      & varname='z_DTempDt', dims=(/'z_Press'/), &
      & longname='tendency', units='K s-1', xtype='double', &
      & history=GTHistOut &
      & )
    call HistoryPut('z_DTempDt',z_DTempDt,history=GTHistOut)


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

  end subroutine OutputFluxTendency

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

end program main
