subroutine ECCM_Dry( a_MolFrIni, Humidity, z_Temp, z_Press, 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)             :: z_MolWtMean(DimZMin:DimZMax) 
                                                   !平均分子量
    real(8), intent(out):: za_MolFr(DimZMin:DimZMax, 1:SpcNum) 
                                                   !モル分率
    real(8)             :: MolWtMeanDry            !乾燥気塊の分子量の作業変数
    real(8)             :: MolWtMeanWet            !湿潤気塊の分子量の作業変数
    real(8)             :: SatPress                !飽和蒸気圧
    real(8)             :: VapPress                !蒸気圧
    real(8)             :: DelMolFr
    integer             :: k, s

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

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

      !(1)乾燥断熱線に沿って k+1 での温度を計算
      z_Temp(k+1) = z_Temp(k) - Grav * MolWtDry / CpDryMol * DelZ
      
      !念為
      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
!        write(*,*) k, s, z_Temp(k), za_MolFr(k,IdxCG(s)), VapPress, SatPress
      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_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, 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_Stab( xz_PotTemp, xz_Exner, xza_MixRt )
!    &                   xz_Stab, xz_StabTemp, xz_StabMolWt )

    use gridset,  only: DimXMin, DimXMax, DimZMin, DimZMax, 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)    :: xz_Stab(DimXMin:DimXMax,   DimZMin:DimZMax)
    real(8)    :: xz_StabTemp(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)    :: 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:GasNum), xza_MolFrAll(i,k,1:GasNum) )
      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 = xz_StabTemp + xz_StabMolWt
!    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 )   

    call StoreStabTemp( xz_StabTemp ) 
    call StoreStabMolWt( xz_StabMolWt ) 

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

  end subroutine ECCM_Stab

  subroutine ECCM_Wet( a_MolFrIni, Humidity, z_Temp, z_Press, 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)             :: z_MolWtMean(DimZMin:DimZMax) 
                                                   !平均分子量
    real(8), intent(out):: za_MolFr(DimZMin:DimZMax, 1:SpcNum) 
                                                   !モル分率
    real(8)             :: MolWtMeanDry            !乾燥気塊の分子量の作業変数
    real(8)             :: MolWtMeanWet            !湿潤気塊の分子量の作業変数
    real(8)             :: SatPress                !飽和蒸気圧
    real(8)             :: VapPress                !蒸気圧
    real(8)             :: DelMolFr
    integer             :: k, s, j
    
    real(8)             :: TempA, PressA, a_MolFrA(SpcNum)
    real(8)             :: TempN, PressN, a_MolFrN(SpcNum)
    real(8)             :: DZ
    real(8)             :: A, B
    real(8)             :: Heat(SpcNum)
    real(8)             :: ReactHeat
    logical             :: MoistFlag(SpcNum), ReactFlag
    
    !-------------------------------------------------------------
    ! 配列の初期化
    !-------------------------------------------------------------
    !地表面の分子量を決める
    za_MolFr(RegZMin+1, 1:SpcNum)   = a_MolFrIni(1:SpcNum) 
    
    !地表面での平均分子量を決める
    MolWtMeanDry = MolWtDry * (1.0d0 - sum(a_MolFrIni))
    MolWtMeanWet = dot_product(MolWtWet, a_MolFrIni) 
    z_MolWtMean(RegZMin+1) = MolWtMeanDry + MolWtMeanWet
  
    !地表面での温度(RegZMin+1 は, 高度 DelZ / 2 に相当)
    z_Temp          = 1.0d-60
    z_Temp(RegZMin+1) = TempSfc - Grav * MolWtDry / CpDryMol * ( DelZ * 5.0d-1 )
    
    !地表面での圧力(RegZMin+1 は, 高度 DelZ / 2 に相当)
    z_Press             = 1.0d-60
    z_Press(RegZMin+1)  = PressSfc *((TempSfc / z_Temp(RegZMin+1)) ** (- CpDryMol /  GasRUniv))

    !-----------------------------------------------------------
    ! 断熱減率 dT/dz の計算. 
    !-----------------------------------------------------------    
    DtDz: do k = RegZMin+1, DimZMax-1
      TempN  = z_Temp(k)
      PressN = z_Press(k)
      a_MolFrN = za_MolFr(k,:)
      
      DZ = DelZ/100
      do j = 1, 100
        MoistFlag = .false.
        Heat = 0.0d0   
        a_MolFrA = a_MolFrN  !初期化
        
        !乾燥断熱線に沿って k+1 での温度を計算
        TempA = TempN - Grav * MolWtDry / CpDryMol * DZ
!        write(*,*) "TempA", TempA
        
        !念為
        if (TempA <= 0.0d0 ) TempA = TempN
        
        !圧力を静水圧平衡から計算
        PressA = PressN * ((TempN / TempA) ** (- CpDryMol / GasRUniv))
        
        ! 凝結が生じるか判定
        do s = 1, CondNum      
          !飽和蒸気圧
          SatPress = SvapPress( SpcWetID(IdxCC(s)), TempA )
          
          !元々のモル分率を用いて現在の蒸気圧を計算
          VapPress = a_MolFrN(IdxCG(s)) * PressA
          
          !飽和蒸気圧から凝結の有無を決める
          if ( VapPress < SatPress ) then         
            !凝結していないので潜熱なし.
            a_MolFrA(IdxCG(s)) = a_MolFrN(IdxCG(s))
            !凝結しないので潜熱なし
            Heat(IdxCG(s)) = 0.0d0
          else
            !潜熱. 現在の高度での値
            Heat(IdxCG(s)) = LatentHeatPerMol( SpcWetID(IdxCC(s)), TempN )
            !凝結のスイッチ
            MoistFlag(s) = .true.
          end if
        end do

        !NH4SH の平衡条件
        if ( IdxNH3 /= 0 ) then      
          DelMolFr = max ( DelMolFrNH4SH( TempA, PressA, a_MolFrN(IdxNH3), a_MolFrN(IdxH2S), Humidity ), 0.0d0 )
          ! 反応熱
          ReactHeat = ReactHeatNH4SHPerMol * DelMolFr
          ! 反応のスイッチ
          if (DelMolFr > 0) ReactFlag = .true. 
        end if

!        write(*,*) "MoistFlag", MoistFlag
!        write(*,*) "ReactFlag", ReactFlag
        
        if (count(MoistFlag) > 0 .or. ReactFlag) then 
          !係数. 高度 k での値で評価
          A = dot_product( Heat(1:SpcNum), a_MolFrN(1:SpcNum)) / ( GasRUniv * TempN )
          B = dot_product(( Heat(1:SpcNum) ** 2.0d0), a_MolFrN(1:SpcNum)) / ( CpDryMol * GasRUniv * ( TempN ** 2.0d0 ) )
          
          !断熱温度減率
          TempA = TempN - Grav * MolWtDry * (1.0d0 + A) / (CpDryMol * (1.0d0 + B)) * DZ + ReactHeat / (CpDryMol * DelZ)
!          write(*,*) "Moist TempA",   TempA
!          write(*,*) "Moist Heat ",   Heat
!          write(*,*) "Moist MolFr ",  a_MolFrN
!          write(*,*) "Moist DTempDZ", - Grav * MolWtDry * (1.0d0 + A) / (CpDryMol * (1.0d0 + B))
          
          !念為
          if (TempA <= 0.0d0 ) TempA = TempN
          
          !圧力を静水圧平衡から計算
          PressA = PressN * ((TempN / TempA) ** (- CpDryMol / GasRUniv)) 
          
          do s = 1, CondNum              
            if (MoistFlag(s)) then 
              !飽和蒸気圧
              SatPress = SvapPress( SpcWetID(IdxCC(s)), TempA )
              ! 飽和モル比
              a_MolFrA(IdxCG(s)) = max(SatPress / PressA, 1.0d-16)
            else
              ! 前のモル比と同じ
              a_MolFrA(IdxCG(s)) = a_MolFrN(IdxCG(s)) 
            end if
          end do

          if ( ReactFlag ) then 
            DelMolFr = max ( DelMolFrNH4SH( TempA, PressA, a_MolFrA(IdxNH3), a_MolFrA(IdxH2S), Humidity ), 0.0d0 )
            a_MolFrA(IdxNH3) = a_MolFrA(IdxNH3) - DelMolFr
            a_MolFrA(IdxH2S) = a_MolFrA(IdxH2S) - DelMolFr
          end if
        end if
        
        TempN  = TempA
        PressN = PressA
        a_MolFrN = a_MolFrA        
      end do
        
      z_Temp(k+1)     = TempA
      z_Press(k+1)    = PressA
      za_MolFr(k+1,:) = a_MolFrA(:)
      
      !平均分子量を決める
      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_Wet2( Idx, Humidity, z_Temp, z_Press, za_MolFr)
    !
    !== 概要
    !  * 湿潤断熱温度減率に沿った温度、圧力を求め、それらに対して指定された相対湿度となるように凝結成分のモル比を決める
    !  * 比熱は乾燥気塊のもので代表させる
    !    * 流体の方程式において比熱は乾燥成分のもので代表させているため
    !  * 大気の平均分子量には湿潤成分の分子量を効果
    !    * 流体の方程式において, 湿潤成分の分子量は考慮しているため
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    integer             :: Idx
    real(8), intent(inout) :: z_Temp(DimZMin:DimZMax)    !下部境界でのモル比
    real(8), intent(inout) :: z_Press(DimZMin:DimZMax) !下部境界でのモル比
    real(8), intent(inout) :: za_MolFr(DimZMin:DimZMax, 1:SpcNum)    !下部境界でのモル比
    real(8), intent(in) :: Humidity                !相対湿度 ( Humidity <= 1.0 )
    real(8)             :: SatPress                !飽和蒸気圧
    real(8)             :: VapPress                !蒸気圧
    real(8)             :: DelMolFr
    integer             :: k, s, j
    
    real(8)             :: TempA, PressA, a_MolFrA(SpcNum)
    real(8)             :: TempN, PressN, a_MolFrN(SpcNum)
    real(8)             :: DZ
    real(8)             :: A, B
    real(8)             :: Heat(SpcNum)
    real(8)             :: ReactHeat
    logical             :: MoistFlag(SpcNum), ReactFlag
    

    !-----------------------------------------------------------
    ! 断熱減率 dT/dz の計算. 
    !-----------------------------------------------------------    
    DtDz: do k = Idx, DimZMax-1
      TempN  = z_Temp(k)
      PressN = z_Press(k)
      a_MolFrN = za_MolFr(k,:)
      
      DZ = DelZ/100
      do j = 1, 100
        MoistFlag = .false.
        Heat = 0.0d0   
        a_MolFrA = a_MolFrN  !初期化
        
        !乾燥断熱線に沿って k+1 での温度を計算
        TempA = TempN - Grav * MolWtDry / CpDryMol * DZ
!        write(*,*) "TempA", TempA
        
        !念為
        if (TempA <= 0.0d0 ) TempA = TempN
        
        !圧力を静水圧平衡から計算
        PressA = PressN * ((TempN / TempA) ** (- CpDryMol / GasRUniv))
        
        ! 凝結が生じるか判定
        do s = 1, CondNum      
          !飽和蒸気圧
          SatPress = SvapPress( SpcWetID(IdxCC(s)), TempA )
          
          !元々のモル分率を用いて現在の蒸気圧を計算
          VapPress = a_MolFrN(IdxCG(s)) * PressA
          
          !飽和蒸気圧から凝結の有無を決める
          if ( VapPress < SatPress ) then         
            !凝結していないので潜熱なし.
            a_MolFrA(IdxCG(s)) = a_MolFrN(IdxCG(s))
            !凝結しないので潜熱なし
            Heat(IdxCG(s)) = 0.0d0
          else
            !潜熱. 現在の高度での値
            Heat(IdxCG(s)) = LatentHeatPerMol( SpcWetID(IdxCC(s)), TempN )
            !凝結のスイッチ
            MoistFlag(s) = .true.
          end if
        end do

        !NH4SH の平衡条件
        if ( IdxNH3 /= 0 ) then      
          DelMolFr = max ( DelMolFrNH4SH( TempA, PressA, a_MolFrN(IdxNH3), a_MolFrN(IdxH2S), Humidity ), 0.0d0 )
          ! 反応熱
          ReactHeat = ReactHeatNH4SHPerMol * DelMolFr
          ! 反応のスイッチ
          if (DelMolFr > 0) ReactFlag = .true. 
        end if

!        write(*,*) "MoistFlag", MoistFlag
!        write(*,*) "ReactFlag", ReactFlag
        
        if (count(MoistFlag) > 0 .or. ReactFlag) then 
          !係数. 高度 k での値で評価
          A = dot_product( Heat(1:SpcNum), a_MolFrN(1:SpcNum)) / ( GasRUniv * TempN )
          B = dot_product(( Heat(1:SpcNum) ** 2.0d0), a_MolFrN(1:SpcNum)) / ( CpDryMol * GasRUniv * ( TempN ** 2.0d0 ) )
          
          !断熱温度減率
          TempA = TempN - Grav * MolWtDry * (1.0d0 + A) / (CpDryMol * (1.0d0 + B)) * DZ + ReactHeat / (CpDryMol * DZ)
!          write(*,*) "Moist TempA",   TempA
!          write(*,*) "Moist Heat ",   Heat
!          write(*,*) "Moist MolFr ",  a_MolFrN
!          write(*,*) "Moist DTempDZ", - Grav * MolWtDry * (1.0d0 + A) / (CpDryMol * (1.0d0 + B))
          
          !念為
          if (TempA <= 0.0d0 ) TempA = TempN
          
          !圧力を静水圧平衡から計算
          PressA = PressN * ((TempN / TempA) ** (- CpDryMol / GasRUniv)) 
          
          do s = 1, CondNum              
            if (MoistFlag(s)) then 
              !飽和蒸気圧
              SatPress = SvapPress( SpcWetID(IdxCC(s)), TempA )
              ! 飽和モル比
              a_MolFrA(IdxCG(s)) = max(SatPress / PressA, 1.0d-16)
            else
              ! 前のモル比と同じ
              a_MolFrA(IdxCG(s)) = a_MolFrN(IdxCG(s)) 
            end if
          end do

          if ( ReactFlag ) then 
            DelMolFr = max ( DelMolFrNH4SH( TempA, PressA, a_MolFrA(IdxNH3), a_MolFrA(IdxH2S), Humidity ), 0.0d0 )
            a_MolFrA(IdxNH3) = a_MolFrA(IdxNH3) - DelMolFr
            a_MolFrA(IdxH2S) = a_MolFrA(IdxH2S) - DelMolFr
          end if
       end if
        
        TempN  = TempA
        PressN = PressA
        a_MolFrN = a_MolFrA        
      end do
        
      z_Temp(k+1)     = TempA
      z_Press(k+1)    = PressA
      za_MolFr(k+1,:) = a_MolFrA(:)
      
   end do DtDz
    
  end subroutine ECCM_Wet2