!= Module ECCM_3d
!
! Authors::   SUGIYAMA Koichiro, ODAKA Masatsugu
! Version::   $Id: eccm_3d.f90,v 1.3 2008-06-19 17:05:40 odakker Exp $
! Tag Name::  $Name: arare4-20080627 $
! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview 
!
!断熱的に上昇する気塊の温度減率を計算し, 静水圧平衡から圧力を決める
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!  * 比熱は定数, 平均分子量は配列
!  * オイラースキームだと精度が足りないので, ルンゲクッタスキームを利用. 
!
!== Future Plans
!

module ECCM_3d

  !モジュール読み込み
  use dc_types, only : DP
  use dc_trace, only : DbgMessage
  use dc_message, only: MessageNotify

  use gridset_3d,only: DimZMin,       &! 配列の Z 方向の下限
    &                 DimZMax,       &! 配列の Z 方向の上限 
    &                 RegZMin,       &!
    &                 SpcNum,        &!
    &                 z_dz, r_dz      !
  use basicset_3d,only:  MolWtDry,      &!
    &                 MolWtWet,      &!
    &                 CpDryMol,      &!
    &                 SpcWetID,      &!
    &                 TempSfc,       &!
    &                 PressSfc,      &!
    &                 Grav            !
  use chemcalc_3d, only: SvapPress,        &!
    &                 LatentHeatPerMol, &!
    &                 ReactHeatNH4SHPerMol
  use moistset, only: CondNum,          &!凝結過程の数
    &                 IdxCG,            &!凝結過程(蒸気)の配列添え字
    &                 IdxCC,            &!凝結過程(雲)の配列添え字
    &                 GasNum,           &!気体の数
    &                 IdxNH3,           &!NH3(蒸気)の配列添え字
    &                 IdxH2S             !H2S(蒸気)の配列添え字
  use ChemData, only: GasRUniv        
  use MoistFunc_3d,only: DelMolFrNH4SH 
  use StoreStab_3d,only: StoreStabTemp, StoreStabMolWt

  !暗黙の型宣言禁止
  implicit none

  !属性の指定
  private

  !関数の公開
  public ECCM_MolFr
  public ECCM_Stab
  public ECCM_Dry
  public ECCM_Wet

contains

!!!------------------------------------------------------------------------------!!!
  subroutine ECCM_Dry( a_MolFrIni, Humidity, z_Temp, z_Press, z_MolWtMean, za_MolFr )
    !
    !== 概要
    !  * 乾燥断熱温度減率に沿った温度、圧力を求め、それらに対して指定された相対湿度となるように凝結成分のモル比を決める
    !  * 比熱は乾燥気塊のもので代表させる
    !    * 流体の方程式において比熱は乾燥成分のもので代表させているため
    !  * 大気の平均分子量には湿潤成分の分子量を効果
    !    * 流体の方程式において, 湿潤成分の分子量は考慮しているため
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    real(DP), intent(in) :: a_MolFrIni(1:SpcNum)    !下部境界でのモル比
    real(DP), intent(in) :: Humidity                !相対湿度 ( Humidity <= 1.0 )
    real(DP), intent(out):: z_Temp(DimZMin:DimZMax) !温度
    real(DP), intent(out):: z_Press(DimZMin:DimZMax)!圧力
    real(DP), intent(out):: z_MolWtMean(DimZMin:DimZMax) 
                                                   !平均分子量
    real(DP), intent(out):: za_MolFr(DimZMin:DimZMax, 1:SpcNum) 
                                                   !モル分率
    real(DP)             :: MolWtMeanDry            !乾燥気塊の分子量の作業変数
    real(DP)             :: MolWtMeanWet            !湿潤気塊の分子量の作業変数
    real(DP)             :: SatPress                !飽和蒸気圧
    real(DP)             :: VapPress                !蒸気圧
    real(DP)             :: DelMolFr
    integer             :: k, s

    !-------------------------------------------------------------
    ! 配列の初期化
    !-------------------------------------------------------------
    !地表面の分子量を決める
    za_MolFr(RegZMin, 1:SpcNum)   = a_MolFrIni(1:SpcNum) 

    !地表面での平均分子量を決める
    MolWtMeanDry = MolWtDry * (1.0d0 - sum(a_MolFrIni))
    MolWtMeanWet = dot_product(MolWtWet, a_MolFrIni) 
    z_MolWtMean(RegZMin) = MolWtMeanDry + MolWtMeanWet
  
    !地表面での温度(RegZMin は, 高度 DelZ / 2 に相当)
    z_Temp          = 1.0d-60
    z_Temp(RegZMin) = TempSfc - Grav * z_MolWtMean(RegZMin) &
      &               / CpDryMol * ( z_dz(RegZMin) * 5.0d-1 )
    
    !地表面での圧力(RegZMin は, 高度 DelZ / 2 に相当)
    z_Press           = 1.0d-60
    z_Press(RegZMin)  = &
      & PressSfc *((TempSfc / z_Temp(RegZMin)) ** (- CpDryMol /  GasRUniv))
    
    !-----------------------------------------------------------
    ! (1) 乾燥断熱線に沿った温度を決める
    ! (2) 静水圧平衡から圧力を求める
    ! (3) (1),(2) の温度圧力に対して, とある相対湿度となるモル比を決める
    !-----------------------------------------------------------    
    DtDz: do k = RegZMin, DimZMax-1

      !(1)乾燥断熱線に沿って k+1 での温度を計算
      z_Temp(k+1) = z_Temp(k) - Grav * z_MolWtMean(k) / CpDryMol * r_dz(k)
      
      !念為
      if (z_Temp(k+1) <= 0.0d0 ) z_Temp(k+1) = z_Temp(k) 
      
      !(2)圧力を静水圧平衡から計算
      z_Press(k+1) =                                                  &
        &  z_Press(k) * ((z_Temp(k) / z_Temp(k+1)) ** (- CpDryMol / GasRUniv)) 

      !(3)モル比の計算
      !  まずはモル比は変化しないものとしてモル比を与える
      !  飽和蒸気圧と平衡定数との平衡条件の前に適用しておく
      za_MolFr(k+1,:) = za_MolFr(k,:)
      
      do s = 1, CondNum      
        !飽和蒸気圧
        SatPress = SvapPress( SpcWetID(IdxCC(s)), z_Temp(k+1) )        
        
        !元々のモル分率を用いて現在の蒸気圧を計算
        VapPress = za_MolFr(k,IdxCG(s)) * z_Press(k+1)
        
        !飽和蒸気圧と圧力から現在のモル比を計算
        if ( VapPress > SatPress ) then         
          za_MolFr(k+1,IdxCG(s)) = max(SatPress * Humidity / z_Press(k+1), 1.0d-16)
        end if
      end do
      
      !NH4SH の平衡条件
      if ( IdxNH3 /= 0 ) then 
        DelMolFr =                                              &
          & max (                                               &
          &    DelMolFrNH4SH(                                   &
          &         z_Temp(k+1), z_Press(k+1),                  &
          &         za_MolFr(k+1,IdxNH3), za_MolFr(k+1,IdxH2S),       &
          &         Humidity                                    &
          &      ),                                             &
          &    0.0d0                                            &
          &  )
        za_MolFr(k+1,IdxNH3) = za_MolFr(k+1,IdxNH3) - DelMolFr
        za_MolFr(k+1,IdxH2S) = za_MolFr(k+1,IdxH2S) - DelMolFr
      end if
      
      !------------------------------------------------------------
      !温度勾配を計算
      !------------------------------------------------------------
      !地表面での平均分子量を決める
      MolWtMeanDry = MolWtDry * (1.0d0 - sum(za_MolFr(k+1, 1:SpcNum)))
      MolWtMeanWet = dot_product(MolWtWet(1:SpcNum), za_MolFr(k+1, 1:SpcNum))
      z_MolWtMean(k+1) = MolWtMeanDry + MolWtMeanWet

    end do DtDz
    
  end subroutine ECCM_Dry


!!------------------------------------------------------------------------------!!!
  subroutine ECCM_Wet( a_MolFrIni, Humidity, z_Temp, z_Press, z_MolWtMean, za_MolFr )
    !
    !== 概要
    !  * 湿潤断熱温度減率に沿った温度、圧力を求め、それらに対して指定された相対湿度となるように凝結成分のモル比を決める
    !  * 比熱は乾燥気塊のもので代表させる
    !    * 流体の方程式において比熱は乾燥成分のもので代表させているため
    !  * 大気の平均分子量には湿潤成分の分子量を効果
    !    * 流体の方程式において, 湿潤成分の分子量は考慮しているため
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    real(DP), intent(in) :: a_MolFrIni(1:SpcNum)    !下部境界でのモル比
    real(DP), intent(in) :: Humidity                !相対湿度 ( Humidity <= 1.0 )
    real(DP), intent(out):: z_Temp(DimZMin:DimZMax) !温度
    real(DP), intent(out):: z_Press(DimZMin:DimZMax)!圧力
    real(DP), intent(out):: z_MolWtMean(DimZMin:DimZMax) 
                                                   !平均分子量
    real(DP), intent(out):: za_MolFr(DimZMin:DimZMax, 1:SpcNum) 
                                                   !モル分率
    real(DP)             :: MolWtMeanDry            !乾燥気塊の分子量の作業変数
    real(DP)             :: MolWtMeanWet            !湿潤気塊の分子量の作業変数
    real(DP)             :: SatPress                !飽和蒸気圧
    real(DP)             :: VapPress                !蒸気圧
    real(DP)             :: DelMolFr
    real(DP)             :: a_MolFr(SpcNum)         !モル比の作業配列
    integer             :: k, s

    real(DP)             :: Temp1, Press1, DTempDZ1
    real(DP)             :: Temp2, Press2, DTempDZ2
    real(DP)             :: Temp3, Press3, DTempDZ3
    real(DP)             :: Temp4, Press4, DTempDZ4
    real(DP)             :: DTempDZ
    
    !-------------------------------------------------------------
    ! 配列の初期化
    !-------------------------------------------------------------
    !地表面の分子量を決める
    za_MolFr(RegZMin, 1:SpcNum)   = a_MolFrIni(1:SpcNum) 
    
    !地表面での平均分子量を決める
    MolWtMeanDry = MolWtDry * (1.0d0 - sum(a_MolFrIni))
    MolWtMeanWet = dot_product(MolWtWet, a_MolFrIni) 
    z_MolWtMean(RegZMin) = MolWtMeanDry + MolWtMeanWet
  
    !地表面での温度(RegZMin は, 高度 DelZ / 2 に相当)
    z_Temp          = 1.0d-60
    z_Temp(RegZMin) = TempSfc - Grav * z_MolWtMean(RegZMin) &
      &               / CpDryMol * ( z_dz(RegZMin) * 5.0d-1 )
    
    !地表面での圧力(RegZMin は, 高度 DelZ / 2 に相当)
    z_Press           = 1.0d-60
    z_Press(RegZMin)  = &
      & PressSfc *((TempSfc / z_Temp(RegZMin)) ** (- CpDryMol /  GasRUniv))
    
    !-----------------------------------------------------------
    ! 断熱減率 dT/dz の計算. 
    !-----------------------------------------------------------    
    DtDz: do k = RegZMin, DimZMax-1
      
      !初期化
      za_MolFr(k+1,:) = za_MolFr(k,:)
      
      !----------------------------------------------------
      !湿潤断熱減率をルンゲクッタ法を用いて計算
      !----------------------------------------------------
      ! (0) 高度 k での値を作業配列に保管
      Temp1  = z_Temp(k)
      Press1 = z_Press(k)
      a_MolFr  = za_MolFr(k,:)

      ! (1) 高度 k での値を用いて温度変化を計算
      call ECCM_DTempDZ( Temp1, Press1, r_dz(k), a_MolFr, DTempDZ1 )

      ! (2) (1) で求めた値を用いて, 高度 k + Δk/2 での値を用いて温度変化を計算
      !     このとき, 分子量は変化しないものとする.
      Temp2  = Temp1 + DTempDZ1 * r_dz(k) * 5.0d-1
      Press2 =                                     &
        & Press1 * ((Temp1 / Temp3) ** (Grav * MolWtDry / (GasRUniv * DTempDZ2))) 
      call ECCM_DTempDZ( Temp2, Press2, r_dz(k), a_MolFr, DTempDZ2 )

      ! (3) (2) で求めた値を用いて, 高度 k + Δk/2 での値を用いて温度変化を計算
      !     このとき, 分子量は変化しないものとする.
      Temp3  = Temp1 + DTempDZ2 * r_dz(k) * 5.0d-1
      Press3 =                                                                   &
        & Press1 * ((Temp1 / Temp3) ** (Grav * MolWtDry / (GasRUniv * DTempDZ2)))
      call ECCM_DTempDZ( Temp3, Press3, r_dz(k), a_MolFr, DTempDZ3 )
      
      ! (4) (3) で求めた値を用いて, 高度 k + Δk での値を用いて温度変化を計算
      !     このとき, 分子量は変化しないものとする.
      Temp4  = Temp1 + DTempDZ3 * r_dz(k)
      Press4 =                                               &
        & Press1 * ((Temp1 / Temp4) ** (Grav * MolWtDry / (GasRUniv * DTempDZ3)))
      call ECCM_DTempDZ( Temp4, Press4, r_dz(k), a_MolFr, DTempDZ4 )
      
      ! (5) 最終的な傾きを求める
      DTempDZ = (DTempDZ1 + DTempDZ2 * 2.0d0 + DTempDZ3 * 2.0d0 + DTempDZ4) / 6.0d0

      !----------------------------------------------------
      !得られた温度減率より温度と圧力を決める
      !----------------------------------------------------
      !温度を計算
      z_Temp(k+1) = z_Temp(k) + DTempDz * r_dz(k)

      !為念
      if(z_Temp(k+1) < 0.0d0) z_Temp(k+1) = z_Temp(k) 
      
      !圧力を静水圧平衡から計算
      z_Press(k+1) =                                                  &
        &  z_Press(k) * ( ( z_Temp(k) / z_Temp(k+1))                  &
        &    ** (Grav * MolWtDry / ( DTempDZ * GasRUniv ) ) )
      
      !----------------------------------------------------
      !モル比の計算
      !----------------------------------------------------
      do s = 1, CondNum      
        !飽和蒸気圧
        SatPress = SvapPress( SpcWetID(IdxCC(s)), z_Temp(k+1) )
        
        !元々のモル分率を用いて現在の蒸気圧を計算
        VapPress = za_MolFr(k,IdxCG(s)) * z_Press(k+1)
        
        !飽和蒸気圧と圧力から現在のモル比を計算
        if ( VapPress > SatPress ) then         
          za_MolFr(k+1,IdxCG(s)) = max(SatPress * Humidity / z_Press(k+1), 1.0d-16)
        end if
      end do
      
      !NH4SH の平衡条件
      if ( IdxNH3 /= 0 ) then 
        DelMolFr =                                              &
          & max (                                               &
          &    DelMolFrNH4SH(                                   &
          &         z_Temp(k+1), z_Press(k+1),                  &
          &         za_MolFr(k+1,IdxNH3), za_MolFr(k+1,IdxH2S),       &
          &         Humidity                                    &
          &      ),                                             &
          &    0.0d0                                            &
          &  )
        za_MolFr(k+1,IdxNH3) = za_MolFr(k+1,IdxNH3) - DelMolFr
        za_MolFr(k+1,IdxH2S) = za_MolFr(k+1,IdxH2S) - DelMolFr
      end if
      
      !------------------------------------------------------------
      !地表面での平均分子量を決める
      !------------------------------------------------------------
      MolWtMeanDry = MolWtDry * (1.0d0 - sum(za_MolFr(k+1, 1:SpcNum)))
      MolWtMeanWet = dot_product(MolWtWet(1:SpcNum), za_MolFr(k+1, 1:SpcNum))
      z_MolWtMean(k+1) = MolWtMeanDry + MolWtMeanWet

    end do DtDz
    
  end subroutine ECCM_Wet

!!!------------------------------------------------------------------------------!!!
  subroutine ECCM_MolFr( a_MolFrIni, Humidity, z_Temp, z_Press, za_MolFr )
    !
    ! 与えられた温度に対し, 気塊が断熱的に上昇した時に実現される
    ! モル比のプロファイルを求める
    !

    
    !暗黙の型宣言禁止
    implicit none
    
    real(DP), intent(in) :: a_MolFrIni(1:SpcNum)
    real(DP), intent(in) :: Humidity
    real(DP), intent(in) :: z_Temp(DimZMin:DimZMax)
    real(DP), intent(in) :: z_Press(DimZMin:DimZMax)
    real(DP), intent(out):: za_MolFr(DimZMin:DimZMax, 1:SpcNum)
    
    real(DP)             :: DelMolFr
    integer             :: k, s
    

    !-----------------------------------------------------------
    ! 配列の初期化
    !-----------------------------------------------------------
    do s = 1, SpcNum
      za_MolFr(:,s) = a_MolFrIni(s) 
    end do

    !-----------------------------------------------------------
    ! 断熱減率 dT/dz の計算. 
    !-----------------------------------------------------------
    do k = RegZMin, DimZMax

      za_MolFr(k,:) = za_MolFr(k-1,:)
      
      !------------------------------------------------------------
      !NH4SH 以外の化学種の平衡条件
      !------------------------------------------------------------
      do s = 1, CondNum

        !モル比を求める
        !モル比は前のステップでのモル比を超えることはない
        za_MolFr(k,IdxCG(s)) =                                 &
          & min(                                                &
          &       za_MolFr(k-1,IdxCG(s)),                      &
          &       SvapPress( SpcWetID(IdxCC(s)), z_Temp(k) ) &
          &        * Humidity / z_Press(k)                      &
          &      )
        
      end do

      !------------------------------------------------------------
      !NH4SH の平衡条件
      !------------------------------------------------------------
      if ( IdxNH3 /= 0 ) then 
        
        !モル比の変化. 
        !とりあえず NH4SH に対する飽和比は 1.0 とする(手抜き...).
        DelMolFr =                                            &
          & max (                                             &
          &    DelMolFrNH4SH(                                 &
          &      z_Temp(k), z_Press(k),                       &
          &      za_MolFr(k,IdxNH3), za_MolFr(k,IdxH2S), Humidity   &
          &     ),                                            &
          &    0.0d0                                          &
          &  )
        
        za_MolFr(k,IdxNH3) = za_MolFr(k,IdxNH3) - DelMolFr 
        za_MolFr(k,IdxH2S) = za_MolFr(k,IdxH2S) - DelMolFr
      end if
      
    end do
  end subroutine ECCM_MolFr

!!!------------------------------------------------------------------------------!!!
  subroutine ECCM_DTempDZ( Temp, Press, DelZ, MolFr, DTempDZ )
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP), intent(in) :: Temp
    real(DP), intent(in) :: Press
    real(DP), intent(in) :: DelZ
    real(DP), intent(inout) :: MolFr(0:SpcNum)    !モル分率
    real(DP), intent(out):: DTempDZ
    real(DP)            :: ReactHeat
    real(DP)            :: Heat(SpcNum)
    real(DP)            :: DelMolFr
    real(DP)            :: SatPress
    real(DP)            :: VapPress
    real(DP)            :: Humidity
    real(DP)            :: A, B
    integer             :: s

    !初期化
    DTempDZ      = 0.0d0
    ReactHeat    = 0.0d0
    Heat         = 0.0d0
    DelMolFr     = 0.0d0
    SatPress     = 0.0d0
    VapPress     = 0.0d0

    !------------------------------------------------------------
    !NH4SH 以外の化学種の平衡条件
    !------------------------------------------------------------
    do s = 1, CondNum      
      
      !飽和蒸気圧
      SatPress = SvapPress( SpcWetID(IdxCC(s)), Temp )
      
      !潜熱. 
      Heat(IdxCG(s)) = LatentHeatPerMol( SpcWetID(IdxCC(s)), Temp )
      
      !元々のモル分率を用いて現在の蒸気圧を計算
      VapPress = MolFr(IdxCG(s)) * Press
      
      !飽和蒸気圧から凝結の有無を決める
      if ( VapPress < SatPress ) then         
        !凝結していないので潜熱なし.
        Heat(IdxCG(s)) = 0.0d0          

      else      

        !飽和蒸気圧と圧力から現在のモル比を計算
        MolFr(IdxCG(s)) = max(SatPress / Press, 1.0d-16)

      end if
    end do
    
    !------------------------------------------------------------
    !NH4SH の平衡条件
    !------------------------------------------------------------
    if ( IdxNH3 /= 0 ) then 
      
      Humidity = 1.0d0
      DelMolFr =                                            &
        & max (                                             &
        &    DelMolFrNH4SH(                                 &
        &         Temp, Press, MolFr(IdxNH3), MolFr(IdxH2S),&
        &         Humidity                                  &
        &      ),                                           &
        &    0.0d0                                          &
        &  )
      MolFr(IdxNH3) = MolFr(IdxNH3) - DelMolFr
      MolFr(IdxH2S) = MolFr(IdxH2S) - DelMolFr

      ReactHeat = ReactHeatNH4SHPerMol * DelMolFr
    end if
    
    !------------------------------------------------------------
    !温度勾配を計算
    !------------------------------------------------------------
    !係数. 温度 Temp(i) で評価
    A = dot_product( Heat(1:SpcNum), MolFr(1:SpcNum)) &
      &  / ( GasRUniv * Temp )
    B = dot_product(( Heat(1:SpcNum) ** 2.0d0), MolFr(1:SpcNum)) &
      &  / ( CpDryMol * GasRUniv * ( Temp ** 2.0d0 ) )
    
    !断熱温度減率
    DTempDZ = - Grav * MolWtDry * (1.0d0 + A) / (CpDryMol * (1.0d0 + B))  &
      &       + ReactHeat / (CpDryMol * DelZ)
    
  end subroutine ECCM_DTempDZ

!!!------------------------------------------------------------------------------!!!
  subroutine ECCM_Stab( xyz_PotTemp, xyz_Exner, xyza_MixRt )

    use gridset_3d,only: DimXMin,       &! 配列の X 方向の下限
      &                 DimXMax,       &! 配列の X 方向の上限 
      &                 DimYMin,       &! 配列の Y 方向の下限
      &                 DimYMax,       &! 配列の Y 方向の上限 
      &                 DimZMin,       &! 配列の Z 方向の下限
      &                 DimZMax,       &! 配列の Z 方向の上限 
      &                 SpcNum          !
    use basicset_3d,only:  MolWtDry,      &!
      &                 MolWtWet,      &!
      &                 CpDry,         &!
      &                 Grav,          &!
      &                 xyz_ExnerBasicZ,   &!
      &                 xyz_PotTempBasicZ, &! 
      &                 xyz_EffMolWtBasicZ, &!
      &                 xyza_MixRtBasicZ

    use xyz_base_module, only : xyz_avr_xyr
    use xyz_deriv_module,only : xyr_dz_xyz
    
    implicit none

    real(DP), intent(in)  :: xyz_PotTemp(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
    real(DP), intent(in)  :: xyz_Exner(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
    real(DP), intent(in)  :: xyza_MixRt(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum)

    real(DP) :: xyz_Stab(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
    real(DP) :: xyz_StabTemp(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
    real(DP) :: xyz_StabMolWt(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)

    real(DP)   :: xyza_MolFrAll(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax,SpcNum)
    real(DP)   :: xyz_TempAll(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
    real(DP)   :: xyz_MolWtWet(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
    integer    :: i, j, k, s

    xyz_TempAll = (xyz_PotTemp + xyz_PotTempBasicZ) * (xyz_Exner + xyz_ExnerBasicZ)
    do s = 1, SpcNum
      xyza_MolFrAll(:,:,:,s) =                             &
        &   (xyza_MixRt(:,:,:,s) + xyza_MixRtBasicZ(:,:,:,s)) &
        &   * MolWtDry / MolWtWet(s) 
    end do
    
    do k = DimZMin, DimZMax
    do j = DimYMin, DimYMax
    do i = DimXMin, DimXMax
      xyz_MolWtWet(i,j,k) = &
        &     dot_product( MolWtWet(1:GasNum), xyza_MolFrAll(i,j,k,1:GasNum) )
    end do
    end do
    end do
    
    xyz_StabTemp =                                               &
      &         Grav / xyz_TempAll                               &
      &           * (   xyz_avr_xyr( xyr_dz_xyz( xyz_TempAll ) ) &
      &               + Grav * xyz_EffMolWtBasicZ / CpDry ) 
    xyz_StabMolWt =                                              &
      &       - Grav * xyz_avr_xyr( xyr_dz_xyz( xyz_MolWtWet ) ) &
      &         / ( MolWtDry * xyz_EffMolWtBasicZ )   
    xyz_Stab = xyz_StabTemp + xyz_StabMolWt

    call StoreStabTemp( xyz_StabTemp ) 
    call StoreStabMolWt( xyz_StabMolWt ) 

    where (xyz_Stab < 1.0d-7) 
      xyz_Stab = 1.0d-7
    end where

  end subroutine ECCM_Stab

end module ECCM_3d