!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2005, 2006. All rights reserved.
!---------------------------------------------------------------------
!= Module MoistFunc
!
!   * Developer: SUGIYAMA Ko-ichiro
!   * Version: $Id: $
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview
!
!便利な計算ツール群
!   * 蒸気圧から混合比へ変換
!   * モル比から混合比へ変換
!   * 混合比から蒸気圧へ変換
!   * 混合比からモル比へ変換
! 
!== Error Handling
!
!== Bugs
!
!== Note
!
!== Future Plans
!
!

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

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

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

  !関数を public に指定
  public Vap2MixRt, xz_Vap2MixRt 
  public MixRt2Vap, xz_MixRt2Vap
  public MolFr2MixRt, xz_MolFr2MixRt
  public MixRt2MolFr
  public Exner2Press, xz_Exner2Press
  public DMixRtSatDPotTemp, xz_DMixRtSatDPotTemp
  public DelMixRtNH4SH, DelMolFrNH4SH
  public xz_DelMixRtNH4SH
contains

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

!!!-----------------------------------------------------------------------!!!
  function Vap2MixRt(MolWt, VapAll, ExnerAll)
    !
    !蒸気圧から比湿を求める. 
    !凝縮成分は少ないため, 混合比は比湿にほぼ等しい
    !
    
    !暗黙の型宣言禁止
    implicit none 
    
    !入出力変数
    real(8), intent(in) :: MolWt            !化学種の分子量
    real(8), intent(in) :: VapAll           !蒸気圧(擾乱 + 基本場)
    real(8), intent(in) :: ExnerAll         !エクスナー関数(擾乱 + 基本場)
    real(8)             :: Vap2MixRt        !蒸気圧から比湿への変換
    
    Vap2MixRt = 0.0d0
    Vap2MixRt = VapAll * MolWt / (MolWtDry * Exner2Press(ExnerAll))
    
  end function Vap2MixRt

  
  function xz_Vap2MixRt(MolWt, xz_VapAll, xz_ExnerAll)
    !
    !蒸気圧から比湿を求める. 
    !凝縮成分は少ないため, 混合比は比湿にほぼ等しい
    !
    
    !暗黙の型宣言禁止
    implicit none 
    
    !入出力変数
    real(8), intent(in) :: MolWt   !化学種の分子量
    real(8), intent(in) :: xz_VapAll(DimXMin:DimXMax, DimZMin:DimZMax) 
                                   !蒸気圧(擾乱成分と基本場の和)   [Pa]
    real(8), intent(in) :: xz_ExnerAll(DimXMin:DimXMax, DimZMin:DimZMax) 
                                   !エクスナー関数(擾乱成分と基本場の和)
    real(8)             :: xz_Vap2MixRt(DimXMin:DimXMax, DimZMin:DimZMax)
                                   !蒸気圧から混合比への変換
    integer             :: i, k

    do k = DimZMin, DimZMax
      do i = DimXMin, DimXMax

        xz_Vap2MixRt(i,k) = Vap2MixRt(MolWt, xz_VapAll(i,k), xz_ExnerAll(i,k))
      end do
    end do

!    write(*,*) minval(xz_VapAll(95,:)), maxval(xz_VapAll(95,:))
!    write(*,*) minval(xz_Vap2MixRt(95,:)), maxval(xz_Vap2MixRt(95,:))

  end function xz_Vap2MixRt
  
  
!!!-----------------------------------------------------------------------!!!
  function MolFr2MixRt(MolWt, MolFr)
    !
    !モル比から比湿を求める. 
    !
    
    !暗黙の型宣言禁止
    implicit none 
    
    !入出力変数
    real(8), intent(in) :: MolWt            !分子量
    real(8), intent(in) :: MolFr            !モル比(擾乱 + 基本場)
    real(8)             :: MolFr2MixRt      !モル比から混合比への変換

    MolFr2MixRt = 0.0d0
    MolFr2MixRt = MolFr * MolWt / MolWtDry     

  end function MolFr2MixRt


  function xz_MolFr2MixRt(MolWt, xz_MolFr)
    !
    !モル比から比湿を求める. 
    !
    
    !暗黙の型宣言禁止
    implicit none 
    
    !入出力変数
    real(8), intent(in) :: MolWt            !分子量
    real(8), intent(in) :: xz_MolFr(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !モル比(擾乱 + 基本場)
    real(8)             :: xz_MolFr2MixRt(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !モル比から混合比への変換
    integer             :: i, k

    do k = DimZMin, DimZMax
      do i = DimXMin, DimXMax
        xz_MolFr2MixRt(i,k) = MolFr2MixRt(MolWt, xz_MolFr(i,k)) 
      end do
    end do
    
  end function xz_MolFr2MixRt


!!!-----------------------------------------------------------------------!!!
  function MixRt2Vap(MolWt, MixRtV, ExnerAll)
    !
    !比湿から蒸気圧を求める. 
    !
    
    !暗黙の型宣言禁止
    implicit none 
    
    !入出力変数
    real(8), intent(in) :: MolWt            !分子量
    real(8), intent(in) :: MixRtV           !混合比(擾乱 + 基本場)
    real(8), intent(in) :: ExnerAll         !エクスナー関数(擾乱 + 基本場)
    real(8)             :: MixRt2Vap        !混合比から蒸気圧への変換
    
    MixRt2Vap = 0.0d0
    MixRt2Vap = MixRtV * Exner2Press(ExnerAll) * MolWtDry / MolWt 
    
  end function MixRt2Vap
  

  function xz_MixRt2Vap(MolWt, xz_MixRtV, xz_ExnerAll)
    !
    !比湿から蒸気圧を求める. 
    !
    
    !暗黙の型宣言禁止
    implicit none 
    
    !入出力変数
    real(8), intent(in) :: MolWt            !分子量
    real(8), intent(in) :: xz_MixRtV(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !混合比(擾乱 + 基本場)
    real(8), intent(in) :: xz_ExnerAll(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !エクスナー関数(擾乱 + 基本場)
    real(8)             :: xz_MixRt2Vap(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !混合比から蒸気圧への変換
    integer             :: i, k

    do k = DimZMin, DimZMax
      do i = DimXMin, DimXMax

        xz_MixRt2Vap(i,k) = MixRt2Vap(MolWt, xz_MixRtV(i,k), xz_ExnerAll(i,k))
      end do
    end do
        
  end function xz_MixRt2Vap
  

!!!-----------------------------------------------------------------------!!!
  function MixRt2MolFr(MolWt, MixRtV)
    !
    !比湿から蒸気圧を求める. 
    !
    
    !暗黙の型宣言禁止
    implicit none 
    
    !入出力変数
    real(8), intent(in) :: MolWt            !分子量
    real(8), intent(in) :: MixRtV           !混合比(擾乱 + 基本場)
    real(8)             :: MixRt2MolFr      !混合比から蒸気圧への変換

    MixRt2MolFr = 0.0d0
    MixRt2MolFr = MixRtV * MolWtDry / MolWt
    
  end function MixRt2MolFr


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

    !暗黙の型宣言禁止
    implicit none 
    
    !入出力変数
    integer, intent(in) :: SpcID
    real(8), intent(in) :: MolWt
    real(8), intent(in) :: TempAll          !温度(擾乱 + 基本場)
    real(8), intent(in) :: ExnerAll         !エクスナー関数(擾乱 + 基本場)
    real(8)             :: DMixRtSatDSvapPress 
    real(8)             :: DMixRtSatDPotTemp   
                           
    DMixRtSatDSvapPress = 0.0d0
    DMixRtSatDPotTemp   = 0.0d0
    
    DMixRtSatDSvapPress =                                   &
      & Exner2Press(ExnerAll) * (MolWt / MolWtDry)          & 
      & / ((Exner2Press(ExnerAll)                           &
      &       - SvapPress(SpcID, TempAll)) ** 2.0d0)
    
    DMixRtSatDPotTemp =                         &
      &   DMixRtSatDSvapPress                   &
      & * DSvapPressDTemp(SpcID, TempAll)       &
      & * ExnerAll
    
  end function DMixRtSatDPotTemp


  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_DMixRtSatDPotTemp(DimXMin:DimXMax, DimZMin:DimZMax)
                           
    integer             :: i, k

    do k = DimZMin, DimZMax
      do i = DimXMin, DimXMax
        xz_DMixRtSatDPotTemp(i,k) = &
          & DMixRtSatDPotTemp(SpcID, MolWt, xz_TempAll(i,k), xz_ExnerAll(i,k))
      end do
    end do
    
  end function xz_DMixRtSatDPotTemp


!!!-----------------------------------------------------------------------!!!
  function Exner2Press(ExnerAll)
    !
    ! エクスナー関数を圧力に換算
    !

    !暗黙の型宣言禁止
    implicit none 
    
    !入出力変数
    real(8), intent(in) :: ExnerAll         !エクスナー関数(擾乱 + 基本場)
    real(8)             :: Exner2Press
    
    Exner2Press = 0.0d0
    Exner2Press = PressSfc * (ExnerAll ** (CpDry / GasRDry))
    
  end function Exner2Press
  

  function xz_Exner2Press(xz_ExnerAll)
    !
    ! エクスナー関数を圧力に換算
    !

    !暗黙の型宣言禁止
    implicit none 
    
    !入出力変数
    real(8), intent(in) :: xz_ExnerAll(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !エクスナー関数(擾乱 + 基本場)
    real(8)             :: xz_Exner2Press(DimXMin:DimXMax, DimZMin:DimZMax)
    integer             :: i, k

    do k = DimZMin, DimZMax
      do i = DimXMin, DimXMax
        xz_Exner2Press(i,k) = PressSfc * (xz_ExnerAll(i,k) ** (CpDry / GasRDry))
      end do
    end do
  end function xz_Exner2Press
  

!!!-----------------------------------------------------------------------!!!
  function DelMixRtNH4SH(TempAll, PressAll, MixRtNH3, MixRtH2S, MixRtNH4SH, &
    &                    MolWtNH3, MolWtH2S, MolWtNH4SH)
    !
    ! NH4SH 生成反応に伴う, NH4SH の生成量(混合比)を求める
    !
    
    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(8), intent(in) :: TempAll       !温度
    real(8), intent(in) :: PressAll      !圧力
    real(8), intent(in) :: MixRtNH3      !NH3 の混合比
    real(8), intent(in) :: MixRtH2S      !H2S の混合比
    real(8), intent(in) :: MixRtNH4SH    !NH4SH の混合比
    real(8), intent(in) :: MolWtNH3      !NH3 の分子量
    real(8), intent(in) :: MolWtH2S      !H2S の分子量
    real(8), intent(in) :: MolWtNH4SH    !NH4SH の分子量
    real(8)             :: DelMixRtNH4SH !NH4SH の混合比
    real(8)             :: EquivConst    !圧平衡定数
    real(8)             :: DelPressNH4SH !NH4SH を蒸発させた時の圧力変化
    real(8)             :: PPress(2)     !作業配列(分圧)
    real(8)             :: Solution(2)   !作業配列(方程式の解)
    
    !アンモニアと硫化水素の分圧
    PPress(1) = MixRt2MolFr( MolWtNH3, MixRtNH3 ) * PressAll
    PPress(2) = MixRt2MolFr( MolWtH2S, MixRtH2S ) * PressAll

    !NH4SH が全て蒸発した場合に生じる圧力変化
    DelPressNH4SH = MixRt2MolFr( MolWtNH4SH, MixRtNH4SH ) * PressAll

    !圧平衡定数
    EquivConst = 61.781d0 - 10834.0d0 / TempAll - dlog(1.0d2)
    
    !気圧変化を二次方程式の解として求める. 
    Solution(1) = 5.0d-1 * (sum(PPress) &
      &         + dsqrt( (PPress(1) - PPress(2))**2.0d0 &
      &                    + 4.0d0 * dexp(EquivConst)) )
    Solution(2) = 5.0d-1 * (sum(PPress) &
      &         - dsqrt( (PPress(1) - PPress(2))**2.0d0 &
      &                    + 4.0d0 * dexp(EquivConst)) )

    !成立条件のチェック
    !  平衡の場合, 圧力変化は以下の条件を満足する必要がある
    !    * NH4SH が全て蒸発した場合に生じる圧力変化よりも大きい
    !    * H2S と NH3 の分圧よりも小さい
    !  満足出来ない場合には, 圧力変化はゼロとなる. 
    if ( minval(Solution, 1, Solution > - DelPressNH4SH) > minval(PPress, 1) ) then 
      DelMixRtNH4SH = 0.0d0
    else      
      DelMixRtNH4SH = MolFr2MixRt( MolWtNH4SH, minval(Solution, 1, Solution > - DelPressNH4SH) / PressAll )
    end if
   
!    write(*,*) 'DelMixRtNH4SH',  EquivConst
!    write(*,*) 'DelMixRtNH4SH',  DelPressNH4SH
!    write(*,*) 'DelMixRtNH4SH',  PPress
!    write(*,*) 'DelMixRtNH4SH',  Solution
!    write(*,*) 'DelMixRtNH4SH',  DelMixRtNH4SH

  end function DelMixRtNH4SH


!!!-----------------------------------------------------------------------!!!
  function xz_DelMixRtNH4SH(xz_TempAll, xz_PressAll, xz_MixRtNH3, xz_MixRtH2S, xz_MixRtNH4SH, &
    &                    MolWtNH3, MolWtH2S, MolWtNH4SH)
    !
    ! 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) :: xz_MixRtNH4SH(DimXMin:DimXMax, DimZMin:DimZMax)
                                         !NH4SH の混合比
    real(8), intent(in) :: MolWtNH3      !NH3 の分子量
    real(8), intent(in) :: MolWtH2S      !H2S の分子量
    real(8), intent(in) :: MolWtNH4SH    !NH4SH の分子量
    real(8)             :: xz_DelMixRtNH4SH(DimXMin:DimXMax, DimZMin:DimZMax)
                                         !NH4SH の混合比
    real(8)             :: EquivConst    !圧平衡定数
    real(8)             :: DelPressNH4SH !NH4SH を蒸発させた時の圧力変化
    real(8)             :: PPress(2)     !作業配列(分圧)
    real(8)             :: Solution(2)   !作業配列(方程式の解)
    integer             :: i, k

    do k = DimZMin, DimZMax
      do i = DimXMin, DimXMax

        !アンモニアと硫化水素の分圧
        PPress(1) = MixRt2MolFr( MolWtNH3, xz_MixRtNH3(i,k) ) * xz_PressAll(i,k)
        PPress(2) = MixRt2MolFr( MolWtH2S, xz_MixRtH2S(i,k) ) * xz_PressAll(i,k)

        !NH4SH が全て蒸発した場合に生じる圧力変化
        DelPressNH4SH = MixRt2MolFr( MolWtNH4SH, xz_MixRtNH4SH(i,k) ) * xz_PressAll(i,k)
        
        !圧平衡定数
        EquivConst = 61.781d0 - 10834.0d0 / xz_TempAll(i,k) - dlog(1.0d2)
        
        !気圧変化を二次方程式の解として求める. 
        Solution(1) = 5.0d-1 * (sum(PPress) &
          &         + dsqrt( (PPress(1) - PPress(2))**2.0d0 &
          &                    + 4.0d0 * dexp(EquivConst)) )
        Solution(2) = 5.0d-1 * (sum(PPress) &
          &         - dsqrt( (PPress(1) - PPress(2))**2.0d0 &
          &                    + 4.0d0 * dexp(EquivConst)) )
        
        !成立条件のチェック
        !  平衡の場合, 圧力変化は以下の条件を満足する必要がある
        !    * NH4SH が全て蒸発した場合に生じる圧力変化よりも大きい
        !    * H2S と NH3 の分圧よりも小さい
        !  満足出来ない場合には, 圧力変化はゼロとなる.
        if ( minval(Solution, 1, Solution > - DelPressNH4SH) > minval(PPress, 1) ) then 
          xz_DelMixRtNH4SH(i,k) = 0.0d0
        else      
          xz_DelMixRtNH4SH(i,k) = MolFr2MixRt( MolWtNH4SH, minval(Solution, 1, Solution > - DelPressNH4SH) / xz_PressAll(i,k) )
        end if
        
      end do
    end do
  end function xz_DelMixRtNH4SH


!!!-----------------------------------------------------------------------!!!
  function DelMolFrNH4SH(TempAll, PressAll, MolFrNH3, MolFrH2S)
    !
    ! 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)             :: DelMolFrNH4SH !NH4SH 生成に伴うモル比変化
    real(8)             :: EquivConst    !圧平衡定数
    real(8)             :: PPress(2)     !作業配列(分圧)
    real(8)             :: Solution(2)   !作業配列(方程式の解)

    !------------------------------------------------------------
    !NH4SH の平衡条件
    !------------------------------------------------------------
    !アンモニアと硫化水素の分圧
    PPress(1) = MolFrNH3 * PressAll
    PPress(2) = MolFrH2S * PressAll

    !圧平衡定数
    EquivConst = 61.781d0 - 10834.0d0 / TempAll - dlog(1.0d2)
    
    !気圧変化を二次方程式の解として求める. 
    Solution(1) = 5.0d-1 * (sum(PPress) &
      &         + dsqrt( (PPress(1) - PPress(2))**2.0d0 &
      &                    + 4.0d0 * dexp(EquivConst)) )
    Solution(2) = 5.0d-1 * (sum(PPress) &
      &         - dsqrt( (PPress(1) - PPress(2))**2.0d0 &
      &                    + 4.0d0 * dexp(EquivConst)) )
    
    !成立条件のチェック
    if ( minval(Solution, 1, Solution > 0.0d0) > minval(PPress, 1) ) then 
      DelMolFrNH4SH = 0.0d0
    else      
      DelMolFrNH4SH = minval(Solution, 1, Solution > 0.0d0) / PressAll
    end if

  end function DelMolFrNH4SH

end module MoistFunc