!= Module MoistFunc
!
! Authors::   SUGIYAMA Ko-ichiro
! Version::   $Id: moistfunc.f90,v 1.13 2007/04/19 10:46:03 sugiyama Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview
!
!便利な計算ツール群
!   * 蒸気圧から混合比へ変換
!   * モル比から混合比へ変換
!   * 混合比から蒸気圧へ変換
!   * 混合比からモル比へ変換
! 
!== Error Handling
!
!== Bugs
!
!== Note
!
!== Future Plans
!
!

module MoistFunc
  !
  !単位の換算を行うためのパッケージ型モジュール
  !   * 蒸気圧から混合比へ変換
  !   * モル比から混合比へ変換
  !   * 混合比から蒸気圧へ変換
  !   * 混合比からモル比へ変換
  ! 

  !モジュール読み込み
  use gridset, only:  DimXMin,            & ! x 方向の配列の下限
    &                 DimXMax,            & ! x 方向の配列の上限
    &                 DimZMin,            & ! z 方向の配列の下限
    &                 DimZMax               ! z 方向の配列の上限
  use basicset, only: PressBasis,         & !温位の標準圧力         [Pa]
    &                 CpDry,              & !乾燥成分の平均定圧比熱 [J/K kg]
    &                 MolWtDry,           & !乾燥成分の平均分子量   [kg/mol]
    &                 GasRDry               !乾燥成分の気体定数     [J/K kg]
  use ChemCalc, only: xz_DSvapPressDTemp    !飽和蒸気圧の温度微分

  !暗黙の型宣言禁止
  implicit none 
  
  !属性の指定
  private

  !関数を public に指定
  public xz_DMixRtSatDPotTemp
  public xz_DelMixRtNH4SH
  public DelMolFrNH4SH

contains

!!!
!!!基本関数
!!!

!!!-----------------------------------------------------------------------!!!
  function xz_DMixRtSatDPotTemp(SpcID, MolWt, xz_TempAll, xz_ExnerAll)
    !
    !飽和蒸気圧の θ 微分を行う
    !実際には, dq/dp * dp/dT * dT/dθ を実行. (但し p は飽和蒸気圧)
    !

    !暗黙の型宣言禁止
    implicit none 
    
    !入出力変数
    integer, intent(in) :: SpcID
    real(8), intent(in) :: MolWt
    real(8), intent(in) :: xz_TempAll(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !温度(擾乱 + 基本場)
    real(8), intent(in) :: xz_ExnerAll(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !エクスナー関数(擾乱 + 基本場)
    real(8)             :: xz_PressAll(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !圧力(擾乱 + 基本場)
    real(8)             :: xz_DMixRtSatDPotTemp(DimXMin:DimXMax, DimZMin:DimZMax)
                           
    xz_DMixRtSatDPotTemp   = 0.0d0
    xz_PressAll = PressBasis * (xz_ExnerAll ** (CpDry / GasRDry))
    
    xz_DMixRtSatDPotTemp =                            &
      &   MolWt / ( MolWtDry * xz_PressAll )          &
      & * xz_DSvapPressDTemp(SpcID, xz_TempAll)       &
      & * xz_ExnerAll
    
  end function xz_DMixRtSatDPotTemp


!!!-----------------------------------------------------------------------!!!
  function xz_DelMixRtNH4SH(xz_TempAll, xz_PressAll, xz_MixRtNH3, xz_MixRtH2S, &
    &                    MolWtNH3, MolWtH2S)
    !
    ! NH4SH 生成反応に伴う, NH4SH の生成量(混合比)を求める
    !
    
    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(8), intent(in) :: xz_TempAll(DimXMin:DimXMax, DimZMin:DimZMax)
                                         !温度
    real(8), intent(in) :: xz_PressAll(DimXMin:DimXMax, DimZMin:DimZMax)
                                         !圧力
    real(8), intent(in) :: xz_MixRtNH3(DimXMin:DimXMax, DimZMin:DimZMax)
                                         !NH3 の混合比
    real(8), intent(in) :: xz_MixRtH2S(DimXMin:DimXMax, DimZMin:DimZMax)
                                         !H2S の混合比
    real(8), intent(in) :: MolWtNH3      !NH3 の分子量
    real(8), intent(in) :: MolWtH2S      !H2S の分子量
    real(8)             :: xz_DelMixRtNH4SH(DimXMin:DimXMax, DimZMin:DimZMax)
                                         !NH4SH の混合比
    real(8)             :: xz_EquivConst(DimXMin:DimXMax, DimZMin:DimZMax)
                                         !圧平衡定数
    real(8)             :: xza_PartialPress(DimXMin:DimXMax, DimZMin:DimZMax,2)
                                         !作業配列(分圧)
    real(8)             :: xz_Solution(DimXMin:DimXMax, DimZMin:DimZMax)
                                         !作業配列(方程式の解)

    !初期化
    xz_DelMixRtNH4SH = 0.0d0

    !アンモニアと硫化水素の分圧. 
    xza_PartialPress(:,:,1) = xz_MixRtNH3 * xz_PressAll * MolWtDry / MolWtNH3 
    xza_PartialPress(:,:,2) = xz_MixRtH2S * xz_PressAll * MolWtDry / MolWtH2S      
    
    !圧平衡定数
    xz_EquivConst = 61.781d0 - 10834.0d0 / xz_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 を
    !  考慮すると, 解の公式のうちが負の方が選択される.
    xz_Solution  =                                                   &
      & (                                                            &
      &     sum(xza_PartialPress, 3)                                        &
      &   - dsqrt( (xza_PartialPress(:,:,1) - xza_PartialPress(:,:,2)) ** 2.0d0 &
      &            + 4.0d0 * dexp( min( 700.0d0, xz_EquivConst ) ) ) &
      &  ) * 5.0d-1

    !生成量を求める
    xz_DelMixRtNH4SH = xz_Solution * ( MolWtNH3 + MolWtH2S ) / ( xz_PressAll * MolWtDry )

  end function xz_DelMixRtNH4SH
  

!!!-----------------------------------------------------------------------!!!
  function DelMolFrNH4SH(TempAll, PressAll, MolFrNH3, MolFrH2S, Humidity)
    !
    ! NH4SH 生成反応に伴う H2S と NH3 のモル比の減少分を求める
    !
    
    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(8), intent(in) :: TempAll       !温度
    real(8), intent(in) :: PressAll      !圧力
    real(8), intent(in) :: MolFrNH3      !NH3 のモル比
    real(8), intent(in) :: MolFrH2S      !H2S のモル比
    real(8), intent(in) :: Humidity      !飽和比
    real(8)             :: DelMolFrNH4SH !NH4SH 生成に伴うモル比変化
    real(8)             :: EquivConst    !圧平衡定数
    real(8)             :: PPress(2)     !作業配列(分圧)
    real(8)             :: 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 MoistFunc
