program main

  use vtype_module
  use fi_module

  use constants0

  use rad_planetconst, only : &
    & RadPlanetConstInit, &
    & Grav,            &
    & Cp,              &
    & MolWtNonRadAct

  use gauss_quad, only : GauLeg

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

  use rad_output, only : RadOutputFluxTendency

  use rad_kd_utils, only : RadKDUtilsCalcTendency

  use hitranconst, only : m_HITRANMolWt

  use io_profile, only : &
    & IOProfileInqDimLen, &
    & IOProfileRead, &
    & IOProfileReadPtcl

  ! 宣言文 ; Declaration statements
  !
  implicit none


  character(stdstr) :: CtlFN0 = 'planetconst.conf'
  real(DP)          :: GravSet
  real(DP)          :: CpSet
  real(DP)          :: MolWtNonRadActSet


  character(stdstr) :: CtlFN1 = 'calc_rte.conf'
  character(stdstr) :: CtlFN2 = 'calc_rte_file.conf'

  integer :: FileUnit
  integer :: iostat



  character(extstr*2) :: InOptPropNcFN
  logical             :: FlagOutSpeFlux
  character(extstr*2) :: OutSpeFluxNcFN
  character(extstr*2) :: OutFluxNcFN
  character(extstr*2) :: OutTendencyNcFN

  character(extstr*2) :: InProfNcFn

  integer           :: InNcID
  character(extstr) :: Mode

  integer, allocatable :: m_MolNum(:)
  real(8), allocatable :: r_Press (:)
  real(8), allocatable :: r_Temp  (:)
  real(8), allocatable :: rm_VMR  (:,:)

  integer                            :: NPtcl
  integer              , allocatable :: a_PtclNum    (:)
  character(len=stdstr), allocatable :: a_PtclName   (:)
  real(8)              , allocatable :: a_PtclDens   (:)
  real(8)              , allocatable :: za_PtclEffRad(:,:)
  real(8)              , allocatable :: za_PtclMMR   (:,:)
  real(8)              , allocatable :: za_DelPtclNum(:,:)

  character(len=extstr) :: PtclOptPropNcFN
  real(DP)              :: StrFluxTOA
  character(len=extstr) :: StrSpeNcFN

  real(DP) :: AlbedoSW
  real(DP) :: AlbedoLW
  real(DP) :: WaveNumAlbedoSwitch

  integer               :: NGQ
  real(DP), allocatable :: a_GQMu    (:)
  real(DP), allocatable :: a_GQWeight(:)

  integer               :: NGQLoop
  real(DP)              :: SZADegEach

  real(DP) :: SZADeg

  real(DP) :: WNIntegStart
  real(DP) :: WNIntegEnd

  real(8), allocatable :: r_MeanMolWt(:)
  real(8), allocatable :: r_VMRRadActSum(:)
  real(8), allocatable :: r_VMRNonRadAct(:)

  real(8)              :: SurfTemp


  real(DP), allocatable :: z_Press(:)
  real(DP), allocatable :: z_DTempDt      (:)
                          ! tendency
  real(DP), allocatable :: r_RadUwFlux    (:)
                          ! upward flux
  real(DP), allocatable :: r_RadDwFlux    (:)
                          ! downward flux
  real(DP), allocatable :: ra_DelRadUwFlux(:,:)
  real(DP), allocatable :: ra_DelRadDwFlux(:,:)

  real(DP), allocatable :: r_RadUwFluxLW    (:)
  real(DP), allocatable :: r_RadDwFluxLW    (:)
  real(DP), allocatable :: r_RadUwFluxSW    (:)
  real(DP), allocatable :: r_RadDwFluxSW    (:)
  real(DP), allocatable :: ra_DelRadUwFluxLW(:,:)
  real(DP), allocatable :: ra_DelRadDwFluxLW(:,:)
  real(DP), allocatable :: ra_DelRadUwFluxSW(:,:)
  real(DP), allocatable :: ra_DelRadDwFluxSW(:,:)
  real(DP), allocatable :: z_DTempDtLW      (:)
  real(DP), allocatable :: z_DTempDtSW      (:)

  real(DP), allocatable :: r_RadUwFluxEach      (:)
  real(DP), allocatable :: r_RadDwFluxEach      (:)
  real(DP), allocatable :: ra_DelRadUwFluxEach  (:,:)
  real(DP), allocatable :: ra_DelRadDwFluxEach  (:,:)
  real(DP), allocatable :: r_RadUwFluxLWEach    (:)
  real(DP), allocatable :: r_RadDwFluxLWEach    (:)
  real(DP), allocatable :: r_RadUwFluxSWEach    (:)
  real(DP), allocatable :: r_RadDwFluxSWEach    (:)
  real(DP), allocatable :: ra_DelRadUwFluxLWEach(:,:)
  real(DP), allocatable :: ra_DelRadDwFluxLWEach(:,:)
  real(DP), allocatable :: ra_DelRadUwFluxSWEach(:,:)
  real(DP), allocatable :: ra_DelRadDwFluxSWEach(:,:)


  integer :: kmax
  integer :: NMol


  character(extstr) :: Method


  integer :: k
  integer :: m
  integer :: n


  ! namelist file の設定

  namelist /planetconst_nml/ &
    & GravSet, CpSet, MolWtNonRadActSet

  ! setting of namelist file
  !
  namelist /rte_nml/ &
    & StrFluxTOA, &
    & StrSpeNcFN, &
    & AlbedoSW, &
    & AlbedoLW, &
    & WaveNumAlbedoSwitch, &
    & NGQ, &
    & SZADeg, &
    & WNIntegStart, &
    & WNIntegEnd

  namelist /rte_file_nml/ &
    & Method, &
    & InProfNcFn, &
    & InOptPropNcFN, &
    & PtclOptPropNcFN, &
    & FlagOutSpeFlux,  &
    & OutSpeFluxNcFN,  &
    & OutFluxNcFN,     &
    & OutTendencyNcFN

  ! for Earth
  !
  GravSet           = 9.80665d0
  CpSet             = 1004.6d0
  MolWtNonRadActSet = 28.96d-3

  Method          = 'LBL'
  InProfNcFn      = 'prof.nc'
  InOptPropNcFN   = 'ac.nc'
  PtclOptPropNcFN = ''
  OutFluxNcFN     = 'flux.nc'
  OutTendencyNcFN = 'tendency.nc'
  FlagOutSpeFlux  = .false.
  OutSpeFluxNcFN  = 'flux_spectrum.nc'

!!$  StrFluxTOA          = 1370.0d0
  StrFluxTOA          =    0.0d0
  StrSpeNcFN          = "strspe.nc"
  AlbedoSW            = 0.2d0
  AlbedoLW            = 0.0d0
  WaveNumAlbedoSwitch = 3250.0d2
  NGQ                 =   -1
  SZADeg              =    0.0d0

  WNIntegStart        = -1.0d0
  WNIntegEnd          = -1.0d0




  ! Read a namelist file

  call fi_open( CtlFN0, "read", FileUnit )
  !
  rewind( FileUnit )
  read( FileUnit, nml = planetconst_nml, iostat = iostat )
  if ( iostat > 0 ) then
    stop 'Namelist, planetconst_nml, is something wrong.'
  end if
  write( 6, nml = planetconst_nml )
  !
  close( FileUnit )

  call fi_open( CtlFN1, "read", FileUnit )
  !
  rewind( FileUnit )
  read( FileUnit, nml = rte_nml, iostat = iostat )
  if ( iostat > 0 ) then
    stop 'Namelist in rte_nml is something wrong.'
  end if
  write( 6, nml = rte_nml )
  !
  close( FileUnit )

  call fi_open( CtlFN2, "read", FileUnit )
  !
  rewind( FileUnit )
  read( FileUnit, nml = rte_file_nml, iostat = iostat )
  if ( iostat > 0 ) then
    stop 'Namelist in calc_rte_file_nml is something wrong.'
  end if
  write( 6, nml = rte_file_nml )
  !
  close( FileUnit )



  ! Initialize planetary constants
  !
  call RadPlanetConstInit( &
    & GravSet           = GravSet,           &
    & CpSet             = CpSet,             &
    & MolWtNonRadActSet = MolWtNonRadActSet  &
    & )


  call IOProfileInqDimLen(      &
    & InProfNcFN, "Press", &
    & kmax &
    & )
  kmax = kmax - 1
  call IOProfileInqDimLen(      &
    & InProfNcFN, "MolNum", &
    & NMol &
    & )
  allocate( m_MolNum   (1:NMol) )
  allocate( r_Press    (0:kmax) )
  allocate( r_Temp     (0:kmax) )
  allocate( rm_VMR     (0:kmax,NMol) )
  call IOProfileRead(      &
    & InProfNcFN, "Press", "MolNum", "Temp", "Ts", "VMR", &
    & kmax+1, NMol, &
    & r_Press, m_MolNum, r_Temp, SurfTemp, rm_VMR    &
    & )

!!$  Mode = "read"
!!$  call ni3_open( InProfNcFN, Mode, InNcID )
!!$  call ni3_inq_dimlen( InNcID, "Press" , kmax )
!!$  kmax = kmax - 1
!!$  call ni3_inq_dimlen( InNcID, "MolNum", NMol )
!!$  allocate( m_MolNum   (1:NMol) )
!!$  allocate( r_Press    (0:kmax) )
!!$  allocate( r_Temp     (0:kmax) )
!!$  allocate( rm_VMR     (0:kmax,NMol) )
!!$  call ni3_get_var( InNcID, "MolNum", m_MolNum )
!!$  call ni3_get_var( InNcID, "Press" , r_Press  )
!!$  call ni3_get_var( InNcID, "Temp"  , r_Temp   )
!!$  call ni3_get_var( InNcID, "VMR"   , rm_VMR   )
!!$  call ni3_get_var( InNcID, 'Ts'    , SurfTemp )


!!$  call ni3_inq_dimlen( InNcID, "PtclNum", NPtcl )
  call IOProfileInqDimLen(      &
    & InProfNcFN, "PtclNum", &
    & NPtcl &
    & )
  !
  allocate( z_Press   (1:kmax) )
  allocate( a_PtclNum (NPtcl) )
  allocate( a_PtclName(NPtcl) )
  allocate( a_PtclDens(NPtcl) )
  allocate( za_PtclEffRad(kmax,NPtcl) )
  allocate( za_PtclMMR   (kmax,NPtcl) )
  allocate( za_DelPtclNum(kmax,NPtcl) )
  !
  write( 6, * ) NPtcl
  if ( NPtcl > 0 ) then
!!$    call ni3_get_var( InNcID, 'PtclName'  , a_PtclName )
!!$    call ni3_get_var( InNcID, 'PtclDens'  , a_PtclDens )
!!$    call ni3_get_var( InNcID, 'PtclEffRad', za_PtclEffRad )
!!$    call ni3_get_var( InNcID, 'PtclMMR'   , za_PtclMMR    )
    call IOProfileReadPtcl(      &
      & InProfNcFN, &
      & kmax, NPtcl, &
      & z_Press, a_PtclNum, a_PtclName, a_PtclDens, za_PtclEffRad, za_PtclMMR  &
      & )
    do m = 1, NPtcl
      do k = 1, kmax
        za_DelPtclNum(k,m) = za_PtclMMR(k,m) &
          & * ( r_Press(k-1) - r_Press(k) ) / Grav &
          & / ( 4.0d0 / 3.0d0 * PI * za_PtclEffRad(k,m)**3 * a_PtclDens(m) )
      end do
    end do
  end if
  !
!!$  a_PtclName(1)      = "water"
!!$  a_PtclName(2)      = "ice"
!!$  za_PtclEffRad(:,1) = 10.0d-6
!!$  za_PtclEffRad(:,2) = 50.0d-6
!!$  za_DelPtclNum      = 0.0d0

  ! 10 g m-2 cloud water
!!$  m = 1
!!$  k = 10
!!$  za_DelPtclNum(k,m) = &
!!$    & 10.0d-3 & 
!!$    & / ( 4.0d0 / 3.0d0 * PI * za_PtclEffRad(k,m)**3 * 1.0d3 )


!!$  call ni3_close( InNcID )


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

  allocate( r_VMRRadActSum(0:kmax) )
  allocate( r_VMRNonRadAct(0:kmax) )
  !
  r_VMRRadActSum = 0.0d0
  do m = 1, NMol
    r_VMRRadActSum = r_VMRRadActSum + rm_VMR(:,m)
  end do
  !
  r_VMRNonRadAct = 1.0d0 - r_VMRRadActSum
  r_VMRNonRadAct = max( r_VMRNonRadAct, 0.0d0 )
  !
  r_MeanMolWt = 0.0d0
  r_MeanMolWt = r_MeanMolWt + r_VMRNonRadAct * MolWtNonRadAct
  do m = 1, NMol
    r_MeanMolWt = r_MeanMolWt + rm_VMR(:,m) * m_HITRANMolWt(m_MolNum(m))
  end do
  r_MeanMolWt = r_MeanMolWt / ( r_VMRRadActSum + r_VMRNonRadAct )
  !
  deallocate( r_VMRRadActSum )
  deallocate( r_VMRNonRadAct )
  !
  write( 6, * ) 'Mean molecular weight'
  do k = 0, kmax
    write( 6, * ) k, r_MeanMolWt(k)
  end do


  allocate( z_DTempDt  (1:kmax) )
  allocate( r_RadUwFlux(0:kmax) )
  allocate( r_RadDwFlux(0:kmax) )
  allocate( ra_DelRadUwFlux(0:kmax,-1:1) )
  allocate( ra_DelRadDwFlux(0:kmax,-1:1) )

  allocate( r_RadUwFluxLW    (0:kmax) )
  allocate( r_RadDwFluxLW    (0:kmax) )
  allocate( r_RadUwFluxSW    (0:kmax) )
  allocate( r_RadDwFluxSW    (0:kmax) )
  allocate( ra_DelRadUwFluxLW(0:kmax,-1:1) )
  allocate( ra_DelRadDwFluxLW(0:kmax,-1:1) )
  allocate( ra_DelRadUwFluxSW(0:kmax,-1:1) )
  allocate( ra_DelRadDwFluxSW(0:kmax,-1:1) )

  allocate( z_DTempDtLW(1:kmax) )
  allocate( z_DTempDtSW(1:kmax) )

  do k = 1, kmax
    z_Press(k) = ( r_Press(k-1) + r_Press(k) ) / 2.0d0
  end do




  call RadFluxInit( &
    & Grav, &
    & WaveNumAlbedoSwitch, &
    & InOptPropNcFN, PtclOptPropNcFN, &
    & StrSpeNcFN, &
    & Method &
    & )
  if ( FlagOutSpeFlux ) then
    call RadFluxOutputInit( &
      & InOptPropNcFN, &
      & OutSpeFluxNcFN, &
      & kmax, r_Press &
      & )
  end if
  allocate( r_RadUwFluxEach      (0:kmax) )
  allocate( r_RadDwFluxEach      (0:kmax) )
  allocate( ra_DelRadUwFluxEach  (0:kmax,-1:1) )
  allocate( ra_DelRadDwFluxEach  (0:kmax,-1:1) )
  allocate( r_RadUwFluxLWEach    (0:kmax) )
  allocate( r_RadDwFluxLWEach    (0:kmax) )
  allocate( r_RadUwFluxSWEach    (0:kmax) )
  allocate( r_RadDwFluxSWEach    (0:kmax) )
  allocate( ra_DelRadUwFluxLWEach(0:kmax,-1:1) )
  allocate( ra_DelRadDwFluxLWEach(0:kmax,-1:1) )
  allocate( ra_DelRadUwFluxSWEach(0:kmax,-1:1) )
  allocate( ra_DelRadDwFluxSWEach(0:kmax,-1:1) )
  if ( NGQ > 0 ) then
    NGQLoop = NGQ + 1
    allocate( a_GQMu    ( NGQ+1 ) )
    allocate( a_GQWeight( NGQ+1 ) )
    call GauLeg( 0.0d0, 1.0d0, NGQ, a_GQMu(1:NGQ), a_GQWeight(1:NGQ) )
    ! divid weights by 2 in order to consider night side
    a_GQWeight = a_GQWeight / 2.0d0
    ! night side
    a_GQMu    (NGQ+1) = -1.0d0
    a_GQWeight(NGQ+1) =  0.5d0
  else
    NGQLoop = 1
    allocate( a_GQMu    ( 1 ) )
    allocate( a_GQWeight( 1 ) )
    n = 1
    a_GQMu    (n) = cos( SZADeg * PI / 180.0d0 )
    a_GQWeight(n) = 1.0d0
  end if
  r_RadUwFlux       = 0.0d0
  r_RadDwFlux       = 0.0d0
  ra_DelRadUwFlux   = 0.0d0
  ra_DelRadDwFlux   = 0.0d0
  r_RadUwFluxLW     = 0.0d0
  r_RadDwFluxLW     = 0.0d0
  r_RadUwFluxSW     = 0.0d0
  r_RadDwFluxSW     = 0.0d0
  ra_DelRadUwFluxLW = 0.0d0
  ra_DelRadDwFluxLW = 0.0d0
  ra_DelRadUwFluxSW = 0.0d0
  ra_DelRadDwFluxSW = 0.0d0
  if ( NGQ > 0 ) then
    write( 6, '(i3,a)' ) NGQ, '-point Gaussian quadrature for dayside + nightside:'
  end if
  do n = 1, NGQLoop
    SZADegEach = acos( a_GQMu(n) ) * 180.0d0/PI
    if ( NGQ > 0 ) then
      write( 6, * ) '   ', n, SZADegEach
    end if
    call RadFluxCalcFlux( &
      & kmax, NMol, &
      & m_MolNum, &
      & r_Press, &
      & SurfTemp, r_Temp, rm_VMR, &
      & NPtcl, a_PtclName, za_PtclEffRad, za_DelPtclNum, &
      & AlbedoSW, AlbedoLW, SZADegEach, &
      & StrFluxTOA, &
      & r_MeanMolWt, &
      & r_RadUwFluxEach, r_RadDwFluxEach, &
      & ra_DelRadUwFluxEach, ra_DelRadDwFluxEach, &
      & r_RadUwFluxLWEach, r_RadDwFluxLWEach, &
      & r_RadUwFluxSWEach, r_RadDwFluxSWEach, &
      & ra_DelRadUwFluxLWEach, ra_DelRadDwFluxLWEach, &
      & ra_DelRadUwFluxSWEach, ra_DelRadDwFluxSWEach, &
      & WNIntegStart, WNIntegEnd &
      & )
    r_RadUwFlux       = r_RadUwFlux       + r_RadUwFluxEach       * a_GQWeight(n)
    r_RadDwFlux       = r_RadDwFlux       + r_RadDwFluxEach       * a_GQWeight(n)
    ra_DelRadUwFlux   = ra_DelRadUwFlux   + ra_DelRadUwFluxEach   * a_GQWeight(n)
    ra_DelRadDwFlux   = ra_DelRadDwFlux   + ra_DelRadDwFluxEach   * a_GQWeight(n)
    r_RadUwFluxLW     = r_RadUwFluxLW     + r_RadUwFluxLWEach     * a_GQWeight(n)
    r_RadDwFluxLW     = r_RadDwFluxLW     + r_RadDwFluxLWEach     * a_GQWeight(n)
    r_RadUwFluxSW     = r_RadUwFluxSW     + r_RadUwFluxSWEach     * a_GQWeight(n)
    r_RadDwFluxSW     = r_RadDwFluxSW     + r_RadDwFluxSWEach     * a_GQWeight(n)
    ra_DelRadUwFluxLW = ra_DelRadUwFluxLW + ra_DelRadUwFluxLWEach * a_GQWeight(n)
    ra_DelRadDwFluxLW = ra_DelRadDwFluxLW + ra_DelRadDwFluxLWEach * a_GQWeight(n)
    ra_DelRadUwFluxSW = ra_DelRadUwFluxSW + ra_DelRadUwFluxSWEach * a_GQWeight(n)
    ra_DelRadDwFluxSW = ra_DelRadDwFluxSW + ra_DelRadDwFluxSWEach * a_GQWeight(n)
  end do
  deallocate( a_GQMu     )
  deallocate( a_GQWeight )
  if ( FlagOutSpeFlux ) then
    call RadFluxOutputFinalize
  end if
  call RadFluxFinalize

  call RadKDUtilsCalcTendency( &
    & Grav, Cp, &
    & kmax, r_Press, &
    & r_RadUwFlux, r_RadDwFlux, z_DTempDt &
    & )
  call RadKDUtilsCalcTendency( &
    & Grav, Cp, &
    & kmax, r_Press, &
    & r_RadUwFluxLW, r_RadDwFluxLW, z_DTempDtLW &
    & )
  call RadKDUtilsCalcTendency( &
    & Grav, Cp, &
    & kmax, r_Press, &
    & r_RadUwFluxSW, r_RadDwFluxSW, z_DTempDtSW &
    & )
  call RadOutputFluxTendency( &
    & OutFluxNcFN, OutTendencyNcFN,   &
    & kmax, r_Press, z_Press, &
    & r_RadUwFlux, r_RadDwFlux, z_DTempDt, &
    & r_RadUwFluxLW, r_RadDwFluxLW, z_DTempDtLW, &
    & r_RadUwFluxSW, r_RadDwFluxSW, z_DTempDtSW &
    & )




  deallocate( z_DTempDt   )
  deallocate( r_RadUwFlux )
  deallocate( r_RadDwFlux )
  deallocate( ra_DelRadUwFlux )
  deallocate( ra_DelRadDwFlux )

  deallocate( r_RadUwFluxLW     )
  deallocate( r_RadDwFluxLW     )
  deallocate( r_RadUwFluxSW     )
  deallocate( r_RadDwFluxSW     )
  deallocate( ra_DelRadUwFluxLW )
  deallocate( ra_DelRadDwFluxLW )
  deallocate( ra_DelRadUwFluxSW )
  deallocate( ra_DelRadDwFluxSW )

  deallocate( z_DTempDtLW )
  deallocate( z_DTempDtSW )

  deallocate( r_RadUwFluxEach       )
  deallocate( r_RadDwFluxEach       )
  deallocate( ra_DelRadUwFluxEach   )
  deallocate( ra_DelRadDwFluxEach   )
  deallocate( r_RadUwFluxLWEach     )
  deallocate( r_RadDwFluxLWEach     )
  deallocate( r_RadUwFluxSWEach     )
  deallocate( r_RadDwFluxSWEach     )
  deallocate( ra_DelRadUwFluxLWEach )
  deallocate( ra_DelRadDwFluxLWEach )
  deallocate( ra_DelRadUwFluxSWEach )
  deallocate( ra_DelRadDwFluxSWEach )


  deallocate( m_MolNum    )
  deallocate( r_Press     )
  deallocate( r_Temp      )
  deallocate( rm_VMR      )

  deallocate( a_PtclName )
  deallocate( a_PtclDens )
  deallocate( za_PtclEffRad )
  deallocate( za_DelPtclNum )

  write( 6, * ) "Finish"

end program main
