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

module ChemCalc
  !
  !化学関連の諸量を計算するためのモジュール. 
  !
  
  !暗黙の型宣言禁止
  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 xz_SvapPressAMP, SvapPressAMP
                                                  !飽和蒸気圧 [Pa]
  public xz_DSvapPressDTemp, DSvapPressDTemp 
                                                  !飽和蒸気圧の温度微分 [Pa/K]
  public xz_LatentHeatPerMass, LatentHeatPerMass
                                                  !潜熱 [J/K kg]
  public xz_LatentHeatPerMassAMP, LatentHeatPerMassAMP
                                                  !潜熱 [J/K kg]
  public xz_LatentHeatPerMol, LatentHeatPerMol
                                                  !潜熱 [J/K mol]
  public xz_LatentHeatPerMolAMP, 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  !化学種の数. 値はダミー

  integer, allocatable      :: SpcID(:)     !系に含まれる化学種. 数字で判別
  character(3), allocatable :: Phase(:)     !相
  
  save ReactHeatNH4SHPerMass, ReactHeatNH4SHPerMol
  save XMin, XMax, ZMin, ZMax, SpcNum
  save SpcID, Phase
  
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
    
  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(OneSpcID, Temp)
    !
    !引数で与えられた化学種に対して, 標準化学ポテンシャルを計算する.
    !
    
    !モジュール呼び出し
    use ChemData, only: ChemData_CpPerMolRef,&!標準状態での単位モル当たりの比熱
      &                 ChemData_Entropy,    &!標準状態でのエントロピー
      &                 TempRef               !基準温度
      
    !暗黙の型宣言禁止
    implicit none
    
    !入出力変数  
    real(8)             :: EntropyRef         !標準化学ポテンシャル
    real(8), intent(in) :: Temp               !温度
    integer, intent(in) :: OneSpcID           !化学種の ID
    real(8)             :: CpPerMol           !比熱 
    real(8)             :: Entropy            !エントロピー
    real(8)             :: sign               !積分のための符号
    
    !データベースから情報取得
    CpPerMol = ChemData_CpPerMolRef(OneSpcID)
    Entropy  = ChemData_Entropy(OneSpcID)
    
    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(OneSpcID, 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) :: OneSpcID           !化学種の ID
    
    Entropy = 0.0d0
    
    Entropy = &
      & EntropyRef(OneSpcID, 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(OneSpcID, Temp)
    !
    !引数で与えられた化学種に対して, 標準化学ポテンシャルを計算する.
    !
    
    !モジュール呼び出し
    use ChemData, only: ChemData_CpPerMolRef,&!標準状態での単位モル当たりの比熱
      &                 ChemData_Enthalpy,   &!標準状態でのエンタルピー
      &                 ChemData_Entropy,    &!標準状態でのエントロピー
      &                 TempRef               !基準温度
      
    !暗黙の型宣言禁止
    implicit none
    
    !入出力変数  
    real(8)             :: ChemPotRef         !標準化学ポテンシャル
    real(8), intent(in) :: Temp               !温度
    integer, intent(in) :: OneSpcID           !化学種の ID
    real(8)             :: CpPerMol           !比熱 
    real(8)             :: Enthalpy           !エンタルピー
    real(8)             :: Entropy            !エントロピー
    real(8)             :: sign               !積分のための符号
    
    !データベースから情報取得
    CpPerMol = ChemData_CpPerMolRef(OneSpcID)
    Enthalpy = ChemData_Enthalpy(OneSpcID)
    Entropy  = ChemData_Entropy(OneSpcID)
    
    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(OneSpcID, 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) :: OneSpcID           !化学種の ID
    
    ChemPot = 0.0d0
    
    ChemPot = &
      & ChemPotRef(OneSpcID, Temp) &
      & + GasRUniv * Temp * dlog( PPress / PressRef)
    
  end function ChemPot
  

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

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

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

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

  end function CpPerMassRef


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

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

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

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

  end function CpPerMolRef


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

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

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

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

  end function CvPerMassRef


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

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

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

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

  end function MolWt


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

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

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

  end function GasR


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

    !モジュール呼び出し
    use ChemData, only: Temp0C,                      & !O度での温度[K] 
      &                 ChemData_SvapPress_AntoineA, & !Antoine 式の A 係数
      &                 ChemData_SvapPress_AntoineB, & !Antoine 式の B 係数
      &                 ChemData_SvapPress_AntoineC, & !Antoine 式の C 係数
      &                 ChemData_SvapPress_AntoineUnit !単位変換用係数

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

    !内部変数
    real(8)             :: A, B, C
    real(8)             :: Unit
    real(8)             :: LogSvapPress
    
    !初期化
    SvapPress = 0.0d0
    A = ChemData_SvapPress_AntoineA(OneSpcID)
    B = ChemData_SvapPress_AntoineB(OneSpcID)
    C = ChemData_SvapPress_AntoineC(OneSpcID)
    Unit = ChemData_SvapPress_AntoineUnit(OneSpcID)

    !飽和蒸気圧の対数を計算
    LogSvapPress = (A - B / (C + Temp - Temp0C)) * dlog(10.0d0) + Unit 

!    write(*,*) 'chemcalc: ', Temp, LogSvapPress 

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

!    write(*,*) 'chemcalc: ', SvapPress, LogSvapPress 

  end function SvapPress


!!!==========================================================================
  function SvapPressAMP(OneSpcID, 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) :: OneSpcID              !化学種の ID
    real(8), intent(in) :: Temp                  !温度 [K]

    !内部変数
    real(8)             :: A, B, C, D, E
    real(8)             :: LogSvapPress
    
    !初期化
    SvapPressAMP = 0.0d0
    A = ChemData_SvapPress_AMPA(OneSpcID)
    B = ChemData_SvapPress_AMPB(OneSpcID)
    C = ChemData_SvapPress_AMPC(OneSpcID)
    D = ChemData_SvapPress_AMPD(OneSpcID)
    E = ChemData_SvapPress_AMPE(OneSpcID)
    
    !飽和蒸気圧の対数を計算
    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 DLogSvapPressDTemp(OneSpcID, Temp)
    !
    !引数で与えられた化学種と温度に対して, 飽和蒸気圧の log の温度微分を計算
    !
    
    !モジュール呼び出し
    use ChemData, only: Temp0C,                      & !O度での温度[K] 
      &                 ChemData_SvapPress_AntoineB, & !Antoine 式の B 係数
      &                 ChemData_SvapPress_AntoineC    !Antoine 式の C 係数

    !暗黙の型宣言禁止
    implicit none
  
    !入出力変数    
    real(8)             :: DLogSvapPressDTemp!飽和蒸気圧の log の温度微分[Pa/K]
    integer, intent(in) :: OneSpcID          !化学種の ID
    real(8), intent(in) :: Temp              !温度 [K]
    
    !内部変数
    real(8)             :: B, C
    
    !初期化
    DLogSvapPressDTemp = 0.0d0
    B = ChemData_SvapPress_AntoineB(OneSpcID)
    C = ChemData_SvapPress_AntoineC(OneSpcID)
    
    !飽和蒸気圧の温度微分
    DLogSvapPressDTemp =                                &
      &  B * dlog(10.0d0) / ( (C + Temp - Temp0C) ** 2.0d0 )
    
  end function DLogSvapPressDTemp


!!!==========================================================================
  function DLogSvapPressAMPDTemp(OneSpcID, 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) :: OneSpcID         !化学種の ID
    real(8), intent(in) :: Temp             !温度 [K]
    
    !内部変数
    real(8)             :: A, C, D, E
    
    !初期化
    DLogSvapPressAMPDTemp = 0.0d0
    A = ChemData_SvapPress_AMPA(OneSpcID)
    C = ChemData_SvapPress_AMPC(OneSpcID)
    D = ChemData_SvapPress_AMPD(OneSpcID)
    E = ChemData_SvapPress_AMPE(OneSpcID)
    
    !飽和蒸気圧の温度微分
    DLogSvapPressAMPDTemp =                                 &
      & - A / (Temp ** 2.0d0) + C / Temp + D + (2.0d0 * E * Temp)
    
  end function DLogSvapPressAMPDTemp


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


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

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

    !暗黙の型宣言禁止
    implicit none

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

    !内部変数
    real(8), parameter  :: GasRUniv = 8.314d0  !普遍気体定数
    
    !初期化
    LatentHeatPerMass  = 0.0d0
    
    !潜熱の計算
    LatentHeatPerMass =                                                        &
      & DLogSvapPressDTemp(OneSpcID, Temp) * GasRUniv * (Temp ** 2.0d0) &
      & / MolWt(OneSpcID)

  end function LatentHeatPerMass


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

    !暗黙の型宣言禁止
    implicit none

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

    !内部変数
    real(8), parameter  :: GasRUniv = 8.314d0  !普遍気体定数
    
    !初期化
    LatentHeatPerMol  = 0.0d0
    
    !潜熱の計算
    LatentHeatPerMol =                                                  &
      & DLogSvapPressDTemp(OneSpcID, Temp) * GasRUniv * (Temp ** 2.0d0) 

  end function LatentHeatPerMol


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

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

    !暗黙の型宣言禁止
    implicit none

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

  end function LatentHeatPerMassAMP


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

     !暗黙の型宣言禁止
    implicit none

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

  end function LatentHeatPerMolAMP


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

!!!==========================================================================
  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
  


!!!
!!! 基本関数のラッパー
!!!

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

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

    !初期化
    xz_SvapPress = 0.0d0 

    !基本関数を呼び出すことで値を決める
    do k = ZMin, ZMax
      do i = XMin, XMax        
        xz_SvapPress(i,k) = SvapPress(OneSpcID, xz_Temp(i,k))

!        write(*,*) 'chemcalc2: ', xz_Temp(i,k), xz_SvapPress(i,k)
      end do
    end do
    
  end function xz_SvapPress


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

    !暗黙の型宣言禁止
    implicit none
    
    ! 入出力変数  
    real(8)             :: xz_SvapPressAMP(XMin:XMax, ZMin:ZMax) 
                                                             !飽和蒸気圧
    integer, intent(in) :: OneSpcID                          !化学種の ID
    real(8), intent(in) :: xz_Temp(XMin:XMax, ZMin:ZMax)     !温度
  
    !内部変数
    integer             :: i, k

    !初期化
    xz_SvapPressAMP = 0.0d0 

    !基本関数を呼び出すことで値を決める
    do k = ZMin, ZMax
      do i = XMin, XMax        
        xz_SvapPressAMP(i,k) = SvapPressAMP(OneSpcID, xz_Temp(i,k))
      end do
    end do
    
  end function xz_SvapPressAMP


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

    !初期化
    xz_DSvapPressDTemp = 0.0d0

    !基本関数を呼び出すことで値を決める
    do k = ZMin, ZMax
      do i = XMin, XMax        
        xz_DSvapPressDTemp(i,k) = DSvapPressDTemp(OneSpcID, xz_Temp(i,k))
      end do
    end do
    
  end function xz_DSvapPressDTemp


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

    !初期化
    xz_DLogSvapPressDTemp = 0.0d0

    !基本関数を呼び出すことで値を決める
    do k = ZMin, ZMax
      do i = XMin, XMax        
        xz_DLogSvapPressDTemp(i,k) = &
          & DLogSvapPressDTemp(OneSpcID, xz_Temp(i,k))
      end do
    end do
    
  end function xz_DLogSvapPressDTemp


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

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

    !内部関数
    integer             :: i, k

    !初期化
    xz_LatentHeatPerMass = 0.0d0

    !基本関数を呼び出すことで値を決める
    do k = ZMin, ZMax
      do i = XMin, XMax        
        xz_LatentHeatPerMass(i,k) = LatentHeatPerMass(OneSpcID, xz_Temp(i,k))
      end do
    end do

  end function xz_LatentHeatPerMass


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

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

    !内部関数
    integer             :: i, k

    !初期化
    xz_LatentHeatPerMol = 0.0d0

    !基本関数を呼び出すことで値を決める
    do k = ZMin, ZMax
      do i = XMin, XMax        
        xz_LatentHeatPerMol(i,k) = LatentHeatPerMol(OneSpcID, xz_Temp(i,k))
      end do
    end do

  end function xz_LatentHeatPerMol


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

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

    !内部関数
    integer             :: i, k

    !初期化
    xz_LatentHeatPerMassAMP = 0.0d0

    !基本関数を呼び出すことで値を決める
    do k = ZMin, ZMax
      do i = XMin, XMax        
        xz_LatentHeatPerMassAMP(i,k) = LatentHeatPerMassAMP(OneSpcID, xz_Temp(i,k))
      end do
    end do

  end function xz_LatentHeatPerMassAMP


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

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

    !内部関数
    integer             :: i, k

    !初期化
    xz_LatentHeatPerMolAMP = 0.0d0

    !基本関数を呼び出すことで値を決める
    do k = ZMin, ZMax
      do i = XMin, XMax        
        xz_LatentHeatPerMolAMP(i,k) = LatentHeatPerMolAMP(OneSpcID, xz_Temp(i,k))
      end do
    end do

  end function xz_LatentHeatPerMolAMP



end module ChemCalc

