!= dcrtm flux 計算 主プログラム
!
!= dcrtm flux calculation main program
!
! Authors::   Masanori Onishi
! Version::   $Id: main_dcrtm.f90,v 1.00 onishi$
! Tag Name::  $Name: dcrtm-20170201 $
! Copyright:: Copyright (C) GFD Dennou Club, 2008-2016. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]

program main_flux
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !

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

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

  ! 格子点設定
  ! grid point setting
  !
  use gridset, only: GridsetSet,  &
    &                GridsetSet2, &
    &                imax,        &
    &                jmax,        &
    &                kmax,        &
    &                imax2,       &
    &                jmax2,       &
    &                kmax2

  ! 2方向近似放射モデル
  ! two stream approximation (Toon et al. (1989))
  !
  use rad_rte_two_stream_app_sf,                & ! 2017/02/08
    &          only: RadRTETwoStreamAppInit_sf, &
    &                RadRTETwoStreamAppLW_sf,   &
    &                RadRTETwoStreamAppSW

  ! 2方向近似放射モデル
  ! two stream approximation for planetary radiation using Two-Stream technique
  !
  use rad_rte_two_stream_app_ts,                & ! 2017/02/01
    &          only: RadRTETwoStreamAppInit_TS, & ! 2017/02/01
    &                RadRTETwoStreamAppLW_TS      ! 2017/02/01

  ! 2方向近似放射モデル
  ! two stream approximation using Source Function Iteration Method (planetary radiation)
  !
  use rad_rte_two_stream_app_itg2,                & ! 2016/12/27
    &          only: RadRTETwoStreamAppInit_itg2, & ! 2016/12/27
    &                RadRTETwoStreamAppLW_itg2      ! 2016/12/27

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

  ! 物理定数設定2
  ! physical constants setteing 2
  !
  use planetconst, only: Grav, &
    &                    AvogadroNum

  ! プランク関数計算
  ! planck function calculation
  !
  use planckf, only: integral_planck

  ! 吸収散乱係数から光学的厚さ, 一次散乱アルベドへの変換
  ! convert absorption and scattering coefficient
  ! to optical depth & single scattering albedo
  !
  use convert, only: KtoTauSSA


  use interpol, only: func_BoundGridNum, &
    &                 ArrInterPol

  use mod_OptDepItr, only: DivLayCount


  ! 宣言文 ; Declaration statements
  !
  implicit none

  ! 設定変数, 他
  ! setteing parameter, values
  !
  integer :: Flag_InputFormat = 1 !(default)
                          ! 入力物理量の指定フラグ
                          ! flag of input value 
                          ! 1: optical depth and single scattering albedo
                          ! 2: absorption, scattering coefficient
  character(99) :: InFileOpt  = 'OptDep_dcrtm.nc' !(default)
                          ! 入力 ファイル名
                          ! input file name
  character(99) :: InFileProf
                          ! 入力プロファイル(Press, Temp, MixRatio, MolWt)
                          ! file name of input profile
  character(99) :: OutFileFlux = 'Flux_dcrtm.nc' !(default)
                          ! 出力ファイル名
                          ! output file name
  real(8) :: SolTemp           = 5800.0_8 !(default)
                          ! 中心星表面温度 [K]
                          ! surface temperature of central star [K]
  real(8) :: SolConst          = 1366.0_8 !(default)
                          ! 太陽定数 [W m-2]
                          ! solar constant [W m-2]
  real(8) :: InAngle           = 2.0_8    !(default)
                          ! sec (太陽光入射角) [1]
                          ! sec (angle of incidence) [1]
  real(8) :: SurfAlbedo        = 0.0_8    !(default)
                          ! 地表面アルベド [1]
                          ! surface albedo [1]
  integer :: Flag_Emit         = 1
                          ! 射出大気の計算フラグ
                          ! flag of emitting atmosphere
                          ! 0: no emitting atmosphere
                          ! 1: emitting atmosphere
  integer :: Flag_Particle     = 0
                          ! 粒子計算フラグ
                          ! flag of particle contained atmosphere
                          ! 0: no particles
                          ! 1: particles are contained
  integer :: Flag_FluxModel    = 1
                          ! フラックス計算手法のフラグ
                          ! flag of radiative transfer model
                          ! 1: Toon et al. (1989)
                          !   for solar radiation: Two-Stream Equation (Eddington App.)
                          !   for planetary radiation: Two-Stream Source Function Method (Hemispheric Mean)
                          ! 2: Toon et al. (1989) + Source Function Iteration Method
                          !   for solar radiation: Two-Stream Equation (Eddington App.)
                          !   for planetary radiation: Two-Stream Source Function Iteration Method (Hemispheric Mean)
  character(99) :: InFileParticle  = 'OptDep_Particle.nc' !(default)
                          ! 粒子の光学特性入力ファイル
                          ! input file of optical properties of particles
  real(8), allocatable :: r_Press(:)
                          ! 気圧 (半整数レベル) [Pa]
                          ! pressure (half level) [Pa]
  real(8), allocatable :: r_Temp(:)
                          ! 気温 (半整数レベル) [K]
                          ! temperature (half level) [K]
  real(8), allocatable :: rm_MixRatio(:,:)
                          ! 体積混合比 (半整数レベル) [1]
                          ! volume mixing ratio (half level) [1]
  real(8), allocatable :: m_MolWt(:)
                          ! 分子量 [kg mol-1]
                          ! molecular weight [kg mol-1]
  real(8) :: MolWt0                 = 28.96d-3 !(default)
                          ! 非吸収分子の平均分子量 [kg mol-1]
                          ! molecular weight of non-absorption gas [kg mol-1]
  real(8), allocatable :: w_WaveNum(:)
                          ! 波数 [m-1]
                          ! wavenumber [m-1]
  real(8), allocatable :: rw_OptDepTOA(:,:)
                          ! 大気上端からの光学的厚さ (半整数レベル) [1]
                          ! optical depth from top of atmosphere (half level) [1]
  real(8), allocatable :: zw_SingleScatA(:,:)
                          ! 1次散乱アルベド (full level) [1]
                          ! single scattering albedo (full level) [1]
  real(8), allocatable :: rmw_LAbsorpK(:,:,:)
                          ! 線吸収の吸収係数 [m2 kg-1]
                          ! absorption coefficient for line absorption [m2 kg-1]
  real(8), allocatable :: rmw_CAbsorpK(:,:,:)
                          ! 連続吸収の吸収係数 [m2 kg-1]
                          ! absorption coefficient for continuum absotption [m2 kg-1]
  real(8), allocatable :: rmw_RayScatK(:,:,:)
                          ! Rayleigh 散乱の吸収係数 [m2 kg-1]
                          ! scattering coefficient for Rayleigh scattering [m2 kg-1]

  ! 計算される物理量
  ! calculated values
  !
  real(8), allocatable :: rw_RadUwFlux(:,:)
                          ! 惑星大気の射出に由来する上向きフラックス [W m-2]
                          ! upward flux originated from emittance of planetary atmospheire [W m-2]
  real(8), allocatable :: rw_RadDwFlux(:,:)
                          ! 惑星大気の射出に由来する下向きフラックス [W m-2]
                          ! downward flux originated from emittance of planetary atmospheire [W m-2]
  real(8), allocatable :: rw_SolUwFlux(:,:)
                          ! 中心星入射に由来する上向きフラックス [W m-2]
                          ! upward flux originated from emittance of central star [W m-2]
  real(8), allocatable :: rw_SolDwFlux(:,:)
                          ! 中心星入射に由来する下向きフラックス [W m-2]
                          ! downward flux originated from emittance of central star [W m-2]

  integer :: fnml = 12
                          ! 装置番号(name list file)
                          ! device number : namelist file
  integer :: wmax
                          ! 波数レベル最大値
                          ! wavenumber level maximum
  integer :: mmax
                          ! 分子数最大値
                          ! molecular number maxim
  integer :: k_w, k_r, k_m
                          ! 作業変数
                          ! work variables
  real(8) :: RatioSolarConst
                          ! 中心星の黒体放射に対する太陽定数の割合
                          ! ratio of solar constant to black body radiation
  real(8), pointer :: a_Axis(:) => null()
                          ! 配列のサイズを調べるための配列
                          ! dummy array
  real(8) :: MaxMixRatio
                          ! 最大混合比
                          ! maximum volume mixing ratio
  integer :: vmax
  integer :: bmax
  integer :: pmax
  integer :: k_v, k_p
  real(8), allocatable :: v_WaveNum      (:)
  real(8), allocatable :: zpv_DelOptDep  (:,:,:)
  real(8), allocatable :: zpv_SingleScatA(:,:,:)
  real(8), allocatable :: zpv_AsymFact   (:,:,:)
  real(8), allocatable :: zv_DelOptDepP  (:,:)
  real(8), allocatable :: zv_SsaTauP     (:,:)
  real(8), allocatable :: zv_SsaTauAFP   (:,:)
  real(8), allocatable :: zw_DelOptDepP  (:,:) ! Particle
  real(8), allocatable :: zw_SsaTauP     (:,:)
  real(8), allocatable :: zw_SsaTauAFP   (:,:)
  real(8), allocatable :: zw_DelOptDepG  (:,:) ! Gas (absorption + rayleigh scat.)
  real(8), allocatable :: zw_DelOptDepR  (:,:) ! Rayleigh scat.
  real(8), allocatable :: zw_AsymFact    (:,:)

  integer :: NLayMax
  integer, allocatable :: z_NLay         (:)
  integer, allocatable :: r_NLaySum      (:)
  real(8), allocatable :: xyb_SSA        (:,:,:)
  real(8), allocatable :: xyb_AF         (:,:,:)
  real(8), allocatable :: xyb_OptDep     (:,:,:)
  real(8), allocatable :: xyb_PFInted    (:,:,:)
  real(8), allocatable :: xyb_RadUwFlux  (:,:,:)
  real(8), allocatable :: xyb_RadDwFlux  (:,:,:)
  real(8), allocatable :: xyba_DelRadUwFlux (:,:,:,:)
  real(8), allocatable :: xyba_DelRadDwFlux (:,:,:,:)

  integer :: k_i, k_j, k_b
  real(8) :: eps_divlay = 1d-3

  ! for dcpam5 radiation code
  !
  !integer :: imax
  !integer :: jmax
  !integer :: kmax
  real(8), allocatable :: xyz_SSA          (:, :, :) 
                          ! single scattering albedo
  real(8), allocatable :: xyz_AF           (:, :, :) 
                          ! asynmetory factor
  real(8), allocatable :: xyr_OptDep       (:, :, :) 
                          ! optical depth from top of atmosphere
  real(8), allocatable :: xy_SurfAlbedo    (:, :) 
                          ! surface albedo
  real(8), allocatable :: xyr_PFInted      (:, :, :) 
                          ! planck function integral
  real(8), allocatable :: xy_SurfPFInted   (:, :) 
                          ! planck function in land surface
  real(8), allocatable :: xy_SurfDPFDTInted(:, :) 
                          ! planck function の温度微分
  real(8), allocatable :: xyr_RadUwFlux    (:, :, :) 
                          ! upward flux
  real(8), allocatable :: xyr_RadDwFlux    (:, :, :) 
                          ! downward flux
  real(8), allocatable :: xyra_DelRadUwFlux(:, :, :, :) 
                          ! 
  real(8), allocatable :: xyra_DelRadDwFlux(:, :, :, :)
                          !
  real(8)              :: SolarFluxTOA
                          ! solar flux of top of atmosphere
  real(8), allocatable :: xy_InAngle       (:, :)
                          ! sec (angle of incidence)


  ! namelist file の設定
  ! setting of namelist file
  !
  namelist /flux_nml/ &
    & Flag_InputFormat,&
    & InFileOpt,       &
    & InFileProf,      &
    & MolWt0,          &
    & OutFileFlux,     &
    & SurfAlbedo,      & 
    & SolTemp,         & 
    & SolConst,        & 
    & InAngle,         &
    & Flag_FluxModel,  &
    & Flag_Emit,       &
    & Flag_Particle,   &
    & InFileParticle


  ! namelist file の読み込み
  ! read namelist file
  !
  open(fnml, file='dcrtm.nml')
  read(fnml, nml=flux_nml)
  close(fnml)

  ! 設定値のチェック
  ! check setting values
  !
  !print *, '==== ==== ==== ===='
  print *, '== Setting Values =='
  print *, 'Flag_InputFormat:', Flag_InputFormat
    if(Flag_InputFormat == 1) then
      print *, '* flux is calculated by Optical Depth and Single Scattering Albedo'
      print *, 'InFileOpt:', InFileOpt
    elseif(Flag_InputFormat == 2) then
      print *, '* flux is calculated by Absorption and Scattering Coefficient'
      print *, 'InFileOpt:', InFileOpt
      print *, 'InFileProf:', InFileProf
      print *, 'MolWt0:', MolWt0
    else
      print *, '* invalid value (Flag_InputFormat = 1 or 2): stop'
      stop
    end if

  print *, 'OutFileFlux:', OutFileFlux
  print *, 'SurfAlbedo:', SurfAlbedo
    if(SurfAlbedo < 0.0_8) then
      print *, '* invalid value (SurfAlbedo >= 0.0): stop'
      stop
    end if

  print *, 'SolTemp:', SolTemp
    if(SolTemp < 0.0_8) then
      print *, '* invalid value (SolTemp >= 0.0): stop'
      stop
    end if

  print *, 'SolConst:', SolConst
    if(SolConst < 0.0_8) then
      print *, '* invalid value (SolConst >= 0.0): stop'
      stop
    end if

  print *, 'InAngle:', InAngle
    if(InAngle >= 0.0) then
      print *, '* day side calculation'
    else
      print *, '* night side calculation (solar radiation calculation is skip)'
    end if

  print *, 'Flag_FluxModel:', Flag_FluxModel
    if    (Flag_FluxModel == 1) then
      print *, '* Toon et al. (1989)'
      print *, '  * for solar radiation: Two-Stream Equation (Eddington App.)'
      print *, '  * for planetary radiation: Two-Stream Source Function Method (Hemispheric Mean)'
    elseif(Flag_FluxModel == 2) then
      print *, '* Toon et al. (1989), Two-Stream Equation (for solar and planetary radiation)'
      print *, '  * for solar radiation: Two-Stream Equation (Eddington App.)'
      print *, '  * for planetary radiation: Two-Stream Equation (Hemispheric Mean)'
    elseif(Flag_FluxModel == 3) then
      print *, '* Toon et al. (1989) + Source Function Iteration Method'
      print *, '  * for solar radiation: Two-Stream Equation (Eddington App.)'
      print *, '  * for planetary radiation: Two-Stream Source Function Iteration Method (Hemispheric Mean)'
    else
      print *, '* invalid value (Flag_FluxModel = 1, 2 or 3): stop'
      stop
    end if

  print *, 'Flag_Emit:', Flag_Emit
    if    (Flag_Emit == 0) then
      print *, '* No emitting atmosphere'
    elseif(Flag_Emit == 1) then
      print *, '* emitting atmospehere'
    else
      print *, '* invalid value (Flag_Emit = 0 or 1): stop'
      stop
    end if

  print *, 'Flag_Particle:', Flag_Particle
    if    (Flag_Particle == 0) then
      print *, '* No particle atmosphere'
    elseif(Flag_Particle == 1) then
      print *, '* particles are contained'
      print *, 'InFileParticle:', InFileParticle

    else
      print *, '* invalid value (Flag_Particle = 0 or 1): stop'
      stop
    end if

  !print *, '====E====N====D===='
  print *, '=setting check end='

  ! 配列サイズを調べる
  ! check array size
  !
  call HistoryGetPointer(''//InFileOpt//'', 'r_Press', a_Axis )
  kmax = size( a_Axis ) - 1
  Nullify( a_Axis )
  !print *, kmax

  call HistoryGetPointer(''//InFileOpt//'', 'w_WaveNum', a_Axis )
  wmax = size( a_Axis ) - 1
  !ResWN = a_Axis(2) - a_Axis(1)

  Nullify( a_Axis )

  !print *, wmax
  !print *, ResWN
  !wmax = abs((MaxWN - MinWN)/ResWN)


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


  call RadRTETwoStreamAppInit_sf

  if( Flag_FluxModel == 2) then
    call RadRTETwoStreamAppInit_TS
  elseif( Flag_FluxModel == 3) then
    !call RadRTETwoStreamAppInit_it
    call RadRTETwoStreamAppInit_itg2
  end if

  !== set1: grid
  imax = 1
  jmax = 1
  !kmax = 32

  call GridsetSet(imax, jmax, kmax)


  allocate (xyz_SSA          (0:imax-1, 1:jmax, 1:kmax) )
  allocate (xyz_AF           (0:imax-1, 1:jmax, 1:kmax) )
  allocate (xyr_OptDep       (0:imax-1, 1:jmax, 0:kmax) )
  allocate (xy_SurfAlbedo    (0:imax-1, 1:jmax)         )
  allocate (xyr_PFInted      (0:imax-1, 1:jmax, 0:kmax) )
  allocate (xy_SurfPFInted   (0:imax-1, 1:jmax)         )
  allocate (xy_SurfDPFDTInted(0:imax-1, 1:jmax)         )
  allocate (xyr_RadUwFlux    (0:imax-1, 1:jmax, 0:kmax) )
  allocate (xyr_RadDwFlux    (0:imax-1, 1:jmax, 0:kmax) )
  allocate (xyra_DelRadUwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1) )
  allocate (xyra_DelRadDwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1) )
  allocate (xy_InAngle       (0:imax-1, 1:jmax)         ) 


  allocate( r_Press       (0:kmax        ) )
  allocate( r_Temp        (0:kmax        ) )
  allocate( w_WaveNum     (        0:wmax) )
  allocate( rw_RadUwFlux  (0:kmax, 0:wmax) )
  allocate( rw_RadDwFlux  (0:kmax, 0:wmax) )
  allocate( rw_SolUwFlux  (0:kmax, 0:wmax) )
  allocate( rw_SolDwFlux  (0:kmax, 0:wmax) )
  allocate( rw_OptDepTOA  (0:kmax, 0:wmax) )
  allocate( zw_SingleScatA(1:kmax, 0:wmax) )
  allocate( zw_AsymFact   (1:kmax, 0:wmax) )
  allocate( z_NLay        (1:kmax        ) )
  allocate( r_NLaySum     (0:kmax        ) )

  ! 鉛直プロファイルの読み込み
  ! read profile data
  !
  if(Flag_InputFormat == 1) then

    call HistoryGet(''//InFileOpt//'', 'r_Press', r_Press)
    call HistoryGet(''//InFileOpt//'', 'w_WaveNum', w_WaveNum)
    call HistoryGet(''//InFileOpt//'', 'r_Temp', r_Temp)
    call HistoryGet(''//InFileOpt//'', 'rw_OptDep', rw_OptDepTOA)
    call HistoryGet(''//InFileOpt//'', 'zw_SingleScatA', zw_SingleScatA)

  elseif(Flag_InputFormat == 2) then

    ! 配列サイズを調べる
    ! check array size
    !
    call HistoryGetPointer(''//InFileOpt//'', 'm_MolWt', a_Axis)
    mmax = size( a_Axis ) - 1
    Nullify( a_Axis )
    !print *, mmax

    allocate( m_MolWt       (        0:mmax        ) )
    allocate( rm_MixRatio   (0:kmax, 0:mmax        ) )
    allocate( rmw_LAbsorpK  (0:kmax, 0:mmax, 0:wmax) )
    allocate( rmw_CAbsorpK  (0:kmax, 0:mmax, 0:wmax) )
    allocate( rmw_RayScatK  (0:kmax, 0:mmax, 0:wmax) )

    rm_MixRatio(:,:) = 0.0_8

    ! 鉛直プロファイルの読み込み
    ! read profile data
    !
    call HistoryGet(''//InFileProf//'', 'r_Press'    , r_Press              )
    call HistoryGet(''//InFileProf//'', 'r_Temp'     , r_Temp               )
    call HistoryGet(''//InFileProf//'', 'm_MolWt'    , m_MolWt    (  1:mmax))
    call HistoryGet(''//InFileProf//'', 'rm_MixRatio', rm_MixRatio(:,1:mmax))
    !call HistoryGet(''//InFileOpt//'', 'r_Press', r_Press)
    !call HistoryGet(''//InFileOpt//'', 'r_Temp', r_Temp)
    !call HistoryGet(''//InFileOpt//'', 'm_MolWt', m_MolWt)
    !call HistoryGet(''//InFileOpt//'', 'rm_MixRatio', rm_MixRatio)
    call HistoryGet(''//InFileOpt//'', 'w_WaveNum', w_WaveNum)
    call HistoryGet(''//InFileOpt//'', 'rmw_LAbsorpK', rmw_LAbsorpK)
    call HistoryGet(''//InFileOpt//'', 'rmw_CAbsorpK', rmw_CAbsorpK)
    call HistoryGet(''//InFileOpt//'', 'rmw_RayScatK', rmw_RayScatK)


    ! 非吸収分子の体積混合比の計算
    ! calculate volume mixing ratio of non-absrption gas
    !
    m_MolWt(0) = MolWt0

    rm_MixRatio(:,0) = 1.0_8
    do k_r = 0, kmax

      do k_m = 1, mmax
        rm_MixRatio(k_r,0) = rm_MixRatio(k_r,0) - rm_Mixratio(k_r,k_m)
      end do

      if(rm_MixRatio(k_r,0) < 0.0_8) then
        print *, 'sum of vlume mixing ratio over 1.0:', 1.0_8 - rm_MixRatio(k_r,0), k_r

        MaxMixRatio = rm_MixRatio(k_r,0)
        do k_m = 1, mmax
          if(rm_MixRatio(k_r,k_m) > MaxMixRatio) then
            MaxMixRatio = rm_MixRatio(k_r,k_m)
          end if
        end do

        do k_m = 1, mmax
          if(rm_MixRatio(k_r,k_m) == MaxMixRatio) then
            exit
          end if
        end do
        rm_MixRatio(k_r,k_m) = rm_MixRatio(k_r,k_m) + rm_MixRatio(k_r,0)
        rm_MixRatio(k_r,0) = 0.0_8

        print *, 'modify volume mixing ratio, k_r=', k_r
        do k_m = 1, mmax
          print *, k_m, rm_MixRatio(k_r,k_m)
        end do
      end if

    end do ! k_r

    call KtoTauSSA (                                          & ! 2016/08/05
      & kmax, mmax, wmax,                                     &
      & r_Press, r_Temp, rm_MixRatio(0:kmax,0:mmax), m_MolWt, &
      & rmw_RayScatK, rmw_LAbsorpK, rmw_CAbsorpK,             &
      & rw_OptDepTOA, zw_SingleScatA                          )

  else

    print *, 'invalid value (Flag_InputFormat):', Flag_InputFormat
    print *, 'stop'
    stop
 
  end if


  zw_AsymFact = 0.0_8

  ! Particle opacities
  !
  if(Flag_Particle == 1) then

    call HistoryGetPointer(''//InFileParticle//'', 'r_Press', a_Axis )
    bmax = size( a_Axis ) - 1
    Nullify( a_Axis )
    !print *, bmax

    if(bmax .ne. kmax) then
      print *, 'vertical dimension is invalid'
      print *, 'gas:', kmax
      print *, 'particles', bmax
      stop
    end if

    call HistoryGetPointer(''//InFileParticle//'', 'w_WaveNum', a_Axis )
    vmax = size( a_Axis ) - 1
    Nullify( a_Axis )
    !print *, vmax

    call HistoryGetPointer(''//InFileParticle//'', 'p_ParNum', a_Axis )
    pmax = size( a_Axis )
    Nullify( a_Axis )
    !print *, pmax

    !allocate ( b_Press        (0:bmax                ) )
    allocate ( v_WaveNum      (                0:vmax) )
    allocate ( zpv_DelOptDep  (1:kmax, 1:pmax, 0:vmax) )
    allocate ( zpv_SingleScatA(1:kmax, 1:pmax, 0:vmax) )
    allocate ( zpv_AsymFact   (1:kmax, 1:pmax, 0:vmax) )
    allocate ( zv_DelOptDepP  (1:kmax,         0:vmax) ) 
    allocate ( zv_SsaTauP     (1:kmax,         0:vmax) ) 
    allocate ( zv_SsaTauAFP   (1:kmax,         0:vmax) )
    allocate ( zw_DelOptDepP  (1:kmax,         0:wmax) ) 
    allocate ( zw_SsaTauP     (1:kmax,         0:wmax) ) 
    allocate ( zw_SsaTauAFP   (1:kmax,         0:wmax) )
    allocate ( zw_DelOptDepG  (1:kmax,         0:wmax) )
    allocate ( zw_DelOptDepR  (1:kmax,         0:wmax) )

    !call HistoryGet(''//InFileParticle//'', 'r_Press', b_Press)
    call HistoryGet(''//InFileParticle//'', 'w_WaveNum', v_WaveNum)
    call HistoryGet(''//InFileParticle//'', 'zpw_DelOptDep', zpv_DelOptDep)
    call HistoryGet(''//InFileParticle//'', 'zpw_SingleScatA', zpv_SingleScatA)
    call HistoryGet(''//InFileParticle//'', 'zpw_AsymFact', zpv_AsymFact)

    zv_DelOptDepP = 0.0_8
    zv_SsaTauP = 0.0_8
    zv_SsaTauAFP = 0.0_8

    do k_r = 1, kmax
      do k_v = 1, vmax
        do k_p = 1, pmax
          zv_DelOptDepP(k_r,k_v)= zv_DelOptDepP(k_r,k_v) + zpv_DelOptDep(k_r,k_p,k_v)
          zv_SsaTauP(k_r,k_v)   = zv_SsaTauP(k_r,k_v) + &
            &                 zpv_SingleScatA(k_r,k_p,k_v)*zpv_DelOptDep(k_r,k_p,k_v)
          zv_SsaTauAFP(k_r,k_v) = zv_SsaTauAFP(k_r,k_v) + &
            &                 zpv_SingleScatA(k_r,k_p,k_v)*zpv_DelOptDep(k_r,k_p,k_v)*zpv_AsymFact(k_r,k_p,k_v)
        end do
      end do
    end do

    do k_r = 1, kmax

      call ArrInterPol(                                 & ! 2016/11/10
        & 1,                                            & ! (in )
        & vmax+1, v_WaveNum, zv_DelOptDepP(k_r,0:vmax), & ! (in )
        & wmax+1, w_WaveNum,                            & ! (in )
        &                    zw_DelOptDepP(k_r,0:wmax)  ) ! (out)

      call ArrInterPol(                              & ! 2016/11/10
        & 1,                                         & ! (in )
        & vmax+1, v_WaveNum, zv_SsaTauP(k_r,0:vmax), & ! (in )
        & wmax+1, w_WaveNum,                         & ! (in )
        &                    zw_SsaTauP(k_r,0:wmax)  ) ! (out)

      call ArrInterPol(                                & ! 2016/11/10
        & 1,                                           & ! (in )
        & vmax+1, v_WaveNum, zv_SsaTauAFP(k_r,0:vmax), & ! (in )
        & wmax+1, w_WaveNum,                           & ! (in )
        &                    zw_SsaTauAFP(k_r,0:wmax)  ) ! (out)

    end do

    do k_r = 1, kmax
      zw_DelOptDepG(k_r,:) = rw_OptDepTOA(k_r-1,:) - rw_OptDepTOA(k_r,:)
      zw_DelOptDepR(k_r,:) = zw_DelOptDepG(k_r,:)* zw_SingleScatA(k_r,:)
    end do

    rw_OptDepTOA = 0.0_8
    do k_r = kmax-1, 0, -1
      rw_OptDepTOA(k_r,:) = rw_OptDepTOA(k_r+1,:) + &
        &                   zw_DelOptDepG(k_r+1,:) + zw_DelOptDepP(k_r+1,:)
    end do

    zw_SingleScatA = 0.0_8
    do k_r = 1, kmax

      do k_w = 0, wmax

        if( (zw_DelOptDepG(k_r,k_w) + zw_DelOptDepP(k_r,k_w)) .ne. 0.0_8 ) then
          zw_SingleScatA(k_r,k_w) = (zw_DelOptDepR(k_r,k_w) + zw_SsaTauP   (k_r,k_w))/&
            &                       (zw_DelOptDepG(k_r,k_w) + zw_DelOptDepP(k_r,k_w))
        else
          zw_SingleScatA(k_r,k_w) = 0.0_8
        end if

        if( ( zw_DelOptDepR(k_r,k_w) + zw_SsaTauP(k_r,k_w)).ne. 0.0_8 ) then
          zw_AsymFact(k_r,k_w) = zw_SsaTauAFP(k_r,k_w)/ &
            &                ( zw_DelOptDepR(k_r,k_w) + zw_SsaTauP(k_r,k_w))
        else
          zw_AsymFact(k_r,k_w) = 0.0_8
        end if

      end do ! k_w in k_r
    end do ! k_r

  end if


  ! フラックスの計算
  ! calculate flux
  !

  ! dcpam5
  !xyz_SSA          (:, :, :) 
  !xyz_AF           (:, :, :) = 0.0_8
  !xyr_OptDep       (:, :, :) 
  xy_SurfAlbedo    (:, :)    = SurfAlbedo
  !xyr_PFInted      (:, :, :)
  !xy_SurfPFInted   (:, :)
  xy_SurfDPFDTInted(:, :)    = 0.0_8
  xy_InAngle       (:, :)    = InAngle


  ! dcpam5 radiation code
  !
  do k_w = 0, wmax
    if( w_WaveNum(k_w) == 0.0_8 ) then
      rw_SolUwFlux(0:kmax,k_w) = 0.0_8
      rw_SolDwFlux(0:kmax,k_w) = 0.0_8
      rw_RadUwFlux(0:kmax,k_w) = 0.0_8
      rw_RadDwFlux(0:kmax,k_w) = 0.0_8
    else    

!      if( w_WaveNum(k_w) <= 300000.0_8 ) then
!        xy_SurfAlbedo(:,:) = 0.0_8
!      else
!      xy_SurfAlbedo(:,:) = SurfAlbedo
!      endif

      do k_r = 1, kmax

        xyz_AF (:,:,k_r) = zw_AsymFact   (k_r,k_w)
        xyz_SSA(:,:,k_r) = zw_SingleScatA(k_r,k_w)

        if(xyz_SSA(0,1,k_r) >= 1.0_8) then ! 2014/04/21
          xyz_SSA(:,:,k_r) = 1.0_8 - 10.0**(-16) !2014/04/22 !0.999999_8
        end if

      end do

      do k_r = 0, kmax
        xyr_OptDep(:,:,k_r) = rw_OptDepTOA(k_r, k_w)
      end do

      xy_SurfPFInted(:, :) = PI * integral_planck(w_WaveNum(k_w) - 0.5_8, w_WaveNum(k_w) + 0.5_8, r_Temp(0))

      do k_r = 0, kmax
        xyr_PFInted(:, :,k_r) = PI * integral_planck(w_WaveNum(k_w) - 0.5_8, w_WaveNum(k_w) + 0.5_8, r_Temp(k_r))
      end do

      SolarFluxTOA = PI * RatioSolarConst * integral_planck(w_WaveNum(k_w) - 0.5_8, w_WaveNum(k_w) + 0.5_8, SolTemp)

      ! dcpam5 radiation code

      if( InAngle >= 0.0_8 ) then

        call RadRTETwoStreamAppSW(            &
          & xyz_SSA, xyz_AF,                  & ! (in)
          & xyr_OptDep,                       & ! (in)
          & xy_SurfAlbedo,                    & ! (in)
          & SolarFluxTOA, xy_InAngle,         & ! (in)
          & xyr_RadUwFlux, xyr_RadDwFlux      & ! (out)
          & )

        rw_SolUwFlux(0:kmax,k_w) = xyr_RadUwFlux(0,1,0:kmax)
        rw_SolDwFlux(0:kmax,k_w) = xyr_RadDwFlux(0,1,0:kmax)

      else

        rw_SolUwFlux(0:kmax,k_w) = 0.0_8
        rw_SolDwFlux(0:kmax,k_w) = 0.0_8

      end if

      if( Flag_Emit == 1 ) then

        if( Flag_FluxModel == 1) then

          call RadRTETwoStreamAppLW_sf(                        &
            & xyz_SSA, xyz_AF,                                 & ! (in)
            & xyr_OptDep,                                      & ! (in)
            & xy_SurfAlbedo,                                   & ! (in)
            & xyr_PFInted, xy_SurfPFInted, xy_SurfDPFDTInted,  & ! (in)
            & xyr_RadUwFlux, xyr_RadDwFlux,                    & ! (out)
            & xyra_DelRadUwFlux, xyra_DelRadDwFlux             & ! (out)
            & )

        elseif( Flag_FluxModel == 2) then

          call RadRTETwoStreamAppLW_TS(                        &
            & xyz_SSA, xyz_AF,                                 & ! (in)
            & xyr_OptDep,                                      & ! (in)
            & xy_SurfAlbedo,                                   & ! (in)
            & xyr_PFInted, xy_SurfPFInted, xy_SurfDPFDTInted,  & ! (in)
            & xyr_RadUwFlux, xyr_RadDwFlux,                    & ! (out)
            & xyra_DelRadUwFlux, xyra_DelRadDwFlux             & ! (out)
            & )

        elseif( Flag_FluxModel == 3) then

          !print *, 'Flag_FluxModel == 2 is now under constraction!'
          !stop

          !== set1: grid
          imax2 = 1
          jmax2 = 1

          NLayMax = 0
          r_NLaySum(0) = 0
          do k_r = 1, kmax

            if(.true.) then
              call DivLayCount (                                &
                & eps_divlay,                                   &
                & xyr_OptDep(0,1,k_r-1)-xyr_OptDep(0,1,k_r),    &
                & xyz_SSA(0,1,k_r), xyz_AF(0,1,k_r),            &
                & xyr_PFInted(0,1,k_r-1), xyr_PFInted(0,1,k_r), &
                & 0.0_8,                                        & ! surface albedo
                & z_NLay(k_r)                                   )
            else
              z_NLay(k_r) = 1 !1 !=========--
            end if
            !print *, 'NLAY:', z_NLay(k_r)

            NLayMax = NLayMax + z_NLay(k_r)
            r_NLaySum(k_r) = NLayMax
            !print *, z_NLaySum(k_r)
          end do ! k_r

          !stop
          !if(.true.) then !======================
          kmax2 = NLayMax

          if(kmax2 > 2**24) then
            print *, 'kmax2 is too big: kmax2 = ', kmax2
            stop
          end if

          call GridsetSet2(imax2, jmax2, kmax2)

          allocate (xyb_SSA          (0:imax2-1, 1:jmax2, 1:kmax2) )
          allocate (xyb_AF           (0:imax2-1, 1:jmax2, 1:kmax2) )
          allocate (xyb_OptDep       (0:imax2-1, 1:jmax2, 0:kmax2) )
          allocate (xyb_PFInted      (0:imax2-1, 1:jmax2, 0:kmax2) )
          allocate (xyb_RadUwFlux    (0:imax2-1, 1:jmax2, 0:kmax2) )
          allocate (xyb_RadDwFlux    (0:imax2-1, 1:jmax2, 0:kmax2) )
          allocate (xyba_DelRadUwFlux(0:imax2-1, 1:jmax2, 0:kmax2, 0:1) )
          allocate (xyba_DelRadDwFlux(0:imax2-1, 1:jmax2, 0:kmax2, 0:1) )


          do k_i = 0, imax2-1
            do k_j = 1, jmax2

              k_r = 1
              xyb_SSA(k_i,k_j,1:r_NLaySum(k_r)) = xyz_SSA(k_i,k_j,k_r)
              xyb_AF (k_i,k_j,1:r_NLaySum(k_r)) = xyz_AF (k_i,k_j,k_r)

              do k_r = 2, kmax
                xyb_SSA(k_i,k_j,r_NLaySum(k_r-1)+1:r_NLaySum(k_r)) = xyz_SSA(k_i,k_j,k_r)
                xyb_AF (k_i,k_j,r_NLaySum(k_r-1)+1:r_NLaySum(k_r)) = xyz_AF (k_i,k_j,k_r)
              end do


             ! k_r = 1
             ! do k_b = 0, r_NLaySum(k_r)-1
             !   xyb_PFInted(k_i,k_j,k_b) = xyr_PFInted(k_i,k_j,k_r-1) - &
             !     & (xyr_PFInted(k_i,k_j,k_r-1) - xyr_PFInted(k_i,k_j,k_r))* &
             !     & k_b/z_NLay(k_r)
             ! end do

              do k_r = 1, kmax
                do k_b = r_NLaySum(k_r-1), r_NLaySum(k_r)-1
                  xyb_PFInted(k_i,k_j,k_b) = xyr_PFInted(k_i,k_j,k_r-1) - &
                    & (xyr_PFInted(k_i,k_j,k_r-1) - xyr_PFInted(k_i,k_j,k_r))* & 
                    & (k_b - r_NLaySum(k_r-1))/z_NLay(k_r)
                end do
              end do

              xyb_PFInted(k_i,k_j,r_NLaySum(kmax)) = xyr_PFInted(k_i,k_j,kmax)

             ! k_r = 1
             ! do k_b = 0, r_NLaySum(k_r)-1
             !   xyb_OptDep(k_i,k_j,k_b) = xyr_OptDep(k_i,k_j,k_r-1) - &
             !     & (xyr_OptDep(k_i,k_j,k_r-1) - xyr_OptDep(k_i,k_j,k_r))* &
             !     & k_b/z_NLay(k_r)
             ! end do

              do k_r = 1, kmax
                do k_b = r_NLaySum(k_r-1), r_NLaySum(k_r)-1
                  xyb_OptDep(k_i,k_j,k_b) = xyr_OptDep(k_i,k_j,k_r-1) - &
                    & (xyr_OptDep(k_i,k_j,k_r-1) - xyr_OptDep(k_i,k_j,k_r))* & 
                    & (k_b - r_NLaySum(k_r-1))/z_NLay(k_r)
                end do
              end do

              xyb_OptDep(k_i,k_j,r_NLaySum(kmax)) = 0.0_8

            end do
          end do

          !print *, '==ORIGINAL=='
          !do k_r = 0, kmax
          !  print *, xyr_OptDep(0,1,k_r), xyr_PFInted(0,1,k_r)
          !end do

          !print *, '==DIVIDED=='
          !do k_b = 0, kmax2
          !  print *, xyb_OptDep(0,1,k_b), xyb_PFInted(0,1,k_b)
          !end do

          call RadRTETwoStreamAppLW_itg2(                      &
            & xyb_SSA, xyb_AF,                                 & ! (in)
            & xyb_OptDep,                                      & ! (in)
            & xy_SurfAlbedo,                                   & ! (in)
            & xyb_PFInted, xy_SurfPFInted, xy_SurfDPFDTInted,  & ! (in)
            & xyb_RadUwFlux, xyb_RadDwFlux,                    & ! (out)
            & xyba_DelRadUwFlux, xyba_DelRadDwFlux             & ! (out)
            & )


          !if(.false.) then
          !  call RadRTETwoStreamAppLW_it(                        &
          !    & xyz_SSA, xyz_AF,                                 & ! (in)
          !    & xyr_OptDep,                                      & ! (in)
          !    & xy_SurfAlbedo,                                   & ! (in)
          !    & xyr_PFInted, xy_SurfPFInted, xy_SurfDPFDTInted,  & ! (in)
          !    & xyr_RadUwFlux, xyr_RadDwFlux,                    & ! (out)
          !    & xyra_DelRadUwFlux, xyra_DelRadDwFlux             & ! (out)
          !    & )
          !end if

          do k_r = 0, kmax
            xyr_RadUwFlux(:,:,k_r) = xyb_RadUwFlux(:,:,r_NLaySum(k_r))
            xyr_RadDwFlux(:,:,k_r) = xyb_RadDwFlux(:,:,r_NLaySum(k_r))
          end do


          deallocate (xyb_SSA           )
          deallocate (xyb_AF            )
          deallocate (xyb_OptDep        )
          deallocate (xyb_PFInted       )
          deallocate (xyb_RadUwFlux     )
          deallocate (xyb_RadDwFlux     )
          deallocate (xyba_DelRadUwFlux )
          deallocate (xyba_DelRadDwFlux )

          !end if !=====================

        end if

        rw_RadUwFlux(0:kmax,k_w) = xyr_RadUwFlux(0,1,0:kmax)
        rw_RadDwFlux(0:kmax,k_w) = xyr_RadDwFlux(0,1,0:kmax)

      else

        rw_RadUwFlux(0:kmax,k_w) = 0.0_8
        rw_RadDwFlux(0:kmax,k_w) = 0.0_8

      end if

    end if

    !print *, 'WN:', w_WaveNum(k_w)
    if(.false.) then
      print *, 'RadUwFlux'
      print *, '=w:', w_WaveNum(k_w)
      do k_r = 0, kmax
        print *, rw_RadUwFlux(k_r, k_w)
      end do

      print *, 'RadDwFlux'
      print *, '=w:', w_WaveNum(k_w)
      do k_r = 0, kmax
        print *, rw_RadDwFlux(k_r, k_w)
      end do
    end if

  end do ! k_w


  ! 出力ファイルの設定
  ! setting of output file
  !
  call HistoryCreate(                                                &
    & file=''//OutFileFlux//'',                                      &
    & title='Flux; '//OutFileFlux//'',                               &
    & source='main_flux.f90;'//InFileOpt//'',                       &
    & institution='DCRTM',                                           &
    & dims=(/'r_Press  ', 'w_WaveNum'/), dimsizes=(/kmax+1,wmax+1/), &
    & longnames=(/'pressure  ', 'wavenumber'/),                      &
    & units=(/'Pa ', 'm-1'/)                                         )

  call HistoryPut('r_Press  ',r_Press)
  call HistoryPut('w_WaveNum',w_WaveNum)

  call HistoryAddVariable( &
    & varname='r_Temp', dims=(/'r_Press'/), &
    & longname='temperature', units='K', xtype='double')
  call HistoryPut('r_Temp',r_Temp)

  call HistoryAddVariable( &
    & varname='rw_RadFluxUp', dims=(/'r_Press  ', 'w_WaveNum'/), &
    & longname='upward flux(planet)', units='W m-2 m', xtype='double')
  call HistoryPut('rw_RadFluxUp',rw_RadUwFlux)

  call HistoryAddVariable( &
    & varname='rw_RadFluxDn', dims=(/'r_Press  ', 'w_WaveNum'/), &
    & longname='downward flux(planet)', units='W m-2 m', xtype='double')
  call HistoryPut('rw_RadFluxDn',rw_RadDwFlux)

  call HistoryAddVariable( &
    & varname='rw_SolFluxUp', dims=(/'r_Press  ', 'w_WaveNum'/), &
    & longname='upward flux(solar)', units='W m-2 m', xtype='double')
  call HistoryPut('rw_SolFluxUp',rw_SolUwFlux)

  call HistoryAddVariable( &
    & varname='rw_SolFluxDn', dims=(/'r_Press  ', 'w_WaveNum'/), &
    & longname='downward flux(solar)', units='W m-2 m', xtype='double')
  call HistoryPut('rw_SolFluxDn',rw_SolDwFlux)

  call HistoryClose


  deallocate( r_Press           )
  deallocate( r_Temp            )
  deallocate( w_WaveNum         )
  deallocate( rw_RadUwFlux      )
  deallocate( rw_RadDwFlux      )
  deallocate( rw_SolUwFlux      )
  deallocate( rw_SolDwFlux      )
  deallocate( rw_OptDepTOA      )
  deallocate( zw_SingleScatA    )

  deallocate (xyz_SSA           )
  deallocate (xyz_AF            )
  deallocate (xyr_OptDep        )
  deallocate (xy_SurfAlbedo     )
  deallocate (xyr_PFInted       )
  deallocate (xy_SurfPFInted    )
  deallocate (xy_SurfDPFDTInted )
  deallocate (xyr_RadUwFlux     )
  deallocate (xyr_RadDwFlux     )
  deallocate (xyra_DelRadUwFlux )
  deallocate (xyra_DelRadDwFlux )
  deallocate (xy_InAngle        )

end program main_flux
