!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2005, 2006. All rights reserved.
!---------------------------------------------------------------------
!= Module ChemCalc
!
!   * Developer: SUGIYAMA Ko-ichiro
!   * Version: $Id: chemcalc.f90,v 1.1.1.1 2006/04/25 03:43:58 deepconv Exp $
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview
!
!化学関連の諸量を計算するためのモジュール. 
!
!Antoine の飽和蒸気圧式を用いて以下を求める. 
!  * 飽和蒸気圧
!  * 飽和蒸気圧の温度微分
!  * 潜熱
! 
!== Error Handling
!
!== Bugs
!
!== Note
!
!== Future Plans
!
!

module ChemCalc
  !
  !化学関連の諸量を計算するためのモジュール. 
  !

  !モジュール呼び出し
  use ChemData, only: ChemData_SpcNum,             & !データベース上の化学種数
    &                 ChemData_MolWt,              & !分子量
    &                 ChemData_SvapPress_AntoineA, & !Antoine 式の A 係数
    &                 ChemData_SvapPress_AntoineB, & !Antoine 式の B 係数
    &                 ChemData_SvapPress_AntoineC, & !Antoine 式の C 係数
    &                 ChemData_SvapPress_AntoineUnit !単位変換用係数
  
  !暗黙の型宣言禁止
  implicit none

  !全体に対して private 属性を付ける. 
  private
  
  !公開するものには public 属性を付ける
  public ChemCalc_Init                            !初期化ルーチン
  public ChemCalcSpc_Init, ChemCalcDim_Init       !初期化ルーチン
  public CheckSvapPress                           !飽和蒸気圧チェックルーチン
  public MolWt                                    !分子量
  public GasR                                     !気体定数
  public Phase                                    !相
  public ChemPotRef, ChemPot, ChemPotNH4SH  
                                                  !化学ポテンシャル
  public EntropyRef, Entropy, EntropyNH4SH  
                                                  !化学ポテンシャル
  public CpPerMassRef, CpPerMolRef, CvPerMassRef  !定圧比熱, 定積比熱
  public xz_SvapPress, SvapPress
                                                  !飽和蒸気圧 [Pa]
  public SvapPressAMP
                                                  !飽和蒸気圧 [Pa]
  public xz_DSvapPressDTemp, DSvapPressDTemp 
                                                  !飽和蒸気圧の温度微分 [Pa/K]
  public xz_LatentHeatPerMass, LatentHeatPerMass
                                                  !潜熱 [J/K kg]
  public LatentHeatPerMassAMP
                                                  !潜熱 [J/K kg]
  public xz_LatentHeatPerMol, LatentHeatPerMol
                                                  !潜熱 [J/K mol]
  public LatentHeatPerMolAMP
                                                  !潜熱 [J/K mol]
  public ReactHeatNH4SHPerMass, ReactHeatNH4SHPerMol
                                                  !NH4SH 反応熱

  real(8)             :: ReactHeatNH4SHPerMass    !NH4SH 生成反応熱 [J/K kg]
  real(8)             :: ReactHeatNH4SHPerMol     !NH4SH 生成反応熱 [J/K mol]
  integer             :: XMin = 10   !X 座標範囲(下限). 値はダミー
  integer             :: XMax = 10   !X 座標範囲(上限). 値はダミー
  integer             :: ZMin = 10   !Z 座標範囲(下限). 値はダミー
  integer             :: ZMax = 10   !Z 座標範囲(上限). 値はダミー
  integer             :: SpcNum = 1  !化学種の数. 値はダミー
  real(8)             :: a_AntA(ChemData_SpcNum)   !Antoine の蒸気圧式の A 係数
  real(8)             :: a_AntB(ChemData_SpcNum)   !Antoine の蒸気圧式の B 係数
  real(8)             :: a_AntC(ChemData_SpcNum)   !Antoine の蒸気圧式の C 係数
  real(8)             :: a_AntU(ChemData_SpcNum)   !Antoine の蒸気圧式の単位換算のための係数
  real(8)             :: a_MolWt(ChemData_SpcNum)  !分子量
  integer, allocatable      :: SpcID(:)     !系に含まれる化学種. 数字で判別
  character(3), allocatable :: Phase(:)     !相
 
  save ReactHeatNH4SHPerMass, ReactHeatNH4SHPerMol
  save XMin, XMax, ZMin, ZMax, SpcNum
  save SpcID, Phase
  save a_AntA, a_AntB, a_AntC, a_AntU, a_MolWt

contains
  
!!!
!!! 初期化ルーチン
!!!
!!!==========================================================================
  subroutine ChemCalc_Init(IMin, IMax, KMin, KMax, SNum, SID)
    !
    !初期化ルーチンのラッパー. 
    !  ChemCalcSpc_init と ChemCalcDim_init を実行
    !

    !モジュール呼び出し
    use ChemData, only: GasRUniv,      &!気体定数
      &                 ChemData_SpcID  !化学種の ID 検索
      
    !暗黙の型宣言禁止
    implicit none

    !入出力変数
    integer,intent(in) :: IMin      !X 方向の配列の下限
    integer,intent(in) :: IMax      !X 方向の配列の上限
    integer,intent(in) :: KMin      !Z 方向の配列の下限
    integer,intent(in) :: KMax      !Z 方向の配列の上限
    integer,intent(in) :: SNum      !化学種の個数
    integer,intent(in) :: SID(SNum) !系に存在する化学種の数
    character(20)      :: Name(1)
    integer            :: dim
    integer            :: id(1)

    !初期化ルーチンを 2 つ呼び出すだけ
    call ChemCalcSpc_Init(SNum, SID)
    call ChemCalcDim_Init(IMin, IMax, KMin, KMax)


    !NH4SH の反応熱の初期化
    !  NH4SH 1kg に対する反応熱にする.
    dim     = 1
    Name(1) = 'NH4SH-s'
    id      = ChemData_SpcID(dim, Name) 
    
    ReactHeatNH4SHPerMol  = GasRUniv * 10834.0d0
    ReactHeatNH4SHPerMass = GasRUniv * 10834.0d0 / MolWt( id(1) )

  end subroutine ChemCalc_Init


!!!==========================================================================
  subroutine ChemCalcSpc_Init(SNum, SID)
    !
    !初期化ルーチン(その 1). 
    !  利用する化学種の確定と, それぞれの物性値を決定する.
    !  配列計算をしない場合には, この初期化をすれば十分である. 
    !

    !モジュール呼び出し
    use ChemData, only: ChemData_Phase    !相

    !暗黙の型宣言禁止
    implicit none

    !入出力変数
    integer,intent(in) :: SNum            !化学種の個数
    integer,intent(in) :: SID(SNum)       !系に存在する化学種の数
    
    !内部変数
    integer            :: s

    !-----------------------------------------------------------
    ! 初期化
    !-----------------------------------------------------------
    SpcNum  = SNum
!    write(*,*) "ChemCalc: SpcNum ", SpcNum
    
    !配列割り当て
    if (allocated(SpcID)) deallocate( SpcID, Phase)
    allocate( SpcID(SpcNum),  Phase(SpcNum) )

    !-----------------------------------------------------------    
    ! データベースから必要なデータだけ取り出す
    !-----------------------------------------------------------
    !系に存在する化学種. 入力された値を取り込む
    SpcID = SID
    
    !相
    do s = 1, SpcNum
      Phase(s) = ChemData_Phase(SpcID(s)) 
    end do
    
    !Antoine の飽和蒸気圧式の係数
    a_AntA = ChemData_SvapPress_AntoineA
    a_AntB = ChemData_SvapPress_AntoineB
    a_AntC = ChemData_SvapPress_AntoineC
    a_AntU = ChemData_SvapPress_AntoineUnit

    !分子量
    a_MolWt = ChemData_MolWt
    
  end subroutine ChemCalcSpc_Init


!!!==========================================================================
  subroutine ChemCalcDim_Init(IMin, IMax, KMin, KMax)
    !
    !初期化ルーチン(その 2)
    !  配列サイズの確定を行う. 配列関数を利用する場合には, 
    !  このサブルーチンで初期化する必要がある. 
    !  先に ChemCalcSpc_init を実行しておく必要がある. 
    !

    !暗黙の型宣言禁止
    implicit none

    !入出力変数
    integer,intent(in) :: IMin  !X 方向の配列の下限
    integer,intent(in) :: IMax  !X 方向の配列の上限
    integer,intent(in) :: KMin  !Z 方向の配列の下限
    integer,intent(in) :: KMax  !Z 方向の配列の上限
    
    !-----------------------------------------------------------
    ! 初期化
    !-----------------------------------------------------------
    XMin = IMin;    XMax = IMax
    ZMin = KMin;    ZMax = KMax
    
  end subroutine ChemCalcDim_Init


!!!
!!! 飽和蒸気圧, 潜熱, etc. の基本関数. 
!!! 化学種の ID と温度に対して値を返す
!!!  

!!!==========================================================================
  function EntropyNH4SH(ID_NH3, ID_H2S, Temp)
    !
    !NH4SH の標準化学ポテンシャルを計算する. 
    !
    
    !モジュール呼び出し
    use ChemData, only: PressRef,            &!基準圧力
      &                 GasRUniv              !気体定数
      
    !暗黙の型宣言禁止
    implicit none
    
    !入出力変数  
    real(8)             :: EntropyNH4SH       !標準化学ポテンシャル
    real(8), intent(in) :: Temp               !温度
    integer, intent(in) :: ID_NH3             !アンモニアの ID
    integer, intent(in) :: ID_H2S             !硫化水素の ID

    EntropyNH4SH = 0.0d0
    EntropyNH4SH = &
      &   EntropyRef(ID_NH3, Temp) &
      & + EntropyRef(ID_H2S, Temp) &
      & - GasRUniv * 61.781d0            &
      & + GasRUniv * dlog( PressRef * PressRef * 1.0d2) 

  end function EntropyNH4SH


!!!==========================================================================
  function EntropyRef(ID, Temp)
    !
    !引数で与えられた化学種に対して, 標準化学ポテンシャルを計算する.
    !
    
    !モジュール呼び出し
    use ChemData, only: ChemData_CpPerMolRef,&!標準状態での単位モル当たりの比熱
      &                 ChemData_Entropy,    &!標準状態でのエントロピー
      &                 TempRef               !基準温度
      
    !暗黙の型宣言禁止
    implicit none
    
    !入出力変数  
    real(8)             :: EntropyRef         !標準化学ポテンシャル
    real(8), intent(in) :: Temp               !温度
    integer, intent(in) :: ID           !化学種の ID
    real(8)             :: CpPerMol           !比熱 
    real(8)             :: Entropy            !エントロピー
    real(8)             :: sign               !積分のための符号
    
    !データベースから情報取得
    CpPerMol = ChemData_CpPerMolRef(ID)
    Entropy  = ChemData_Entropy(ID)
    
    EntropyRef = 0.0d0
    sign = 1.0d0
    if (Temp < TempRef) sign = -1.0d0

    EntropyRef =                &
      &   Temp * Entropy              &
      & + sign * CpPerMol * dlog( (Temp / TempRef) ** sign )
    
  end function EntropyRef


!!!==========================================================================
  function Entropy(ID, Temp, PPress)
    !
    !引数で与えられた化学種に対して, 化学ポテンシャルを計算する.
    !
    
    !モジュール呼び出し
    use ChemData, only: PressRef,           & !基準圧力
      &                 GasRUniv              !気体定数
    
    !暗黙の型宣言禁止
    implicit none
    
    !入出力変数  
    real(8)             :: Entropy            !標準化学ポテンシャル
    real(8), intent(in) :: Temp               !温度
    real(8), intent(in) :: PPress             !圧力
    integer, intent(in) :: ID           !化学種の ID
    
    Entropy = 0.0d0
    
    Entropy = &
      & EntropyRef(ID, Temp) - GasRUniv * dlog(PPress / PressRef)
    
  end function Entropy
  

!!!==========================================================================
  function ChemPotNH4SH(ID_NH3, ID_H2S, Temp)
    !
    !NH4SH の標準化学ポテンシャルを計算する. 
    !
    
    !モジュール呼び出し
    use ChemData, only: PressRef,            &!基準圧力
      &                 GasRUniv              !気体定数
      
    !暗黙の型宣言禁止
    implicit none
    
    !入出力変数  
    real(8)             :: ChemPotNH4SH       !標準化学ポテンシャル
    real(8), intent(in) :: Temp               !温度
    integer, intent(in) :: ID_NH3             !アンモニアの ID
    integer, intent(in) :: ID_H2S             !硫化水素の ID

    ChemPotNH4SH = 0.0d0
    ChemPotNH4SH = &
      &   ChemPotRef(ID_NH3, Temp) &
      & + ChemPotRef(ID_H2S, Temp) &
      & + GasRUniv * Temp * 61.781d0 - GasRUniv * 10834.0d0  &
      & - GasRUniv * Temp * dlog( PressRef * PressRef * 1.0d2) 

  end function ChemPotNH4SH


!!!==========================================================================
  function ChemPotRef(ID, Temp)
    !
    !引数で与えられた化学種に対して, 標準化学ポテンシャルを計算する.
    !
    
    !モジュール呼び出し
    use ChemData, only: ChemData_CpPerMolRef,&!標準状態での単位モル当たりの比熱
      &                 ChemData_Enthalpy,   &!標準状態でのエンタルピー
      &                 ChemData_Entropy,    &!標準状態でのエントロピー
      &                 TempRef               !基準温度
      
    !暗黙の型宣言禁止
    implicit none
    
    !入出力変数  
    real(8)             :: ChemPotRef         !標準化学ポテンシャル
    real(8), intent(in) :: Temp               !温度
    integer, intent(in) :: ID           !化学種の ID
    real(8)             :: CpPerMol           !比熱 
    real(8)             :: Enthalpy           !エンタルピー
    real(8)             :: Entropy            !エントロピー
    real(8)             :: sign               !積分のための符号
    
    !データベースから情報取得
    CpPerMol = ChemData_CpPerMolRef(ID)
    Enthalpy = ChemData_Enthalpy(ID)
    Entropy  = ChemData_Entropy(ID)
    
    ChemPotRef = 0.0d0
    sign = 1.0d0
    if (Temp < TempRef) sign = -1.0d0

    ChemPotRef =                &
      &   Enthalpy                    &
      & - Temp * Entropy              &
      & + CpPerMol * (Temp - TempRef) &
      & - sign * CpPerMol * Temp * dlog( (Temp / TempRef) ** sign )
    
  end function ChemPotRef


!!!==========================================================================
  function ChemPot(ID, Temp, PPress)
    !
    !引数で与えられた化学種に対して, 化学ポテンシャルを計算する.
    !
    
    !モジュール呼び出し
    use ChemData, only: PressRef,           & !基準圧力
      &                 GasRUniv              !気体定数
    
    !暗黙の型宣言禁止
    implicit none
    
    !入出力変数  
    real(8)             :: ChemPot            !標準化学ポテンシャル
    real(8), intent(in) :: Temp               !温度
    real(8), intent(in) :: PPress             !化学種の分圧
    integer, intent(in) :: ID           !化学種の ID
    
    ChemPot = 0.0d0
    
    ChemPot = &
      & ChemPotRef(ID, Temp) &
      & + GasRUniv * Temp * dlog( PPress / PressRef)
    
  end function ChemPot
  

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

    !モジュール呼び出し
    use ChemData, only: ChemData_CpPerMassRef !標準状態での単位質量当たりの比熱

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

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

  end function CpPerMassRef


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

    !モジュール呼び出し
    use ChemData, only: ChemData_CpPerMolRef !標準状態での単位質量当たりの比熱

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

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

  end function CpPerMolRef


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

    !モジュール呼び出し
    use ChemData, only: ChemData_CvPerMassRef !標準状態での単位質量当たりの比熱

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

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

  end function CvPerMassRef


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

    !モジュール呼び出し
    use ChemData, only: ChemData_MolWt   !分子量

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

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

  end function MolWt


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

    !モジュール呼び出し
    use ChemData, only: ChemData_GasR    !気体定数 [J/K kg]

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

  end function GasR


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

    !モジュール呼び出し
    use ChemData, only:                          &
      &                 ChemData_SvapPress_AMPA, & !AMP 式の A 係数
      &                 ChemData_SvapPress_AMPB, & !AMP 式の B 係数
      &                 ChemData_SvapPress_AMPC, & !AMP 式の C 係数
      &                 ChemData_SvapPress_AMPD, & !AMP 式の D 係数
      &                 ChemData_SvapPress_AMPE    !AMP 式の E 係数

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

    !内部変数
    real(8)             :: A, B, C, D, E
    real(8)             :: LogSvapPress
    
    !初期化
    SvapPressAMP = 0.0d0
    A = ChemData_SvapPress_AMPA(ID)
    B = ChemData_SvapPress_AMPB(ID)
    C = ChemData_SvapPress_AMPC(ID)
    D = ChemData_SvapPress_AMPD(ID)
    E = ChemData_SvapPress_AMPE(ID)
    
    !飽和蒸気圧の対数を計算
    LogSvapPress =  A / Temp + B + C * dlog(Temp) + D * Temp &
      &             + E * (temp ** 2) + dlog(1.0d-1) 
    
    !対数が大きくなりすぎないようにする. 
    ! Fujitsu Fortran Compiler では 700 より大きい数の exp を取ると警告が出る.
    LogSvapPress = min(7.0d2, LogSvapPress)

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


!!!==========================================================================
  function DLogSvapPressAMPDTemp(ID, Temp)
    !
    !引数で与えられた化学種と温度に対して, 飽和蒸気圧の log の温度微分を計算
    !AMP 式を利用
    !
    
    !モジュール呼び出し
    use ChemData, only:                          & !
      &                 ChemData_SvapPress_AMPA, & !Antoine 式の B 係数
      &                 ChemData_SvapPress_AMPC, & !Antoine 式の B 係数
      &                 ChemData_SvapPress_AMPD, & !Antoine 式の B 係数
      &                 ChemData_SvapPress_AMPE    !Antoine 式の C 係数
    
    !暗黙の型宣言禁止
    implicit none
  
    !入出力変数    
    real(8)             :: DLogSvapPressAMPDTemp 
                                            !飽和蒸気圧の log の温度微分 [Pa/K]
    integer, intent(in) :: ID         !化学種の ID
    real(8), intent(in) :: Temp             !温度 [K]
    
    !内部変数
    real(8)             :: A, C, D, E
    
    !初期化
    DLogSvapPressAMPDTemp = 0.0d0
    A = ChemData_SvapPress_AMPA(ID)
    C = ChemData_SvapPress_AMPC(ID)
    D = ChemData_SvapPress_AMPD(ID)
    E = ChemData_SvapPress_AMPE(ID)
    
    !飽和蒸気圧の温度微分
    DLogSvapPressAMPDTemp =                                 &
      & - A / (Temp ** 2.0d0) + C / Temp + D + (2.0d0 * E * Temp)
    
  end function DLogSvapPressAMPDTemp


!!!==========================================================================
  function DSvapPressAMPDTemp(ID, Temp)
    !
    !引数で与えられた化学種と温度に対して, 飽和蒸気圧の温度微分を計算
    !
    
    !暗黙の型宣言禁止
    implicit none
  
    !入出力変数    
    real(8)             :: DSvapPressAMPDTemp  
                                                  !飽和蒸気圧の温度微分 [Pa/K]
    integer, intent(in) :: ID               !化学種名
    real(8), intent(in) :: Temp                   !温度 [K]
    
    !初期化
    DSvapPressAMPDTemp = 0.0d0
    
    !飽和蒸気圧の温度微分
    DSvapPressAMPDTemp =                          &
      & SvapPressAMP(ID, Temp)              &
      & * DLogSvapPressAMPDTemp(ID, Temp)
    
  end function DSvapPressAMPDTemp
  

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

    !モジュール呼び出し
    use ChemData, only: ChemData_MolWt         !分子量

    !暗黙の型宣言禁止
    implicit none

    !入出力変数
    real(8)             :: LatentHeatPerMassAMP !潜熱
    integer, intent(in) :: ID            !化学種名
    real(8), intent(in) :: Temp                !温度
    
    !内部変数
    real(8), parameter  :: GasRUniv = 8.314d0  !普遍気体定数
    
    !初期化
    LatentHeatPerMassAMP  = 0.0d0
    
    !潜熱の計算
    LatentHeatPerMassAMP =                           &
      & DLogSvapPressAMPDTemp(ID, Temp) &
      & * GasRUniv * (Temp ** 2.0d0)                &
      & / ChemData_MolWt(ID)

  end function LatentHeatPerMassAMP


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

     !暗黙の型宣言禁止
    implicit none

    !入出力変数
    real(8)             :: LatentHeatPerMolAMP !潜熱
    integer, intent(in) :: ID            !化学種名
    real(8), intent(in) :: Temp                !温度
    
    !内部変数
    real(8), parameter  :: GasRUniv = 8.314d0  !普遍気体定数
    
    !初期化
    LatentHeatPerMolAMP  = 0.0d0
    
    !潜熱の計算
    LatentHeatPerMolAMP =                           &
      & DLogSvapPressAMPDTemp(ID, Temp)       &
      & * GasRUniv * (Temp ** 2.0d0)                

  end function LatentHeatPerMolAMP



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

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

    !内部変数
    real(8)             :: LogSvapPress
    real(8), parameter  :: Temp0C = 273.15d0
    
    !飽和蒸気圧の対数を計算
    LogSvapPress =                                     &
      & (                                              &
      &    a_AntA(ID)                                  &
      &    - a_AntB(ID) / (a_AntC(ID) + Temp - Temp0C) &
      &  ) * dlog(10.0d0)                              &
      & + a_AntU(ID)

    !対数が大きくなりすぎないようにする. 
    ! Fujitsu Fortran Compiler では 700 より大きい数の exp を取ると警告が出る.
    LogSvapPress = min(7.0d2, LogSvapPress)
    
    !飽和蒸気圧を計算
    SvapPress =  dexp( LogSvapPress )

  end function SvapPress


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

    !内部変数
    real(8)             :: LogSvapPress
    real(8)             :: DLogSvapPressDTemp
    real(8), parameter  :: Temp0C = 273.15d0
    

    !飽和蒸気圧の対数を計算
    LogSvapPress =                                     &
      & (                                              &
      &    a_AntA(ID)                                  &
      &    - a_AntB(ID) / (a_AntC(ID) + Temp - Temp0C) &
      &  ) * dlog(10.0d0)                              &
      & + a_AntU(ID)

    DLogSvapPressDTemp =                               &
      &  a_AntB(ID) * dlog(10.0d0)                     &
      &   / ( (a_AntC(ID) + Temp - Temp0C) ** 2.0d0 )
    
    !対数が大きくなりすぎないようにする. 
    ! Fujitsu Fortran Compiler では 700 より大きい数の exp を取ると警告が出る.
    LogSvapPress = min(7.0d2, LogSvapPress)
    
    !飽和蒸気圧の温度微分
    DSvapPressDTemp = dexp( LogSvapPress ) * DLogSvapPressDTemp
    
  end function DSvapPressDTemp


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

    !暗黙の型宣言禁止
    implicit none

    !入出力変数
    real(8)             :: LatentHeatPerMass   !潜熱
    integer, intent(in) :: ID                  !化学種名
    real(8), intent(in) :: Temp                !温度

    !内部関数
    real(8)             :: DLogSvapPressDTemp
    real(8), parameter  :: GasRUniv = 8.314d0  !普遍気体定数
    real(8), parameter  :: Temp0C = 273.15d0

    
    !潜熱の計算
    DLogSvapPressDTemp =                               &
      &  a_AntB(ID) * dlog(10.0d0)                     &
      &   / ( (a_AntC(ID) + Temp - Temp0C) ** 2.0d0 )

    LatentHeatPerMass =                                 &
      & DLogSvapPressDTemp * GasRUniv * (Temp ** 2.0d0) &
      & / MolWt(ID)

  end function LatentHeatPerMass


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

    !暗黙の型宣言禁止
    implicit none

    !入出力変数
    real(8)             :: LatentHeatPerMol    !潜熱
    integer, intent(in) :: ID                  !化学種名
    real(8), intent(in) :: Temp                !温度

    !内部関数
    real(8)             :: DLogSvapPressDTemp
    real(8), parameter  :: GasRUniv = 8.314d0  !普遍気体定数
    real(8), parameter  :: Temp0C = 273.15d0

    
    !潜熱の計算
    DLogSvapPressDTemp =                               &
      &  a_AntB(ID) * dlog(10.0d0)                     &
      &   / ( (a_AntC(ID) + Temp - Temp0C) ** 2.0d0 )

    LatentHeatPerMol =                                &
      & DLogSvapPressDTemp * GasRUniv * (Temp ** 2.0d0) 

  end function LatentHeatPerMol


!!!
!!! 飽和蒸気圧チェックルーチン
!!!

!!!==========================================================================
  subroutine CheckSvapPress()
    !
    !飽和蒸気圧のチェックサブルーチン. 
    !蒸気圧の測定値と計算値の差を出力する.
    !

    !モジュール読み込み
    use ChemData, only : ChemData_DataNum,         & !データ点数 
      &                  ChemData_SVapPress,       & !飽和蒸気圧
      &                  ChemData_SVapPress_Temp     !温度
    
    !暗黙の型宣言禁止
    implicit none
    
    !内部変数
    integer             :: SNum
    integer             :: DNum
    real(8)             :: SVapPressData(SpcNum, ChemData_DataNum)
    real(8)             :: Temp(SpcNum, ChemData_DataNum)
    real(8)             :: SVapPressCalc

    !データベースから情報を取得する
    do SNum = 1, SpcNum
      SvapPressData(SNum,:) = ChemData_SvapPress(SpcID(SNum),:) 
      Temp(SNum,:)          = ChemData_SvapPress_Temp(SpcID(SNum),:) 
    end do

    do SNum = 1, SpcNum
      if (Phase(SNum) == 'Gas') cycle   !気相の場合
      
      write(*,*) "SpcName:      ", SpcID(SNum)
      
      do DNum = 1, ChemData_DataNum
        if ( Temp(SNum,DNum) == 0.0d0) cycle
        SvapPressCalc = SvapPress( SpcID(SNum), Temp(SNum,DNum) )
        
        write(*,*) "Temp:         ", Temp(SNum,DNum)
        write(*,*) "SVapPress:    ", SVapPressData(SNum,DNum), SVapPressCalc
        write(*,*) "DelSVapPress: ", abs(SVapPressData(SNum,DNum) - SVapPressCalc)
      end do
      SVapPressCalc = 0.0d0 !初期化
    end do
    
  end subroutine CheckSvapPress
  

!!!
!!! 空間 2 次元の関数群
!!!

!!!==========================================================================
  function xz_SvapPress( ID, xz_Temp )
    !
    !飽和蒸気圧の計算. Antoine の式を利用
    !

    !暗黙の型宣言禁止
    implicit none
    
    ! 入出力変数  
    real(8)             :: xz_SvapPress(XMin:XMax, ZMin:ZMax) 
                                                          !飽和蒸気圧
    real(8), intent(in) :: xz_Temp(XMin:XMax, ZMin:ZMax)  !温度
    integer, intent(in) :: ID                             !化学種の ID
  
    !内部変数
    real(8)             :: xz_LogSvapPress(XMin:XMax, ZMin:ZMax) 
    real(8), parameter  :: Temp0C = 273.15d0


    !飽和蒸気圧の log を計算
    xz_LogSvapPress =                                     &
      & (                                                 &
      &    a_AntA(ID)                                     &
      &    - a_AntB(ID) / (a_AntC(ID) + xz_Temp - Temp0C) &
      &  ) * dlog(10.0d0)                                 &
      & + a_AntU(ID)
    
    !対数が大きくなりすぎないようにする. 
    ! Fujitsu Fortran Compiler では 700 より大きい数の exp を取ると警告が出る.
    where ( xz_LogSvapPress > 7.0d2 )
      xz_LogSvapPress = 7.0d2
    end where
    
    !飽和蒸気圧を計算
    xz_SvapPress =  dexp( xz_LogSvapPress )

  end function xz_SvapPress


!!!==========================================================================
  function xz_DSvapPressDTemp(ID, xz_Temp)
    !
    !飽和蒸気圧の温度微分の計算. Antoine の式を利用
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    !入出力変数
    real(8)             :: xz_DSvapPressDTemp(XMin:XMax, ZMin:ZMax)
                                                  !飽和蒸気圧の温度微分 [Pa/K]
    real(8), intent(in) :: xz_Temp(XMin:XMax, ZMin:ZMax)           
                                                  !温度 [K]
    integer, intent(in) :: ID                     !化学種の ID
  
    !内部変数
    real(8)             :: xz_LogSvapPress(XMin:XMax, ZMin:ZMax)
    real(8)             :: xz_DLogSvapPressDTemp(XMin:XMax, ZMin:ZMax)
    real(8), parameter  :: Temp0C = 273.15d0


    !飽和蒸気圧の log を計算
    xz_LogSvapPress =                                      &
      & (                                                  &
      &     a_AntA(ID)                                     &
      &     - a_AntB(ID) / (a_AntC(ID) + xz_Temp - Temp0C) &
      &  ) * dlog(10.0d0)                                  &
      & + a_AntU(ID)
    
    xz_DLogSvapPressDTemp =                                &
      & a_AntB(ID) * dlog(10.0d0)                          &
      &  / ( (a_AntC(ID) + xz_Temp - Temp0C) ** 2.0d0 )
    
    !対数が大きくなりすぎないようにする. 
    ! Fujitsu Fortran Compiler では 700 より大きい数の exp を取ると警告が出る.
    where ( xz_LogSvapPress > 7.0d2 )
      xz_LogSvapPress = 7.0d2
    end where
    
    !飽和蒸気圧の温度微分を計算
    xz_DSvapPressDTemp = dexp( xz_LogSvapPress ) * xz_DLogSvapPressDTemp
    
  end function xz_DSvapPressDTemp
  

!!!==========================================================================
  function xz_LatentHeatPerMol(ID, xz_Temp)
    !
    !飽和蒸気圧より潜熱を計算する. Antoine の式を利用
    !
    
    !暗黙の型宣言禁止
    implicit none

    !入出力変数
    real(8)             :: xz_LatentHeatPerMol(XMin:XMax, ZMin:ZMax)
                                                            !潜熱[J/K kg]
    real(8), intent(in) :: xz_Temp(XMin:XMax, ZMin:ZMax)    !温度[K]
    integer, intent(in) :: ID                               !化学種の ID

    !内部関数
    real(8)             :: xz_DLogSvapPressDTemp(XMin:XMax, ZMin:ZMax)
    real(8), parameter  :: GasRUniv = 8.314d0
    real(8), parameter  :: Temp0C = 273.15d0


    !飽和蒸気圧の log を計算した後, 潜熱を計算する. 
    xz_DLogSvapPressDTemp =                                    &
      & a_AntB(ID) * dlog(10.0d0)                              &
      &  / ( (a_AntC(ID) + xz_Temp - Temp0C) ** 2.0d0 )
    
    xz_LatentHeatPerMol =                                      &
      & xz_DLogSvapPressDTemp * GasRUniv * (xz_Temp ** 2.0d0) 
    
  end function xz_LatentHeatPerMol


!!!==========================================================================
  function xz_LatentHeatPerMass(ID, xz_Temp)
    !
    !飽和蒸気圧より潜熱を計算する. Antoine の式を利用
    !
    
    !暗黙の型宣言禁止
    implicit none

    !入出力変数
    real(8)             :: xz_LatentHeatPerMass(XMin:XMax, ZMin:ZMax)
                                                            !潜熱[J/K kg]
    real(8), intent(in) :: xz_Temp(XMin:XMax, ZMin:ZMax)    !温度[K]
    integer, intent(in) :: ID                       !化学種の ID

    !内部関数
    real(8)             :: xz_DLogSvapPressDTemp(XMin:XMax, ZMin:ZMax)
    real(8), parameter  :: GasRUniv = 8.314d0
    real(8), parameter  :: Temp0C = 273.15d0

    !初期化
    xz_DLogSvapPressDTemp = 0.0d0
    xz_LatentHeatPerMass  = 0.0d0

    !飽和蒸気圧の log を計算した後, 潜熱を計算する. 
    xz_DLogSvapPressDTemp =                                   &
      & a_AntB(ID) * dlog(10.0d0)                             &
      &  / ( (a_AntC(ID) + xz_Temp - Temp0C) ** 2.0d0 )
    
    xz_LatentHeatPerMass =                                    &
      & xz_DLogSvapPressDTemp * GasRUniv * (xz_Temp ** 2.0d0) &
      &  / a_MolWt(ID)
    
  end function xz_LatentHeatPerMass

end module ChemCalc

