!= Module ChemCalc
!
! Authors::   SUGIYAMA Ko-ichiro
! Version::   $Id: chemcalc.f90,v 1.12 2006/09/30 05:53:42 odakker Exp $
! Tag Name::  $Name: arare4-20061019 $
! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview
!
!化学関連の諸量を計算するためのモジュール. 
!AMP と Antoine の飽和蒸気圧式を用いて以下を求める. 
!デフォルトでは AMP 式を使うようにしてある. 
!  * 飽和蒸気圧
!  * 飽和蒸気圧の温度微分
!  * 潜熱
! 
!== Error Handling
!
!== Bugs
!
!== Note
!
!== Future Plans
!
!

module ChemCalc
  !
  !化学関連の諸量を計算するためのモジュール. 
  !AMP と Antoine の飽和蒸気圧式を用いて以下を求める. 
  !デフォルトでは AMP 式を使うようにしてある. 
  !  * 飽和蒸気圧
  !  * 飽和蒸気圧の温度微分
  !  * 潜熱
  !

  !モジュール呼び出し
  use dc_message, only: MessageNotify

  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,& !単位変換用係数
    &                 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

  !全体に対して private 属性を付ける. 
  private
  
  !公開するものには public 属性を付ける
  public ChemCalc_Init                            !初期化ルーチン
  public ChemCalcSpc_Init, ChemCalcDim_Init       !初期化ルーチン
  public MolWt                                    !分子量
  public GasR                                     !気体定数
  public Phase                                    !相
  public CpRef, CpPerMolRef, CvRef                !定圧比熱, 定積比熱
  public SvapPress, xz_SvapPress                  !飽和蒸気圧 [Pa]
  public xz_DSvapPressDTemp                       !飽和蒸気圧の温度微分 [Pa/K]
  public xz_LatentHeat                            !潜熱 [J/K kg]
  public LatentHeatPerMol                         !潜熱 [J/K mol]
  public ReactHeatNH4SH, ReactHeatNH4SHPerMol     !NH4SH 反応熱

  real(8)             :: ReactHeatNH4SH           !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_ampA(ChemData_SpcNum)   !AMP 式の蒸気圧式の A 係数
  real(8)             :: a_ampB(ChemData_SpcNum)   !AMP 式の蒸気圧式の B 係数
  real(8)             :: a_ampC(ChemData_SpcNum)   !AMP 式の蒸気圧式の C 係数
  real(8)             :: a_ampD(ChemData_SpcNum)   !AMP 式の蒸気圧式の D 係数
  real(8)             :: a_ampE(ChemData_SpcNum)   !AMP 式の蒸気圧式の E 係数
  real(8)             :: a_MolWt(ChemData_SpcNum)  !分子量
  integer, allocatable      :: SpcID(:)     !系に含まれる化学種. 数字で判別
  character(3), allocatable :: Phase(:)     !相
 
  save ReactHeatNH4SH, ReactHeatNH4SHPerMol
  save XMin, XMax, ZMin, ZMax, SpcNum
  save SpcID, Phase
  save a_antA, a_antB, a_antC, a_antU, a_MolWt
  save a_ampA, a_ampB, a_ampC, a_ampD, a_ampE

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

    !モジュール呼び出し
    use ChemData, only: GasRUniv,      &   !気体定数
      &                 ChemData_OneSpcID  !化学種の 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
    integer            :: id

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

    !NH4SH の反応熱の初期化
    !  NH4SH 1kg に対する反応熱にする.
    Name = 'NH4SH-s'
    id   = ChemData_OneSpcID( Name )  
    
    ReactHeatNH4SHPerMol  = GasRUniv * 10834.0d0
    ReactHeatNH4SH = GasRUniv * 10834.0d0 / MolWt( id )
!    write(*,*)     ReactHeatNH4SH, 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


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

    call MessageNotify( "M", &
      & "ChemCalcSpc_Init", "SpcNum = %d", i=(/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

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

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

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

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

  end function CpRef


!!!==========================================================================
  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 CvRef(ID)
    !
    !引数で与えられた化学種に対して, 標準状態での単位質量当たりの定圧比熱を計算
    !

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

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

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

  end function CvRef


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



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

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

    !暗黙の型宣言禁止
    implicit none
    
    ! 入出力変数  
    real(8)             :: xz_SvapPressAnt(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 を計算
    !対数が大きくなりすぎないようにする. 
    ! Fujitsu Fortran Compiler では 700 より大きい数の exp を取ると警告が出る.
    xz_LogSvapPress =                                          &
      & min(                                                   & 
      &       (                                                &
      &         a_antA(ID)                                     &
      &         - a_antB(ID) / (a_antC(ID) + xz_Temp - Temp0C) &
      &        ) * dlog(10.0d0)                                &
      &        + a_antU(ID),                                   &
      &      7.0d2                                             &
      &    )
    
    !飽和蒸気圧を計算
    xz_SvapPressAnt =  dexp( xz_LogSvapPress )

  end function xz_SvapPressAnt


!!!==========================================================================
  function xz_DSvapPressDTempAnt(ID, xz_Temp)
    !
    !飽和蒸気圧の温度微分の計算. Antoine の式を利用
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    !入出力変数
    real(8)             :: xz_DSvapPressDTempAnt(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 を計算
    !対数が大きくなりすぎないようにする. 
    ! Fujitsu Fortran Compiler では 700 より大きい数の exp を取ると警告が出る.
    xz_LogSvapPress =                                          &
      & min(                                                   &
      &      (                                                 &
      &          a_antA(ID)                                    &
      &        - a_antB(ID) / (a_antC(ID) + xz_Temp - Temp0C)  &
      &       ) * dlog(10.0d0)                                 &
      &       + a_antU(ID),                                    &
      &     7.0d2                                              &
      &    )
    
    xz_DLogSvapPressDTemp =                                &
      & a_antB(ID) * dlog(10.0d0)                          &
      &  / ( (a_antC(ID) + xz_Temp - Temp0C) ** 2.0d0 )
    
    !飽和蒸気圧の温度微分を計算
    xz_DSvapPressDTempAnt = xz_DLogSvapPressDTemp * dexp( xz_LogSvapPress ) 
    
  end function xz_DSvapPressDTempAnt
  

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

    !入出力変数
    real(8)             :: xz_LatentHeatAnt(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_LatentHeatAnt =                                        &
      & xz_DLogSvapPressDTemp * GasRUniv * (xz_Temp ** 2.0d0) &
      &  / a_MolWt(ID)

  end function xz_LatentHeatAnt


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

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

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

    !飽和蒸気圧を計算
    xz_SvapPress =  dexp( xz_LogSvapPress )
    
  end function xz_SvapPress


!!!==========================================================================
  function xz_DSvapPressDTemp(ID, xz_Temp)
    !
    !引数で与えられた化学種と温度に対して, 飽和蒸気圧の温度微分を計算
    !
    
    !暗黙の型宣言禁止
    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                            !化学種名

    !内部変数
    real(8)             :: xz_LogSvapPress(XMin:XMax, ZMin:ZMax)
    real(8)             :: xz_DLogSvapPressDTemp(XMin:XMax, ZMin:ZMax)
    

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

    !飽和蒸気圧の温度微分
    xz_DLogSvapPressDTemp =                   &
      & - a_ampA(ID) / (xz_Temp ** 2.0d0)        &
      & + a_ampC(ID) / xz_Temp                   &
      & + a_ampD(ID)                             &
      & + a_ampE(ID) * 2.0d0 * xz_Temp
    
    !飽和蒸気圧の温度微分
    xz_DSvapPressDTemp =                          &
      & xz_DLogSvapPressDTemp * dexp( xz_LogSvapPress ) 
    
  end function xz_DSvapPressDTemp
  

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

    !暗黙の型宣言禁止
    implicit none

    !入出力変数
    real(8)             :: xz_LatentHeat(XMin:XMax, ZMin:ZMax)
                                                         !潜熱
    real(8), intent(in) :: xz_Temp(XMin:XMax, ZMin:ZMax) !温度
    integer, intent(in) :: ID                            !化学種名
    
    !内部変数
    real(8)             :: xz_DLogSvapPressDTemp(XMin:XMax, ZMin:ZMax)
    real(8), parameter  :: GasRUniv = 8.314d0  !普遍気体定数
    

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

  end function xz_LatentHeat


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

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

    !内部変数
    real(8)             :: 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 LatentHeatPerMol(ID, Temp)
    !
    !引数で与えられた化学種と温度に対して, 潜熱を計算
    !

    !暗黙の型宣言禁止
    implicit none

    !入出力変数
    real(8)             :: LatentHeatPerMol   !潜熱
    real(8), intent(in) :: Temp            !温度
    integer, intent(in) :: ID                 !化学種名
    
    !内部変数
    real(8)             :: DLogSvapPressDTemp
    real(8), 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


end module ChemCalc

