program main

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

  use dc_types

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

  use gridset, only : GridSetInit

  use rad_kd_utils, only : InqDimLength
  use rad_interpolate, only : Interpolate1D

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

  use planetconst, only : Grav, Cp, MeanMolWt, H2OMolWt


  ! 宣言文 ; Declaration statements
  !
  implicit none


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

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

!!$  integer, parameter :: MaxNmlArySize = 1024


  ! 設定変数, 他
  ! setteing parameter, values
  !
  character(128) :: InKDTblNcFN
                          ! 
                          ! input file name for KD table
  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 :: r_MeanMolWt(:)

  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


  character(128) :: Method

  integer :: k


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


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


  namelist /file_nml/ &
    & InProfNcFn, &
!!$    & InNcFN, OutNcFN
    & InKDTblNcFN, &
    & FlagOutSpeFlux, &
    & OutSpeFluxNcFN, &
    & OutFluxNcFN,    &
    & OutOptPropNcFN


  InProfNcFn = ''

  InNcFN  = ''
  OutOptPropNcFN = ''


  ! 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 )
  !
  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                         &
    & )

  allocate( r_MeanMolWt(0:kmax) )
  r_MeanMolWt = MeanMolWt


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


  call GridSetInit( &
    & kmax &
    & )



  Method = 'KD'
  call RadFluxInit( &
    & Grav, H2OMolWt, &
    & InKDTblNcFN, &
    & 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, &
    & r_MeanMolWt, &
    & -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 &
    & )



  deallocate( z_DTempDt   )
  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( 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
