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

module ChemEquiv
  !暗黙の型宣言禁止
  implicit none

  !全体に対して private 属性を付ける. 
  private
  
  !公開するものには public 属性を付ける
  public ChemEquiv_Init                           !初期化ルーチン
  public ss_SvapPress_ss,  SvapPress              !飽和蒸気圧 [Pa]
  public ss_DSvapPressDTemp_ss, DSvapPressDTemp   !飽和蒸気圧の温度微分 [Pa/K]
  public ss_LatentHeat_ss, LatentHeat             !潜熱 [J/K kg]

  integer       :: XMin = 10   !X 座標範囲(下限). 値はダミー
  integer       :: XMax = 10   !X 座標範囲(上限). 値はダミー
  integer       :: ZMin = 10   !Z 座標範囲(下限). 値はダミー
  integer       :: ZMax = 10   !Z 座標範囲(上限). 値はダミー
  integer       :: SpcNum = 1  !化学種の数. 値はダミー

  integer, allocatable :: SpcName(:)   !系に含まれる化学種. 数字で判別
  real(8), allocatable :: MolWt(:)     !分子量
  real(8), allocatable :: ACoeff(:,:)  !Antoine の飽和蒸気圧式の係数
  
  save XMin, XMax, ZMin, ZMax, SpcNum
  save SpcName, MolWt, ACoeff
  
contains

  subroutine ChemEquiv_Init(IMin, IMax, KMin, KMax, SNum, SName)
    !
    !初期化ルーチン. 
    !配列サイズの確定と, Antoine の飽和蒸気圧式の係数を決定する
    !

    !モジュール呼び出し
    use ChemDataBase, only: ChemData_SvapPress_Antoine, &  !Antoine 式の係数
      &                     ChemData_MolWt                 !分子量

    !暗黙の型宣言禁止
    implicit none

    !入出力変数
    integer,intent(in) :: IMin, IMax      !水平方向の添え字の下限・上限
    integer,intent(in) :: KMin, KMax      !鉛直方向の添え字の下限・上限
    integer,intent(in) :: SNum            !化学種の個数
    integer,intent(in) :: SName(SpcNum)   !系に存在する化学種の数

    !内部変数
    integer            :: s

    !-----------------------------------------------------------
    ! 初期化
    !-----------------------------------------------------------
    XMin = IMin;    XMax = IMax
    ZMin = KMin;    ZMax = KMax
    SpcNum  = SNum
    
    !配列割り当て
    allocate( SpcName(SpcNum), MolWt(SpcNum), ACoeff(SpcNum, 3) )

    !-----------------------------------------------------------    
    ! データベースから必要なデータだけ取り出す
    !-----------------------------------------------------------
    !系に存在する化学種. 入力された値を取り込む
    SpcName = SName

!    !単位質量当たりの比熱
!    do s = 1, SpcNum
!      Cp(s) = ChemData_CpPerMass(SpcName(s)) 
!    end do

    !分子量
    do s = 1, SpcNum
      MolWt(s) = ChemData_MolWt(SpcName(s)) 
    end do

    !Antoine の係数
    do s = 1, SpcNum
      ACoeff(s,:) = ChemData_SvapPress_Antoine(SpcName(s),:) 
    end do
    
  end subroutine ChemEquiv_Init
  

  function ss_SvapPress_ss( ss_Temp )
    !
    !飽和蒸気圧の計算. Antoine の式を利用
    !

    !暗黙の型宣言禁止
    implicit none
    
    ! 入出力変数  
    real(8)             :: ss_SvapPress_ss(XMin:XMax, ZMin:ZMax, SpcNum) 
                                                             !飽和蒸気圧
    real(8), intent(in) :: ss_Temp(XMin:XMax, ZMin:ZMax)     !温度
  
    !内部変数
    integer             :: s

    !初期化
    ss_SvapPress_ss = 0.0d0 
    
    do s = 1, SpcNum
      if (ACoeff(s,1) == 0) cycle   !気相の場合は係数がゼロ
      
      ss_SvapPress_ss(:,:,s) = &
        & dexp( ACoeff(s,1) - ACoeff(s,2) / (ACoeff(s,3) + ss_Temp(:,:)) )
    end do
    
  end function ss_SvapPress_ss


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

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

    !内部変数
    integer             :: s

    !初期化
    SvapPress = 0.0d0 
    s = minloc(SpcName, 1, SpcName == Name)   !入力名に対応する配列添え字を得る

    !飽和蒸気圧の計算
    SvapPress = dexp( ACoeff(s,1) - ACoeff(s,2) / (ACoeff(s,3) + Temp) ) 
    
  end function SvapPress


  function ss_DSvapPressDTemp_ss(ss_Temp)
    !
    !飽和蒸気圧の温度微分の計算. Antoine の式を利用
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    !入出力変数
    real(8)             :: ss_DSvapPressDTemp_ss(XMin:XMax, ZMin:ZMax, SpcNum)
                                                  !飽和蒸気圧の温度微分 [Pa/K]
    real(8), intent(in) :: ss_Temp(XMin:XMax, ZMin:ZMax)             !温度 [K]
  
    !内部変数
    integer             :: s

    !初期化
    ss_DSvapPressDTemp_ss = 0.0d0

    do s = 1, SpcNum
      if (ACoeff(s,1) == 0) cycle   !気相の場合は係数がゼロ

      ss_DSvapPressDTemp_ss(:,:,s) = &
        & dexp( ACoeff(s,1) - ACoeff(s,2) / (ACoeff(s,3) + ss_Temp(:,:)) ) &
        & * ( ACoeff(s,2) / ((ACoeff(s,3) + ss_Temp(:,:)) ** 2.0d0) )
    end do
    
  end function ss_DSvapPressDTemp_ss
  

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

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

    !内部変数
    integer             :: s

    !初期化
    DSvapPressDTemp = 0.0d0
    s = minloc(SpcName, 1, SpcName == Name)   !入力名に対応する配列添え字を得る

    !飽和蒸気圧の温度微分
    DSvapPressDTemp = &
      & dexp( ACoeff(s,1) - ACoeff(s,2) / (ACoeff(s,3) + Temp) ) &
      & * ( ACoeff(s,2) / ((ACoeff(s,3) + Temp) ** 2.0d0) )
    
  end function DSvapPressDTemp


  function ss_LatentHeat_ss(ss_Temp)
    !
    !飽和蒸気圧より潜熱を計算する. Antoine の式を利用
    !

    !暗黙の型宣言禁止
    implicit none

    !入出力変数
    real(8)             :: ss_LatentHeat_ss(XMin:XMax, ZMin:ZMax, SpcNum)
                                                            !潜熱[J/K kg]
    real(8), intent(in) :: ss_Temp(XMin:XMax, ZMin:ZMax)    !温度[K]

    !内部変数
    integer             :: s
    real(8)             :: ss_DSvapPressLogDTemp(XMin:XMax, ZMin:ZMax, SpcNum)
                                               !飽和蒸気圧の対数の温度微分 
    real(8), parameter  :: GasRUniv = 8.314d0  !普遍気体定数

    !初期化
    ss_DSvapPressLogDTemp = 0.0d0
    ss_LatentHeat_ss = 0.0d0

    do s = 1, SpcNum
      if (ACoeff(s,1) == 0) cycle    !気相の場合は係数がゼロ
      
      ss_DSvapPressLogDTemp(:,:,s) = &
        & ACoeff(s,2) / ((ACoeff(s,3) + ss_Temp(:,:)) ** 2.0d0)
      
      ss_LatentHeat_ss(:,:,s) = &
        & ss_DSvapPressLogDTemp(:,:,s) * GasRUniv * (ss_Temp(:,:) ** 2.0d0) &
        & / MolWt(s)
    end do

  end function ss_LatentHeat_ss


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

    !暗黙の型宣言禁止
    implicit none

    !入出力変数
    real(8), intent(in) :: Name                !化学種名
    real(8), intent(in) :: Temp                !温度
    real(8)             :: LatentHeat          !潜熱
    
    !内部変数
    integer             :: s
    real(8)             :: DSvapPressLogDTemp  !飽和蒸気圧の対数の温度微分 
    real(8), parameter  :: GasRUniv = 8.314d0  !普遍気体定数
    
    !初期化
    DSvapPressLogDTemp = 0.0d0
    LatentHeat         = 0.0d0
    s = minloc(SpcName, 1, SpcName == Name)   !入力名に対応する配列添え字を得る

    !潜熱の計算
    DSvapPressLogDTemp = ACoeff(s,2) / ((ACoeff(s,3) + Temp) ** 2.0d0)
    LatentHeat = DSvapPressLogDTemp * GasRUniv * (Temp ** 2.0d0) / MolWt(s)
  
  end function LatentHeat
  
end module ChemEquiv
