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

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

    !地表面での温度(nx は, 高度 DelZ / 2 に相当)
    z_Temp    = 1.0d-60
    z_Temp(1) = TempSfc - Grav * MolWtDry / CpDryMol * ( z_dz(1) * 5.0d-1 )

!    write(*,*) z_Temp(1)
    
    !地表面での圧力(1 は, 高度 DelZ / 2 に相当)
    z_Press    = 1.0d-60
    z_Press(1) = PressSfc *((TempSfc / z_Temp(1)) ** (- CpDryMol /  GasRUniv))

!    write(*,*) z_Press(1)
    
    !-----------------------------------------------------------
    ! (1) 乾燥断熱線に沿った温度を決める
    ! (2) 静水圧平衡から圧力を求める
    ! (3) (1),(2) の温度圧力に対して, とある相対湿度となるモル比を決める
    !-----------------------------------------------------------    
    DtDz: do k = 1, kmax-1

      !(1)乾燥断熱線に沿って k+1 での温度を計算
      z_Temp(k+1) = z_Temp(k) - Grav * MolWtDry / 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
      
    end do DtDz
    
  end subroutine ECCM_Dry

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

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

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

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

      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( xyz_PTemp, xyz_Exner, xyzf_QMix )

    use gridset,only: imin, imax, jmin, jmax, kmin, kmax, ncmax          !
    use basicset,only: xyz_ExnerBZ, xyz_PTempBZ, xyzf_QMixBZ
    use constants,only: MolWtDry, CpDry, Grav          
    use composition,only:  MolWtWet
    use axesset, only : xyz_avr_xyr
    use xyz_deriv_module,only : xyr_dz_xyz
    
    implicit none

    real(DP), intent(in)  :: xyz_PTemp(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: xyz_Exner(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: xyzf_QMix(imin:imax,jmin:jmax,kmin:kmax, ncmax)

    real(DP) :: xyz_Stab(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xyz_StabTemp(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xyz_StabMolWt(imin:imax,jmin:jmax,kmin:kmax)

    real(DP)   :: xyzf_MolFrAll(imin:imax,jmin:jmax,kmin:kmax,ncmax)
    real(DP)   :: xyz_TempAll(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)   :: xyz_MolWtWet(imin:imax,jmin:jmax,kmin:kmax)
    integer    :: i, j, k, s

    xyz_TempAll = (xyz_PTemp + xyz_PTempBZ) * (xyz_Exner + xyz_ExnerBZ)
    do s = 1, ncmax
      xyzf_MolFrAll(:,:,:,s) = (xyzf_QMix(:,:,:,s) + xyzf_QMixBZ(:,:,:,s)) * MolWtDry / MolWtWet(s) 
    end do
    
    do k = kmin, kmax
    do j = jmin, jmax
    do i = imin, imax
      xyz_MolWtWet(i,j,k) = dot_product( MolWtWet(1:GasNum), xyzf_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 / CpDry ) 
    xyz_StabMolWt = - Grav * xyz_avr_xyr( xyr_dz_xyz( xyz_MolWtWet ) ) / MolWtDry 
    xyz_Stab = xyz_StabTemp + xyz_StabMolWt

    where (xyz_Stab < 1.0d-7) 
      xyz_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(DP), intent(in) :: a_MolFrIni(1:ncmax)    !下部境界でのモル比
    real(DP), intent(in) :: Humidity                !相対湿度 ( Humidity <= 1.0 )
    real(DP), intent(out):: z_Temp(kmin:kmax) !温度
    real(DP), intent(out):: z_Press(kmin:kmax)!圧力
                                                   !平均分子量
    real(DP), intent(out):: za_MolFr(kmin:kmax, 1:ncmax) 
                                                   !モル分率
    real(DP)             :: SatPress                !飽和蒸気圧
    real(DP)             :: VapPress                !蒸気圧
    real(DP)             :: DelMolFr
    real(DP)             :: a_MolFr(ncmax)         !モル比の作業配列
    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  = 1.0d-60
    za_MolFr(1, 1:ncmax)   = a_MolFrIni(1:ncmax) 
    
    !地表面での温度(1 は, 高度 DelZ / 2 に相当)
    z_Temp          = 1.0d-60
    z_Temp(1) = TempSfc - Grav * MolWtDry / CpDryMol * ( z_dz(1) * 5.0d-1 )
    
    !地表面での圧力(1 は, 高度 DelZ / 2 に相当)
    z_Press           = 1.0d-60
    z_Press(1)  = PressSfc *((TempSfc / z_Temp(1)) ** (- CpDryMol /  GasRUniv))
    
    !-----------------------------------------------------------
    ! 断熱減率 dT/dz の計算. 
    !-----------------------------------------------------------    
    DtDz: do k = 1, kmax-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 / Temp2) ** (Grav * MolWtDry / (GasRUniv * DTempDZ1))) 
      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
      
    end do DtDz
    
  end subroutine ECCM_Wet