!2012/01/25(linebyline.f90から引き継ぎ)
!2014/05/21
!2014/06/03 ! add continuum
!2014/06/06

!line by line で吸収スペクトルを計算

! 吸収係数、光学的厚さの計算の定式化は Rothman et al., 1998
! 但し、line shift は温度依存させる(2013/05/01)
! 鉛直分布の重みづけについて、温度、混合比の依存性を追加(2013/06/07)
! physical constants: CODATA 2010
! voigt profile test: Kuntz 1997 and Ruyten 2004, and Humlicek 1982 (2013/07/01)
! add line shape function (2013/07/08)

module mod_LineByLine

  implicit none

  ! 物理定数, 計算条件
  ! Pysical constants, condition of calculation
  !
  real(8) :: PI = 3.1415926535897932_8
                              ! $ \pi $ .
                              ! 円周率.  Circular constant
  real(8) :: LightSpeed = 2.99792458e+8_8 
                              ! $ \c $ [m s-1].
                              ! 光速. Light Speed
  real(8) :: PlanckConst = 6.62606957e-34_8 !CODATA 2010
                              ! $ \h $ [J s].
                              ! プランク定数.
  real(8) :: BoltzmannConst = 1.3806488e-23_8 !CODATA 2010
                              ! $ \k_b $ [J K-1]
                              ! ボルツマン定数
  real(8) :: AvogadroNum = 6.02214129e+23_8 !CODATA 2010
                              ! $ \N_A $ [mol-1]
                              ! アボガドロ数
  real(8) :: SecRadConst = 6.62606957e-34_8 * 2.99792458e+8_8 / 1.3806488e-23_8 !CODATA 2010
                              ! $ \c_2 $ [m K-1]
                              ! Second Radiation Constant (hc/k)
  real(8) :: Grav = 9.8_8
                              ! $ g $ [m s-2]. 
                              ! 重力加速度. 
                              ! Gravitational acceleration

  real(8) :: Press_Ref = 101325.0_8
                              ! $ p_0 $ [Pa].
                              ! 圧力
                              ! pressure
  real(8) :: Temp_Ref = 296.0_8
                              ! $ T $ [K].
                              ! 温度
                              ! temperature

  real(8) :: MolVal(1:50, 1:10, 1:3)


contains 

  subroutine MolecularWeight(Press, Temp, MixRatioH, MolAve)

    real(8), intent(in ) :: Press(1:2)
    real(8), intent(in ) :: Temp(1:2)
    real(8), intent(in ) :: MixRatioH(1:2, 1:50)
    real(8), intent(out) :: MolAve
    integer :: k_mol, k_iso
    real(8) :: MolWeight

    MolVal = 0.0_8
    call SetMolVal(2)

    MolAve = 0.0_8 ! 2013/07/18
    
    do k_mol = 1, 50
      if( (MixRatioH(1, k_mol) .ne. 0.0_8) .and. &
        & (MixRatioH(2, k_mol) .ne. 0.0_8) ) then
        MolWeight = 0.0_8
        do k_iso = 1, 10
          MolWeight = MolWeight + MolVal(k_mol, k_iso, 1) * MolVal(k_mol, k_iso, 3)        
        end do
        MolAve = MolAve + MolWeight * &
          & (MixRatioH(1, k_mol) * Press(1)/Temp(1) + MixRatioH(2, k_mol) * Press(2)/Temp(2))/&
          & (Press(1)/Temp(1) + Press(2)/Temp(2))
      end if
    end do  

  end subroutine MolecularWeight

  subroutine cont_opticaldepth(kmax, Press, &
    & Temp, contVMR, MinWaveNum, ResWaveNum, MolAve, &
    & OptDepth)

    integer, intent(in ) :: kmax

    real(8), intent(in ) :: Press(1:2)
                              ! $ p_0 $ [Pa].
                              ! 圧力
                              ! pressure
    real(8), intent(in ) :: Temp(1:2)
                              ! $ T $ [K].
                              ! 温度
                              ! temperature
    real(8), intent(in ) :: contVMR(1:2, 1:5)

    real(8), intent(in ) :: MinWaveNum

    real(8), intent(in ) :: ResWaveNum

    real(8), intent(out) :: OptDepth(0:kmax)

    real(8), intent(in ) :: MolAve

    real(8) :: Absorp(0:kmax)

    !integer :: sk

    call cont_absorption(kmax, Press, &
      & Temp, contVMR, MinWaveNum, ResWaveNum, &
      & Absorp) 

    !do sk = 0, kmax
    !  print *, Absorp(sk)
    !end do

    OptDepth(0:kmax) = abs(Press(1) - Press(2)) * AvogadroNum/&
      & (Grav * MolAve) * Absorp(0:kmax)

  end subroutine cont_opticaldepth

  subroutine cont_absorption(kmax, Press, &
    & Temp, contVMR, MinWaveNum, ResWaveNum, &
    & Absorp) 

    integer, intent(in ) :: kmax

    real(8), intent(in ) :: Press(1:2)
                              ! $ p_0 $ [Pa].
                              ! 圧力
                              ! pressure
    real(8), intent(in ) :: Temp(1:2)
                              ! $ T $ [K].
                              ! 温度
                              ! temperature
    real(8), intent(in ) :: contVMR(1:2, 1:5)

    real(8), intent(in ) :: MinWaveNum

    real(8), intent(in ) :: ResWaveNum

    real(8), intent(out) :: Absorp(0:kmax)

    real(8) :: AbsorpG(1:2, 0:kmax)
    real(8) :: NumDensity(1:2)
    real(8) :: WaveNum(0:kmax)
    real(8) :: ArrVMR_all(1:5)
    real(8) :: InWNmin ![m-1]
    real(8) :: InWNmax ![m-1]
    real(8) :: InWNres ![m-1]
    
    InWNmin = MinWaveNum
    InWNmax = MinWaveNum + kmax * ResWaveNum
    InWNres = ResWaveNum

    NumDensity(1:2) = Press(1:2)/ Temp(1:2)

    call call_mtckd(Press(1), Temp(1), contVMR(1,1:5), &
      & InWNmin, InWNmax, INWNres, &
      & WaveNum, AbsorpG(1, 0:kmax))

    call call_mtckd(Press(2), Temp(2), contVMR(2,1:5), &
      & InWNmin, InWNmax, INWNres, &
      & WaveNum, AbsorpG(2, 0:kmax))

    Absorp(0:kmax) = (AbsorpG(1, 0:kmax) * NumDensity(1) + AbsorpG(2, 0:kmax) * NumDensity(2))/&
      & (NumDensity(1) + NumDensity(2))

  end subroutine cont_absorption

  subroutine call_mtckd(InPress, InTemp, ArrVMR_all, &
    & InWNmin, InWNmax, INWNres, &
    & WaveNum, AbsorpCS)
    real(8), intent(in ) :: InPress ![Pa]
    real(8), intent(in ) :: InTemp  ![K]
    real(8), intent(in ) :: ArrVMR_all(1:5)
    real(8), intent(in ) :: InWNmin ![m-1]
    real(8), intent(in ) :: InWNmax ![m-1]
    real(8), intent(in ) :: InWNres ![m-1]
    real(8), intent(out) :: WaveNum (0:int((InWNmax-InWNmin)/InWNres)+1)
    real(8), intent(out) :: AbsorpCS(0:int((InWNmax-InWNmin)/InWNres)+1)
    real(8) :: Press
    real(8) :: WNmin
    real(8) :: WNmax
    real(8) :: WNres
    real(8) :: AbsorpCS_m(1:2, 1:5001)
    real(8) :: ArrVMR_dry(1:5)
    real(8) :: SumVMR
    real(8) :: N_wn
    integer :: Lmax
    integer :: k


    ! change VMR: VMR_all to VMR_dry
    !
!    SumVMR = 0.0_8
!    do k = 1, 5
!      SumVMR = SumVMR + ArrVMR_all(k)
!    end do

!    if( SumVMR > 1.0_8 ) then
!      print *, 'mixing ratio is over the range: SumVMR:', SumVMR
!      print *, 'stop (call_mtckd)'
!      stop
!    end if

    ! mt_ckd trick??
    !
    ! we need to chane the VMR of N2 for "correct" continuum calculation.
    ! for example, vmr(H2O) = 0.01, any other molecules are no absorption 
    ! i.e., ArrVMR_all(2:5) = 0, MT_CKD does not work (output is bad value).
    ! this case, MT_CKD needs ArrVMR(5) = 1. - vmr(H2O)
    !

    if( ArrVMR_all(1) == 1.0_8 ) then ! All modeclules are H2O
      ArrVMR_dry(1) = 1.0_8
      ArrVMR_dry(2:5) = 0.0_8
    elseif( ArrVMR_all(1) == 0.5_8 ) then !
      !ArrVMR_all(:) = ArrVMR_all(:) * 10.0_8
      ArrVMR_dry(:) = ArrVMR_all(:)*10.0_8/(1.0_8 - ArrVMR_all(1)*10.0_8)
    else
      ArrVMR_dry(:) = ArrVMR_all(:)/(1.0_8 - ArrVMR_all(1))
    end if

    !print *, ArrVMR_dry

    ! change scale: SI to MT_CKD(mbr, cm-1)
    !
    Press = InPress * 1e-2_8
    WNmin = InWNmin * 1e-2_8
    WNmax = InWNmax * 1e-2_8
    WNres = InWNres * 1e-2_8

    N_wn = 1.0_8 + (WNmax - WNmin)/WNres
    Lmax = int( (N_wn - 1.0_8)/5000.0_8 ) + 1
    !print *, 'N_wn:', N_wn, 'Lmax:', Lmax


    do k = 1, Lmax

      if( k == Lmax ) then

        call MT_CKD(Press, InTemp, ArrVMR_dry, &
          & WNmin + (k-1)*5000*WNres, WNmax, &
          & WNres, AbsorpCS_m)

        WaveNum(5000*(Lmax-1):int((WNmax-WNmin)/WNres)+1) = &
          & AbsorpCS_m(1,1:int(N_wn) - 5000*(Lmax-1)) * 100.0_8 ![m-1]

        if( ArrVMR_all(1) == 0.5_8 ) then
          AbsorpCS(5000*(Lmax-1):int((WNmax-WNmin)/WNres)+1) = &
            & AbsorpCS_m(2,1:int(N_wn) - 5000*(Lmax-1)) * 1e-4_8 *0.10_8![m2/molecule]
        else
          AbsorpCS(5000*(Lmax-1):int((WNmax-WNmin)/WNres)+1) = &
            & AbsorpCS_m(2,1:int(N_wn) - 5000*(Lmax-1)) * 1e-4_8 ![m2/molecule]
        end if

      else

        call MT_CKD(Press, InTemp, ArrVMR_dry, &
          & WNmin + (k-1)*5000*WNres, WNmin + (5000*k-1)*WNres, &
          & WNres, AbsorpCS_m)

        WaveNum ((k-1)*5000:5000*k-1) = AbsorpCS_m(1,1:5000) * 100.0_8 ![m-1]

        if( ArrVMR_all(1) == 0.5_8 ) then
          AbsorpCS((k-1)*5000:5000*k-1) = AbsorpCS_m(2,1:5000) * 1e-4_8 *0.10_8![m2/molecule]
        else
          AbsorpCS((k-1)*5000:5000*k-1) = AbsorpCS_m(2,1:5000) * 1e-4_8 ![m2/molecule]
        end if

      end if

    end do

  end subroutine call_mtckd


  subroutine mol_opticaldepth(DataBase, kmax, k_mol, Press, &
    & Temp, MixRatioH, MinWaveNum, ResWaveNum, MolAve, &
    & OptDepth)

    character(10), intent(in ) :: DataBase

    integer, intent(in ) :: kmax

    real(8), intent(in ) :: Press(1:2)
                              ! $ p_0 $ [Pa].
                              ! 圧力
                              ! pressure
    real(8), intent(in ) :: Temp(1:2)
                              ! $ T $ [K].
                              ! 温度
                              ! temperature
    real(8), intent(in ) :: MixRatioH(1:2)
                          ! Mix Ratio of molecules in Hitran database
    real(8), intent(in ) :: MinWaveNum

    real(8), intent(in ) :: ResWaveNum

    real(8), intent(out) :: OptDepth(0:kmax)

    integer, intent(in ) :: k_mol

    real(8), intent(in ) :: MolAve

    real(8) :: Absorp(0:kmax)

    real(8) :: VolMixRatio

    call mol_absorption(DataBase, kmax, k_mol, Press, Temp, &
      & MixRatioH, MinWaveNum, ResWaveNum, &
      & Absorp)

    VolMixRatio = (MixRatioH(1) * Press(1)/Temp(1) + MixRatioH(2) * Press(2)/Temp(2))/&
      & ( Press(1)/Temp(1) + Press(2)/Temp(2) )
    
    OptDepth(0:kmax) = VolMixRatio * abs(Press(1) - Press(2)) * AvogadroNum/&
      & (Grav * MolAve) * Absorp(0:kmax)

    ! 2013/07/17, 2014/05/20
    !print *, 'k_mol=', k_mol
    !print *, 'VolMixRatio=', VolMixRatio
    !print *, 'Molecular Amount=', VolMixRatio * abs(Press(1) - Press(2)) * AvogadroNum/(Grav * MolAve)

  end subroutine mol_opticaldepth

  subroutine mol_absorption(DataBase, kmax, k_mol, Press, &
    & Temp, MixRatioH, MinWaveNum, ResWaveNum, &
    & Absorp) 

    character(10), intent(in ) :: DataBase

    integer, intent(in ) :: kmax

    real(8), intent(in ) :: Press(1:2)
                              ! $ p_0 $ [Pa].
                              ! 圧力
                              ! pressure
    real(8), intent(in ) :: Temp(1:2)
                              ! $ T $ [K].
                              ! 温度
                              ! temperature
    real(8), intent(in ) :: MixRatioH(1:2)
                          ! Mix Ratio of molecules in Hitran database
    real(8), intent(in ) :: MinWaveNum

    real(8), intent(in ) :: ResWaveNum

    real(8), intent(out) :: Absorp(0:kmax)

    integer, intent(in ) :: k_mol

    real(8) :: AbsorpG(1:2, 0:kmax)

    real(8) :: NumDensity(1:2)

    NumDensity(1:2) = MixRatioH(1:2) * Press(1:2)/ Temp(1:2)

    call linebyline(DataBase, kmax, k_mol, Press(1), Temp(1), &
      & MixRatioH(1), MinWaveNum, ResWaveNum, &
      & AbsorpG(1, 0:kmax))

    call linebyline(DataBase, kmax, k_mol, Press(2), Temp(2), &
      & MixRatioH(2), MinWaveNum, ResWaveNum, &
      & AbsorpG(2, 0:kmax))

    Absorp(0:kmax) = (AbsorpG(1, 0:kmax) * NumDensity(1) + AbsorpG(2, 0:kmax) * NumDensity(2))/&
      & (NumDensity(1) + NumDensity(2))

  end subroutine mol_absorption

  subroutine linebyline(DataBase, kmax, k_mol, Press, & 
    & Temp, MixRatioH, MinWaveNum, ResWaveNum, &
    & AbsorpG)

    character(10), intent(in ) :: DataBase

    integer, intent(in ) :: kmax

    real(8), intent(in ) :: Press
                              ! $ p_0 $ [Pa].
                              ! 圧力
                              ! pressure
    real(8), intent(in ) :: Temp
                              ! $ T $ [K].
                              ! 温度
                              ! temperature
    real(8), intent(in ) :: MixRatioH
                          ! Mix Ratio of molecules in Hitran database
    real(8), intent(in ) :: MinWaveNum

    real(8), intent(in ) :: ResWaveNum

    real(8), intent(out) :: AbsorpG(0:kmax)

    integer, intent(in ) :: k_mol

    integer, allocatable :: MolNum(:)
                              ! 
                              ! Molecule number (input HITRAN)
    integer, allocatable :: IsoNum(:)
                              !
                              ! Isotopologue number (input HITRAN)
    real(8), allocatable :: WaveNum0(:)
                              ! $ \v_0 $ [m-1]
                              ! 波数. Vacuume wavenumber (input HITRAN)
    real(8), allocatable :: Int_Ref(:)
                              ! $ \S_0 $ [m molecule-1]
                              ! 吸収線強度. Intensity at 296[K] (input HITRAN)
    real(8), allocatable :: EinsA(:)
                              ! $ \A $ [s-1]
                              ! Einstein A-coefficient (input HITRAN)
    real(8), allocatable :: AirBroadHarf(:)
                              ! $ \gamma_{air} $ [m-1 Pa-1]
                              ! Air-broadened half-width at 296[K] (input HITRAN)
    real(8), allocatable :: SelfBroadHarf(:)
                              ! $ \gamma_{self} $ [m-1 Pa-1]
                              ! Self-broadened hafl-width at 296[K] (input HITRAN)
    real(8), allocatable :: LowStatE(:)
                              ! $ \E^" $ [m-1]
                              ! Lower-State energy (input HITRAN)
    real(8), allocatable :: N_air(:)
                              ! $ n_{air} $ [1]
                              ! Temperature dependence exponent for gamma_air (input HITRAN)
    real(8), allocatable :: delta_air(:)
                              ! $ \delta_{air} $ [m-1 Pa-1]
                              ! Air pressure-induced line shift (input HITRAN)
    real(8) :: Intensity
                              ! $ \S $ [m molecule-1]
                              ! Intensity
    real(8) :: Alpha_D
                              ! $ \alpha_D $ [m-1]
                              ! Harf-width of Doppler broadening
    real(8) :: Alpha_L
                              ! $ \alpha_L $ [m-1]
                              ! Harf-width of Pressure broadening
    real(8) :: WaveNumS
                              ! $ \v $ [m-1]
                              ! 波数. pressure shift wavenumber
    integer :: FileEnd(1:50)

    integer :: ParSumMol(1:50, 1:10)

    real(8) :: Qval(0:108, 1:2)

    real(8) :: Q_TT(1:50, 1:10) !Q_Tからこちらに変更する

    integer :: cutoff = 1 !1: 2500[m-1], 2: HarfWidth dependent

    integer :: k_iso

    integer :: k_min, k_max

    integer :: k

    integer :: k_ln !, k_g

    real(8) :: BroadHarf

    real(8) :: CalcWidth

!    real(8) :: OffSet = 2500.0_8

    ! Hitran data 読み込み用
    integer :: k_Hi
    integer :: fHi = 13
    character(31) :: filenameHi
    character( 1) :: SnumHi
    character( 2) :: DnumHi


    ! HITRAN dataファイルの行数: FileEndの設定
    if( DataBase == 'HITRAN1996') then   ! 2014/05/21
       FileEnd(1:50) = (/ 49444,  60802, 275133,  26174,    4477, & !  1: 5
        &                 48032,   6292,  15331,  38853,  100680, & !  6:10
        &                 11152, 165426,   8676,    107,     533, & ! 11:15
        &                   576,    237,   7230,    858,    2702, & ! 16:20
        &                 15565,    120,    772,   9355,    5444, & ! 21:25
        &                  1668,   4749,   2886,  54866,   11520, & ! 26:30
        &                  7151,   3388,  26963,      2,   32199, & ! 31:35
        &                     0,      0,      0,      0,       0, & ! 36:40
        &                     0,      0,      0,      0,       0, & ! 41:45
        &                     0,      0,      0,      0,       0 /) ! 46:50

    else if( DataBase == 'HITRAN2008') then
      FileEnd(1:50) = (/  69201, 314919, 409686,  47843,    4477, & !  1: 5
        &                290091,   6428, 105079,  58250,  104223, & !  6:10
        &                 29084, 487254,  31976,    107,     613, & ! 11:15
        &                  1293,    806,  11501,  29242,   37050, & ! 16:20
        &                 16276,    120,   4253, 196171,  126983, & ! 21:25
        &                 11340,  22402,  20099,  70601, 2889065, & ! 26:30
        &                 20788,  62684,  38804,      2,   32199, & ! 31:35
        &                  1206,   4358,  18378,  19897,   36911, & ! 36:40
        &                  3572,  60033,      0,      0,       0, & ! 41:45
        &                     0,      0,      0,      0,       0 /) ! 46:50

    else if( DataBase == 'HITEMP2010' ) then
      FileEnd( 1) = 114241164
      FileEnd( 2) =  11193608
      FileEnd( 5) =    113631
      FileEnd( 8) =    115610
      FileEnd(13) =     41557

    else if( DataBase == 'HITRAN2012' ) then
      FileEnd(1:50) = (/  224515, 471847, 422116,  47843,    4606, & !  1: 5
        &                 468013,  13975, 105079,  95121,  104223, & !  6:10
        &                  46392, 961962,  31979,  34376,   83691, & ! 11:15
        &                   8980,   4751,  11501,  29361,   44601, & ! 16:20
        &                  16276,   1268,   4253, 212496,  126983, & ! 21:25
        &                  20410,  49629,  22189, 184104, 2889065, & ! 26:30
        &                  54235,  62684,  38804,      2,   32199, & ! 31:35
        &                   1206,   4358,  18378,  19897,   36911, & ! 36:40
        &                   3572,  60033, 124126, 180332,    9146, & ! 41:45
        &                   2078,  10881,      0,      0,       0 /) ! 46:50
    else
      print *, 'invalid database, DataBase:', DataBase
      print *, 'stop (linebyline)'
      stop
    end if
 
    ! parsum.datのコラムと分子種の対応マトリックス
    ParSumMol(1:50, 1:10) = 0.0_8

    ParSumMol( 1, 1:6 ) = (/ 1, 2, 3, 4, 5, 6 /) !H2O
    ParSumMol( 2, 1:8 ) = (/ 7, 8, 9, 10, 11, 12, 13, 14 /) !CO2
    ParSumMol( 3, 1:10) = (/ 15, 16, 17, 18, 19, 20, 21, 22, 23, 24 /) !O3 !32まである
    ParSumMol( 4, 1:5 ) = (/ 33, 34, 35, 36, 37 /) !N2O
    ParSumMol( 5, 1:6 ) = (/ 38, 39, 40, 41, 42, 43 /) !CO2
    ParSumMol( 6, 1:3 ) = (/ 44, 45, 46 /) !CH4
    ParSumMol( 7, 1:3 ) = (/ 47, 48, 49 /) !O3
    ParSumMol( 8, 1:3 ) = (/ 50, 51, 52 /) !NO
    ParSumMol( 9, 1:2 ) = (/ 53, 54 /) !SO2
    ParSumMol(10, 1:1 ) = (/ 55 /) !NO2
    ParSumMol(11, 1:2 ) = (/ 56, 57 /) !NH3
    ParSumMol(12, 1:1 ) = (/ 58 /) !HNO3
    ParSumMol(13, 1:3 ) = (/ 59, 60, 61 /) !OH
    ParSumMol(14, 1:1 ) = (/ 62 /) !HF
    ParSumMol(15, 1:2 ) = (/ 63, 64 /) !HCl
    ParSumMol(16, 1:2 ) = (/ 65, 66 /) !HBr
    ParSumMol(17, 1:1 ) = (/ 67 /) !HI
    ParSumMol(18, 1:2 ) = (/ 68, 69 /) !ClO
    ParSumMol(19, 1:5 ) = (/ 70, 71, 72, 73, 74 /) !OCS
    ParSumMol(20, 1:3 ) = (/ 75, 76, 77 /) !H2CO
    ParSumMol(21, 1:2 ) = (/ 78, 79 /) !HOCl
    ParSumMol(22, 1:1 ) = (/ 80 /) !N2
    ParSumMol(23, 1:3 ) = (/ 81, 82, 83 /) !HCN
    ParSumMol(24, 1:2 ) = (/ 84, 85 /) !CH3Cl
    ParSumMol(25, 1:1 ) = (/ 86 /) !H2O2
    ParSumMol(26, 1:2 ) = (/ 87, 88 /) !C2H2
    ParSumMol(27, 1:1 ) = (/ 89 /) !C2H6
    ParSumMol(28, 1:1 ) = (/ 90 /) !PH3
    ParSumMol(29, 1:1 ) = (/ 91 /) !COF2
    ParSumMol(30, 1:1 ) = (/ 92 /) !SF6
    ParSumMol(31, 1:3 ) = (/ 93, 94, 95 /) !H2S
    ParSumMol(32, 1:1 ) = (/ 96 /) !HCOOH
    ParSumMol(33, 1:1 ) = (/ 97 /) !HO2
    ParSumMol(34, 1:1 ) = (/ 98 /) !O
    ParSumMol(35, 1:2 ) = (/ 99, 100 /) !ClONO2
    ParSumMol(36, 1:1 ) = (/ 101 /) !NO+
    ParSumMol(37, 1:2 ) = (/ 102, 103 /) !HOBr
    ParSumMol(38, 1:2 ) = (/ 104, 105 /) !C2H4
    ParSumMol(39, 1:1 ) = (/ 106 /) !CH3OH
    !ParSumMol(40, 1:2 ) = (/ /) !CH3Br
    ParSumMol(41, 1:1 ) = (/ 107 /) !CH3CN
    ParSumMol(42, 1:1 ) = (/ 108 /) !CF4

    ! MolValの設定
    ! 1: Earth, HITRANは地球の同位体比に規格化されている
    ! 2: 同位体比を適当に設定
    !
    call SetMolVal(1)

    ! Hitran dataの読み込み ======================Read HITRAN Start
    !
    allocate( MolNum       (1:FileEnd(k_mol)) )
    allocate( IsoNum       (1:FileEnd(k_mol)) )
    allocate( WaveNum0     (1:FileEnd(k_mol)) )
    allocate( Int_Ref      (1:FileEnd(k_mol)) )
    allocate( EinsA        (1:FileEnd(k_mol)) )
    allocate( AirBroadHarf (1:FileEnd(k_mol)) )
    allocate( SelfBroadHarf(1:FileEnd(k_mol)) )
    allocate( LowStatE     (1:FileEnd(k_mol)) )
    allocate( N_air        (1:FileEnd(k_mol)) )
    allocate( delta_air    (1:FileEnd(k_mol)) )

    if( k_mol < 10 ) then   ! 2014/04/15
      write(SnumHi, '(I1)') k_mol
      DnumHi = '0' // SnumHi
    else
      write(DnumHi, '(I2)') k_mol
    end if

    if( DataBase == 'HITRAN1996') then                     ! 2014/05/21
      open(fHi, file='./../data/'// DnumHi //'_HIT96.PAR')
    else if( DataBase == 'HITRAN2008') then
      open(fHi, file='./../data/'// DnumHi //'_hit08.par')
    else if( DataBase == 'HITEMP2010') then
      open(fHi, file='./../data/'// DnumHi //'_HITEMP2010.par')
    else if( DataBase == 'HITRAN2012') then
      open(fHi, file='./../data/'// DnumHi //'_hit12.par')
    end if
    !open (fHi, file='./data/01_0-50_HITEMP2010.par') ! 25cm-1 check 20131225
    !print *, 'start read HITRAN'
    do k_Hi = 1, FileEnd(k_mol)
      read (fHi, '(I2, I1, F12.6, E10.3, E10.3, F5.4, F5.4, F10.4, F4.2, F8.6)') &
        & MolNum(k_Hi), IsoNum(k_Hi), WaveNum0(k_Hi), Int_Ref(k_Hi), EinsA(k_Hi), &
        & AirBroadHarf(k_Hi), SelfBroadHarf(k_Hi), LowStatE(k_Hi), N_air(k_Hi), delta_air(k_Hi)
      !print *, '====k_Hi:', k_Hi, '===='
      !print *, 'Mol', MolNum(k_Hi)
      !print *, 'Iso', IsoNum(k_Hi)
      !print *, 'WN', WaveNum0(k_Hi)
      !print *, 'Int', Int_Ref(k_Hi)
      !print *, 'Eins', EinsA(k_Hi)
      !print *, 'AirB', AirBroadHarf(k_Hi)
      !print *, 'SelfB', SelfBroadHarf(k_Hi)
      !print *, 'LowE', LowStatE(k_Hi)
      !print *, 'Nair', N_air(k_Hi)
      !print *, 'delta', delta_air(k_Hi)
      ! 単位をSIに変換
      WaveNum0(k_Hi) = WaveNum0(k_Hi) * 100.0_8
      Int_Ref(k_Hi) = Int_Ref(k_Hi) * 0.010_8
      AirBroadHarf(k_Hi) = AirBroadHarf(k_Hi) * 9.869233e-4_8 ! 1 [cm-1 atm-1] = 10^2/101325 [m-1 Pa-1]
      SelfBroadHarf(k_Hi) = SelfBroadHarf(k_Hi) * 9.869233e-4_8
      LowStatE(k_Hi) = LowStatE(k_Hi) * 100.0_8
      delta_air(k_Hi) = delta_air(k_Hi) * 9.869233e-4_8
    end do
    close(fHi) !=================================Read HITRAN end

    !!
    ! Q_Tの計算をこちらに移動
    call ReadParSum( Temp, Qval )

    ! Q_Tの計算
    !
    do k_iso = 1, 10
      if(  ( MolVal(k_mol, k_iso, 1) .ne. 0.0_8 ) .and. &
        &  ( ParSumMol(k_mol, k_iso) .ne. 0.0_8 ) ) then   ! parsumdatに無いものも除く
        Q_TT(k_mol, k_iso) = func_LinearInterPol( &
          & Qval(0, 1), Qval(0, 2), Qval(ParSumMol(k_mol, k_iso), 1), &
          & Qval(ParSumMol(k_mol, k_iso), 2), Temp )
      end if
      !print *, Q_TT(k_mol, k_iso)
    end do ! k_iso

    ! print *, "mixratioH=/0 kmol=", k_mol
    ! 波数範囲にあるデータを読み込む
    !
    AbsorpG(:) = 0.0_8 ! 2014/05/15
    do k_ln = 1, FileEnd(k_mol)

      if ( ( MolVal(k_mol, IsoNum(k_ln), 1) .ne. 0.0_8 ) .and. &
        &  ( ParSumMol(k_mol, IsoNum(k_ln)) .ne. 0.0_8 ) ) then   ! parsumdatに無いものも除く
        !print *, 'call calc_absorption'
        call Calc_absorption( kmax, MinWaveNum, ResWaveNum, k_mol, Press, Temp, MixRatioH, &
          & WaveNum0(k_ln), Int_Ref(k_ln), AirBroadHarf(k_ln), SelfBroadHarf(k_ln), LowStatE(k_ln), &
          & N_air(k_ln), delta_air(k_ln), MolVal(k_mol, IsoNum(k_ln), 1), &
          & MolVal(k_mol, IsoNum(k_ln), 2), MolVal(k_mol, IsoNum(k_ln), 3), Q_TT(k_mol, IsoNum(k_ln)), &
          & AbsorpG(0:kmax))

      end if

    end do ! k_ln

    deallocate( MolNum        )
    deallocate( IsoNum        )
    deallocate( WaveNum0      )
    deallocate( Int_Ref       )
    deallocate( EinsA         )
    deallocate( AirBroadHarf  )
    deallocate( SelfBroadHarf )
    deallocate( LowStatE      )
    deallocate( N_air         )
    deallocate( delta_air     )

  end subroutine linebyline

  subroutine Calc_absorption( kmax, MinWaveNum, ResWaveNum, &
    & k_mol, Press, Temp, MixRatioH, &
    & WaveNum0, Int_Ref, AirBroadHarf, SelfBroadHarf, LowStatE, &
    & N_air, delta_air, MolVal_1, MolVal_2, MolVal_3, Q_TT, &
    & AbsorpG)

    integer, intent(in ) :: kmax
    real(8), intent(in ) :: MinWaveNum
    real(8), intent(in ) :: ResWaveNum
    integer, intent(in ) :: k_mol
    real(8), intent(in ) :: Press
    real(8), intent(in ) :: Temp
    real(8), intent(in ) :: MixRatioH
    real(8), intent(in ) :: WaveNum0
    real(8), intent(in ) :: Int_Ref
    real(8), intent(in ) :: AirBroadHarf
    real(8), intent(in ) :: SelfBroadHarf
    real(8), intent(in ) :: LowStatE
    real(8), intent(in ) :: N_air
    real(8), intent(in ) :: delta_air
    real(8), intent(in ) :: MolVal_1, MolVal_2, MolVal_3
    real(8), intent(in ) :: Q_TT
    real(8), intent(inout) :: AbsorpG(0:kmax)
    real(8) :: WaveNumS
    real(8) :: Alpha_D, Alpha_L
    real(8) :: BroadHarf
    real(8) :: Intensity
    real(8) :: CalcWidth
    integer :: cutoff = 1
    integer :: k_min, k_max
    integer :: k_w

    integer :: VOIGT = 2 ! SET voigt profile 2013/07/01
                         ! Humlicek を使う 2014/04/17

    ! pressure shift of wacuum wavenumber
    !
    WaveNumS = WaveNum0 + delta_air * Press * (Temp_Ref/ Temp) ! 温度依存性の導入 2013/05/01
    !print *, 'WNS', WaveNumS

    if(WaveNumS < 0.0_8) then           ! 2013/12/25
      AbsorpG(0:kmax) = AbsorpG(0:kmax) 
    else

      ! 半値幅の計算
      !
      Alpha_D = ( WaveNumS / LightSpeed ) * &
        & ( 2.0_8 * BoltzmannConst * Temp * AvogadroNum/ MolVal_3 )**0.50_8
        !print *, "AlphaD", Alpha_D

      BroadHarf = SelfBroadHarf * MixRatioH * MolVal_1 + &
        & AirBroadHarf * (1.0_8 - MixRatioH * MolVal_1)
      Alpha_L = BroadHarf * Press * ( Temp_Ref / Temp )**N_air
      !print *, "AlphaL", Alpha_L

      ! 計算する波数範囲をcalcWidthとする
      if( cutoff == 1 ) then
        CalcWidth = 2500.0 ! 2013/08/26

      elseif( cutoff == 2 ) then
        if( Alpha_D .gt. Alpha_L ) then
          CalcWidth = Alpha_D * 5.0_8
        else
          CalcWidth = Alpha_L * 100.0_8
        end if
      end if
 
      ! Intensityの計算
      !
      Intensity = func_IntensityTempVar( &
        & Temp, MolVal_2, Q_TT, &
        & Int_Ref, WaveNum0, LowStatE ) !WaveNum0
      !print *, "Temp", Temp, "LowStatE", LowStatE, "Intensity", Intensity, "N", N_air, "Q", Q_T, &
      !  & Intensity * func_Voigt( 0.0_8 , ( Alpha_L / Alpha_D )  ) / ( Alpha_D * PI**0.50_8 )
      ! 2013/06/19
   !   print *, "Intensity: ", Intensity
  
      ! スペクトルの計算 ! バグ修正 2013/05/01
      !
      if( nint((WaveNumS - CalcWidth - MinWaveNum) /ResWaveNum ) .le. 0 ) then
        k_min = 0
      else
        k_min = nint((WaveNumS - CalcWidth - MinWaveNum) /ResWaveNum )
      end if
      ! print *, "WaveNumS - CalcWidth - MinWaveNum", WaveNumS - CalcWidth - MinWaveNum, "k_min", k_min
      if( nint((WaveNumS + CalcWidth - MinWaveNum) /ResWaveNum ) .ge. kmax ) then
        k_max = kmax
      elseif( nint((WaveNumS + CalcWidth - MinWaveNum) /ResWaveNum ) .le. 0 ) then
        k_max = -1
      else
        k_max = nint((WaveNumS + CalcWidth - MinWaveNum) /ResWaveNum )
      end if
      ! print *, "WaveNumS + CalcWidth - MinWaveNum", WaveNumS + CalcWidth - MinWaveNum, "k_max", k_max

      ! AbsorpG(:) = 0.0_8 ! 2013/07/19 : これは間違い
      do k_w = k_min, k_max ! in kk, k_mol

        if( VOIGT == 1 ) then ! add func_LineShape 2013/07/08
          AbsorpG(k_w) = AbsorpG(k_w) + &
            & Intensity * func_LineShape( k_mol, MinWaveNum + ResWavenum * k_w - WaveNumS ) * &
            & func_Voigt( ( MinWaveNum + ResWavenum * k_w - WaveNumS )/ Alpha_D, &
            & ( Alpha_L / Alpha_D )  ) / ( Alpha_D * PI**0.50_8 )
        elseif( VOIGT == 2 ) then
          AbsorpG(k_w) = AbsorpG(k_w) + &
            & Intensity * func_LineShape( k_mol, MinWaveNum + ResWavenum * k_w - WaveNumS ) * &
            & real(Humlicek( ( MinWaveNum + ResWavenum * k_w - WaveNumS )/ Alpha_D, &
            & ( Alpha_L / Alpha_D )  ) / ( Alpha_D * PI**0.50_8 ))
        end if

      end do ! k_w in kk, k_mol

    end if

  end subroutine Calc_absorption

  function func_LineShape( k_mol, delWaveNum ) result(a) ! 2013/07/08

    integer :: k_mol
    real(8) :: delWaveNum, a
    integer :: LineShape ! 1: LineShape function is 1.0, 2: line tale is 0
    real(8) :: WaveNumCutOff = 2500.0_8
    real(8) :: Coeff

    Coeff = delWaveNum/WaveNumCutOff

    LineShape = 1 ! SET

    if(LineShape == 1) then
      a = 1.0
    elseif(LineShape == 2) then
      if(k_mol == 2) then
        a = 1.0_8 - Coeff * Coeff * ( Coeff * Coeff - 2.0_8 )
      else
        a = 1.0_8 - Coeff * Coeff
      end if
    end if

  end function func_LineShape

  function func_IntensityTempVar( Temp, Q_Ref, Q_T, Int_Ref, WaveNum0, LowStatE ) result(a)

    real(8) :: Temp, Q_Ref, Q_T, Int_Ref, WaveNum0, LowStatE, a

    a = Int_Ref * exp( ( SecRadConst * LowStatE ) * ( 1.0_8 / Temp_Ref - 1.0_8 / Temp ) ) &
      & * ((1.0_8 -exp(-SecRadConst * WaveNum0 / Temp     )) / &
      &    (1.0_8 -exp(-SecRadConst * WaveNum0 / Temp_Ref )) ) &
      & * ( Q_Ref / Q_T )

  end function func_IntensityTempVar

  function func_LinearInterPol( x1, x2, y1, y2, x0 ) result(a)

    real(8) :: x1, x2, y1, y2, x0, a

    if ( (x2-x1) .ne. 0.0_8 ) then
      a = ( y2 - y1 ) * ( x0 - x1 ) / ( x2 - x1 ) + y1
    else
      print *, "値が不適切です"
    end if

  end function func_LinearInterPol

  !==voigt function
  function Humlicek(X, Y) result(a)  ! 2013/07/01
    real(8) :: X, Y
    real(8) :: S
    complex :: T, U
    complex :: a

    T = cmplx(Y, -X)
    S = ABS(X) + Y
    if(S .LT. 15. ) GOTO 1

      a = T*.5641896/(.5 + T*T)
      return
1   if(S .LT. 5.5) GOTO 2
      U = T*T
      a = T*(1.410474 + U*.5641896)/(.75+U*(3.+U))
      return
2   if(Y .LT. .195*ABS(X) -.176) GOTO 3
      a = (16.4955 + T*(20.20933 + T*(11.96482 + T*(3.778987 + T*.5642236))))/&
      & (16.4955 + T*(38.82363 + T*(39.27121 + T*(21.69274 + T*(6.699398+T)))))
      return
3   U = T*T
      a = CEXP(U) - T*(36183.31 - U*(3321.9905-U*(1540.787-U*(219.0313-U*(35.76683-U*(1.320522-U*.56419))))))/&
      & (32066.6-U*(24322.84-U*(9022.228-U*(2186.181-U*(364.2191-U*(61.57037-U*(1.841439-U)))))))
      return
  end function Humlicek  ! 2013/07/01

  function func_Voigt( x, y ) result(a)
    !
    real(8) :: x, y, a
    real(8) :: a1, b1, a2, b2
    real(8) :: a3, b3, c3, d3, a4, b4, c4, d4
    real(8) :: a5, b5, c5, d5, e5, a6, b6, c6, d6, e6
    real(8) :: a7, b7, c7, d7, e7, f7, g7, h7, o7, p7, q7, r7, s7, t7
    real(8) :: a8, b8, c8, d8, e8, f8, g8, h8, o8, p8, q8, r8, s8, t8 

    if ( ( abs(x) + y ) .ge. 15.0_8 ) then         !>=確認
      a1 = 0.2820948_8 * y + 0.5641896 * y**3
      b1 = 0.5641896_8 * y
      a2 = 0.5_8 + y**2 + y**4
      b2 = -1.0_8 + 2.0_8 * y**2

      a = ( a1 + b1 * x**2) / ( a2 + b2 * x**2 + x**4 )
      
    elseif ( ( ( abs(x) + y ) .ge. 5.5_8 ) .and. &      !>=確認
      &  ( ( abs(x) + y ) .lt. 15.0_8)       ) then
      a3 =   1.05786_8 * y + 4.65456_8 * y**3 + 3.10304_8 * y**5 + 0.56419_8 * y**7
      b3 =     2.962_8 * y + 0.56419_8 * y**3 + 1.69257_8 * y**5
      c3 = - 2.53885_8 * y + 1.69257_8 * y**3
      d3 =   0.56419_8 * y
      a4 = 0.5625_8 + 4.5_8 * y**2 + 10.5_8 * y**4 + 6.0_8 * y**6 + y**8
      b4 =   -4.5_8 + 9.0_8 * y**2 +  6.0_8 * y**4 + 4.0_8 * y**6
      c4 =   10.5_8 - 6.0_8 * y**2 +  6.0_8 * y**4
      d4 =   -6.0_8 + 4.0_8 * y**2

      a = ( a3 + b3 * x**2 + c3 * x**4 + d3 * x**6 ) / &
        & ( a4 + b4 * x**2 + c4 * x**4 + d4 * x**6 + x**8 )

    elseif ( ( ( abs(x) + y ) .lt. 5.5_8 ) .and. &
      &  ( y .gt. ( 0.195_8 * abs(x) - 0.176_8 ) ) ) then
      a5 = 272.102_8 + 973.778_8 * y + 1629.76_8 * y**2 + 1678.33_8 * y**3 + 1174.8_8 * y**4 + &
        & 581.746 * y**5 + 204.501_8 * y**6 + 49.5213_8 * y**7 + 7.55895_8 * y**8 + 0.564224_8 * y**9
      b5 = -60.5644_8 - 2.34403_8 * y + 220.843 * y**2 + 336.364_8 * y**3 + 247.198_8 * y**4 + &
        & 100.705_8 * y**5 + 22.6778_8 * y**6 + 2.25689_8 * y**7
      c5 = 4.58029_8 + 18.546_8 * y + 42.5683_8 * y**2 + 52.8454_8 * y**3 + 22.6798_8 * y**4 + &
        & 3.38534_8 * y**5
      d5 = -0.128922_8 + 1.66203_8 * y + 7.56186_8 * y**2 + 2.25689_8 * y**3
      e5 = 0.000971457_8 + 0.564224_8 * y
      a6 = 272.102_8 + 1280.83_8 * y + 2802.87_8 * y**2 + 3764.97 * y**3 + 3447.63_8 * y**4 + &
        & 2256.98_8 * y**5 + 1074.41 * y**6 + 369.199_8 * y**7 + 88.2674_8 * y**8 + 13.3988_8 * y**9 + y**10
      b6 = 211.678_8 + 902.306_8 * y + 1758.34_8 * y**2 + 2037.31_8 * y**3 + 1549.68_8 * y**4 + &
        & 793.427_8 * y**5 + 266.299_8 * y**6 + 53.5952_8 * y**7 + 5.0_8 * y**8
      c6 = 78.866_8 + 308.186_8 * y + 497.302_8 * y**2 + 479.258_8 * y**3 + 269.292_8 * y**4 + &
        & 80.3928_8 * y**5 + 10.0_8 * y**6
      d6 = 22.0353_8 + 55.0293_8 * y + 92.7568_8 * y**2 + 53.5952_8 * y**3 + 10.0_8 * y**4
      e6 = 1.49645_8 + 13.3988_8 * y + 5.0_8 * y**2

      a = ( a5 + b5 * x**2 + c5 * x**4 + d5 * x**6 + e5 * x**8 ) / &
        & ( a6 + b6 * x**2 + c6 * x**4 + d6 * x**6 + e6 * x**8 + x**10 )

    elseif ( ( ( abs(x) + y ) .lt. 5.5_8 ) .and. &
      &  ( y .lt. ( 0.195_8 * abs(x) - 0.176_8 ) ) ) then
      a7 = 1.16028e+9_8 * y - 9.86604e+8_8 * y**3 + 4.56662e+8_8 * y**5 - 1.53575e+8_8 * y**7 + &
        & 4.08168e+7_8 * y**9 - 9.69463e+6_8 * y**11 + 1.6841e+6_8 * y**13 - 320772.0_8 * y**15 + &
        & 40649.2_8 * y**17 - 5860.68_8 * y**19 + 571.687_8 * y**21 - 72.9359_8 * y**23 + &
        & 2.35944_8 * y**25 - 0.56419_8 * y**27
      b7 = -5.60505e+8_8 * y - 9.85386e+8_8 * y**3 + 8.06985e+8_8 * y**5 - 2.91876e+8_8 * y**7 + &
        & 8.64829e+7_8 * y**9 - 7.72359e+6_8 * y**11 + 3.59915e+6_8 * y**13 - 234417.0_8 * y**15 + &
        & 45251.3_8 * y**17 - 2269.19_8 * y**19 - 234.143_8 * y**21 + 23.0312_8 * y**23 &
        & -7.33447_8 * y**25
      c7 = -6.51523e+8_8 * y + 2.47157e+8_8 + y**3 + 2.94262e+8_8 * y**5 - 2.04467e+8_8 * y**7 + &
        & 2.29302e+7_8 * y**9 - 2.3818e+7_8 * y**11 + 576054.0_8 * y**13 + 98079.1_8 * y**15 &
        & -25338.3_8 * y**17 + 1097.77_8 * y**19 + &
        & 97.6203_8 * y**21 - 44.0068_8 * y**23
      d7 = -2.63894e+8_8 * y + 2.70167e+8_8 * y**3 - 9.96224e+7_8 * y**5 - 4.15013e+7_8 * y**7 + &
        & 3.83112e+7_8 * y**9 + 2.2404e+6_8 * y**11 - 303569.0_8 * y**13 - 66431.2_8 * y**15 + &
        & 8381.97_8 * y**17 + 228.563_8 * y**19 - 161.358_8 * y**21
      e7 = -6.31771e+7_8 * y + 1.40677e+8_8 * y**3 + 5.56965e+6_8 * y**5 + 2.46201e+7_8 * y**7 + &
        & 468142.0_8 * y**9 - 1.003e+6_8 * y**11 - 66212.1_8 * y**13 + 23507.6_8 * y**15 + &
        & 296.38_8 * y**17 - 403.396_8 * y**19
      f7 = -1.69846e+7 * y + 4.07382e+6_8 * y**3 - 3.32896e+7_8 * y**5 - 1.93114e+6_8 * y**7 &
        & -934717.0_8 * y**9 + 8820.94_8 * y**11 + 37544.8_8 * y**13 + 125.591_8 * y**15 &
        & -726.113_8 * y**17
      g7 = -1.23165e+6_8 * y + 7.52883e+6_8 * y**3 - 900010.0_8 * y**5 - 186682.0_8 * y**7 + 79902.5_8 * y**9 + &
        & 37371.9_8 * y**11 - 260.198_8 * y**13 - 968.15_8 * y**15
      h7 = -610622.0_8 * y + 86407.6_8 * y**3 + 153468.0_8 * y**5 + 72520.9_8 * y**7 + 23137.1_8 * y**9 &
        & -571.645_8 * y**11 - 968.15_8 * y**13
      o7 = -23586.5 * y + 49883.8_8 * y**3 + 26538.5_8 * y**5 + 8073.15_8 * y**7 - 575.164_8 * y**9 &
        & -726.113 * y**11
      p7 = -8009.1_8 * y + 2198.86_8 * y**3 + 953.655_8 * y**5 - 352.467_8 * y**7 - 403.396_8 * y**9
      q7 = -622.056_8 * y - 271.202_8 * y**3 - 134.792_8 * y**5 - 161.358_8 * y**7
      r7 = -77.0535_8 * y - 29.7896_8 * y**3 -44.0068_8 * y**5
      s7 = -2.92264_8 * y - 7.33447_8 * y**3
      t7 = -0.56419_8 * y

      a8 = 1.02827e+9_8 - 1.5599e+9_8 * y**2 + 1.17022e+9_8 * y**4 -5.79099e+8_8 * y**6 + 2.11107e+8_8 * y**8 &
        & -6.11148e+7_8 * y**10 + 1.44647e+7_8 * y**12 - 2.85721e+6_8 * y**14 + 483737.0_8 * y**16 - 70946.1_8 * y**18 + &
        & 9504.65_8 * y**20 - 955.194_8 * y**22 + 126.532_8 * y**24 - 3.68288_8 * y**26 + 1.0_8 * y**28
      b8 = 1.5599e+9_8 - 2.28855e+9_8 * y**2 + 1.66421e+9_8 * y**4 - 7.53828e+8_8 * y**6 + 2.89676e+8_8 * y**8 &
        & -7.01358e+7_8 * y**10 + 1.39465e+7_8 * y**12 - 2.84954e+6_8 * y**14 + 498334.0_8 * y**16 - 55600.0_8 * y**18 + &
        & 3058.26_8 * y**20 + 533.254_8 * y**22 - 40.5117_8 * y**24 + 14.0_8 * y**26
      c8 = 1.17022e+9_8 - 1.66421e+9_8 * y**2 + 1.06002e+9_8 * y**4 - 6.60078e+8_8 * y**6 + 6.33496e+7_8 * y**8 &
        & - 4.60396e+7_8 * y**10 + 1.4841e+7_8 * y**12 - 1.06352e+6_8 * y**14 - 217801.0_8 * y**16 + 48153.3_8 * y**18 &
        & - 1500.17_8 * y**20 - 198.876_8 * y**22 + 91.0_8 * y**24
      d8 = 5.79099e+8_8 - 7.53828e+8_8 * y**2 + 6.60078e+8_8 * y**4 + 5.40367e+7_8 * y**6 + 1.99846e+8_8 * y**8 &
        & - 6.87656e+6_8 * y**10 - 6.89002e+6_8 * y**12 + 280428.0_8 * y**14 + 161461.0_8 * y**16 - 16493.7_8 * y**18 &
        & - 567.164_8 * y**20 + 364.0_8 * y**22
      e8 = 2.11107e+8_8 - 2.89676e+8_8 * y**2 + 6.33496e+7_8 * y**4 - 1.99846e+8_8 * y**6 - 5.01017e+7_8 * y**8 &
        & - 5.25722e+6_8 * y**10 + 1.9547e+6_8 * y**12 + 240373.0_8 * y**14 - 55582.0_8 * y**16 - 1012.79_8 * y**18 &
        & + 1001.0_8 * y**20
      f8 = 6.11148e+7_8 - 7.01358e+7_8 * y**2 + 4.60396e+7_8 * y**4 - 6.87656e+6_8 * y**6 + 5.25722e+6_8 * y**8 &
        & + 3.04316e+6_8 * y**10 + 123052.0_8 * y**12 - 106663.0_8 * y**14 - 1093.82_8 * y**16 + 2002.0_8 * y**18
      g8 = 1.44647e+7_8 - 1.39465e+7_8 * y**2 + 1.4841e+7_8 * y**4 + 6.89002e+6_8 * y**6 + 1.9547e+6 * y**8 &
        & - 123052.0_8 * y**10 - 131337.0_8 * y**12 - 486.14_8 * y**14 + 3003.0_8 * y**16
      h8 = 2.85721e+6_8 - 2.84954e+6_8 * y**2 + 1.06352e+6_8 * y**4 + 280428.0_8 * y**6 - 240373.0_8 * y**8 &
        & - 106663.0_8 * y**10 + 486.14_8 * y**12 + 3432.0_8 * y**14
      o8 = 483737.0_8 - 498334.0_8 * y**2 - 217801.0_8 * y**4 - 161461.0_8 * y**6 - 55582.0_8 * y**8 &
        & + 1093.82_8 * y**10 + 3003.0_8 * y**12
      p8 = 70946.1_8 - 55600.0_8 * y**2 - 48153.3_8 * y**4 - 16493.7_8 * y**6 + 1012.79_8 * y**8 &
        & + 2002.0_8 * y**10
      q8 = 9504.65_8 - 3058.26_8 * y**2 - 1500.17_8 * y**4 + 567.164_8 * y**6 + 1001.0_8 * y**8
      r8 = 955.194_8 + 533.254_8 * y**2 + 198.876_8 * y**4 + 364.0_8 * y**6
      s8 = 126.532_8 + 40.5117_8 * y**2 + 91.0_8 * y**4
      t8 = 3.68288_8 + 14.0_8 * y**2

      a = exp( y**2 - x**2 ) * cos( 2.0_8 * x * y ) - &
        & (a7 + b7 * x**2 + c7 * x**4 + d7 * x**6 + e7 * x**8 + f7 * x**10 + g7 * x**12 + h7 * x**14 + &
        &  o7 * x**16 + p7 * x**18 + q7 * x**20 + r7 * x**22 + s7 * x**24 + t7 * x**26 ) / &
        & (a8 + b8 * x**2 + c8 * x**4 + d8 * x**6 + e8 * x**8 + f8 * x**10 + g8 * x**12 + h8 * x**14 + &
        &  o8 * x**16 + p8 * x**18 + q8 * x**20 + r8 * x**22 + s8 * x**24 + t8 * x**26 + x**28 )

    end if

  end function func_Voigt

  subroutine ReadParSum( Temp, Qval )

    real(8), intent(in ) :: Temp
    real(8), intent(out) :: Qval(0:108, 1:2)
    integer :: count, ksr
    integer :: fi = 14

    !print *, 'Temp', Temp

    if( (Temp < 70.0_8) .or. (Temp > 3000.0_8) ) then ! 2014/04/15
      print *, 'Temperature is out of range(70-3000[K]), T:', Temp
      print *, 'stop (ReadParSum)'
      stop
    end if

    count = int(Temp - 70.0_8)

    Qval = 0.0_8

    open(fi, file = './../data/parsum.dat') ! 2013/08/05

    do ksr = 1, count + 1
      read(fi, *)
      !read(fi, '()')  !2014/05/20 No Good !
    end do

    if( Temp == 3000.0_8 ) then ! 2014/04/15
      read(fi, *) Qval(:, 1)
      Qval(:, 2) = 0.0_8 
    else
      read(fi, *) Qval
    end if

    close(fi)

  end subroutine ReadParSum

  subroutine SetMolVal(setnum)
    integer :: setnum
  ! 分子の存在比, Q_296[K], 分子の質量: MolValの設定

    if( setnum == 1 ) then
  ! H2O (1)
  MolVal(1,1,1:3) = (/ 1.0_8, 1.7464E+02_8, 18.010565E-3_8/)
  MolVal(1,2,1:3) = (/ 1.0_8, 1.7511E+02_8, 20.014811E-3_8/)
  MolVal(1,3,1:3) = (/ 1.0_8, 1.0479E+03_8, 19.014780E-3_8/)
  MolVal(1,4,1:3) = (/ 1.0_8, 8.5901E+02_8, 19.016740E-3_8/)
  MolVal(1,5,1:3) = (/ 1.0_8, 8.7519E+02_8, 21.020985E-3_8/)
  MolVal(1,6,1:3) = (/ 1.0_8, 5.2204E+03_8, 20.020956E-3_8/)
  ! CO2 (2)
  MolVal(2,1,1:3) = (/ 1.0_8, 2.8694E+02_8, 43.989830E-3_8/)
  MolVal(2,2,1:3) = (/ 1.0_8, 5.7841E+02_8, 44.993185E-3_8/)
  MolVal(2,3,1:3) = (/ 1.0_8, 6.0948E+02_8, 45.994076E-3_8/)
  MolVal(2,4,1:3) = (/ 1.0_8, 3.5527E+03_8, 44.994045E-3_8/)
  MolVal(2,5,1:3) = (/ 1.0_8, 1.2291E+03_8, 46.997431E-3_8/)
  MolVal(2,6,1:3) = (/ 1.0_8, 7.1629E+03_8, 45.997400E-3_8/)
  MolVal(2,7,1:3) = (/ 1.0_8, 3.2421E+02_8, 47.998322E-3_8/)
  MolVal(2,8,1:3) = (/ 1.0_8, 3.7764E+03_8, 46.998291E-3_8/)
  MolVal(2,9,1:3) = (/ 1.0_8, 6.5350E+02_8, 49.001675E-3_8/)
  MolVal(2,10,1:3) = (/ 1.0_8, 7.4229E+03_8, 48.001646E-3_8/)
  ! O3 (3)
  MolVal(3,1,1:3) = (/ 1.0_8, 3.4838E+03_8, 47.984745E-3_8/)
  MolVal(3,2,1:3) = (/ 1.0_8, 7.4657E+03_8, 49.988991E-3_8/)
  MolVal(3,3,1:3) = (/ 1.0_8, 3.6471E+03_8, 49.988991E-3_8/)
  MolVal(3,4,1:3) = (/ 1.0_8, 4.3331E+04_8, 48.988960E-3_8/)
  MolVal(3,5,1:3) = (/ 1.0_8, 2.1405E+04_8, 48.988960E-3_8/)
  ! N2O (4)
  MolVal(4,1,1:3) = (/ 1.0_8, 5.0018E+03_8, 44.001062E-3_8/)
  MolVal(4,2,1:3) = (/ 1.0_8, 3.3619E+03_8, 44.998096E-3_8/)
  MolVal(4,3,1:3) = (/ 1.0_8, 3.4586E+03_8, 44.998096E-3_8/)
  MolVal(4,4,1:3) = (/ 1.0_8, 5.3147E+03_8, 46.005308E-3_8/)
  MolVal(4,5,1:3) = (/ 1.0_8, 3.0971E+04_8, 45.005278E-3_8/)
  ! CO (5)
  MolVal(5,1,1:3) = (/ 1.0_8, 1.0712E+02_8, 27.994915E-3_8/)
  MolVal(5,2,1:3) = (/ 1.0_8, 2.2408E+02_8, 28.998270E-3_8/)
  MolVal(5,3,1:3) = (/ 1.0_8, 1.1247E+02_8, 29.999161E-3_8/)
  MolVal(5,4,1:3) = (/ 1.0_8, 6.5934E+02_8, 28.999130E-3_8/)
  MolVal(5,5,1:3) = (/ 1.0_8, 2.3582E+02_8, 31.002516E-3_8/)
  MolVal(5,6,1:3) = (/ 1.0_8, 1.3809E+03_8, 30.002485E-3_8/)
  ! CH4 (6)
  MolVal(6,1,1:3) = (/ 1.0_8, 5.9045E+02_8, 16.031300E-3_8/)
  MolVal(6,2,1:3) = (/ 1.0_8, 1.1808E+03_8, 17.034655E-3_8/)
  MolVal(6,3,1:3) = (/ 1.0_8, 4.7750E+03_8, 17.037475E-3_8/)
  MolVal(6,4,1:3) = (/ 1.0_8, 1.5980E+03_8, 18.040830E-3_8/)
  ! O2 (7)
  MolVal(7,1,1:3) = (/ 1.0_8, 2.1577E+02_8, 31.989830E-3_8/)
  MolVal(7,2,1:3) = (/ 1.0_8, 4.5230E+02_8, 33.994076E-3_8/)
  MolVal(7,3,1:3) = (/ 1.0_8, 2.6406E+03_8, 32.994045E-3_8/)
  ! NO (8)
  MolVal(8,1,1:3) = (/ 1.0_8, 1.1421E+03_8, 29.997989E-3_8/)
  MolVal(8,2,1:3) = (/ 1.0_8, 7.8926E+02_8, 30.995023E-3_8/)
  MolVal(8,3,1:3) = (/ 1.0_8, 1.2045E+03_8, 32.002234E-3_8/)
  ! SO2 (9)
  MolVal(9,1,1:3) = (/ 1.0_8, 6.3403E+03_8, 63.961901E-3_8/)
  MolVal(9,2,1:3) = (/ 1.0_8, 6.3689E+03_8, 65.957695E-3_8/)
  ! NO2 (10)
  MolVal(10,1,1:3) = (/ 1.0_8, 1.3578E+04_8, 45.992904E-3_8/)
  ! NH3 (11)
  MolVal(11,1,1:3) = (/ 1.0_8, 1.7252E+03_8, 17.026549E-3_8/)
  MolVal(11,2,1:3) = (/ 1.0_8, 1.1527E+03_8, 18.023583E-3_8/)
  ! HNO3 (12)
  MolVal(12,1,1:3) = (/ 1.0_8, 2.1412E+05_8, 62.995644E-3_8/)
  ! OH (13)
  MolVal(13,1,1:3) = (/ 1.0_8, 8.0362E+01_8, 17.002740E-3_8/)
  MolVal(13,2,1:3) = (/ 1.0_8, 8.0882E+01_8, 19.006986E-3_8/)
  MolVal(13,3,1:3) = (/ 1.0_8, 2.0931E+02_8, 18.008915E-3_8/)
  ! HF (14)
  MolVal(14,1,1:3) = (/ 1.0_8, 4.1466E+01_8, 20.006229E-3_8/)
  ! HCl (15)
  MolVal(15,1,1:3) = (/ 1.0_8, 1.6066E+02_8, 35.976678E-3_8/)
  MolVal(15,2,1:3) = (/ 1.0_8, 1.6089E+02_8, 37.973729E-3_8/)
  ! HBr (16)
  MolVal(16,1,1:3) = (/ 1.0_8, 2.0018E+02_8, 79.926160E-3_8/)
  MolVal(16,2,1:3) = (/ 1.0_8, 2.0024E+02_8, 81.924115E-3_8/)
  ! HI (17)
  MolVal(17,1,1:3) = (/ 1.0_8, 3.8900E+02_8, 127.912297E-3_8/)
  ! ClO (18)
  MolVal(18,1,1:3) = (/ 1.0_8, 3.2746E+03_8, 50.963768E-3_8/)
  MolVal(18,2,1:3) = (/ 1.0_8, 3.3323E+03_8, 52.960819E-3_8/)
  ! OCS (19)
  MolVal(19,1,1:3) = (/ 1.0_8, 1.2210E+03_8, 59.966986E-3_8/)
  MolVal(19,2,1:3) = (/ 1.0_8, 1.2535E+03_8, 61.962780E-3_8/)
  MolVal(19,3,1:3) = (/ 1.0_8, 2.4842E+03_8, 60.970341E-3_8/)
  MolVal(19,4,1:3) = (/ 1.0_8, 4.9501E+03_8, 60.966371E-3_8/)
  MolVal(19,5,1:3) = (/ 1.0_8, 1.3137E+03_8, 61.971231E-3_8/)
  ! H2CO (20)
  MolVal(20,1,1:3) = (/ 1.0_8, 2.8467E+03_8, 30.010565E-3_8/)
  MolVal(20,2,1:3) = (/ 1.0_8, 5.8376E+03_8, 31.013920E-3_8/)
  MolVal(20,3,1:3) = (/ 1.0_8, 2.9864E+03_8, 32.014811E-3_8/)
  ! HOCl (21)
  MolVal(21,1,1:3) = (/ 1.0_8, 1.9274E+04_8, 51.971593E-3_8/)
  MolVal(21,2,1:3) = (/ 1.0_8, 1.9616E+04_8, 53.968644E-3_8/)
  ! N2 (22)
  MolVal(22,1,1:3) = (/ 1.0_8, 4.6598E+02_8, 28.006147E-3_8/)
  ! HCN (23)
  MolVal(23,1,1:3) = (/ 1.0_8, 8.9529E+02_8, 27.010899E-3_8/)
  MolVal(23,2,1:3) = (/ 1.0_8, 1.8403E+03_8, 28.014254E-3_8/)
  MolVal(23,3,1:3) = (/ 1.0_8, 6.2141E+02_8, 28.007933E-3_8/)
  ! CH3Cl (24)
  MolVal(24,1,1:3) = (/ 1.0_8, 1.1583E+05_8, 49.992328E-3_8/)
  MolVal(24,2,1:3) = (/ 1.0_8, 1.1767E+05_8, 51.989379E-3_8/)
  ! H2O2 (25)
  MolVal(25,1,1:3) = (/ 1.0_8, 9.8198E+03_8, 34.005480E-3_8/)
  ! C2H2 (26)
  MolVal(26,1,1:3) = (/ 1.0_8, 4.1403E+02_8, 26.015650E-3_8/)
  MolVal(26,2,1:3) = (/ 1.0_8, 1.6562E+03_8, 27.019005E-3_8/)
  MolVal(26,3,1:3) = (/ 1.0_8, 1.5818E+03_8, 27.021825E-3_8/)
  ! C2H6 (27)
  MolVal(27,1,1:3) = (/ 1.0_8, 7.0881E+04_8, 30.046950E-3_8/)
  MolVal(27,2,1:3) = (/ 1.0_8, 8.0001E+04_8, 31.050305E-3_8/)
  ! PH3 (28)
  MolVal(28,1,1:3) = (/ 1.0_8, 3.2486E+03_8, 33.997238E-3_8/)
  ! COF2 (29)
  MolVal(29,1,1:3) = (/ 1.0_8, 7.0044E+04_8, 65.991722E-3_8/)
  ! SF6 (30)
  MolVal(30,1,1:3) = (/ 1.0_8, 1.6233E+06_8, 145.962492E-3_8/)
  ! H2S (31)
  MolVal(31,1,1:3) = (/ 1.0_8, 5.0307E+02_8, 33.987721E-3_8/)
  MolVal(31,2,1:3) = (/ 1.0_8, 5.0435E+02_8, 35.983515E-3_8/)
  MolVal(31,3,1:3) = (/ 1.0_8, 2.0149E+03_8, 34.987105E-3_8/)
  ! HCOOH (32)
  MolVal(32,1,1:3) = (/ 1.0_8, 3.9133E+04_8, 46.005480E-3_8/)
  ! HO2 (33)
  MolVal(33,1,1:3) = (/ 1.0_8, 4.3004E+03_8, 32.997655E-3_8/)
  ! O (34)
  MolVal(34,1,1:3) = (/ 1.0_8, 6.7212E+00_8, 15.994915E-3_8/)
  ! ClONO2 (35)
  MolVal(35,1,1:3) = (/ 1.0_8, 4.7884E+06_8, 96.956672E-3_8/)
  MolVal(35,2,1:3) = (/ 1.0_8, 4.9102E+06_8, 98.953723E-3_8/)
  ! NO+ (36)
  MolVal(36,1,1:3) = (/ 1.0_8, 3.1168E+02_8, 29.997989E-3_8/)
  ! HOBr (37)
  MolVal(37,1,1:3) = (/ 1.0_8, 2.8339E+04_8, 95.921076E-3_8/)
  MolVal(37,2,1:3) = (/ 1.0_8, 2.8238E+04_8, 97.919027E-3_8/)
  ! C2H4 (38)
  MolVal(38,1,1:3) = (/ 1.0_8, 1.1041E+04_8, 28.031300E-3_8/)
  MolVal(38,2,1:3) = (/ 1.0_8, 4.5197E+04_8, 29.034655E-3_8/)
  ! CH3OH (39)
  MolVal(39,1,1:3) = (/ 1.0_8, 3.5314E+04_8, 32.026215E-3_8/)
  ! CH3Br (40)
  MolVal(40,1,1:3) = (/ 1.0_8, 2.0762E+04_8, 93.941811E-3_8/)
  MolVal(40,2,1:3) = (/ 1.0_8, 2.0849E+04_8, 95.939764E-3_8/)
  ! CH3CN (41)
  MolVal(41,1,1:3) = (/ 1.0_8, 8.8670E+04_8, 41.026549E-3_8/)
  ! CF4 (42)
  MolVal(42,1,1:3) = (/ 1.0_8, 1.2127E+05_8, 87.993616E-3_8/)
  ! CS (46)
  MolVal(46,1,1:3) = (/ 1.0_8, 2.5653E+02_8, 43.971036E-3_8/)
  MolVal(46,2,1:3) = (/ 1.0_8, 2.6069E+02_8, 45.966787E-3_8/)
  MolVal(46,3,1:3) = (/ 1.0_8, 5.4345E+02_8, 44.974368E-3_8/)
  MolVal(46,4,1:3) = (/ 1.0_8, 1.0346E+03_8, 44.970399E-3_8/)
  end if

    if( setnum == 2 ) then
  ! H2O (1)
  MolVal(1,1,1:3) = (/0.997317E+00_8, 1.7464E+02_8, 18.010565E-3_8/)
  MolVal(1,2,1:3) = (/ 1.99983E-03_8, 1.7511E+02_8, 20.014811E-3_8/)
  MolVal(1,3,1:3) = (/ 3.71884E-04_8, 1.0479E+03_8, 19.014780E-3_8/)
  MolVal(1,4,1:3) = (/ 3.10693E-04_8, 8.5901E+02_8, 19.016740E-3_8/)
  MolVal(1,5,1:3) = (/ 6.23003E-07_8, 8.7519E+02_8, 21.020985E-3_8/)
  MolVal(1,6,1:3) = (/ 1.15853E-07_8, 5.2204E+03_8, 20.020956E-3_8/)
  ! CO2 (2)
  MolVal(2,1,1:3) = (/ .984204E+00_8, 2.8694E+02_8, 43.989830E-3_8/)
  MolVal(2,2,1:3) = (/ 1.10574E-02_8, 5.7841E+02_8, 44.993185E-3_8/)
  MolVal(2,3,1:3) = (/ 3.94707E-03_8, 6.0948E+02_8, 45.994076E-3_8/)
  MolVal(2,4,1:3) = (/ 7.33989E-04_8, 3.5527E+03_8, 44.994045E-3_8/)
  MolVal(2,5,1:3) = (/ 4.43446E-05_8, 1.2291E+03_8, 46.997431E-3_8/)
  MolVal(2,6,1:3) = (/ 8.24623E-06_8, 7.1629E+03_8, 45.997400E-3_8/)
  MolVal(2,7,1:3) = (/ 3.95734E-06_8, 3.2421E+02_8, 47.998322E-3_8/)
  MolVal(2,8,1:3) = (/ 1.47180E-06_8, 3.7764E+03_8, 46.998291E-3_8/)
  MolVal(2,9,1:3) = (/ 4.44600E-08_8, 6.5350E+02_8, 49.001675E-3_8/)
  MolVal(2,10,1:3) = (/ 1.65354E-08_8, 7.4229E+03_8, 48.001646E-3_8/)
  ! O3 (3)
  MolVal(3,1,1:3) = (/ .992901E+00_8, 3.4838E+03_8, 47.984745E-3_8/)
  MolVal(3,2,1:3) = (/ 3.98194E-03_8, 7.4657E+03_8, 49.988991E-3_8/)
  MolVal(3,3,1:3) = (/ 1.99097E-03_8, 3.6471E+03_8, 49.988991E-3_8/)
  MolVal(3,4,1:3) = (/ 7.40475E-04_8, 4.3331E+04_8, 48.988960E-3_8/)
  MolVal(3,5,1:3) = (/ 3.70237E-04_8, 2.1405E+04_8, 48.988960E-3_8/)
  ! N2O (4)
  MolVal(4,1,1:3) = (/ .990333E+00_8, 5.0018E+03_8, 44.001062E-3_8/)
  MolVal(4,2,1:3) = (/ 3.64093E-03_8, 3.3619E+03_8, 44.998096E-3_8/)
  MolVal(4,3,1:3) = (/ 3.64093E-03_8, 3.4586E+03_8, 44.998096E-3_8/)
  MolVal(4,4,1:3) = (/ 1.98582E-03_8, 5.3147E+03_8, 46.005308E-3_8/)
  MolVal(4,5,1:3) = (/ 3.69280E-04_8, 3.0971E+04_8, 45.005278E-3_8/)
  ! CO (5)
  MolVal(5,1,1:3) = (/ .986544E+00_8, 1.0712E+02_8, 27.994915E-3_8/)
  MolVal(5,2,1:3) = (/ 1.10836E-02_8, 2.2408E+02_8, 28.998270E-3_8/)
  MolVal(5,3,1:3) = (/ 1.97822E-03_8, 1.1247E+02_8, 29.999161E-3_8/)
  MolVal(5,4,1:3) = (/ 3.67867E-04_8, 6.5934E+02_8, 28.999130E-3_8/)
  MolVal(5,5,1:3) = (/ 2.22250E-05_8, 2.3582E+02_8, 31.002516E-3_8/)
  MolVal(5,6,1:3) = (/ 4.13292E-06_8, 1.3809E+03_8, 30.002485E-3_8/)
  ! CH4 (6)
  MolVal(6,1,1:3) = (/ .988274E+00_8, 5.9045E+02_8, 16.031300E-3_8/)
  MolVal(6,2,1:3) = (/ 1.11031E-02_8, 1.1808E+03_8, 17.034655E-3_8/)
  MolVal(6,3,1:3) = (/ 6.15751E-04_8, 4.7750E+03_8, 17.037475E-3_8/)
  MolVal(6,4,1:3) = (/ 6.91785E-06_8, 1.5980E+03_8, 18.040830E-3_8/)
  ! O2 (7)
  MolVal(7,1,1:3) = (/ .995262E+00_8, 2.1577E+02_8, 31.989830E-3_8/)
  MolVal(7,2,1:3) = (/ 3.99141E-03_8, 4.5230E+02_8, 33.994076E-3_8/)
  MolVal(7,3,1:3) = (/ 7.42235E-04_8, 2.6406E+03_8, 32.994045E-3_8/)
  ! NO (8)
  MolVal(8,1,1:3) = (/ .993974E+00_8, 1.1421E+03_8, 29.997989E-3_8/)
  MolVal(8,2,1:3) = (/ 3.65431E-03_8, 7.8926E+02_8, 30.995023E-3_8/)
  MolVal(8,3,1:3) = (/ 1.99312E-03_8, 1.2045E+03_8, 32.002234E-3_8/)
  ! SO2 (9)
  MolVal(9,1,1:3) = (/ .945678E+00_8, 6.3403E+03_8, 63.961901E-3_8/)
  MolVal(9,2,1:3) = (/ 4.19503E-02_8, 6.3689E+03_8, 65.957695E-3_8/)
  ! NO2 (10)
  MolVal(10,1,1:3) = (/ .991616E+00_8, 1.3578E+04_8, 45.992904E-3_8/)
  ! NH3 (11)
  MolVal(11,1,1:3) = (/ .995872E+00_8, 1.7252E+03_8, 17.026549E-3_8/)
  MolVal(11,2,1:3) = (/ 3.66129E-03_8, 1.1527E+03_8, 18.023583E-3_8/)
  ! HNO3 (12)
  MolVal(12,1,1:3) = (/ .989110E+00_8, 2.1412E+05_8, 62.995644E-3_8/)
  ! OH (13)
  MolVal(13,1,1:3) = (/ .997473E+00_8, 8.0362E+01_8, 17.002740E-3_8/)
  MolVal(13,2,1:3) = (/ 2.00014E-03_8, 8.0882E+01_8, 19.006986E-3_8/)
  MolVal(13,3,1:3) = (/ 1.55371E-04_8, 2.0931E+02_8, 18.008915E-3_8/)
  ! HF (14)
  MolVal(14,1,1:3) = (/ .999844E+00_8, 4.1466E+01_8, 20.006229E-3_8/)
  ! HCl (15)
  MolVal(15,1,1:3) = (/ .757587E+00_8, 1.6066E+02_8, 35.976678E-3_8/)
  MolVal(15,2,1:3) = (/ .242257E+00_8, 1.6089E+02_8, 37.973729E-3_8/)
  ! HBr (16)
  MolVal(16,1,1:3) = (/ .506781E+00_8, 2.0018E+02_8, 79.926160E-3_8/)
  MolVal(16,2,1:3) = (/ .493063E+00_8, 2.0024E+02_8, 81.924115E-3_8/)
  ! HI (17)
  MolVal(17,1,1:3) = (/ .999844E+00_8, 3.8900E+02_8, 127.912297E-3_8/)
  ! ClO (18)
  MolVal(18,1,1:3) = (/ .755908E+00_8, 3.2746E+03_8, 50.963768E-3_8/)
  MolVal(18,2,1:3) = (/ .241720E+00_8, 3.3323E+03_8, 52.960819E-3_8/)
  ! OCS (19)
  MolVal(19,1,1:3) = (/ .937395E+00_8, 1.2210E+03_8, 59.966986E-3_8/)
  MolVal(19,2,1:3) = (/ 4.15828E-02_8, 1.2535E+03_8, 61.962780E-3_8/)
  MolVal(19,3,1:3) = (/ 1.05315E-02_8, 2.4842E+03_8, 60.970341E-3_8/)
  MolVal(19,4,1:3) = (/ 7.39908E-03_8, 4.9501E+03_8, 60.966371E-3_8/)
  MolVal(19,5,1:3) = (/ 1.87967E-03_8, 1.3137E+03_8, 61.971231E-3_8/)
  ! H2CO (20)
  MolVal(20,1,1:3) = (/ .986237E+00_8, 2.8467E+03_8, 30.010565E-3_8/)
  MolVal(20,2,1:3) = (/ 1.10802E-02_8, 5.8376E+03_8, 31.013920E-3_8/)
  MolVal(20,3,1:3) = (/ 1.97761E-03_8, 2.9864E+03_8, 32.014811E-3_8/)
  ! HOCl (21)
  MolVal(21,1,1:3) = (/ .755790E+00_8, 1.9274E+04_8, 51.971593E-3_8/)
  MolVal(21,2,1:3) = (/ .241683E+00_8, 1.9616E+04_8, 53.968644E-3_8/)
  ! N2 (22)
  MolVal(22,1,1:3) = (/ .992687E+00_8, 4.6598E+02_8, 28.006147E-3_8/)
  ! HCN (23)
  MolVal(23,1,1:3) = (/ .985114E+00_8, 8.9529E+02_8, 27.010899E-3_8/)
  MolVal(23,2,1:3) = (/ 1.10676E-02_8, 1.8403E+03_8, 28.014254E-3_8/)
  MolVal(23,3,1:3) = (/ 3.62174E-03_8, 6.2141E+02_8, 28.007933E-3_8/)
  ! CH3Cl (24)
  MolVal(24,1,1:3) = (/ .748937E+00_8, 1.1583E+05_8, 49.992328E-3_8/)
  MolVal(24,2,1:3) = (/ .239491E+00_8, 1.1767E+05_8, 51.989379E-3_8/)
  ! H2O2 (25)
  MolVal(25,1,1:3) = (/ .994952E+00_8, 9.8198E+03_8, 34.005480E-3_8/)
  ! C2H2 (26)
  MolVal(26,1,1:3) = (/ .977599E+00_8, 4.1403E+02_8, 26.015650E-3_8/)
  MolVal(26,2,1:3) = (/ 2.19663E-02_8, 1.6562E+03_8, 27.019005E-3_8/)
  MolVal(26,3,1:3) = (/ 3.04550E-04_8, 1.5818E+03_8, 27.021825E-3_8/)
  ! C2H6 (27)
  MolVal(27,1,1:3) = (/ .976990E+00_8, 7.0881E+04_8, 30.046950E-3_8/)
  MolVal(27,2,1:3) = (/ 2.19526E-02_8, 8.0001E+04_8, 31.050305E-3_8/)
  ! PH3 (28)
  MolVal(28,1,1:3) = (/ .999533E+00_8, 3.2486E+03_8, 33.997238E-3_8/)
  ! COF2 (29)
  MolVal(29,1,1:3) = (/ .986544E+00_8, 7.0044E+04_8, 65.991722E-3_8/)
  ! SF6 (30)
  MolVal(30,1,1:3) = (/ .950180E+00_8, 1.6233E+06_8, 145.962492E-3_8/)
  ! H2S (31)
  MolVal(31,1,1:3) = (/ .949884E+00_8, 5.0307E+02_8, 33.987721E-3_8/)
  MolVal(31,2,1:3) = (/ 4.21369E-02_8, 5.0435E+02_8, 35.983515E-3_8/)
  MolVal(31,3,1:3) = (/ 7.49766E-03_8, 2.0149E+03_8, 34.987105E-3_8/)
  ! HCOOH (32)
  MolVal(32,1,1:3) = (/ .983898E+00_8, 3.9133E+04_8, 46.005480E-3_8/)
  ! HO2 (33)
  MolVal(33,1,1:3) = (/ .995107E+00_8, 4.3004E+03_8, 32.997655E-3_8/)
  ! O (34)
  MolVal(34,1,1:3) = (/ .997628E+00_8, 6.7212E+00_8, 15.994915E-3_8/)
  ! ClONO2 (35)
  MolVal(35,1,1:3) = (/ .749570E+00_8, 4.7884E+06_8, 96.956672E-3_8/)
  MolVal(35,2,1:3) = (/ .239694E+00_8, 4.9102E+06_8, 98.953723E-3_8/)
  ! NO+ (36)
  MolVal(36,1,1:3) = (/ .993974E+00_8, 3.1168E+02_8, 29.997989E-3_8/)
  ! HOBr (37)
  MolVal(37,1,1:3) = (/ .505579E+00_8, 2.8339E+04_8, 95.921076E-3_8/)
  MolVal(37,2,1:3) = (/ .491894E+00_8, 2.8238E+04_8, 97.919027E-3_8/)
  ! C2H4 (38)
  MolVal(38,1,1:3) = (/ .977294E+00_8, 1.1041E+04_8, 28.031300E-3_8/)
  MolVal(38,2,1:3) = (/ 2.19595E-02_8, 4.5197E+04_8, 29.034655E-3_8/)
  ! CH3OH (39)
  MolVal(39,1,1:3) = (/ .985930E+00_8, 3.5314E+04_8, 32.026215E-3_8/)
  ! CH3Br (40)
  MolVal(40,1,1:3) = (/ .500995E+00_8, 2.0762E+04_8, 93.941811E-3_8/)
  MolVal(40,2,1:3) = (/ .487433E+00_8, 2.0849E+04_8, 95.939764E-3_8/)
  ! CH3CN (41)
  MolVal(41,1,1:3) = (/ .973866E+00_8, 8.8670E+04_8, 41.026549E-3_8/)
  ! CF4 (42)
  MolVal(42,1,1:3) = (/ .988890E+00_8, 1.2127E+05_8, 87.993616E-3_8/)
  ! CS (46)
  MolVal(46,1,1:3) = (/ .939624E+00_8, 2.5653E+02_8, 43.971036E-3_8/)
  MolVal(46,2,1:3) = (/ .416817E-01_8, 2.6069E+02_8, 45.966787E-3_8/)
  MolVal(46,3,1:3) = (/ .105565E-01_8, 5.4345E+02_8, 44.974368E-3_8/)
  MolVal(46,4,1:3) = (/ .741668E-02_8, 1.0346E+03_8, 44.970399E-3_8/)
  end if

  end subroutine SetMolVal

end module mod_LineByLine

  Include 'sub_cntnm_progr_MO.f90' ! 2014/06/03
