subroutine ChemCalc_Init( )
    !
    !初期化ルーチン
    !

    !暗黙の型宣言禁止
    implicit none

    !入出力変数
    character(20)      :: Name
    integer            :: id

    !-----------------------------------------------------------
    ! 初期化
    !

    ! データベースの初期化
    call chemdata_init( )

    !Antoine の飽和蒸気圧式の係数
    a_antA = ChemData_SvapPress_AntoineA
    a_antB = ChemData_SvapPress_AntoineB
    a_antC = ChemData_SvapPress_AntoineC
    a_antU = ChemData_SvapPress_AntoineUnit

    !AMP 式の飽和蒸気圧式の係数
    a_ampA = ChemData_SvapPress_AMPA
    a_ampB = ChemData_SvapPress_AMPB
    a_ampC = ChemData_SvapPress_AMPC
    a_ampD = ChemData_SvapPress_AMPD
    a_ampE = ChemData_SvapPress_AMPE

    !分子量
    a_MolWt = ChemData_MolWt
    
    !NH4SH の反応熱の初期化
    !  NH4SH 1kg に対する反応熱にする.
    Name = 'NH4SH-s'
    id   = ChemData_OneSpcID( Name )  
    
    ReactHeatNH4SHPerMol  = GasRUniv * 10834.0d0
    ReactHeatNH4SH = GasRUniv * 10834.0d0 / MolWt( id )

    call MessageNotify( "M", "ChemCalc_Init", "ReactHeatNH4SH = %f", d=(/ReactHeatNH4SH/) )
    call MessageNotify( "M", "ChemCalc_Init", "NH4SH MolWt = %f", d=(/MolWt(id)/) )

  end subroutine ChemCalc_Init

  function CpPerMolRef(ID)
    !
    !引数で与えられた化学種に対して, 標準状態での単位モル当たりの定圧比熱を計算
    !

    !暗黙の型宣言禁止
    implicit none
    
    !入出力変数  
    real(DP)            :: CpPerMolRef  !標準状態での単位モル当たりの比熱
    integer, intent(in) :: ID           !化学種の ID

    
    !データベースから情報取得
    CpPerMolRef = ChemData_CpPerMolRef(ID)

  end function CpPerMolRef

  function CpRef(ID)
    !
    !引数で与えられた化学種に対して, 標準状態での単位質量当たりの定圧比熱を計算
    !

    !暗黙の型宣言禁止
    implicit none
    
    !入出力変数  
    real(DP)            :: CpRef        !標準状態での単位質量当たりの比熱
    integer, intent(in) :: ID           !化学種の ID

    
    !データベースから情報取得
    CpRef = ChemData_CpRef(ID)

  end function CpRef

  function CvRef(ID)
    !
    !引数で与えられた化学種に対して, 標準状態での単位質量当たりの定圧比熱を計算
    !

    !暗黙の型宣言禁止
    implicit none
    
    !入出力変数  
    real(DP)            :: CvRef       !標準状態での単位質量当たりの比熱
    integer, intent(in) :: ID           !化学種の ID

    
    !データベースから情報取得
    CvRef = ChemData_CvRef(ID)

  end function CvRef

  function DelMolFrNH4SH(TempAll, PressAll, MolFrNH3, MolFrH2S, Humidity)
    !
    ! NH4SH 生成反応に伴う H2S と NH3 のモル比の減少分を求める
    !
    
    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(DP),intent(in) :: TempAll       !温度
    real(DP),intent(in) :: PressAll      !圧力
    real(DP),intent(in) :: MolFrNH3      !NH3 のモル比
    real(DP),intent(in) :: MolFrH2S      !H2S のモル比
    real(DP),intent(in) :: Humidity      !飽和比
    real(DP)            :: DelMolFrNH4SH !NH4SH 生成に伴うモル比変化
    real(DP)            :: EquivConst    !圧平衡定数
    real(DP)            :: PPress(2)     !作業配列(分圧)
    real(DP)            :: Solution      !作業配列(方程式の解)

    !------------------------------------------------------------
    !NH4SH の平衡条件
    !------------------------------------------------------------
    !アンモニアと硫化水素の分圧
    PPress(1) = MolFrNH3 * PressAll
    PPress(2) = MolFrH2S * PressAll

    !圧平衡定数
    EquivConst = 61.781d0 - 10834.0d0 / TempAll - dlog(1.0d2) - 2.0d0 * dlog( Humidity )
    
    !気圧変化を二次方程式の解として求める. 
    Solution = 5.0d-1 * (sum(PPress) - dsqrt( (PPress(1) - PPress(2))**2.0d0 + 4.0d0 * dexp( min( 700.0d0, EquivConst ))) )
    
    !NH4SH の生成量. 
    DelMolFrNH4SH = Solution / PressAll

  end function DelMolFrNH4SH

  function GasR(ID)
    !
    !引数で与えられた化学種に対して, 気体定数を計算
    !

    !暗黙の型宣言禁止
    implicit none
    
    !入出力変数  
    real(DP)            :: GasR          !分子量
    integer, intent(in) :: ID      !化学種の ID
    
    
    !データベースから情報取得
    GasR = ChemData_GasR(ID)

  end function GasR

  function LatentHeatPerMol(ID, Temp)
    !
    !引数で与えられた化学種と温度に対して, 潜熱を計算
    !

    !暗黙の型宣言禁止
    implicit none

    !入出力変数
    real(DP)            :: LatentHeatPerMol   !潜熱
    real(DP),intent(in) :: Temp            !温度
    integer, intent(in) :: ID                 !化学種名
    
    !内部変数
    real(DP)            :: DLogSvapPressDTemp
    real(DP),parameter  :: GasRUniv = 8.314d0  !普遍気体定数
    

    !飽和蒸気圧の温度微分
    DLogSvapPressDTemp = - a_ampA(ID) / (Temp ** 2.0d0) + a_ampC(ID) / Temp + a_ampD(ID) + a_ampE(ID) * 2.0d0 * Temp
        
    !潜熱の計算
    LatentHeatPerMol = DLogSvapPressDTemp * GasRUniv * (Temp ** 2.0d0) 

  end function LatentHeatPerMol

  function MolWt(ID)
    !
    !引数で与えられた化学種に対して, 分子量を計算
    !

    !暗黙の型宣言禁止
    implicit none
    
    !入出力変数  
    real(DP)            :: MolWt         !分子量
    integer, intent(in) :: ID      !化学種の ID

    
    !データベースから情報取得
    MolWt = ChemData_MolWt(ID)

  end function MolWt

  function SvapPress(ID, Temp)
    !
    !引数で与えられた化学種と温度に対して, 飽和蒸気圧を計算. AMP 式を利用
    !

    !暗黙の型宣言禁止
    implicit none
    
    !入出力変数  
    real(DP)            :: SvapPress   !飽和蒸気圧
    real(DP),intent(in) :: Temp        !温度 [K]
    integer, intent(in) :: ID          !化学種の ID

    !内部変数
    real(DP)            :: LogSvapPress
    
    !飽和蒸気圧の対数を計算
    !対数が大きくなりすぎないようにする. 
    ! Fujitsu Fortran Compiler では 700 より大きい数の exp を取ると警告が出る.
    LogSvapPress = min( ( a_ampA(ID) / Temp + a_ampB(ID) + a_ampC(ID) * dlog( Temp ) + a_ampD(ID) * Temp + a_ampE(ID) * ( temp ** 2 ) + dlog(1.0d-1) ), 7.0d2 )

    !飽和蒸気圧を計算
    SvapPress =  dexp( LogSvapPress )
    
  end function SvapPress

  function xyz_DQMixSatDPTemp(SpcID, MolWt, xyz_TempAll, xyz_ExnerAll)
    !
    !飽和蒸気圧の θ 微分を行う
    !実際には, dq/dp * dp/dT * dT/dθ を実行. (但し p は飽和蒸気圧)
    !

    !暗黙の型宣言禁止
    implicit none 
    
    !入出力変数
    integer, intent(in) :: SpcID
    real(DP),intent(in) :: MolWt
    real(DP),intent(in) :: xyz_TempAll(imin:imax,jmin:jmax,kmin:kmax)
                                            !温度(擾乱 + 基本場)
    real(DP),intent(in) :: xyz_ExnerAll(imin:imax,jmin:jmax,kmin:kmax)
                                            !エクスナー関数(擾乱 + 基本場)
    real(DP)            :: xyz_PressAll(imin:imax,jmin:jmax,kmin:kmax)
                                            !圧力(擾乱 + 基本場)
    real(DP)            :: xyz_DQMixSatDPTemp(imin:imax,jmin:jmax,kmin:kmax)
                           
!    xyz_DQMixSatDPTemp   = 0.0d0
    xyz_PressAll = PressBasis * (xyz_ExnerAll ** (CpDry / GasRDry))
    
    xyz_DQMixSatDPTemp = MolWt / ( MolWtDry * xyz_PressAll ) * xyz_DSvapPressDTemp(SpcID, xyz_TempAll) * xyz_ExnerAll
    
  end function xyz_DQMixSatDPTemp

  function xyz_DSvapPressDTemp(ID, xyz_Temp)
    !
    !引数で与えられた化学種と温度に対して, 飽和蒸気圧の温度微分を計算
    !
    
    !暗黙の型宣言禁止
    implicit none
  
    !入出力変数    
    real(DP)            :: xyz_DSvapPressDTemp(imin:imax,jmin:jmax,kmin:kmax)
                                                         !飽和蒸気圧の温度微分 [Pa/K]
    real(DP),intent(in) :: xyz_Temp(imin:imax,jmin:jmax,kmin:kmax) 
                                                         !温度 [K]
    integer, intent(in) :: ID                            !化学種名

    !内部変数
    real(DP)            :: xyz_LogSvapPress(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)            :: xyz_DLogSvapPressDTemp(imin:imax,jmin:jmax,kmin:kmax)
    

    !飽和蒸気圧の対数を計算
    !対数が大きくなりすぎないようにする. 
    ! Fujitsu Fortran Compiler では 700 より大きい数の exp を取ると警告が出る.
    xyz_LogSvapPress = min( ( a_ampA(ID) / xyz_Temp + a_ampB(ID) + a_ampC(ID) * dlog( xyz_Temp ) + a_ampD(ID) * xyz_Temp + a_ampE(ID) * ( xyz_temp ** 2 ) + dlog(1.0d-1) ), 7.0d2 )

    !飽和蒸気圧の温度微分
    xyz_DLogSvapPressDTemp = - a_ampA(ID) / (xyz_Temp ** 2.0d0) + a_ampC(ID) / xyz_Temp + a_ampD(ID) + a_ampE(ID) * 2.0d0 * xyz_Temp
    
    !飽和蒸気圧の温度微分
    xyz_DSvapPressDTemp = xyz_DLogSvapPressDTemp * dexp( xyz_LogSvapPress ) 
    
  end function xyz_DSvapPressDTemp

  function xyz_DelQMixNH4SH(xyz_TempAll, xyz_PressAll, xyz_PressDry, xyz_QMixNH3, xyz_QMixH2S, MolWtNH3, MolWtH2S)
    !
    ! NH4SH 生成反応に伴う, NH4SH の生成量(混合比)を求める
    !
    
    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(DP),intent(in) :: xyz_TempAll(imin:imax,jmin:jmax,kmin:kmax)
                                         !温度
    real(DP),intent(in) :: xyz_PressAll(imin:imax,jmin:jmax,kmin:kmax)
                                         !圧力
    real(DP),intent(in) :: xyz_PressDry(imin:imax,jmin:jmax,kmin:kmax)
                                         !圧力
    real(DP),intent(in) :: xyz_QMixNH3(imin:imax,jmin:jmax,kmin:kmax)
                                         !NH3 の混合比
    real(DP),intent(in) :: xyz_QMixH2S(imin:imax,jmin:jmax,kmin:kmax)
                                         !H2S の混合比
    real(DP),intent(in) :: MolWtNH3      !NH3 の分子量
    real(DP),intent(in) :: MolWtH2S      !H2S の分子量

    real(DP) :: xyz_DelQMixNH4SH(imin:imax,jmin:jmax,kmin:kmax)
                                         !NH4SH の混合比
    real(DP) :: xyz_EquivConst(imin:imax,jmin:jmax,kmin:kmax)
                                         !圧平衡定数
    real(DP) :: xyzf_PartialPress(imin:imax,jmin:jmax,kmin:kmax,2)
                                         !作業配列(分圧)
    real(DP) :: xyz_Solution(imin:imax,jmin:jmax,kmin:kmax)
                                         !作業配列(方程式の解)

    !初期化
!    xyz_DelQMixNH4SH = 0.0d0
    
    !アンモニアと硫化水素の分圧. 
    xyzf_PartialPress(:,:,:,1) = xyz_QMixNH3 * xyz_PressAll * MolWtDry / MolWtNH3 
    xyzf_PartialPress(:,:,:,2) = xyz_QMixH2S * xyz_PressAll * MolWtDry / MolWtH2S 

    !圧平衡定数
    xyz_EquivConst = 61.781d0 - 10834.0d0 / xyz_TempAll - dlog(1.0d2)

    !気圧変化を求める. 
    !  (P_NH3 - X) * (P_H2S - X) = exp(Kp)
    !  DelX^2 - (P_NH3 + P_H2S) * DelX + P_NH3 * P_H2S * exp( Kp ) = 0
    !  という二次方程式を求める必要があるが, P_NH3 > P_H2S と X < P_H2S を
    !  考慮すると, 解の公式のうちが負の方が選択される.
    xyz_Solution  = ( sum(xyzf_PartialPress, 4) - dsqrt( (xyzf_PartialPress(:,:,:,1) - xyzf_PartialPress(:,:,:,2)) ** 2.0d0 + 4.0d0 * dexp( min( 700.0d0, xyz_EquivConst ) ) ) ) * 5.0d-1

    !生成量を求める
    xyz_DelQMixNH4SH = xyz_Solution * ( MolWtNH3 + MolWtH2S ) / ( xyz_PressDry * MolWtDry )

  end function xyz_DelQMixNH4SH

  function xyz_LatentHeat(ID, xyz_Temp)
    !
    !引数で与えられた化学種と温度に対して, 潜熱を計算
    !

    !暗黙の型宣言禁止
    implicit none

    !入出力変数
    real(DP)            :: xyz_LatentHeat(imin:imax,jmin:jmax,kmin:kmax)
                                                         !潜熱
    real(DP),intent(in) :: xyz_Temp(imin:imax,jmin:jmax,kmin:kmax) 
                                                         !温度
    integer, intent(in) :: ID                            !化学種名
    
    !内部変数
    real(DP)            :: xyz_DLogSvapPressDTemp(imin:imax,jmin:jmax,kmin:kmax)
    real(DP),parameter  :: GasRUniv = 8.314d0  !普遍気体定数
    

    !飽和蒸気圧の温度微分
    xyz_DLogSvapPressDTemp = - a_ampA(ID) / (xyz_Temp ** 2.0d0) + a_ampC(ID) / xyz_Temp + a_ampD(ID) + a_ampE(ID) * 2.0d0 * xyz_Temp
        
    !潜熱の計算
    xyz_LatentHeat = xyz_DLogSvapPressDTemp * GasRUniv * (xyz_Temp ** 2.0d0) / a_MolWt(ID)

  end function xyz_LatentHeat

  function xyz_SvapPress(ID, xyz_Temp)
    !
    !引数で与えられた化学種と温度に対して, 飽和蒸気圧を計算. AMP 式を利用
    !

    !暗黙の型宣言禁止
    implicit none
    
    !入出力変数  
    real(DP)            :: xyz_SvapPress(imin:imax,jmin:jmax,kmin:kmax)
                                                         !飽和蒸気圧
    real(DP),intent(in) :: xyz_Temp(imin:imax,jmin:jmax,kmin:kmax) 
                                                         !温度 [K]
    integer, intent(in) :: ID                            !化学種の ID

    !内部変数
    real(DP)            :: xyz_LogSvapPress(imin:imax,jmin:jmax,kmin:kmax)
    
    !飽和蒸気圧の対数を計算
    !対数が大きくなりすぎないようにする. 
    ! Fujitsu Fortran Compiler では 700 より大きい数の exp を取ると警告が出る.
    xyz_LogSvapPress = min( ( a_ampA(ID) / xyz_Temp + a_ampB(ID) + a_ampC(ID) * dlog( xyz_Temp ) + a_ampD(ID) * xyz_Temp + a_ampE(ID) * ( xyz_temp ** 2 ) + dlog(1.0d-1) ), 7.0d2 )

    !飽和蒸気圧を計算
    xyz_SvapPress =  dexp( xyz_LogSvapPress )
    
  end function xyz_SvapPress