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

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

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

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

  ! 物理定数設定1
  ! physical constants setteing 1
  !
  use constants0, only: GasRUniv

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

  ! opacity 計算
  ! opacity calculation
  !
  use mod_LinebyLine, only:                & !mol_absorpkOld,&
    &                       mol_absorpk,   &
    &                       ray_scatk,     &
    &                       FileCheck,     &
    &                       LineCount,     & 
    &                       SetIsoMolWt,   & 
    &                       SetIsoQ0Val,   &
    &                       SetIsoMRatioE, &
    &                       ReadHITRAN,    &
    &                       ReadParSumNC,  &
    &                       SetQTm

  ! continuum absorption 計算
  ! continuum absorption calculation
  !
  use mod_contabs, only: cont_absorpk

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

  ! 宣言文 ; Declaration statements
  !
  implicit none

  ! 設定変数, 他
  ! setteing parameter, values
  !
  character(99) :: PATHQ           = './src/optdep/parsum.dat.nc' !(default)
                          ! parsum.dat.nc file のパス
                          ! path of parsum.dat.nc file
  character(10) :: DataBase        = 'HITRAN2012' !(default)
                          ! 線吸収データベース
                          ! database of line absorption 
  character(99) :: InFileName
                          ! 入力プロファイル
                          ! file name of input profile
  character(99) :: OutFileOptDep   = 'OptDep_dcrtm.nc' !(default)
                          ! 出力プロファイル
                          ! file name of output profile (optical depth)
  character(99) :: OutFileOptProp  = 'OptProp_dcrtm.nc' !(default)
                          ! 出力プロファイル
                          ! file name of output profile (absorption coefficient)
  real(8) :: MinWN                 = 0.0_8 !(default)
                          ! 最小波数 [m-1]
                          ! mininum wavenumber [m-1]
  real(8) :: MaxWN                 = 300000.0_8 !(default)
                          ! 最大波数 [m-1]
                          ! maximum wavenumber ]m-1]
  real(8) :: ResWN                 = 100.0_8 !(default)
                          ! 波数分解能 [m-1]
                          ! resolution of wavenumber [m-1]
  integer :: Flag_ISOTOPE          = 0 !(default)
                          ! 同位体比変更フラグ
                          ! flag of isotope abandance
                          ! 0: Earth abandance, 1: read user defined abandance
  integer :: Flag_LINE             = 1 !(default)
                          ! 線吸収計算フラグ
                          ! flag of line absorption
                          ! 0: no line calculation, 1: line calculation
  integer :: Flag_CONT             = 1 !(default)
                          ! 連続吸収計算フラグ
                          ! flag of continuum absorption
                          ! 0: no continuum, 1: continuum (MT_CKD model), 2: continuum (Pollack+1993)
  integer :: Flag_RAYLEIGH         = 1 !(default)
                          ! Rayleigh 散乱計算フラグ
                          ! flag of Rayleigh scattering
                          ! 0: no Rayleigh scattering, 1: Rayleigh scattering
  integer :: Flag_LINESHAPE        = 1 !(default)
                          ! 吸収線形フラグ
                          ! flag of line shape
                          ! 1: voigt (all), 2: voigt & CO2 sub-Lorentzian
                          ! 3: Lorentz (all), 4: Lorentz & CO2 sub-Lorentzian
  integer :: Flag_CUTOFF           = 0 !(default)
                          ! cut off フラグ
                          ! flag of cut off
                          ! 0: 25cm-1 all molecules, 1: read user defined cut off value
  real(8), allocatable :: r_Press(:)
                          ! 気圧 (半整数レベル) [Pa]
                          ! pressure (half level) [Pa]
  real(8), allocatable :: r_Temp(:)
                          ! 気温 (半整数レベル) [K]
                          ! temperature (half level) [K]
  real(8) :: MolWt0                 = 28.96d-3 !(default)
                          ! 非吸収分子の平均分子量 [kg mol-1]
                          ! molecular weight of non-absorption gas [kg mol-1]
  integer, allocatable :: m_MolNum(:)
                          ! 分子の番号 [1]
                          ! Molecule's Number [1]
  real(8), allocatable :: m_MolWt(:)
                          ! 分子量 [kg mol-1]
                          ! molecular weight [kg mol-1]
  real(8), allocatable :: rm_MixRatio(:,:)
                          ! 体積混合比 (半整数レベル) [1]
                          ! volume mixing ratio (half level) [1]
  real(8) :: a_alpha(1:50) 
                          ! coefficient of continuum absorption (Pollack et al. (1993) )
                          ! [kg2 m-4 molecule-1]

  integer :: Flag_OUTPUT                = 3 !(default)
                          ! flag of output
                          ! Flag_OUTPUT = 1: optical depth only, 2: coefficient only, 3: optical depth and coefficient
  real(8), allocatable :: m_CutOffVal(:)
                          ! cut off 波数
                          ! value of cut off wavenumber [m-1]

  ! 計算される物理量
  ! calculated values
  !
  real(8), allocatable :: r_MolAve(:)
                          ! 平均分子量 (半整数レベル) [kg mol-1]
                          ! average molerular weight [kg mol-1]
  real(8), allocatable :: w_WaveNum(:)
                          ! 波数 [m-1]
                          ! wavenumber [m-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_TAbsorpK(:,:,:)
                          ! 線吸収, 連続吸収をあわせた吸収係数 [m2 kg-1]
                          ! absorption coefficient for line & continuum absotption [m2 kg-1]
  real(8), allocatable :: rmw_RayScatK(:,:,:)
                          ! Rayleigh 散乱の吸収係数 [m2 kg-1]
                          ! scattering coefficient for Rayleigh scattering [m2 kg-1]

  real(8), allocatable :: z_Press(:)
                          ! 気圧 (整数レベル) [Pa]
                          ! pressure (full level) [Pa]
  real(8), allocatable :: z_Temp(:)
                          ! 気温 (整数レベル) [K]
                          ! temperature (full level) [K]
  real(8), allocatable :: zm_MixRatio(:,:)
                          ! 体積混合比 (整数レベル) [1]
                          ! volume mixing ratio (full level) [1]
  real(8), allocatable :: z_MolAve(:)
                          ! 平均分子量 (整数レベル) [kg mol-1]
                          ! average molerular weight [kg mol-1]
  real(8), allocatable :: rw_OptDepTOA(:,:)
                          ! 大気上端からの光学的厚さ (半整数レベル) [1]
                          ! optical depth from top of atmosphere (half level) [1]
  real(8), allocatable :: zw_OptDepA(:,:)
                          ! 大気層あたりの光学的厚さ [1]
                          ! optical depth of atmospheric layer [1]
  real(8), allocatable :: zw_OptDepS(:,:)
                          ! 大気層あたりの光学的厚さ [1]
                          ! optical depth of atmospheric layer [1]
  real(8), allocatable :: zw_SingleScatA(:,:)
                          ! 1次散乱アルベド (full level) [1]
                          ! single scattering albedo (full level) [1]

  ! 作業変数, 他
  ! work variables
  !
  integer :: mfo = 11
                          ! 装置番号
                          ! device number
  integer :: fnml = 12
                          ! 装置番号(name list file)
                          ! device number : namelist file
  integer :: kmax
                          ! 鉛直レベル最大値
                          ! vertical level maximum
  integer :: wmax
                          ! 波数レベル最大値
                          ! wavenumber level maximum
  integer :: mmax
                          ! 分子数最大値
                          ! molecular number maximum
  integer :: k_r, k_w , k_m, k_i
                          ! 作業変数
                          ! work variables
  character( 3) :: charData
                          ! database 略称
                          ! database short name
  character( 4) :: Lnum4
                          ! 大気層番号
                          ! layer number for output file
  character( 3) :: charCONT
                          ! 連続吸収フラグ
                          !
  real(8) :: contVMR(1:5)
                          ! 連続吸収計算用体積混合比 [1]
                          ! volume mixing ratio for MT_CKD [1]
  real(8) :: MaxMixRatio  
                          ! 最大混合比
                          ! maximum volume mixing ratio
  real(8), allocatable :: m_SumMixRatio(:)
                          ! 混合比の和
                          ! summation of volume mixing ratio
  real(8), pointer :: a_Axis(:) => null()
                          ! 配列のサイズを調べるための配列
                          ! dummy array
  real(8) :: a_CutOffVal(1:50)
                          ! namelist から cut off 読むための配列
                          ! array for namelist of cut off value
  real(8) :: ia_IsoMolWt(1:10, 1:50)
                          ! 同位体の分子質量 [kg mol-1]
                          ! molecular weight of isotope
  real(8) :: ia_IsoQ0Val(1:10, 1:50)
                          ! Q value (269K)
                          !
  real(8) :: ia_IsoMRatioE(1:10, 1:50)
                          ! 地球の同位体の存在比
                          ! existence of isotope in Earth
  real(8) :: ia_IsoMRatioP(1:10, 1:50)
                          ! 惑星大気の同位体存在比
                          ! existence of isotope in planetary atmosphere
  real(8) :: ia_IsoARatio(1:10, 1:50)
                          ! 同位体存在比の地球に対するののの比
                          ! Ratio of isotope existence in planetary atmosphere to that in Earth
  integer :: FileEnd
                          ! HITRAN データベースの行数
                          ! Number of lines of HITRAN database
  integer, allocatable :: af_IntHITRAN(:,:)
                          ! HITRAN data for integer 
                          !
  real(8), allocatable :: af_RealHITRAN(:,:)
                          ! HITRAN data for real
                          !
  real(8), allocatable :: imt_Qval(:,:,:)
                          ! parsum.dat.nc のデータ
                          ! data of parsum.dat.nc
  real(8) :: i_Qval(1:10)
                          ! 温度, 分子種を指定した parsum.dat.nc のデータ
                          ! data of parsum.dat.nc for a Temperature and a Molecule

  integer :: dummyN, k
  real(8) :: cputime1, cputime2
  !character(99) :: PATHQold = './src/optdep/parsum.dat'

  ! namelist file の設定
  ! setting of namelist file
  !
  namelist /optdep_nml/ &
    & InFileName,     &
    & DataBase,       &
    & PATHQ,          &
    & MolWt0,         &
    & MinWN,          &
    & MaxWN,          &
    & ResWN,          &
    & Flag_ISOTOPE,   &
    & Flag_LINE,      &
    & Flag_LINESHAPE, &
    & Flag_CUTOFF,    &
    & Flag_CONT,      &
    & Flag_RAYLEIGH,  &
    & Flag_OUTPUT,    &
    & OutFileOptDep,  &
    & OutFileOptProp

  namelist /cutoff_nml/ &
    & a_CutOffVal

  namelist /isomixratio_nml/ &
    & ia_IsoMRatioP

  namelist /rhosquarecs_nml/ &
    & a_alpha


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

  ! 設定値のチェック
  ! check setting values
  !
  !print *, '==== ==== ==== ===='
  print *, '== Setting Values =='
  print *, 'InFileName:', InFileName
  print *, 'DataBase:', DataBase
    select case (DataBase)
      case ('HITRAN1996', 'HITRAN2008', 'HITRAN2012', 'HITEMP2010')
      case ('TestData')
        print *, '* test line data is used'
      case default
        print *, '* invalid value (DataBase): stop'
        stop
    end select

  print *, 'PATHQ:', PATHQ

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

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

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

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

  print *, 'Flag_ISOTOPE:', Flag_ISOTOPE
    if(Flag_ISOTOPE == 0) then
      print *, '* Earth isotope abandance'
    elseif(Flag_ISOTOPE == 1) then
      print *, '* read user defined isotope abandance'
    else
      print *, '* invalid value (Flag_ISOTOPE = 0 or 1): stop'
      stop
    end if

  print *, 'Flag_LINE:', Flag_LINE
    if    (Flag_LINE == 0) then
      print *, '* line absorption is NOT calculated'
    elseif(Flag_LINE == 1) then
      print *, '* line absorption is calculated'

      print *, 'Flag_LINESHAPE:', Flag_LINESHAPE
        if    (Flag_LINESHAPE == 1) then
          print *, '* voigt (all)'
        elseif(Flag_LINESHAPE == 2) then
          print *, '* voigt & CO2 sub-Lorentzian'
        elseif(Flag_LINESHAPE == 3) then
          print *, 'Lorentz (all)'
        elseif(Flag_LINESHAPE == 4) then
          print *, 'Lorentz & CO2 sub-Lorentzian'
        else
          print *, '* invalid value (Flag_LINESHAPE = 1, 2, 3  or 4): stop'
          stop
        end if

        print *, 'Flag_CUTOFF:', Flag_CUTOFF
        if    (Flag_CUTOFF == 0) then
          print *, '* 25cm-1 all molecules'
        elseif(Flag_LINE == 1) then
          print *, '* read user defined cut off value (cutoff_nml)'
        else
          print *, '* invalid value (Flag_CUTOFF = 0 or 1): stop'
          stop
        end if

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

  print *, 'Flag_CONT:', Flag_CONT
    if    (Flag_CONT == 0) then
      print *, '* no continuum'
    elseif(Flag_CONT == 1) then
      print *, '* continuum (MT_CKD model)'
    elseif(Flag_CONT == 2) then
      print *, '* continuum (Pollack+1993)'
    else
      print *, '* invalid value (Flag_LINE = 0, 1 or 2): stop'
      stop
    end if

  print *, 'Flag_RAYLEIGH:', Flag_RAYLEIGH
    if    (Flag_RAYLEIGH == 0) then
      print *, '* Rayleigh scattering is NOT calculated'
    elseif(Flag_RAYLEIGH == 1) then
      print *, '* Rayleigh scattering is calculated'
    else
      print *, '* invalid value (Flag_RAYLEIGH = 0 or 1): stop'
      stop
    end if

  print *, 'Flag_OUTPUT:', Flag_OUTPUT
    if    (Flag_OUTPUT == 1) then
      print *, '* optical depth & single scattering albedo are calculated'
    elseif(Flag_OUTPUT == 2) then
      print *, '* absorption and scattering coefficients are calculated'
    elseif(Flag_OUTPUT == 3) then
      print *, '* optical depth, single scattering albedo and optical coefficient are calculated'
    else
      print *, '* invalid value (Flag_LINE = 1, 2 or 3): stop'
      stop
    end if

  if( (Flag_OUTPUT == 1) .or. (Flag_OUTPUT == 3) ) then
    print *, 'OutFileOptDep:', OutFileOptDep
  end if

  if( (Flag_OUTPUT == 2) .or. (Flag_OUTPUT == 3) ) then
    print *, 'OutFileOptProp:', OutFileOptProp
  end if
!  print *, '====E====N====D===='
  print *, '=setting check end='

  ! 配列サイズを調べる
  ! check array size
  !
  call HistoryGetPointer(''//InFileName//'', 'r_Press', a_Axis )

  kmax = size( a_Axis ) - 1

  Nullify( a_Axis )

  !print *, kmax

  call HistoryGetPointer(''//InFileName//'', 'm_MolWt', a_Axis)

  mmax = size( a_Axis )

  Nullify( a_Axis )

  !print *, mmax

  wmax = Int( ( MaxWN - MinWN ) /ResWN )

  allocate( r_Press       (0:kmax                ) )
  allocate( r_Temp        (0:kmax                ) )
  allocate( r_MolAve      (0:kmax                ) )
  allocate( m_MolNum      (        0:mmax        ) )
  allocate( m_MolWt       (        0:mmax        ) )
  allocate( m_CutOffVal   (        0:mmax        ) )
  allocate( rm_MixRatio   (0:kmax, 0:50          ) ) !mmax) )
  allocate( m_SumMixRatio (        0:mmax        ) )
  allocate( w_WaveNum     (                0:wmax) )
  allocate( rmw_LAbsorpK  (0:kmax, 0:mmax, 0:wmax) )
  allocate( rmw_CAbsorpK  (0:kmax, 0:mmax, 0:wmax) )
  allocate( rmw_TAbsorpK  (0:kmax, 0:mmax, 0:wmax) )
  allocate( rmw_RayScatK  (0:kmax, 0:mmax, 0:wmax) )

  allocate( z_Press       (1:kmax                ) )
  allocate( z_Temp        (1:kmax                ) )
  allocate( zm_MixRatio   (1:kmax, 0:mmax        ) )
  allocate( z_MolAve      (1:kmax                ) )
  allocate( rw_OptDepTOA  (0:kmax,         0:wmax) )
  allocate( zw_OptDepA    (1:kmax,         0:wmax) )
  allocate( zw_OptDepS    (1:kmax,         0:wmax) )
  allocate( zw_SingleScatA(1:kmax,         0:wmax) )


  write(charData, '(A1, A2)') 'H', DataBase(9:10)
  write(charCONT, '(A2, I1)') 'CN', Flag_CONT

  !
  ! 同位体比の設定
  ! set isotope abandance ratio
  !
  ia_IsoARatio = 0.0_8

  call SetIsoMRatioE(ia_IsoMRatioE)

  if(Flag_ISOTOPE == 0) then

    do k_m = 1, mmax
      do k_i = 1, 10
        if( ia_IsoMRatioE(k_i, k_m) .ne. 0.0_8) then
          ia_IsoARatio(k_i,k_m) = 1.0_8
        end if
      end do
    end do

  elseif(Flag_ISOTOPE == 1) then

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

    do k_m = 1, mmax
      do k_i = 1, 10
        if( ia_IsoMRatioE(k_i, k_m) .ne. 0.0_8) then
          ia_IsoARatio(k_i,k_m) = ia_IsoMRatioP(k_i,k_m)/ia_IsoMRatioE(k_i,k_m)
        end if
      end do
    end do

  else

    print *, 'Flag_ISOTOPE is wrong number:', Flag_ISOTOPE
    print *, 'stop'
    stop

  end if


  call SetIsoMolWt(ia_IsoMolWt)

  call SetIsoQ0Val(ia_IsoQ0Val)


  ! cut off の設定
  ! setting of cut off
  ! 
  print *, 'Flag_CUTOFF:', Flag_CUTOFF
  if( Flag_CUTOFF == 0 ) then

    m_CutOffVal = 2500.0_8

    print *, '==cut off value=='
    do k_m = 1, mmax
      print *, k_m, m_CutOffVal(k_m)
    end do
  
  elseif( Flag_CUTOFF == 1 ) then

    open(fnml, file='dcrtm.nml')
    read(fnml, nml=cutoff_nml)
    close(fnml)

    print *, '==cut off value=='
    do k_m = 1, mmax
      m_CutOffVal(k_m) = a_CutOffVal(k_m)
      print *, k_m, m_CutOffVal(k_m)
    end do
  
  else
    print *, 'Flag_CUTOFF is wrong number: ', Flag_CUTOFF
    print *, 'stop'
    stop
  end if


  ! 波数配列の計算
  ! wavenumber calculation
  !
  do k_w = 0, wmax
    w_WaveNum(k_w) = MinWN + ResWN * k_w
    !print *, k_w, w_WaveNum(k_w)
  end do

  ! 
  ! 
  do k_m = 0, mmax
    m_MolNum(k_m) = k_m
  end do

  rm_MixRatio(:,:) = 0.0_8

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

  !! 大気上端の圧力のチェック              ! 2016/12/07
  !! check TOA pressure is NOT zero or not
  !!
  !if ( r_Press(kmax) == 0.0_8 ) then
  !  print *, 'pressure of TOA = 0 is not alloweded: stop'
  !  stop
  !end if  

  ! 非吸収分子の体積混合比の計算
  ! 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


  m_SumMixRatio = 0.0_8
  do k_r = 0, kmax
    m_SumMixRatio(0:mmax) = m_SumMixRatio(0:mmax) + rm_MixRatio(k_r,0:mmax)
  end do


  !
  ! HITRAN data check
  !
  if(Flag_Line == 1) then
    do k_m = 1, mmax
      if( m_SumMixRatio(k_m) .ne. 0.0_8 ) then
        call FileCheck(DataBase, k_m)
      end if
    end do
  end if

  !
  ! parsum.dat.nc 読み込み
  !
  allocate( imt_Qval(1:10,1:mmax,1:2931) )
  call ReadParSumNC( PATHQ, mmax, imt_Qval )


  ! 平均分子量の計算
  ! calculate average molecular weight
  !
  do k_r = 0, kmax
    r_MolAve(k_r) = 0.0_8
    do k_m = 0, mmax
      r_MolAve(k_r) = r_MolAve(k_r) + m_MolWt(k_m)*rm_MixRatio(k_r,k_m)
    end do
  end do


  rmw_LAbsorpK = 0.0_8
  rmw_CAbsorpK = 0.0_8
  rmw_TAbsorpK = 0.0_8
  rmw_RayScatK = 0.0_8


  ! Rayleigh 散乱の計算
  ! Rayleigh Scattering Cross Section 
  !  
  print *, '== Rayleigh Scattering calculation=='
  print *, 'Flag_RAYLEIGH:', Flag_RAYLEIGH
  if(Flag_RAYLEIGH == 1) then
    do k_r = 0, kmax

      do k_m = 0, mmax
        if(m_SumMixRatio(k_m) .ne. 0.0_8) then
          call ray_scatk (                           & ! 2016/06/21
            & k_m, m_MolWt(k_m), wmax, MinWN, ResWN, & ! (in )
            & rmw_RayScatK(k_r,k_m,:)                ) ! (out)
        end if
      end do

    end do ! k_r
  end if


  ! 連続吸収の計算
  ! calculate continuum
  !
  print *, '== Continuum Absorption calculation=='
  print *, 'Flag_CONT:', Flag_CONT

  if(Flag_CONT == 2) then

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

    do k_m = 1, mmax
      print *, k_m, 'alpha=', a_alpha(k_m)
    end do

  else

    a_alpha = 0.0_8

  end if

  call cont_absorpk(                & !
    & Flag_CONT, mmax, kmax,        & !(in )
    & r_Press, r_Temp, rm_MixRatio, & !(in )
    & m_MolWt, MinWN, MaxWN, ResWN, & !(in )
    & a_alpha,                      & !(in )
    & rmw_CAbsorpK                  ) !(out)

  !call cpu_time(cputime1)

  ! 線吸収の計算
  ! calculate line absorption
  !
  print *, '==Line Absorption calculation=='
  print *, 'Flag_LINE:', Flag_LINE
  if(Flag_LINE == 1) then !2015/02/04

    do k_m = 1, mmax !!!!!!!!!!!2016/10/26 k_r, k_m loop transpose

    print *, 'molecule:', k_m

  !if(.true.) then

    if( m_SumMixRatio(k_m) .ne. 0.0_8 ) then

      call LineCount( DataBase, k_m, FileEnd )
      !print *, 'FileEnd:', FileEnd

      allocate( af_IntHITRAN (1:2,1:FileEnd) )
      allocate( af_RealHITRAN(1:8,1:FileEnd) )

      call ReadHITRAN(                       &
        & DataBase, k_m, FileEnd,            & !(in )
        & af_IntHITRAN, af_RealHITRAN        ) !(out)

      do k_r = 0, kmax !kmax, 0, -1

        if( rm_MixRatio(k_r,k_m) .ne. 0.0_8 ) then

          call SetQTm(                        & ! 2016/10/27
            & r_Temp(k_r), imt_Qval(:,k_m,:), & ! (in )
            & i_Qval                          ) ! (out)

          call mol_absorpk(                                                & !2016/10/27
            & Flag_LINESHAPE, m_CutOffVal(k_m), i_Qval,                    & !(in )
            & ia_IsoARatio(:,k_m), ia_IsoQ0Val(:,k_m), ia_IsoMolWt(:,k_m), & !(in )
            & af_IntHITRAN, af_RealHITRAN, FileEnd,                        & !(in )
            & k_m, m_MolWt(k_m), r_Press(k_r), r_Temp(k_r),      & !(in )
            & rm_MixRatio(k_r,k_m), wmax, MinWN, ResWN,                    & !(in )
            & rmw_LAbsorpK(k_r,k_m,:)                                      ) !(out)

        else
          rmw_LAbsorpK(k_r,k_m,:) = 0.0_8
        end if

      end do ! k_r in k_m

      deallocate( af_IntHITRAN   )
      deallocate( af_RealHITRAN  )

    else
      rmw_LAbsorpK(:,k_m,:) = 0.0_8
    end if

  !else

!      do k_r = kmax, 0, -1
!      !print *, 'absorption coefficient: k_r=', k_r
!  !    do k_m = 1, mmax

!        if( rm_MixRatio(k_r, k_m) .ne. 0.0_8 ) then

!          call mol_absorpkOld(                                        & !      2016/06/21
!            & Flag_LINESHAPE, m_CutOffVal(k_m), PATHQold,             & !(in ) 2016/10/18
!            & ia_IsoARatio(:,k_m), ia_IsoQ0Val(:,k_m), ia_IsoMolWt(:,k_m), & ! 2016/10/20
!            & DataBase, k_m, m_MolWt(k_m), r_Press(k_r), r_Temp(k_r), & !(in )
!            & rm_MixRatio(k_r,k_m), wmax, MinWN, ResWN,               & !(in )
!            & rmw_LAbsorpK(k_r,k_m,:)                                 ) !(out)

!        else

!          rmw_LAbsorpK(k_r,k_m,:) = 0.0_8 

!        end if

  !    end do ! k_m in k_r
!      end do !k_r

!  endif

    end do ! k_m

  end if !

  !call cpu_time(cputime2)

  !print *, 'cputime:', cputime2 - cputime1

  !do k_w = 0, wmax
  !  print *, w_WaveNum(k_w)
  !end do

  !do k_w = 0, wmax
  !  do k_r = 0, kmax
  !    print *, rmw_LAbsorpK(k_r,:,k_w)
  !  end do
  !end do

  ! データ出力
  ! data output
  !
  if(Flag_OUTPUT .ne. 1) then

    rmw_TAbsorpK = rmw_LAbsorpK + rmw_CAbsorpK

    ! 出力ファイルの設定
    ! setting of output file
    !
    call HistoryCreate(                                    &
      & file=''//OutFileOptProp//'', &
      & title='absorption scattering cross section;'//DataBase//','//charCONT//'', &
      & source='main_optdep.f90:'//DataBase//','//charCONT//'', &
      & institution='DCRTM',       &
      & dims=(/'r_Press  ','m_MolNum ', 'w_WaveNum'/), dimsizes=(/kmax+1,mmax+1,wmax+1/),      &
      & longnames=(/'boundary pressure',  &
      &             'molecular number ',  &
      &             'wavenumber       '/),&
      & units=(/'Pa ', '1  ', 'm-1'/)     )

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

    call HistoryAddVariable( &
      & varname='m_MolWt', dims=(/'m_MolNum'/), &
      & longname='molecular weight', units='kg mol-1', xtype='double')
    call HistoryPut('m_MolWt',m_MolWt)

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

    call HistoryAddVariable( &
      & varname='rm_MixRatio', dims=(/'r_Press  ','m_MolNum '/), &
      & longname='volume mixing ratio', units='1', xtype='double')
    call HistoryPut('rm_MixRatio',rm_MixRatio)

    call HistoryAddVariable( &
      & varname='rmw_LAbsorpK', dims=(/'r_Press  ','m_MolNum ', 'w_WaveNum'/), &
      & longname='line absorption coefficient', units='m2 kg-1', xtype='double')
    call HistoryPut('rmw_LAbsorpK',rmw_LAbsorpK)

    call HistoryAddVariable( &
      & varname='rmw_CAbsorpK', dims=(/'r_Press  ', 'm_MolNum ', 'w_WaveNum'/), &
      & longname='continuum absorption coefficient', units='m2 kg-1', xtype='double')
    call HistoryPut('rmw_CAbsorpK',rmw_CAbsorpK)

    call HistoryAddVariable( &
      & varname='rmw_TAbsorpK', dims=(/'r_Press  ', 'm_MolNum ', 'w_WaveNum'/), &
      & longname='line & continuum absorption coefficient', units='m2 kg-1', xtype='double')
    call HistoryPut('rmw_TAbsorpK',rmw_TAbsorpK)

    call HistoryAddVariable( &
      & varname='rmw_RayScatK', dims=(/'r_Press  ', 'm_MolNum ', 'w_WaveNum'/), &
      & longname='rayleigh scattering coefficient', units='m2 kg-1', xtype='double')
    call HistoryPut('rmw_RayScatK',rmw_RayScatK)

    call HistoryClose

  end if

  if(Flag_OUTPUT .ne. 2) then

    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                          )

    do k_r = 1, kmax
      z_Press(k_r) = ( &
        & r_Press(k_r-1) * r_Press(k_r-1)/r_Temp(k_r-1) + &
        & r_Press(k_r  ) * r_Press(k_r  )/r_Temp(k_r  ) ) &
        & /(r_Press(k_r-1)/r_Temp(k_r-1) + r_Press(k_r)/r_Temp(k_r))
    end do


    ! 出力ファイルの設定
    ! setting of output file
    !
    call HistoryCreate(                                    &
      & file=''//OutFileOptDep//'', &
      & title='Optical Depth from TOA;'//DataBase//','//charCONT//'', &
      & source='main_optdep.f90:'//DataBase//','//charCONT//'', &
      & institution='DCRTM',       &
      & dims=(/'r_Press  ', 'w_WaveNum', 'z_Press  '/), dimsizes=(/kmax+1,wmax+1,kmax/),      &
      & longnames=(/'boundary pressure ', 'wavenumber        ', 'mid-layer pressure'/),       &
      & units=(/'Pa ', 'm-1', 'Pa '/)                    )

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

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

    call HistoryAddVariable( &
      & varname='rw_OptDep', dims=(/'r_Press  ', 'w_WaveNum'/), &
      & longname='optical depth from TOA', units='1', xtype='double')
    call HistoryPut('rw_OptDep',rw_OptDepTOA)

    call HistoryAddVariable( &
      & varname='zw_SingleScatA', dims=(/'z_Press  ', 'w_WaveNum'/), & 
      & longname='single scattering albedo', units='1', xtype='double')
    call HistoryPut('zw_SingleScatA',zw_SingleScatA)

    call HistoryClose
  end if

end program main_optdep
