!= Module ChemCalc_3d
!
! Authors::   SUGIYAMA Ko-ichiro
! Version::   $Id: chemcalc.f90,v 1.5 2012-01-12 11:53:10 sugiyama Exp $
! Tag Name::  $Name: arare5-20120511 $
! 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_types,       only : DP

  use dc_message, only: MessageNotify

  use ChemData, only: ChemData_init,                 & !初期化
    &                 ChemData_SpcNum,               & !データベース上の化学種数
    &                 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 係数
    &                 GasRUniv,                      & !気体定数
    &                 ChemData_OneSpcID,             & !化学種の ID 検索     
    &                 ChemData_CpRef,                & !標準状態での単位質量当たりの比熱
    &                 ChemData_CpPerMolRef,          & !標準状態での単位モル当たりの比熱
    &                 ChemData_CvRef,                & !標準状態での単位質量当たりの比熱
    &                 ChemData_MolWt,                & !分子量
    &                 ChemData_GasR                    !気体定数 [J/K kg]
  use gridset, only:  imin,                     &! 配列の X 方向の下限
    &                 imax,                        &! 配列の X 方向の上限
    &                 jmin,                        &! 配列の Y 方向の下限
    &                 jmax,                        &! 配列の Y 方向の上限
    &                 kmin,                        &! 配列の Z 方向の下限
    &                 kmax                          ! 配列の Z 方向の上限
  use constants,only: PressBasis,         & !温位の標準圧力         [Pa]
    &                 CpDry,              & !乾燥成分の平均定圧比熱 [J/K kg]
    &                 MolWtDry,           & !乾燥成分の平均分子量   [kg/mol]
    &                 GasRDry               !乾燥成分の気体定数     [J/K kg]

  !暗黙の型宣言禁止
  implicit none

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

  real(DP)            :: ReactHeatNH4SH           !NH4SH 生成反応熱 [J/K kg]
  real(DP)            :: ReactHeatNH4SHPerMol     !NH4SH 生成反応熱 [J/K mol]
  real(DP)            :: a_antA(ChemData_SpcNum)   !Antoine の蒸気圧式の A 係数
  real(DP)            :: a_antB(ChemData_SpcNum)   !Antoine の蒸気圧式の B 係数
  real(DP)            :: a_antC(ChemData_SpcNum)   !Antoine の蒸気圧式の C 係数
  real(DP)            :: a_antU(ChemData_SpcNum)   !Antoine の蒸気圧式の単位換算のための係数
  real(DP)            :: a_ampA(ChemData_SpcNum)   !AMP 式の蒸気圧式の A 係数
  real(DP)            :: a_ampB(ChemData_SpcNum)   !AMP 式の蒸気圧式の B 係数
  real(DP)            :: a_ampC(ChemData_SpcNum)   !AMP 式の蒸気圧式の C 係数
  real(DP)            :: a_ampD(ChemData_SpcNum)   !AMP 式の蒸気圧式の D 係数
  real(DP)            :: a_ampE(ChemData_SpcNum)   !AMP 式の蒸気圧式の E 係数
  real(DP)            :: a_MolWt(ChemData_SpcNum)  !分子量
 
  save ReactHeatNH4SH, ReactHeatNH4SHPerMol
  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( )
    !
    !初期化ルーチン
    !

    !暗黙の型宣言禁止
    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



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

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

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

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

  end function CpRef


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

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

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

  end function CpPerMolRef


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

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

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

  end function CvRef


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

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

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

  end function MolWt


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

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

  end function GasR



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

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

    !暗黙の型宣言禁止
    implicit none
    
    ! 入出力変数  
    real(DP)            :: xyz_SvapPressAnt(imin:imax,jmin:jmax,kmin:kmax) 
                                                          !飽和蒸気圧
    real(DP),intent(in) :: xyz_Temp(imin:imax,jmin:jmax,kmin:kmax)  
                                                          !温度
    integer, intent(in) :: ID                             !化学種の ID
  
    !内部変数
    real(DP)            :: xyz_LogSvapPress(imin:imax,jmin:jmax,kmin:kmax) 
    real(DP),parameter  :: Temp0C = 273.15d0


    !飽和蒸気圧の log を計算
    !対数が大きくなりすぎないようにする. 
    ! Fujitsu Fortran Compiler では 700 より大きい数の exp を取ると警告が出る.
    xyz_LogSvapPress =                                          &
      & min(                                                   & 
      &       (                                                &
      &         a_antA(ID)                                     &
      &         - a_antB(ID) / (a_antC(ID) + xyz_Temp - Temp0C) &
      &        ) * dlog(10.0d0)                                &
      &        + a_antU(ID),                                   &
      &      7.0d2                                             &
      &    )
    
    !飽和蒸気圧を計算
    xyz_SvapPressAnt =  dexp( xyz_LogSvapPress )

  end function xyz_SvapPressAnt


!!!==========================================================================
  function xyz_DSvapPressDTempAnt(ID, xyz_Temp)
    !
    !飽和蒸気圧の温度微分の計算. Antoine の式を利用
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    !入出力変数
    real(DP)            :: xyz_DSvapPressDTempAnt(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                     !化学種の ID
  
    !内部変数
    real(DP)            :: xyz_LogSvapPress(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)            :: xyz_DLogSvapPressDTemp(imin:imax,jmin:jmax,kmin:kmax)
    real(DP),parameter  :: Temp0C = 273.15d0


    !飽和蒸気圧の log を計算
    !対数が大きくなりすぎないようにする. 
    ! Fujitsu Fortran Compiler では 700 より大きい数の exp を取ると警告が出る.
    xyz_LogSvapPress =                                          &
      & min(                                                   &
      &      (                                                 &
      &          a_antA(ID)                                    &
      &        - a_antB(ID) / (a_antC(ID) + xyz_Temp - Temp0C)  &
      &       ) * dlog(10.0d0)                                 &
      &       + a_antU(ID),                                    &
      &     7.0d2                                              &
      &    )
    
    xyz_DLogSvapPressDTemp =                                &
      & a_antB(ID) * dlog(10.0d0)                          &
      &  / ( (a_antC(ID) + xyz_Temp - Temp0C) ** 2.0d0 )
    
    !飽和蒸気圧の温度微分を計算
    xyz_DSvapPressDTempAnt = xyz_DLogSvapPressDTemp * dexp( xyz_LogSvapPress ) 
    
  end function xyz_DSvapPressDTempAnt
  

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

    !入出力変数
    real(DP)            :: xyz_LatentHeatAnt(imin:imax,jmin:jmax,kmin:kmax)
                                                            !潜熱[J/K kg]
    real(DP),intent(in) :: xyz_Temp(imin:imax,jmin:jmax,kmin:kmax)
                                                    !温度[K]
    integer, intent(in) :: ID                       !化学種の ID

    !内部関数
    real(DP)            :: xyz_DLogSvapPressDTemp(imin:imax,jmin:jmax,kmin:kmax)
    real(DP),parameter  :: GasRUniv = 8.314d0
    real(DP),parameter  :: Temp0C = 273.15d0

    !飽和蒸気圧の log を計算した後, 潜熱を計算する. 
    xyz_DLogSvapPressDTemp =                                   &
      & a_antB(ID) * dlog(10.0d0)                             &
      &  / ( (a_antC(ID) + xyz_Temp - Temp0C) ** 2.0d0 )
    
    xyz_LatentHeatAnt =                                        &
      & xyz_DLogSvapPressDTemp * GasRUniv * (xyz_Temp ** 2.0d0) &
      &  / a_MolWt(ID)

  end function xyz_LatentHeatAnt


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


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

end module ChemCalc

