!= Module MoistFunc_3d
!
! Authors::   SUGIYAMA Ko-ichiro ODAKA Masatsugu
! Version::   $Id: moistfunc_3d.f90,v 1.1 2008-06-19 16:53:25 odakker 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_3d
  !
  !単位の換算を行うためのパッケージ型モジュール
  !   * 蒸気圧から混合比へ変換
  !   * モル比から混合比へ変換
  !   * 混合比から蒸気圧へ変換
  !   * 混合比からモル比へ変換
  ! 

  !モジュール読み込み
  use dc_types,       only : STRING, DP

  use gridset, only : DimXMin, DimXMax,   & !x 方向の配列の下限と上限
    &                    DimYMin, DimYMax,   & !y 方向の配列の下限と上限
    &                    DimZMin, DimZMax      !z 方向の配列の下限と上限
  use constants,  only : PressBasis,         & !温位の標準圧力         [Pa]
    &                    CpDry,              & !乾燥成分の平均定圧比熱 [J/K kg]
    &                    MolWtDry,           & !乾燥成分の平均分子量   [kg/mol]
    &                    GasRDry               !乾燥成分の気体定数     [J/K kg]
  use ChemCalc_3d,only : xyz_DSvapPressDTemp   !飽和蒸気圧の温度微分

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

  !関数を public に指定
  public xyz_DMixRtSatDPotTemp
  public xyz_DelMixRtNH4SH
  public DelMolFrNH4SH

contains

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

!!!-----------------------------------------------------------------------!!!
  function xyz_DMixRtSatDPotTemp(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(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                            !温度(擾乱 + 基本場)
    real(DP),intent(in) :: xyz_ExnerAll(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                            !エクスナー関数(擾乱 + 基本場)
    real(DP)            :: xyz_PressAll(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                            !圧力(擾乱 + 基本場)
    real(DP)            :: xyz_DMixRtSatDPotTemp(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                           
    xyz_DMixRtSatDPotTemp   = 0.0d0
    xyz_PressAll = PressBasis * (xyz_ExnerAll ** (CpDry / GasRDry))
    
    xyz_DMixRtSatDPotTemp =                            &
      &   MolWt / ( MolWtDry * xyz_PressAll )          &
      & * xyz_DSvapPressDTemp(SpcID, xyz_TempAll)       &
      & * xyz_ExnerAll
    
  end function xyz_DMixRtSatDPotTemp


!!!-----------------------------------------------------------------------!!!
  function xyz_DelMixRtNH4SH(xyz_TempAll, xyz_PressAll, xyz_MixRtNH3, xyz_MixRtH2S, &
    &                    MolWtNH3, MolWtH2S)
    !
    ! NH4SH 生成反応に伴う, NH4SH の生成量(混合比)を求める
    !
    
    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(DP),intent(in) :: xyz_TempAll(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                         !温度
    real(DP),intent(in) :: xyz_PressAll(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                         !圧力
    real(DP),intent(in) :: xyz_MixRtNH3(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                         !NH3 の混合比
    real(DP),intent(in) :: xyz_MixRtH2S(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                         !H2S の混合比
    real(DP),intent(in) :: MolWtNH3      !NH3 の分子量
    real(DP),intent(in) :: MolWtH2S      !H2S の分子量

    real(DP) :: xyz_DelMixRtNH4SH(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                         !NH4SH の混合比
    real(DP) :: xyz_EquivConst(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                         !圧平衡定数
    real(DP) :: xyza_PartialPress(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax,2)
                                         !作業配列(分圧)
    real(DP) :: xyz_Solution(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                         !作業配列(方程式の解)

    !初期化
    xyz_DelMixRtNH4SH = 0.0d0

    !アンモニアと硫化水素の分圧. 
    xyza_PartialPress(:,:,:,1) = xyz_MixRtNH3 * xyz_PressAll * MolWtDry / MolWtNH3 
    xyza_PartialPress(:,:,:,2) = xyz_MixRtH2S * 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(xyza_PartialPress, 4)                                        &
      &   - dsqrt( (xyza_PartialPress(:,:,:,1) - xyza_PartialPress(:,:,:,2)) ** 2.0d0 &
      &            + 4.0d0 * dexp( min( 700.0d0, xyz_EquivConst ) ) ) &
      &  ) * 5.0d-1

    !生成量を求める
    xyz_DelMixRtNH4SH = xyz_Solution * ( MolWtNH3 + MolWtH2S ) &
      &                   / ( xyz_PressAll * MolWtDry )

  end function xyz_DelMixRtNH4SH
  

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