!= Module ECCM
!
! Authors::   SUGIYAMA Koichiro, ODAKA Masatsugu
! Version::   $Id: eccm.f90,v 1.18 2006/11/22 15:40:55 sugiyama Exp $
! Tag Name::  $Name: arare4-20061224 $
! 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

  !モジュール読み込み
  use dc_message, only: MessageNotify

  use gridset,  only: DimZMin,       &! 配列の Z 方向の下限
    &                 DimZMax,       &! 配列の Z 方向の上限 
    &                 RegZMin,       &!
    &                 SpcNum,        &!
    &                 DelZ            !
  use basicset,only:  MolWtDry,      &!
    &                 MolWtWet,      &!
    &                 CpDry,         &!
    &                 CpDryMol,      &!
    &                 SpcWetID,      &!
    &                 SpcWetSymbol,  &!
    &                 TempSfc,       &!
    &                 PressSfc,      &!
    &                 Grav            !
  use chemcalc, only: SvapPress,     &!
    &                 LatentHeatPerMol, &!
    &                 ReactHeatNH4SHPerMol
  use ChemData, only: GasRUniv,      &!
    &                 ChemData_OneSpcID
  use MoistFunc,only: DelMolFrNH4SH 

  !暗黙の型宣言禁止
  implicit none

  !属性の指定
  private

  !関数の公開
  public ECCM_Init
!  public ECCM_Temp_Press
  public ECCM_MolFr
  public ECCM_Stab
  public ECCM_Dry

  !変数定義
  integer     :: LoopNum = 0
  integer     :: GasNum(10) = 0
  integer     :: CloudNum(10) = 0
  integer     :: NH3Num = 0
  integer     :: H2SNum = 0
  
  save LoopNum, GasNum, CloudNum
  save NH3Num, H2SNum
  
contains

  subroutine ECCM_Init( )

    !暗黙の型宣言禁止
    implicit none

    !変数定義
    integer                  :: s
    integer                  :: n1

    !-----------------------------------------------------------
    ! 雲粒と気体の ID の組を作る
    !-----------------------------------------------------------
    !初期化
    LoopNum = 0

    !化学種の中から雲粒を作るものを選び, その配列添え字と分子量を保管.
    SelectCloud: do s = 1, SpcNum
      
      ! NH4SH については無視
      if ( trim(SpcWetSymbol(s)) == 'NH4SH-s-Cloud' ) then 
        cycle SelectCloud
      end if

      !'Cloud' という文字列が含まれるものの個数を数える
      n1 = index(SpcWetSymbol(s), '-Cloud' )
      if (n1 /= 0) then
        LoopNum          = LoopNum + 1
        CloudNum(LoopNum)= s
        GasNum(LoopNum)  = minloc(SpcWetID, 1, SpcWetID == ChemData_OneSpcID(SpcWetSymbol(s)(1:n1-3) // '-g'))
      end if
    end do SelectCloud
    
    !-----------------------------------------------------------
    ! 硫化アンモニウム, およびアンモニアと硫化水素の ID を取得
    !-----------------------------------------------------------
    NH3Num   = minloc(SpcWetID, 1, SpcWetID == ChemData_OneSpcID('NH3-g'))
    H2SNum   = minloc(SpcWetID, 1, SpcWetID == ChemData_OneSpcID('H2S-g'))

    !-----------------------------------------------------------
    ! 確認
    !-----------------------------------------------------------
    if ( LoopNum == 0 ) then 
      call MessageNotify( "W",           &
        &  "ECCM_init", "CloudNum = 0, please comment out of MoistAdjust" )
!      write(*,*) "ECCM: CloudNum = 0, please comment out of MoistAdjust"
!      stop
    end if
    
    call MessageNotify( "M", "ECCM_Init", "LoopNum = %d",  i=(/LoopNum/)  )
    call MessageNotify( "M", "ECCM_Init", "CloudNum = %d", i=(/CloudNum/) )
    call MessageNotify( "M", "ECCM_Init", "GasNum = %d",   i=(/GasNum/)   ) 
    call MessageNotify( "M", "ECCM_Init", "NH3Num = %d",   i=(/NH3Num/)   )
    call MessageNotify( "M", "ECCM_Init", "H2SNum = %d",   i=(/H2SNum/)   )

  end subroutine ECCM_Init


!!!------------------------------------------------------------------------------!!!
  subroutine ECCM_Dry( a_MolFrIni, Humidity, z_Temp, z_Press, z_MolWtMean, za_MolFr )
    
    !暗黙の型宣言禁止
    implicit none
    
    real(8), intent(in) :: a_MolFrIni(1:SpcNum)    !下部境界でのモル比
    real(8), intent(in) :: Humidity                !相対湿度 ( Humidity <= 1.0 )
    real(8), intent(out):: z_Temp(DimZMin:DimZMax) !温度
    real(8), intent(out):: z_Press(DimZMin:DimZMax)!圧力
    real(8), intent(out):: z_MolWtMean(DimZMin:DimZMax) 
                                                   !平均分子量
    real(8), intent(out):: za_MolFr(DimZMin:DimZMax, 1:SpcNum) 
                                                   !モル分率
    real(8)             :: a_MolWt(0:SpcNum)       !分子量の作業配列
    real(8)             :: a_MolFr(0:SpcNum)       !モル比の作業配列
    real(8)             :: SatPress                !飽和蒸気圧
    real(8)             :: VapPress                !蒸気圧
    real(8)             :: DelMolFr
    integer             :: k, s

    
    
    !-------------------------------------------------------------
    ! 配列の初期化
    !-------------------------------------------------------------
    !分子量の配列を用意
    a_MolWt(0)        = MolWtDry
    a_MolWt(1:SpcNum) = MolWtWet(1:SpcNum)

    !初期モル比
    za_MolFr = 0.0d0
    za_MolFr(RegZMin, 1:SpcNum)   = a_MolFrIni(1:SpcNum) 

    !初期の平均分子量
    a_MolFr(0)        = 1.0d0 - sum(a_MolFrIni)
    a_MolFr(1:SpcNum) = a_MolFrIni
    z_MolWtMean(RegZMin) = dot_product(a_MolWt, a_MolFr)
    
    !地表面での温度(RegZMin は, 高度 DelZ / 2 に相当)
    z_Temp          = 1.0d-60
    z_Temp(RegZMin) = TempSfc - Grav * z_MolWtMean(RegZMin) &
      &               / CpDryMol * ( DelZ * 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 * DelZ
      
!      !念為
!      if (z_Temp(k+1) <= 0.0d0 ) exit DtDz
      
      !(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, LoopNum      
        !飽和蒸気圧
        SatPress = SvapPress( SpcWetID(CloudNum(s)), z_Temp(k+1) )
        
        !元々のモル分率を用いて現在の蒸気圧を計算
        VapPress = za_MolFr(k,GasNum(s)) * z_Press(k+1)
        
        !飽和蒸気圧と圧力から現在のモル比を計算
        if ( VapPress > SatPress ) then         
          za_MolFr(k+1,GasNum(s)) = max(SatPress * Humidity / z_Press(k+1), 1.0d-16)
        end if
      end do
      
      !NH4SH の平衡条件
      if ( NH3Num /= 0 ) then 
        DelMolFr =                                              &
          & max (                                               &
          &    DelMolFrNH4SH(                                   &
          &         z_Temp(k+1), z_Press(k+1),                  &
          &         za_MolFr(k+1,NH3Num), za_MolFr(k+1,H2SNum), &
          &         Humidity                                    &
          &      ),                                             &
          &    0.0d0                                            &
          &  )
        za_MolFr(k+1,NH3Num) = za_MolFr(k+1,NH3Num) - DelMolFr
        za_MolFr(k+1,H2SNum) = za_MolFr(k+1,H2SNum) - DelMolFr
      end if
      
      !------------------------------------------------------------
      !温度勾配を計算
      !------------------------------------------------------------
      !モル比
      a_MolFr(0) = 1.0d0 - sum( za_MolFr(k+1, 1:SpcNum) )
      a_MolFr(1:SpcNum) = za_MolFr(k+1, 1:SpcNum)

      !平均分子量
      z_MolWtMean(k+1) = dot_product(a_MolWt, a_MolFr)
      
    end do DtDz
    
  end subroutine ECCM_Dry



!!!------------------------------------------------------------------------------!!!
!  subroutine ECCM_Temp_Press( MolFrIni, z_Temp, z_Press )
!
! このサブルーチンでは, 分子量の変化が考慮されていない. 編集が必要だが, 
! 時間がないので, とりあえずコメントアウト 
!    
!    !暗黙の型宣言禁止
!    implicit none
!    
!    real(8), intent(in) :: MolFrIni(1:SpcNum)      !下部境界でのモル比
!    real(8), intent(out):: z_Temp(DimZMin:DimZMax) !温度
!    real(8), intent(out):: z_Press(DimZMin:DimZMax)!圧力
!    real(8)             :: z_MolFr(DimZMin:DimZMax, 0:SpcNum) 
!                                                   !モル分率
!    real(8)             :: z_DTempDZ(DimZMin:DimZMax) 
!                                                   !温度減率
!!    real(8)             :: z_MolWtMean(DimZMin:DimZMax) 
!!                                                   !平均分子量
!    real(8)             :: MolFr(0:SpcNum)         !モル比の作業配列
!!    real(8)             :: MolWt(0:SpcNum)         !分子量の作業配列
!    integer             :: k
!
!    real(8)             :: Temp1, Press1, DTempDZ1
!    real(8)             :: Temp2, Press2, DTempDZ2
!    real(8)             :: Temp3, Press3, DTempDZ3
!    real(8)             :: Temp4, Press4, DTempDZ4
!    
!    
!    !-------------------------------------------------------------
!    ! 配列の初期化
!    !-------------------------------------------------------------
!    !初期モル比
!    z_MolFr = 0.0d0
!    z_MolFr(RegZMin, 1:SpcNum)   = MolFrIni(1:SpcNum) 
!    z_MolFr(RegZMin, 0)          = 1.0d0 - sum(MolFrIni)
!    z_MolFr(RegZMin-1, 1:SpcNum) = MolFrIni(1:SpcNum) 
!    z_MolFr(RegZMin-1, 0)        = 1.0d0 - sum(MolFrIni)
!    
!!    !分子量の配列を用意
!!    MolWt(0)        = MolWtDry
!!    MolWt(1:SpcNum) = MolWtWet(1:SpcNum)
!
!    !地表面での温度(RegZMin は, 高度 DelZ / 2 に相当)
!    z_Temp          = 1.0d-60
!    z_Temp(RegZMin) = TempSfc - Grav / CpDry * ( DelZ * 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
!      
!      !初期化
!      z_MolFr(k,:) = z_MolFr(k-1,:)
!      
!      !==============================================================
!      !湿潤断熱減率をルンゲクッタ法を用いて計算
!      Temp1  = z_Temp(k)
!      Press1 = z_Press(k)
!      MolFr  = z_MolFr(k-1,:)
!      call ECCM_DTempDZ( Temp1, Press1, z_MolFr(k,:), DTempDZ1 )
!      
!      Temp2  = Temp1 + DTempDZ1 * DelZ * 5.0d-1
!      if (Temp2 < 0.0d0) exit DtDz
!      Press2 =                                               &
!        & Press1 * ( ( Temp1 / Temp2 )                       &
!        &      ** (Grav * MolWtDry    / (GasRUniv * DTempDZ1) ) )
!      call ECCM_DTempDZ( Temp2, Press2, MolFr, DTempDZ2 )
!
!      Temp3  = Temp1 + DTempDZ2 * DelZ * 5.0d-1
!      if (Temp3 < 0.0d0) exit DtDz
!      Press3 =                                               &
!        & Press1 * ( ( Temp1 / Temp3 )                       &
!        &      ** (Grav * MolWtDry    / (GasRUniv * DTempDZ2) ) )
!      call ECCM_DTempDZ( Temp3, Press3, MolFr, DTempDZ3 )
!      
!      Temp4  = Temp1 + DTempDZ3 * DelZ  
!      if (Temp4 < 0.0d0) exit DtDz
!      Press4 =                                               &
!        & Press1 * ( ( Temp1 / Temp4 )                       &
!        &      ** (Grav * MolWtDry    / (GasRUniv * DTempDZ3) ) )
!      call ECCM_DTempDZ( Temp4, Press4, MolFr, DTempDZ4 )
!      
!      z_DTempDZ(k) = ( + DTempDZ1 + DTempDZ2 * 2.0d0           &
!        &              + DTempDZ3 * 2.0d0 + DTempDZ4 ) / 6.0d0
!
!      !==============================================================
!
!      !温度を計算
!      z_Temp(k+1) = z_Temp(k) + z_DTempDz(k) * DelZ
!      
!      !圧力を静水圧平衡から計算
!      z_Press(k+1) =                                                  &
!        &  z_Press(k) * ( ( z_Temp(k) / z_Temp(k+1))                  &
!        &    ** (Grav * MolWtDry / ( z_DTempDZ(k) * GasRUniv ) ) )
!
!      !モル比 (DimZMax でも値を持つようにするため)
!      z_MolFr(k+1,:) = z_MolFr(k,:)
!      z_MolFr(k+1,0) = 1.0d0 - sum( z_MolFr(k,1:SpcNum) )
!      
!      if (z_Temp(k+1) <= 0.0d0 ) exit DtDz
!
!    end do DtDz
!    
!  end subroutine ECCM_Temp_Press



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

    
    !暗黙の型宣言禁止
    implicit none
    
    real(8), intent(in) :: a_MolFrIni(1:SpcNum)
    real(8), intent(in) :: Humidity
    real(8), intent(in) :: z_Temp(DimZMin:DimZMax)
    real(8), intent(in) :: z_Press(DimZMin:DimZMax)
    real(8), intent(out):: za_MolFr(DimZMin:DimZMax, 1:SpcNum)
    
    real(8)             :: 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, LoopNum      

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

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



!  subroutine ECCM_DTempDZ( Temp, Press, MolFr, DTempDZ )
!    
! このサブルーチンでは, 分子量の変化が考慮されていない. 編集が必要だが, 
! 時間がないので, とりあえずコメントアウト 
!
!    !暗黙の型宣言禁止
!    implicit none
!    
!    !変数定義
!    real(8), intent(in) :: Temp
!    real(8), intent(in) :: Press
!    real(8), intent(inout) :: MolFr(0:SpcNum)    !モル分率
!    real(8), intent(out):: DTempDZ
!    real(8)             :: ReactHeat
!    real(8)             :: Heat(SpcNum)
!    real(8)             :: DelMolFr
!    real(8)             :: SatPress
!    real(8)             :: VapPress
!    real(8)             :: Humidity
!    real(8)             :: 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, LoopNum      
!      
!      !飽和蒸気圧
!      SatPress = SvapPress( SpcWetID(CloudNum(s)), Temp )
!      
!      !潜熱. 
!      Heat(GasNum(s)) = &
!        & LatentHeatPerMol( SpcWetID(CloudNum(s)), Temp )
!        
!      !元々のモル分率を用いて現在の蒸気圧を計算
!      VapPress = MolFr(GasNum(s)) * Press
!      
!      !飽和蒸気圧から凝結の有無を決める
!      if ( VapPress < SatPress ) then         
!        !モル比は変化しない
!        MolFr(GasNum(s)) = MolFr(GasNum(s))
!        
!        !凝結していないので潜熱なし.
!        Heat(GasNum(s)) = 0.0d0          
!      else      
!        !飽和蒸気圧と圧力から現在のモル比を計算
!        MolFr(GasNum(s)) = max(SatPress / Press, 1.0d-16)
!      end if
!    end do
!
!    !------------------------------------------------------------
!    !NH4SH の平衡条件
!    !------------------------------------------------------------
!    if ( NH3Num /= 0 ) then 
!
!      Humidity = 1.0d0
!      DelMolFr =                                            &
!        & max (                                             &
!        &    DelMolFrNH4SH(                                 &
!        &         Temp, Press, MolFr(NH3Num), MolFr(H2SNum),&
!        &         Humidity                                  &
!        &      ),                                           &
!        &    0.0d0                                          &
!        &  )
!      MolFr(NH3Num) = MolFr(NH3Num) - DelMolFr
!      MolFr(H2SNum) = MolFr(H2SNum) - DelMolFr
!
!      ReactHeat = ReactHeatNH4SHPerMol * DelMolFr
!    end if
!
!    !------------------------------------------------------------
!    !温度勾配を計算
!    !------------------------------------------------------------
!    !モル比
!    MolFr(0) = 1.0d0 - sum( MolFr(1:SpcNum) )
!
!    !分子量
!
!    !係数. 温度 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 )
!   
!!    write(*,*) A, B, ReactHeat, DTempDZ 
!!    write(*,*) MolWtDry, CpDryMol 
!!    write(*,*) CpDryMol / MolWtDry, CpDry
!  end subroutine ECCM_DTempDZ



  subroutine ECCM_Stab( xz_PotTemp, xz_Exner, xza_MixRt, &
    &                   xz_Stab, xz_StabTemp, xz_StabMolWt )

    use gridset,  only: DimXMin,       &! 配列の X 方向の下限
      &                 DimXMax,       &! 配列の X 方向の上限 
      &                 DimZMin,       &! 配列の Z 方向の下限
      &                 DimZMax,       &! 配列の Z 方向の上限 
      &                 SpcNum          !
    use basicset,only:  MolWtDry,      &!
      &                 MolWtWet,      &!
      &                 CpDry,         &!
      &                 Grav,          &!
      &                 xz_ExnerBasicZ,   &!
      &                 xz_PotTempBasicZ, &! 
      &                 xz_EffMolWtBasicZ, &!
      &                 xza_MixRtBasicZ
    use average, only:  xz_avr_xr
    use differentiate_center2, only: xr_dz_xz
    
    implicit none

    real(8), intent(in)  :: xz_PotTemp(DimXMin:DimXMax,DimZMin:DimZMax)
    real(8), intent(in)  :: xz_Exner(DimXMin:DimXMax,  DimZMin:DimZMax)
    real(8), intent(in)  :: xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
    real(8), intent(out) :: xz_Stab(DimXMin:DimXMax,   DimZMin:DimZMax)
    real(8), intent(out) :: xz_StabTemp(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8), intent(out) :: xz_StabMolWt(DimXMin:DimXMax, DimZMin:DimZMax)

    real(8)    :: xza_MolFrAll(DimXMin:DimXMax,DimZMin:DimZMax,SpcNum)
    real(8)    :: xz_TempAll(DimXMin:DimXMax,  DimZMin:DimZMax)
    real(8)    :: xz_MolWtWet(DimXMin:DimXMax, DimZMin:DimZMax)
    integer    :: i, k, s

    xz_TempAll = (xz_PotTemp + xz_PotTempBasicZ) * (xz_Exner + xz_ExnerBasicZ)
    do s = 1, SpcNum
      xza_MolFrAll(:,:,s) =                             &
        &   (xza_MixRt(:,:,s) + xza_MixRtBasicZ(:,:,s)) &
        &   * MolWtDry / MolWtWet(s) 
    end do
    
    do k = DimZMin, DimZMax
      do i = DimXMin, DimXMax
        xz_MolWtWet(i,k) = dot_product( MolWtWet(1:3), xza_MolFrAll(i,k,1:3) )
      end do
    end do
    
    xz_StabTemp =                                           &
      &         Grav / xz_TempAll                           &
      &           * (   xz_avr_xr( xr_dz_xz( xz_TempAll ) ) &
      &               + Grav * xz_EffMolWtBasicZ / CpDry ) 
    xz_StabMolWt =                                          &
      &       - Grav * xz_avr_xr( xr_dz_xz( xz_MolWtWet ) ) &
      &         / ( MolWtDry * xz_EffMolWtBasicZ )   
    xz_Stab =                                               &
      &         Grav / xz_TempAll                           &
      &           * ( xz_avr_xr( xr_dz_xz( xz_TempAll ) )   &
      &               + Grav * xz_EffMolWtBasicZ / CpDry )  &
      &       - Grav * xz_avr_xr( xr_dz_xz( xz_MolWtWet ) ) &
      &         / ( MolWtDry * xz_EffMolWtBasicZ )   

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

  end subroutine ECCM_Stab

end module ECCM
