!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2005, 2006. All rights reserved.
!---------------------------------------------------------------------
!= Module MoistFunc
!
!   * Developer: SUGIYAMA Ko-ichiro
!   * Version: $Id: moistfunc.f90,v 1.1.1.1 2006/04/25 03:43:58 deepconv Exp $
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview
!
!$BJXMx$J7W;;%D!<%k72(B
!   * $B>x5$05$+$i:.9gHf$XJQ49(B
!   * $B%b%kHf$+$i:.9gHf$XJQ49(B
!   * $B:.9gHf$+$i>x5$05$XJQ49(B
!   * $B:.9gHf$+$i%b%kHf$XJQ49(B
! 
!== Error Handling
!
!== Bugs
!
!== Note
!
!== Future Plans
!
!

module MoistFunc
  !
  !$BC10L$N49;;$r9T$&$?$a$N%Q%C%1!<%87?%b%8%e!<%k(B
  !   * $B>x5$05$+$i:.9gHf$XJQ49(B
  !   * $B%b%kHf$+$i:.9gHf$XJQ49(B
  !   * $B:.9gHf$+$i>x5$05$XJQ49(B
  !   * $B:.9gHf$+$i%b%kHf$XJQ49(B
  ! 

  !$B%b%8%e!<%kFI$_9~$_(B
  use gridset, only:  DimXMin,            & ! x $BJ}8~$NG[Ns$N2<8B(B
    &                 DimXMax,            & ! x $BJ}8~$NG[Ns$N>e8B(B
    &                 DimZMin,            & ! z $BJ}8~$NG[Ns$N2<8B(B
    &                 DimZMax               ! z $BJ}8~$NG[Ns$N>e8B(B
  use basicset, only: PressSfc,           & !$BI8=`05NO(B               [Pa]
    &                 CpDry,              & !$B4%Ag@.J,$NJ?6QDj05HfG.(B [J/K kg]
    &                 MolWtDry,           & !$B4%Ag@.J,$NJ?6QJ,;RNL(B   [kg/mol]
    &                 GasRDry               !$B4%Ag@.J,$N5$BNDj?t(B     [J/K kg]
  use ChemCalc, only: SvapPress,           &!$BK0OB>x5$05(B
    &                 xz_SvapPress,        &!$BK0OB>x5$05(B
    &                 DSvapPressDTemp,     &!$BK0OB>x5$05$N29EYHyJ,(B
    &                 xz_DSvapPressDTemp,  &!$BK0OB>x5$05$N29EYHyJ,(B
    &                 ReactHeatNH4SHPerMass

  !$B0EL[$N7?@k8@6X;_(B
  implicit none 
  
  !$BB0@-$N;XDj(B
  private

  !$B4X?t$r(B public $B$K;XDj(B
!  public Vap2MixRt, xz_Vap2MixRt 
!  public MixRt2Vap, xza_MixRt2Vap
!  public MolFr2MixRt, xz_MolFr2MixRt
!  public MixRt2MolFr
!  public Exner2Press, xz_Exner2Press
  public DMixRtSatDPotTemp, xz_DMixRtSatDPotTemp
  public DelMixRtNH4SH, DelMolFrNH4SH
  public xz_DelMixRtNH4SH
contains

!!!
!!!$B4pK\4X?t(B
!!!

!!!-----------------------------------------------------------------------!!!
  function Vap2MixRt(MolWt, VapAll, ExnerAll)
    !
    !$B>x5$05$+$iHf<>$r5a$a$k(B. 
    !$B6E=L@.J,$O>/$J$$$?$a(B, $B:.9gHf$OHf<>$K$[$\Ey$7$$(B
    !
    
    !$B0EL[$N7?@k8@6X;_(B
    implicit none 
    
    !$BF~=PNOJQ?t(B
    real(8), intent(in) :: MolWt            !$B2=3X<o$NJ,;RNL(B
    real(8), intent(in) :: VapAll           !$B>x5$05(B($B>qMp(B + $B4pK\>l(B)
    real(8), intent(in) :: ExnerAll         !$B%(%/%9%J!<4X?t(B($B>qMp(B + $B4pK\>l(B)
    real(8)             :: Vap2MixRt        !$B>x5$05$+$iHf<>$X$NJQ49(B
    
    Vap2MixRt = 0.0d0
    Vap2MixRt = VapAll * MolWt / (MolWtDry * Exner2Press(ExnerAll))
    
  end function Vap2MixRt

  
  function xz_Vap2MixRt(MolWt, xz_VapAll, xz_ExnerAll)
    !
    !$B>x5$05$+$iHf<>$r5a$a$k(B. 
    !$B6E=L@.J,$O>/$J$$$?$a(B, $B:.9gHf$OHf<>$K$[$\Ey$7$$(B
    !
    
    !$B0EL[$N7?@k8@6X;_(B
    implicit none 
    
    !$BF~=PNOJQ?t(B
    real(8), intent(in) :: MolWt   !$B2=3X<o$NJ,;RNL(B
    real(8), intent(in) :: xz_VapAll(DimXMin:DimXMax, DimZMin:DimZMax) 
                                   !$B>x5$05(B($B>qMp@.J,$H4pK\>l$NOB(B)   [Pa]
    real(8), intent(in) :: xz_ExnerAll(DimXMin:DimXMax, DimZMin:DimZMax) 
                                   !$B%(%/%9%J!<4X?t(B($B>qMp@.J,$H4pK\>l$NOB(B)
    real(8)             :: xz_Vap2MixRt(DimXMin:DimXMax, DimZMin:DimZMax)
                                   !$B>x5$05$+$i:.9gHf$X$NJQ49(B
    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)
    !
    !$B%b%kHf$+$iHf<>$r5a$a$k(B. 
    !
    
    !$B0EL[$N7?@k8@6X;_(B
    implicit none 
    
    !$BF~=PNOJQ?t(B
    real(8), intent(in) :: MolWt            !$BJ,;RNL(B
    real(8), intent(in) :: MolFr            !$B%b%kHf(B($B>qMp(B + $B4pK\>l(B)
    real(8)             :: MolFr2MixRt      !$B%b%kHf$+$i:.9gHf$X$NJQ49(B

    MolFr2MixRt = 0.0d0
    MolFr2MixRt = MolFr * MolWt / MolWtDry     

  end function MolFr2MixRt


  function xz_MolFr2MixRt(MolWt, xz_MolFr)
    !
    !$B%b%kHf$+$iHf<>$r5a$a$k(B. 
    !
    
    !$B0EL[$N7?@k8@6X;_(B
    implicit none 
    
    !$BF~=PNOJQ?t(B
    real(8), intent(in) :: MolWt            !$BJ,;RNL(B
    real(8), intent(in) :: xz_MolFr(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !$B%b%kHf(B($B>qMp(B + $B4pK\>l(B)
    real(8)             :: xz_MolFr2MixRt(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !$B%b%kHf$+$i:.9gHf$X$NJQ49(B
    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)
    !
    !$BHf<>$+$i>x5$05$r5a$a$k(B. 
    !
    
    !$B0EL[$N7?@k8@6X;_(B
    implicit none 
    
    !$BF~=PNOJQ?t(B
    real(8), intent(in) :: MolWt            !$BJ,;RNL(B
    real(8), intent(in) :: MixRtV           !$B:.9gHf(B($B>qMp(B + $B4pK\>l(B)
    real(8), intent(in) :: ExnerAll         !$B%(%/%9%J!<4X?t(B($B>qMp(B + $B4pK\>l(B)
    real(8)             :: MixRt2Vap        !$B:.9gHf$+$i>x5$05$X$NJQ49(B
    
    MixRt2Vap = 0.0d0
    MixRt2Vap = MixRtV * Exner2Press(ExnerAll) * MolWtDry / MolWt 
    
  end function MixRt2Vap
  

  function xza_MixRt2Vap(MolWt, xza_MixRtV, xz_ExnerAll)
    !
    !$BHf<>$+$i>x5$05$r5a$a$k(B. 
    !
    
    !$B0EL[$N7?@k8@6X;_(B
    implicit none 
    
    !$BF~=PNOJQ?t(B
    real(8), intent(in) :: MolWt            !$BJ,;RNL(B
    real(8), intent(in) :: xza_MixRtV(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !$B:.9gHf(B($B>qMp(B + $B4pK\>l(B)
    real(8), intent(in) :: xz_ExnerAll(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !$B%(%/%9%J!<4X?t(B($B>qMp(B + $B4pK\>l(B)
    real(8)             :: xza_MixRt2Vap(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !$B:.9gHf$+$i>x5$05$X$NJQ49(B
    integer             :: i, k

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

        xza_MixRt2Vap(i,k) = MixRt2Vap(MolWt, xza_MixRtV(i,k), xz_ExnerAll(i,k))
      end do
    end do
        
  end function xza_MixRt2Vap
  

!!!-----------------------------------------------------------------------!!!
  function MixRt2MolFr(MolWt, MixRtV)
    !
    !$BHf<>$+$i>x5$05$r5a$a$k(B. 
    !
    
    !$B0EL[$N7?@k8@6X;_(B
    implicit none 
    
    !$BF~=PNOJQ?t(B
    real(8), intent(in) :: MolWt            !$BJ,;RNL(B
    real(8), intent(in) :: MixRtV           !$B:.9gHf(B($B>qMp(B + $B4pK\>l(B)
    real(8)             :: MixRt2MolFr      !$B:.9gHf$+$i>x5$05$X$NJQ49(B

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


!!!-----------------------------------------------------------------------!!!
  function DMixRtSatDPotTemp(SpcID, MolWt, TempAll, ExnerAll)
    !
    !$BK0OB>x5$05$N(B $B&H(B $BHyJ,$r9T$&(B
    !$B<B:]$K$O(B, dq/dp * dp/dT * dT/d$B&H(B $B$r<B9T(B. ($BC"$7(B p $B$OK0OB>x5$05(B)
    !

    !$B0EL[$N7?@k8@6X;_(B
    implicit none 
    
    !$BF~=PNOJQ?t(B
    integer, intent(in) :: SpcID
    real(8), intent(in) :: MolWt
    real(8), intent(in) :: TempAll          !$B29EY(B($B>qMp(B + $B4pK\>l(B)
    real(8), intent(in) :: ExnerAll         !$B%(%/%9%J!<4X?t(B($B>qMp(B + $B4pK\>l(B)
    real(8)             :: PressAll         !$B05NO(B($B>qMp(B + $B4pK\>l(B)
    real(8)             :: DMixRtSatDSvapPress 
    real(8)             :: DMixRtSatDPotTemp   
                           
    DMixRtSatDSvapPress = 0.0d0
    DMixRtSatDPotTemp   = 0.0d0
    PressAll = PressSfc * (ExnerAll ** (CpDry / GasRDry))
    
    DMixRtSatDSvapPress =                                    &
      & PressAll * ( MolWt / MolWtDry)                       &  
      & / ( (PressAll - SvapPress(SpcID, TempAll)) ** 2.0d0 )
    
    DMixRtSatDPotTemp =                         &
      &   DMixRtSatDSvapPress                   &
      & * DSvapPressDTemp(SpcID, TempAll)       &
      & * ExnerAll
    
  end function DMixRtSatDPotTemp


  function xz_DMixRtSatDPotTemp(SpcID, MolWt, xz_TempAll, xz_ExnerAll)
    !
    !$BK0OB>x5$05$N(B $B&H(B $BHyJ,$r9T$&(B
    !$B<B:]$K$O(B, dq/dp * dp/dT * dT/d$B&H(B $B$r<B9T(B. ($BC"$7(B p $B$OK0OB>x5$05(B)
    !

    !$B0EL[$N7?@k8@6X;_(B
    implicit none 
    
    !$BF~=PNOJQ?t(B
    integer, intent(in) :: SpcID
    real(8), intent(in) :: MolWt
    real(8), intent(in) :: xz_TempAll(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !$B29EY(B($B>qMp(B + $B4pK\>l(B)
    real(8), intent(in) :: xz_ExnerAll(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !$B%(%/%9%J!<4X?t(B($B>qMp(B + $B4pK\>l(B)
    real(8)             :: xz_PressAll(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !$B05NO(B($B>qMp(B + $B4pK\>l(B)
    real(8)             :: xz_DMixRtSatDPotTemp(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)             :: xz_DMixRtSatDSvapPress(DimXMin:DimXMax, DimZMin:DimZMax)
                           
    xz_DMixRtSatDSvapPress = 0.0d0
    xz_DMixRtSatDPotTemp   = 0.0d0
    xz_PressAll = PressSfc * (xz_ExnerAll ** (CpDry / GasRDry))
    
    xz_DMixRtSatDSvapPress =                                   &
      & xz_PressAll * ( MolWt / MolWtDry)                      &  
      & / ( (xz_PressAll - xz_SvapPress(SpcID, xz_TempAll)) ** 2.0d0 )
    
    xz_DMixRtSatDPotTemp =                            &
      &   xz_DMixRtSatDSvapPress                      &
      & * xz_DSvapPressDTemp(SpcID, xz_TempAll)       &
      & * xz_ExnerAll
    
  end function xz_DMixRtSatDPotTemp


!!!-----------------------------------------------------------------------!!!
  function Exner2Press(ExnerAll)
    !
    ! $B%(%/%9%J!<4X?t$r05NO$K49;;(B
    !

    !$B0EL[$N7?@k8@6X;_(B
    implicit none 
    
    !$BF~=PNOJQ?t(B
    real(8), intent(in) :: ExnerAll         !$B%(%/%9%J!<4X?t(B($B>qMp(B + $B4pK\>l(B)
    real(8)             :: Exner2Press
    
    Exner2Press = 0.0d0
    Exner2Press = PressSfc * (ExnerAll ** (CpDry / GasRDry))
    
  end function Exner2Press
  

  function xz_Exner2Press(xz_ExnerAll)
    !
    ! $B%(%/%9%J!<4X?t$r05NO$K49;;(B
    !

    !$B0EL[$N7?@k8@6X;_(B
    implicit none 
    
    !$BF~=PNOJQ?t(B
    real(8), intent(in) :: xz_ExnerAll(DimXMin:DimXMax, DimZMin:DimZMax)
                                            !$B%(%/%9%J!<4X?t(B($B>qMp(B + $B4pK\>l(B)
    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 $B@8@.H?1~$KH<$&(B, NH4SH $B$N@8@.NL(B($B:.9gHf(B)$B$r5a$a$k(B
    !
    
    !$B0EL[$N7?@k8@6X;_(B
    implicit none

    !$BJQ?tDj5A(B
    real(8), intent(in) :: TempAll       !$B29EY(B
    real(8), intent(in) :: PressAll      !$B05NO(B
    real(8), intent(in) :: MixRtNH3      !NH3 $B$N:.9gHf(B
    real(8), intent(in) :: MixRtH2S      !H2S $B$N:.9gHf(B
    real(8), intent(in) :: MixRtNH4SH    !NH4SH $B$N:.9gHf(B
    real(8), intent(in) :: MolWtNH3      !NH3 $B$NJ,;RNL(B
    real(8), intent(in) :: MolWtH2S      !H2S $B$NJ,;RNL(B
    real(8), intent(in) :: MolWtNH4SH    !NH4SH $B$NJ,;RNL(B
    real(8)             :: DelMixRtNH4SH !NH4SH $B$N:.9gHf(B
    real(8)             :: EquivConst    !$B05J?9UDj?t(B
    real(8)             :: DelPressNH4SH !NH4SH $B$r>xH/$5$;$?;~$N05NOJQ2=(B
    real(8)             :: PPress(2)     !$B:n6HG[Ns(B($BJ,05(B)
!    real(8)             :: Solution      !$B:n6HG[Ns(B($BJ}Dx<0$N2r(B)
    real(8)             :: Solution(2)   !$B:n6HG[Ns(B($BJ}Dx<0$N2r(B)
    
    !$B%"%s%b%K%"$HN22=?eAG$NJ,05(B
    PPress(1) = MixRtNH3 * PressAll * MolWtDry / MolWtNH3   
    PPress(2) = MixRtH2S * PressAll * MolWtDry / MolWtH2S

    !NH4SH $B$,A4$F>xH/$7$?>l9g$K@8$8$k05NOJQ2=(B
    DelPressNH4SH = MixRtNH4SH * PressAll * MolWtDry / MolWtNH4SH

    !$B05J?9UDj?t(B
    EquivConst = 61.781d0 - 10834.0d0 / TempAll - dlog(1.0d2)
    
    !$B5$05JQ2=$rFs<!J}Dx<0$N2r$H$7$F5a$a$k(B. 
    Solution(1) =                                   &
      & (                                           &
      &     sum(PPress)                             &
      &   + dsqrt( (PPress(1) - PPress(2)) ** 2.0d0 &
      &   + 4.0d0 * dexp(EquivConst))               &
      &  ) * 5.0d-1
    Solution(2) =                                   &
      & (                                           &
      &     sum(PPress)                             &
      &   - dsqrt( (PPress(1) - PPress(2)) ** 2.0d0 &
      &   + 4.0d0 * dexp(EquivConst))               &
      &  ) * 5.0d-1
    !$B5$05JQ2=$r5a$a$k(B. 
    !  $BK\Mh$O(B DelX^2 - (P_NH3 + P_H2S) * DelX + P_NH3 * P_H2S * exp( Kp ) = 0
    !  $B$H$$$&Fs<!J}Dx<0$r5a$a$kI,MW$,$"$k$,(B, DelX^2 $B$O==J,>.$5$$$H$_$J$7(B, 
    !  (P_NH3 + P_H2S) * DelX - P_NH3 * P_H2S * exp( Kp ) = 0 $B$r2r$/(B. 
!    Solution =                                               &
!      & max(                                                 &
!      &        ( PPress(1) * PPress(2) - dexp(EquivConst) )  &
!      &         / sum(PPress),                               &
!      &        1.0d-40                                       &
!      &     ) 
    
    !$B@.N)>r7o$N%A%'%C%/(B
    !  $BJ?9U$N>l9g(B, $B05NOJQ2=$O0J2<$N>r7o$rK~B-$9$kI,MW$,$"$k(B
    !    * NH4SH $B$,A4$F>xH/$7$?>l9g$K@8$8$k05NOJQ2=$h$j$bBg$-$$(B
    !    * H2S $B$H(B NH3 $B$NJ,05$h$j$b>.$5$$(B
    !  $BK~B-=PMh$J$$>l9g$K$O(B, $B05NOJQ2=$O%<%m$H$J$k(B. 
    if ( minval(Solution, 1, Solution > - DelPressNH4SH) > minval(PPress, 1) ) then
      DelMixRtNH4SH = 0.0d0
    else      
      DelMixRtNH4SH = minval(Solution, 1, Solution > - DelPressNH4SH)&
        &              * MolWtNH4SH / ( PressAll * MolWtDry )
    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 $B@8@.H?1~$KH<$&(B, NH4SH $B$N@8@.NL(B($B:.9gHf(B)$B$r5a$a$k(B
    !
    
    !$B0EL[$N7?@k8@6X;_(B
    implicit none

    !$BJQ?tDj5A(B
    real(8), intent(in) :: xz_TempAll(DimXMin:DimXMax, DimZMin:DimZMax)
                                         !$B29EY(B
    real(8), intent(in) :: xz_PressAll(DimXMin:DimXMax, DimZMin:DimZMax)
                                         !$B05NO(B
    real(8), intent(in) :: xz_MixRtNH3(DimXMin:DimXMax, DimZMin:DimZMax)
                                         !NH3 $B$N:.9gHf(B
    real(8), intent(in) :: xz_MixRtH2S(DimXMin:DimXMax, DimZMin:DimZMax)
                                         !H2S $B$N:.9gHf(B
    real(8), intent(in) :: xz_MixRtNH4SH(DimXMin:DimXMax, DimZMin:DimZMax)
                                         !NH4SH $B$N:.9gHf(B
    real(8), intent(in) :: MolWtNH3      !NH3 $B$NJ,;RNL(B
    real(8), intent(in) :: MolWtH2S      !H2S $B$NJ,;RNL(B
    real(8), intent(in) :: MolWtNH4SH    !NH4SH $B$NJ,;RNL(B
    real(8)             :: xz_DelMixRtNH4SH(DimXMin:DimXMax, DimZMin:DimZMax)
                                         !NH4SH $B$N:.9gHf(B
    real(8)             :: xz_EquivConst(DimXMin:DimXMax, DimZMin:DimZMax)
                                         !$B05J?9UDj?t(B
    real(8)             :: xz_DelPressNH4SH(DimXMin:DimXMax, DimZMin:DimZMax)
                                         !NH4SH $B$r>xH/$5$;$?;~$N05NOJQ2=(B
    real(8)             :: xza_PPress(DimXMin:DimXMax, DimZMin:DimZMax,2)
                                         !$B:n6HG[Ns(B($BJ,05(B)
    real(8)             :: xz_PPress(DimXMin:DimXMax, DimZMin:DimZMax)
                                         !$B:n6HG[Ns(B($BJ,05(B)
    real(8)             :: xz_Zero(DimXMin:DimXMax, DimZMin:DimZMax)
                                         !$B:n6HG[Ns(B($BCM$,%<%m$N9TNs(B)
!    real(8)             :: xz_Solution(DimXMin:DimXMax, DimZMin:DimZMax)
!                                         !$B:n6HG[Ns(B($BJ}Dx<0$N2r(B)
    real(8)             :: xza_Solution(DimXMin:DimXMax, DimZMin:DimZMax,3)
                                         !$B:n6HG[Ns(B($BJ}Dx<0$N2r(B)
    logical             :: xz_SolutionSW(DimXMin:DimXMax, DimZMin:DimZMax)
                                         !$B:n6HG[Ns(B($BJ}Dx<0$N2r(B)
    integer             :: i,k

    !$B=i4|2=(B
    xz_DelMixRtNH4SH = 0.0d0
    xz_Zero          = 0.0d0

    !$B8=;~E@$G$N%"%s%b%K%"$HN22=?eAG$NJ,05(B
    xza_PPress(:,:,1) = xz_MixRtNH3 * xz_PressAll * MolWtDry / MolWtNH3   
    xza_PPress(:,:,2) = xz_MixRtH2S * xz_PressAll * MolWtDry / MolWtH2S

    !NH4SH $B$,A4$F>xH/$7$?>l9g$K@8$8$k05NOJQ2=(B
    xz_DelPressNH4SH = xz_MixRtNH4SH * xz_PressAll * MolWtDry / MolWtNH4SH
    
    !$B05J?9UDj?t(B
    xz_EquivConst = 61.781d0 - 10834.0d0 / xz_TempAll - dlog(1.0d2)
        
    !$B5$05JQ2=$r5a$a$k(B. 
    xza_Solution(:,:,1) =                                           &
      & (                                                           &
      &     sum(xza_PPress, 3)                                      &
      &   + dsqrt( (xza_PPress(:,:,1) - xza_PPress(:,:,2)) ** 2.0d0 &
      &   + 4.0d0 * dexp(xz_EquivConst))                            &
      &  ) * 5.0d-1
    xza_Solution(:,:,2) =                                           &
      & (                                                           &
      &     sum(xza_PPress, 3)                                      &
      &   - dsqrt( (xza_PPress(:,:,1) - xza_PPress(:,:,2)) ** 2.0d0 &
      &   + 4.0d0 * dexp(xz_EquivConst))                            &
      &  ) * 5.0d-1
    
    !  $BK\Mh$O(B DelX^2 - (P_NH3 + P_H2S) * DelX + P_NH3 * P_H2S * exp( Kp ) = 0
    !  $B$H$$$&Fs<!J}Dx<0$r5a$a$kI,MW$,$"$k$,(B, DelX^2 $B$O==J,>.$5$$$H$_$J$7(B, 
    !  (P_NH3 + P_H2S) * DelX - P_NH3 * P_H2S * exp( Kp ) = 0 $B$r2r$/(B. 
!    xz_Solution =                                                        &
!      & ( xza_PPress(:,:,1) * xza_PPress(:,:,2) - dexp(xz_EquivConst) )  &
!      &  / sum(xza_PPress, 3) 

    !$B@.N)>r7o$N%A%'%C%/(B
    !  $BJ?9U$N>l9g(B, $B05NOJQ2=$O0J2<$N>r7o$rK~B-$9$kI,MW$,$"$k(B
    !    * NH4SH $B$,A4$F>xH/$7$?>l9g$K@8$8$k05NOJQ2=$h$j$bBg$-$$(B
    !    * H2S $B$H(B NH3 $B$NJ,05$h$j$b>.$5$$(B
    !  $BK~B-=PMh$J$$>l9g$K$O(B, $B05NOJQ2=$O%<%m$H$J$k(B.
    ! 
    do k = DimZmin, DimZMax
      do i = DimXMin, DimXMax
        
        if ( minval(xza_Solution(i,k,:), 1, xza_Solution(i,k,:) > - xz_DelPressNH4SH(i,k)) > minval(xza_PPress(i,k,:), 1) ) then
          xz_DelMixRtNH4SH(i,k) = 0.0d0
        else      
          xz_DelMixRtNH4SH(i,k) = minval(xza_Solution(i,k,:), 1, xza_Solution(i,k,:) > - xz_DelPressNH4SH(i,k))&
            &              * MolWtNH4SH / ( xz_PressAll(i,k) * MolWtDry )
        end if

      end do
    end do

!    xz_PPress = merge( xza_PPress(:,:,1),  xza_PPress(:,:,2),   &
!      &                xza_PPress(:,:,1) < xza_PPress(:,:,2) ) 
!
!    xz_SolutionSW =                                      &
!      &       xza_Solution(:,:,1) > - xz_DelPressNH4SH   &
!      & .AND. xza_Solution(:,:,1) < xz_PPress  
!
!    xza_Solution(:,:,3) =                                &
!      & merge( xza_Solution(:,:,1), xz_Zero, xz_SolutionSW )  
!
!    xz_SolutionSW =                                      &
!      &       xza_Solution(:,:,2) > - xz_DelPressNH4SH   &
!      & .AND. xza_Solution(:,:,2) < xz_PPress  
!
!    xz_DelMixRtNH4SH =                                                   &
!      & merge( xza_Solution(:,:,2), xza_Solution(:,:,3), xz_SolutionSW ) &  
!      & * MolWtNH4SH / ( MolWtDry * xz_PressAll )

  end function xz_DelMixRtNH4SH


!!!-----------------------------------------------------------------------!!!
  function DelMolFrNH4SH(TempAll, PressAll, MolFrNH3, MolFrH2S)
    !
    ! NH4SH $B@8@.H?1~$KH<$&(B H2S $B$H(B NH3 $B$N%b%kHf$N8:>/J,$r5a$a$k(B
    !
    
    !$B0EL[$N7?@k8@6X;_(B
    implicit none

    !$BJQ?tDj5A(B
    real(8), intent(in) :: TempAll       !$B29EY(B
    real(8), intent(in) :: PressAll      !$B05NO(B
    real(8), intent(in) :: MolFrNH3      !NH3 $B$N%b%kHf(B
    real(8), intent(in) :: MolFrH2S      !H2S $B$N%b%kHf(B
    real(8)             :: DelMolFrNH4SH !NH4SH $B@8@.$KH<$&%b%kHfJQ2=(B
    real(8)             :: EquivConst    !$B05J?9UDj?t(B
    real(8)             :: PPress(2)     !$B:n6HG[Ns(B($BJ,05(B)
!    real(8)             :: Solution      !$B:n6HG[Ns(B($BJ}Dx<0$N2r(B)
    real(8)             :: Solution(2)   !$B:n6HG[Ns(B($BJ}Dx<0$N2r(B)

    !------------------------------------------------------------
    !NH4SH $B$NJ?9U>r7o(B
    !------------------------------------------------------------
    !$B%"%s%b%K%"$HN22=?eAG$NJ,05(B
    PPress(1) = MolFrNH3 * PressAll
    PPress(2) = MolFrH2S * PressAll

    !$B05J?9UDj?t(B
    EquivConst = 61.781d0 - 10834.0d0 / TempAll - dlog(1.0d2)
    
    !$B5$05JQ2=$rFs<!J}Dx<0$N2r$H$7$F5a$a$k(B. 
    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)) )
    !$B5$05JQ2=$r5a$a$k(B. 
    !  $BK\Mh$O(B DelX^2 - (P_NH3 + P_H2S) * DelX + P_NH3 * P_H2S * exp( Kp ) = 0
    !  $B$H$$$&Fs<!J}Dx<0$r5a$a$kI,MW$,$"$k$,(B, DelX^2 $B$O==J,>.$5$$$H$_$J$7(B, 
    !  (P_NH3 + P_H2S) * DelX - P_NH3 * P_H2S * exp( Kp ) = 0 $B$r2r$/(B. 
!    Solution =                                               &
!      & max(                                                 &
!      &        ( PPress(1) * PPress(2) - dexp(EquivConst) )  &
!      &         / sum(PPress),                               &
!      &        1.0d-40                                       &
!      &     ) 
!    write(*,*)  PPress(1), PPress(2), - dexp(EquivConst) 
    
    !$B@.N)>r7o$N%A%'%C%/(B
    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

!    if ( Solution > minval(PPress, 1) ) then 
!      DelMolFrNH4SH = 0.0d0
!    else      
!      DelMolFrNH4SH = Solution / PressAll
!    end if

  end function DelMolFrNH4SH

end module MoistFunc
