!= Module ECCM
!
! Authors::   SUGIYAMA Koichiro, ODAKA Masatsugu
! Version::   $Id: eccm.f90,v 1.29 2011-02-28 12:00:24 sugiyama Exp $
! Tag Name::  $Name:  $
! 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 gtool_history                   !外部にあるnetCDFからデータを読み込むための
                                      !モジュール(for ECC_N1994)
  use dc_message, only: MessageNotify

  use gridset,  only: DimZMin,       &! 配列の Z 方向の下限
    &                 DimZMax,       &! 配列の Z 方向の上限 
    &                 RegZMin,       &!
    &                 SpcNum
  use axesset, only:  s_Z,           &!
    &                 DelZ            !
  use constants,only: MolWtDry,     &!
    &                 CpDryMol,      &!
    &                 TempSfc,       &!
    &                 PressSfc,      &!
    &                 Grav,          &!
    &                 GasRDry,       &!追加(for ECC_N1994)
    &                 CpDry           !追加(for ECC_N1994)
  use chemcalc, only: SvapPress,        &!
    &                 LatentHeatPerMol, &!
    &                 ReactHeatNH4SHPerMol
  use composition,only: MolWtWet,       &!
    &                 SpcWetID,         &!
    &                 CondNum,          &!凝結過程の数
    &                 IdxCG,            &!凝結過程(蒸気)の配列添え字
    &                 IdxCC,            &!凝結過程(雲)の配列添え字
    &                 GasNum,           &!気体の数
    &                 IdxNH3,           &!NH3(蒸気)の配列添え字
    &                 IdxH2S             !H2S(蒸気)の配列添え字
  use ChemData, only: GasRUniv        
  use MoistFunc,only: DelMolFrNH4SH 
!  use StoreStab,only: StoreStabTemp, StoreStabMolWt

  !暗黙の型宣言禁止
  implicit none

  !属性の指定
  private

  !関数の公開
  public ECCM_init           !NAMELISTを使用するための初期化ルーチン 
  public ECCM_MolFr
  public ECCM_Stab
  public ECCM_Dry
  public ECCM_Wet
  public ECCM_Wet2
  public HUM                 !ECCM_N1994で用いる相対湿度の基本場を作成
  public ECCM_N1994            !Yamasaki(1983)の温度と相対湿度で基本場を作成
  public HUM_Takemi2007      !Takemi(2007)の相対湿度の基本場(一部)を作成
  public HUM_Takemi2007_TDRY !Takemi(2007)の相対湿度の基本場(一部)を作成
  public HUM_Takemi2007_MDRY !Takemi(2007)の相対湿度の基本場(一部)を作成
  public ECCM_Takemi2007     !Takemi(2007)の温度と湿度で基本場を作成
  public HUM_Takemi2007ALL   !Takemi(2007)の相対湿度の基本場を作成(for ECCM_MolFr)

  real(8)  :: DryHeight = 0.0d0

  save DryHeight
  
contains

!!!-------------------------------------------------------------------!!!
  subroutine ECCM_init(cfgfile)
    !
    ! ECCM モジュールの初期化ルーチン
    !

    !暗黙の型宣言禁止
    implicit none

    !
    !NAMELIST から放射強制の設定を取得
    !

  character(*), intent(in) :: cfgfile

    ! NAMELIST から情報を取得
    NAMELIST /eccm/ DryHeight

    open (10, FILE=cfgfile)
    read(10, NML=eccm)
    close(10)

    !確認
    call MessageNotify( "M", "ECCM_init", "DryHeight = %f", d=(/DryHeight/) )

  end subroutine ECCM_init
!!!------------------------------------------------------------------------------!!!
!  subroutine ECCM_Dry( a_MolFrIni, Humidity, z_Temp, z_Press, z_MolWtMean, za_MolFr )
  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_Wet( a_MolFrIni, Humidity, z_Temp, z_Press, z_MolWtMean, za_MolFr )
  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 は, 高度 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

!!------------------------------------------------------------------------------!!!
  subroutine HUM(z_humid, kansoku_altitude, kansoku_HumZ)
    !
    !== 概要
    !
    !  * Yamasaki, 1983 の湿度の観測データから基本場を作成する 
    !    * 相対湿度の観測データをNetCDFファイル化したものから値を読み込む
    !    * 読み込んだ値を線形補間して相対湿度の基本場を作成する
    !

    implicit none

    real(8) :: Observed_Humidity(1:100,1:88)        !NetCDF から読み込んだ相対湿度の観測値
                                                    !水平一様な2次元データ
    real(8),intent(out) :: kansoku_HumZ(1:88)       !鉛直1次元にした相対湿度の観測値の代入先
    real(8),intent(out) :: z_humid(DimZMin:DimZMax) !線形補間した相対湿度
    real(8),intent(out) :: kansoku_altitude(1:89)   !観測データに対応する高度
    integer :: k,m

    z_humid = 0.0d0
    kansoku_altitude = 0.0d0
    kansoku_HumZ = 0.0d0

    ! 観測データに対応する高度を与える
    ! これは線形補間の際に使用する
    ! なお, Yamsaki(1983)にある元データはデータ間隔がばらばらであったので
    ! 250m間隔に直したものを使用している
    do m = 1,89
     kansoku_altitude(m) = 25.0d1*(m-1)
    end do

    ! NetCDF ファイルから相対湿度を読み込む
    call HistoryGet('Ob_Humidity.nc', 'Ob_Humidity', array=Observed_Humidity)

    ! 相対湿度の観測値のファイル"Ob_Humidity.nc"は2次元データ(x,z)であるため
    ! 1次元データにする
    do m=1,88
      kansoku_HumZ(m)  = Observed_Humidity(1,m)
    end do

    ! kansoku_HumZ を線形補間してz_humidを作成する
    ! kansoku_altitude(m) < s_Z(k) < kansoku_altitude(m+1),
    ! kansoku_altitude = s_Z(k),
    ! kansoku_altitude < s_Z(k)
    ! の3つで場合分け
    do k = RegZMin+1,DimZMax
     do m = 1,87
      if (kansoku_altitude(m).ne.s_Z(k).and.s_Z(k).gt.kansoku_altitude(m) &
         & .and.s_Z(k).lt.kansoku_altitude(m+1)) then
        z_humid(k) = kansoku_HumZ(m)+((kansoku_HumZ(m+1) - kansoku_HumZ(m)) &
        & /(kansoku_altitude(m+1)-kansoku_altitude(m)))*(s_Z(k)-kansoku_altitude(m))

      else if (kansoku_altitude(m).eq.s_Z(k)) then
        z_humid(k) = kansoku_HumZ(m)

      else if (kansoku_altitude(m).lt.s_Z(k)) then
        z_humid(k) = z_humid(k-1)
      end if
     end do
    end do

  end subroutine HUM

!!------------------------------------------------------------------------------!!!
  subroutine ECCM_N1994( z_Temp, z_Press, za_MolFr )
    !
    !== 概要
    !  * Yamasaki, 1983 の温度と相対湿度の観測値で基本場を作成する
    !    * 温度の基本場
    !      * 観測データをNetCDFファイル化したものから値を読み込む
    !        * 読み込んだ値を線形補間して温度の基本場を作成する
    !    * 湿度の基本場
    !      * 相対湿度は subroutine HUM で作成済み
    !        * ここでは相対湿度をモル比に変換して湿度の基本場を作成する 
    !    * 気圧の基本場
    !      * 湿潤成分と乾燥成分の分子量差を考慮した静水圧平衡から計算する
    !

    implicit none

    real(8), intent(out):: z_Press(DimZMin:DimZMax)!圧力
    real(8), intent(out):: z_Temp(DimZMin:DimZMax) !温度
    real(8), intent(out):: za_MolFr(DimZMin:DimZMax, 1:SpcNum) !モル比
    real(8)             :: Ob_altitude(1:89)       !観測データに対応する高度
    real(8)             :: Ob_HumZ(1:88)           !相対湿度の観測値(1次元)
    real(8)             :: z_HumZ(DimZMin:DimZMax) !相対湿度の基本場の値(%)
    real(8)        :: Observed_Temp(1:100,1:89)!NetCDF から読み込んだ気温の観測値(2次元)
    real(8)        :: Ob_TempZ(1:89)           !NetCDF から読み込んだ気温の観測値(1次元)
    real(8)             :: z_PressSaturated(DimZMin:DimZMax)!Tetens(1930) の飽和蒸気圧
    integer             :: k, s, i, m, b

    z_HumZ = 0.0d0

    ! subroutine HUMから相対湿度の基本場を取得
    call HUM(z_HumZ, Ob_altitude, Ob_HumZ)

    ! 温度データをNetCDFから取得し線形補間する

    ! NetCDF ファイルからの気温の観測値の読み込み
    call HistoryGet('Ob_Temp.nc', 'Ob_Temp', array=Observed_Temp)

    ! 温度の観測値のファイル"Ob_Temp.nc"は2次元データ(x,z)なので
    ! 1次元データ(z)にする
    do m=1,89
      Ob_TempZ(m) = Observed_Temp(1,m)
    end do

    ! Ob_TempZ を線形補間して温度の基本場 z_Temp を作成する
    ! kansoku_altitude(m) < s_Z(k) < kansoku_altitude(m+1),
    ! kansoku_altitude = s_Z(k),
    ! kansoku_altitude < s_Z(k)
    ! の3つで場合分け

    ! 初期化
    z_Temp = 1.0d-60      

    do k = RegZMin+1,DimZMax
     do m = 1,88

      ! s_Z が Ob_altitude のある区間に挟まれるとき,
      ! その区間の Obaltitude(m), Obaltitude(m+1) を結ぶ
      ! 直線で Ob_TempZ を線形補間する
      if (Ob_altitude(m).ne.s_Z(k).and.s_Z(k).gt.Ob_altitude(m) &
          & .and.s_Z(k).lt.Ob_altitude(m+1)) then
        z_Temp(k) = Ob_TempZ(m)+((Ob_TempZ(m+1)-Ob_TempZ(m)) &
        & /(Ob_altitude(m+1)-Ob_altitude(m)))*(s_Z(k)-Ob_altitude(m))

      ! s_Z(k) = Obaltitude(m) の時は, Ob_TempZ(m) = z_Temp(k) である
      else if (Ob_altitude(m).eq.s_Z(k)) then
        z_Temp(k) = Ob_TempZ(m)

      ! s_Z(k) > Ob_altitude では観測データが無いので等温大気にする
      else if (Ob_altitude(m).lt.s_Z(k)) then
        z_Temp(k) = z_Temp(k-1)

      end if
     end do
    end do

    ! 気圧とモル比を計算する

    ! 初期化
    z_Press = 1.0d-60

    ! 地表面気圧と温度から大気最下層の気圧を求める
    ! 静水圧の式 dP/dz = - \rho * g を使用
    ! 簡単のため, \rho = 1 kg/m^3 の乾燥大気を仮定している
    z_Press(RegZMin+1) = PressSfc - (Grav * DelZ)

    ! 大気最下層のモル比を計算する
    za_MolFr(RegZMin+1,1) = SvapPress( 6,z_Temp(RegZMin+1)) * &
                           &(z_HumZ(RegZMin+1)/100)/z_Press(RegZMin+1)

    ! Tetens(1930) の飽和蒸気圧も使ってみたいので以下で計算
    do k=RegZMin+1,DimZMax
      z_PressSaturated(k) = 6.11*(10**(7.5*(z_Temp(k)-273)/(z_Temp(k)-35.7)))*100
    end do

    ! 大気最下層の値を元に, 気圧と湿度の基本場を求める
    ! 気圧は湿潤成分と乾燥成分の分子量差を考慮した静水圧平衡から計算する
    ! 相対湿度はモル比に変換する
    ! 以下のif文のDryHeightは基本場に乾燥した層を入れる場合に必要な値
    ! NAMELISTで鉛直領域以上の値を指定しておけば乾燥した層は入らない

    ! 湿度の基本場の計算
    do k = RegZMin+1,DimZMax-1
     if (s_Z(k+1).lt.DryHeight) then
        za_MolFr(k+1,1) = SvapPress( 6,z_Temp(k)) * (z_HumZ(k)/100)/z_Press(k)

    ! 気圧の基本場の計算
      z_Press(k+1) = &
         & z_Press(k)-(Grav*z_Press(k)*DelZ) &
         & /((GasRUniv/ &
         & (MolWtDry + za_MolFr(k,1) * (MolWtWet(1) - MolWtDry))) &
         & *z_Temp(k))

      !DryHeightより上に領域があるときは上で求めたモル比をゼロにして
      !乾燥大気における静水圧の式から気圧を求める
      else if (s_Z(k+1).ge.DryHeight) then
        za_MolFr(k+1,1) = 0.0d0

        z_Press(k+1) = &
           & z_Press(k)-(Grav*z_Press(k)*DelZ) &
           & /((GasRUniv/MolWtDry)*z_Temp(k))
      end if
    end do


 end subroutine ECCM_N1994

!!!------------------------------------------------------------------------------!!!
  subroutine HUM_Takemi2007(z_RH,Ztr,TR,LL)
    !
    !== 概要
    ! * deepconv の地球用のテスト計算としてTakemi(2007)の再現計算を
    !   するための湿度の基本場を作成する
    !   * Takemi(2007) の熱帯場と中緯度場の相対湿度の"BASE CASE"の
    !     高度 1.5km 以上の湿度を計算している(詳細はTakemi(2007)の図2を参照のこと)
    !     * Takemi(2007)では高度1.5km以上はある関数で与えているが,
    !       高度1.5km以下は一定の混合比で与えているため,
    !       高度1.5km以下の湿度は温度や気圧が求まらないと得られないので後で求める
    !
    implicit none

    real(8),intent(out) :: z_RH(DimZMin:DimZMax) !Relative Humidity(相対湿度)
    real(8),intent(out) :: Ztr                   !対流圏界面の高度
    integer,intent(out) :: TR                    !対流圏界面の高度の格子番号
    integer,intent(out) :: LL                    !高度 1.5 km 直下の格子番号 (Low Level)
    integer :: k

    z_RH = 0.0d0
    Ztr  = 0.0d0

    ! Ztrに対流圏界面高度の値を与える
    ! Takemi(2007)では12kmだが, 必ずしも格子点がちょうど12kmには乗らないため,
    ! 12kmか,12kmの直上の位置を対流圏界面にする
    do k = RegZMin+1,DimZMax

      ! s_Z が ちょうど 12 km なら 12 km が対流圏界面
      if (s_Z(k).eq.12000) then
        Ztr = s_Z(k)

      ! s_Z がちょうど12kmではなかった場合,12kmの直上のs_Zが対流圏界面Ztrになる
      else if (s_Z(k).gt.12000) then
        Ztr = minval( s_Z, 1, s_Z > 12000 )
      end if
    end do

    ! 対流圏界面高度の格子点を求める (Ztr = s_Z(TR))
    do k = RegZMin+1,DimZMax
      if (s_Z(k).ge.12000.and.s_Z(k).le.(12000 + DelZ)) then
        TR = k
      end if
    end do

    ! 高度1.5kmもしくはその直下の格子番号を求める
    do k = RegZMin+1,DimZMax
      if (s_Z(k).le.1500.and.s_Z(k).ge.(1500 - DelZ)) then
        LL = k
      end if
    end do

    ! Takemi(2007)の式を用いて相対湿度を求める
    ! とりあえず高度1.5km以下の湿度はゼロにしておく
    do k = RegZMin+1, DimZMax
      if (k.le.LL) then
        z_RH(k) = 0.0d0

      else if (k.le.TR.and.k.gt.LL) then
        z_RH(k) = 1.0d0 - 0.750d0 * (s_Z(k) / Ztr)**1.25

      else if (k.gt.TR) then
        z_RH(k) = 0.25
      end if
    end do
  end subroutine HUM_Takemi2007

!!!------------------------------------------------------------------------------!!!
  subroutine HUM_Takemi2007_TDRY(z_RH,Ztr,TR,LL)
    !
    !== 概要
    ! * deepconv の地球用のテスト計算としてTakemi(2007)の再現計算を
    !   するための湿度の基本場を作成する
    !   * Takemi(2007) の熱帯場と中緯度場の相対湿度の"DRY CASE"の
    !     高度 1.5km 以上の湿度を計算している(詳細はTakemi(2007)の図2を参照のこと)
    !     * Takemi(2007)では高度1.5km以上はある関数で与えているが,
    !       高度1.5km以下は一定の混合比で与えているため,
    !       高度1.5km以下の湿度は温度や気圧が求まらないと得られないので後で求める
    !     * 中緯度場とは乾燥のさせ方が違う
    !

    implicit none

    real(8),intent(out) :: z_RH(DimZMin:DimZMax) !Relative Humidity(相対湿度)
    real(8),intent(out) :: Ztr                   !対流圏界面の高度
    integer,intent(out) :: TR                    !対流圏界面の高度の格子番号
    integer,intent(out) :: LL                    !高度 1.5 km 直下の格子番号 (Low Level)
    integer :: k

    z_RH = 0.0d0
    Ztr  = 0.0d0

    ! Ztrに対流圏界面高度の値を与える
    ! Takemi(2007)では12kmだが, 必ずしも格子点がちょうど12kmには乗らないため,
    ! 12kmか,12kmの直上の位置を対流圏界面にする
    do k = RegZMin+1,DimZMax

      ! s_Z が ちょうど 12 km なら 12 km が対流圏界面
      if (s_Z(k).eq.12000) then
        Ztr = s_Z(k)

      ! s_Z がちょうど12kmではなかった場合,12kmの直上のs_Zが対流圏界面Ztrになる
      else if (s_Z(k).gt.12000) then
        Ztr = minval( s_Z, 1, s_Z > 12000 )
      end if
    end do

    ! 対流圏界面高度の格子点を求める (Ztr = s_Z(TR))
    do k = RegZMin+1,DimZMax
      if (s_Z(k).ge.12000.and.s_Z(k).le.(12000 + DelZ)) then
        TR = k
      end if
    end do

    ! 高度1.5kmもしくはその直下の格子番号を求める
    do k = RegZMin+1,DimZMax
      if (s_Z(k).le.1500.and.s_Z(k).ge.(1500 - DelZ)) then
        LL = k
      end if
    end do

    ! Takemi(2007)の式を用いて相対湿度を求める
    ! とりあえず高度1.5km以下の湿度はゼロにしておく
    ! さらに, 乾燥大気のエントレインメントを考慮した熱帯の"DRY CASE"では
    ! 高度 2.5 km, 5.0 km, 7.5 km のどこかで相対湿度を 20 % 減らすDRY1,2,3ケースを用意する
    do k = RegZMin+1, DimZMax
      if (k.le.LL) then
        z_RH(k) = 0.0d0

      ! 1.5 km から (2.5 km, 5.0 km, 7.5 km) まではそのまま変化

!     else if (s_Z(k).lt.2500.and.k.gt.LL) then   ! DRY1
!     else if (s_Z(k).lt.5000.and.k.gt.LL) then   ! DRY2
      else if (s_Z(k).lt.7500.and.k.gt.LL) then   ! DRY3
        z_RH(k) = 1.0d0 - 0.750d0 * (s_Z(k) / Ztr)**1.25

      ! (2.5 km, 5.0 km, 7.5 km)で20 % 湿度を減少させる 
!     else if (s_Z(k).ge.2500) then   ! DRY1
!     else if (s_Z(k).ge.5000) then   ! DRY2
      else if (s_Z(k).ge.7500) then   ! DRY3
        z_RH(k) = 1.0d0 - 0.750d0 * (s_Z(k) / Ztr)**1.25 - 0.20d0

      end if
    end do

    ! 上空で相対湿度が 25% 以下になる場合は 25% で固定する
    do k = RegZMin+1, DimZMax
      if (z_RH(k).le.(0.25).and.k.gt.LL) then
        z_RH(k) = 0.250d0
      end if
    end do

  end subroutine HUM_Takemi2007_TDRY

!!!------------------------------------------------------------------------------!!!
  subroutine HUM_Takemi2007_MDRY(z_RH,Ztr,TR,LL)
    !
    !== 概要
    ! * deepconv の地球用のテスト計算としてTakemi(2007)の再現計算を
    !   するための湿度の基本場を作成する
    !   * Takemi(2007) の熱帯場と中緯度場の相対湿度の"DRY CASE"の
    !     高度 1.5km 以上の湿度を計算している(詳細はTakemi(2007)の図2を参照のこと)
    !     * Takemi(2007)では高度1.5km以上はある関数で与えているが,
    !       高度1.5km以下は一定の混合比で与えているため,
    !       高度1.5km以下の湿度は温度や気圧が求まらないと得られないので後で求める
    !     * 熱帯場とは乾燥のさせ方が違う

    implicit none

    real(8),intent(out) :: z_RH(DimZMin:DimZMax) !Relative Humidity(相対湿度)
    real(8),intent(out) :: Ztr                   !対流圏界面の高度
    integer,intent(out) :: TR                    !対流圏界面の高度の格子番号
    integer,intent(out) :: LL                    !高度 1.5 km 直下の格子番号 (Low Level)
    integer :: k

    z_RH = 0.0d0
    Ztr  = 0.0d0

    ! Ztrに対流圏界面高度の値を与える
    ! Takemi(2007)では12kmだが, 必ずしも格子点がちょうど12kmには乗らないため,
    ! 12kmか,12kmの直上の位置を対流圏界面にする
    do k = RegZMin+1,DimZMax

      ! s_Z が ちょうど 12 km なら 12 km が対流圏界面
      if (s_Z(k).eq.12000) then
        Ztr = s_Z(k)

      ! s_Z がちょうど12kmではなかった場合,12kmの直上のs_Zが対流圏海面Ztrになる
      else if (s_Z(k).gt.12000) then
        Ztr = minval( s_Z, 1, s_Z > 12000 )
      end if
    end do

    ! 対流圏界面高度の格子点を求める (Ztr = s_Z(TR))
    do k = RegZMin+1,DimZMax
      if (s_Z(k).ge.12000.and.s_Z(k).le.(12000 + DelZ)) then
        TR = k
      end if
    end do

    ! 高度1.5kmもしくはその直下の格子番号を求める
    do k = RegZMin+1,DimZMax
      if (s_Z(k).le.1500.and.s_Z(k).ge.(1500 - DelZ)) then
        LL = k
      end if
    end do


    ! Takemi(2007)の式を用いて相対湿度を求める
    ! とりあえず高度1.5kmまでの湿度はゼロにしておく
    ! 乾燥大気のエントレインメントを考慮した中緯度場の"DRY CASE"では
    ! 高度 2.5 km で相対湿度を 13 % もしくは 30 % 減らす DRY1,2 ケースを用意する.
    do k = RegZMin+1, DimZMax
      if (k.le.LL) then
        z_RH(k) = 0.0d0

      ! 1.5 km から2.5 km まではそのまま変化
      else if (s_Z(k).lt.2500.and.k.gt.LL) then
        z_RH(k) = 1.0d0 - 0.750d0 * (s_Z(k) / Ztr)**1.25

      ! 2.5 km 以降は 13 % or 30 % 減らす
      else if (s_Z(k).ge.2500) then
!       z_RH(k) = 1.0d0 - 0.750d0 * (s_Z(k) / Ztr)**1.25 - 0.13d0   ! DRY1
        z_RH(k) = 1.0d0 - 0.750d0 * (s_Z(k) / Ztr)**1.25 - 0.30d0   ! DRY2

      end if
    end do

    ! 上空で相対湿度が 25% 以下になる場合は 25% で固定する
    do k = RegZMin+1, DimZMax
      if (z_RH(k).le.(0.25).and.k.gt.LL) then
        z_RH(k) = 0.250d0
      end if
    end do

  end subroutine HUM_Takemi2007_MDRY

!!!------------------------------------------------------------------------------!!!
  subroutine ECCM_Takemi2007( z_Temp, z_Press, za_MolFr, z_LLhumidity )
    !
    !== 概要
    ! * deepconv の地球用のテスト計算としてTakemi(2007)の再現計算を
    !   するための湿度の基本場を作成する
    !   * 基本場の温度の式が温位で与えられているため, 温度に変換する必要がある
    !

    implicit none

    real(8), intent(out):: z_Press(DimZMin:DimZMax)!圧力
    real(8), intent(out):: z_Temp(DimZMin:DimZMax) !温度
    real(8), intent(out):: za_MolFr(DimZMin:DimZMax, 1:SpcNum) !モル比
    real(8), intent(out) :: z_LLhumidity(DimZMin:DimZMax)  ! 1.5 km までの相対湿度
    !
    ! ECCM_Takemiでのみ必要な変数
    !
    real(8) :: z_TakemiRH(DimZMin:DimZMax)         !高度 1.5km 以上の相対湿度
    real(8) :: z_TakemiTheta(DimZMin:DimZMax)      !基本場の温位
    real(8) :: Theta_0                             ! 地表面の温位
    real(8) :: Theta_tr                            ! 対流圏界面の温位
    real(8) :: TakemiZtr                           ! 対流圏界面の高度
    real(8) :: HumSfc                              ! 地表面の相対湿度
    real(8) :: a_MolFrSfc                          ! 地表面のモル比
!!--------------------------------------------------
    integer             :: k, s, i, m, b, T, L

    z_TakemiRH = 0.0d0
    z_TakemiTheta = 0.0d0

    ! 湿度場の選択
    call HUM_Takemi2007(z_TakemiRH,TakemiZtr,T,L)
!   call HUM_Takemi2007_TDRY(z_TakemiRH,TakemiZtr,T,L)
!   call HUM_Takemi2007_MDRY(z_TakemiRH,TakemiZtr,T,L)

    ! 温度・湿度・気圧を求める

    !地表面の温位は300K
    Theta_0 = 300.0d0

    ! 中緯度場と熱帯場の選択
    ! 対流圏界面の温位で区別している
!   Theta_tr = 343.0d0     ! 中緯度
    Theta_tr = 358.0d0     ! 熱帯

    ! 大気最下層の温位
    z_TakemiTheta(RegZMin+1) = Theta_0 + &
                           & (Theta_tr - Theta_0)*(s_Z(RegZMin+1) / TakemiZtr)**1.25

    z_Temp  = 1.0d-60
    z_Press = 1.0d-60
    z_LLhumidity = 0.0d0

    ! 大気最下層の気温・気圧・水蒸気混合比・分子量の計算

    ! 地表面のモル比
    ! 10g/kg〜18g/kgの混合比で与えられているためモル比に直す
    ! なお熱帯場はQ18のみ使用
    a_MolFrSfc = 0.0180d0 * MolWtDry / MolWtWet(1)  ! Q18
!   a_MolFrSfc = 0.0160d0 * MolWtDry / MolWtWet(1)  ! Q16
!   a_MolFrSfc = 0.0140d0 * MolWtDry / MolWtWet(1)  ! Q14
!   a_MolFrSfc = 0.0120d0 * MolWtDry / MolWtWet(1)  ! Q12
!   a_MolFrSfc = 0.0100d0 * MolWtDry / MolWtWet(1)  ! Q10 

    ! 地表面の相対湿度
    ! モル比から逆算. 気圧を求めるのに必要
    HumSfc = a_MolFrSfc * PressSfc / SvapPress(6,TempSfc)

    ! 大気最下層の気圧
    z_Press(RegZMin+1) = &
           & PressSfc - (Grav * PressSfc * DelZ) &
           & /((GasRUniv/ &
           & (MolWtDry + a_MolFrSfc * (MolWtWet(1) - MolWtDry))) &
           & *TempSfc)

    ! 大気最下層の温度.　温位を元に計算
    ! 地表面気圧は 1000hPa と仮定する
    z_Temp(RegZMin+1) = z_TakemiTheta(RegZMin+1) &
                        &* (z_Press(RegZMin+1) / 100000.0d0)**(GasRDry / CpDry)

    ! 大気最下層のモル比
    ! 熱帯場はQ18のみ使用
    za_MolFr(RegZMin+1,1) = 0.0180d0 * MolWtDry / MolWtWet(1)  ! Q18
!   za_MolFr(RegZMin+1,1) = 0.0160d0 * MolWtDry / MolWtWet(1)  ! Q16
!   za_MolFr(RegZMin+1,1) = 0.0140d0 * MolWtDry / MolWtWet(1)  ! Q14
!   za_MolFr(RegZMin+1,1) = 0.0120d0 * MolWtDry / MolWtWet(1)  ! Q12
!   za_MolFr(RegZMin+1,1) = 0.0100d0 * MolWtDry / MolWtWet(1)  ! Q10
 
    ! 大気最下層の相対湿度
    ! モル比から逆算. 気圧を求めるのに必要
    z_LLhumidity(RegZMin+1) = za_MolFr(RegZMin+1,1) &
                              &* z_Press(RegZMin+1) / SvapPress(6,z_Temp(RegZMin+1))

    ! 大気最下層の種々の物理量をもとに温度・湿度気圧の基本場を作成する
    ! 高度 0 - 1.5km と 1.5km - で湿度の与え方が異なり,
    ! さらに, 対流圏界面以上では温度の式が異なるので, 以下の 3パターンに分けて
    ! 基本場を求める
    ! 高度 0 - 1.5km, 1.5km - 対流圏界面(Ztr), Ztr以上 
    do k = RegZMin+1, DimZMax

    !
    ! 基本場(I)
    ! 大気最下層から高度 1.5 km 以下(k が L-1 以下)の基本場
    !
     if (k.le.L-1) then

      ! モル比
      ! 高度 1.5km まではモル比が任意の値で固定されている
      ! ただし, 熱帯場は 18 g/kg のみ
      za_MolFr(k+1,1) = 0.0180d0 * MolWtDry / MolWtWet(1)
!     za_MolFr(k+1,1) = 0.0160d0 * MolWtDry / MolWtWet(1)
!     za_MolFr(k+1,1) = 0.0140d0 * MolWtDry / MolWtWet(1)
!     za_MolFr(k+1,1) = 0.0120d0 * MolWtDry / MolWtWet(1)
!     za_MolFr(k+1,1) = 0.0100d0 * MolWtDry / MolWtWet(1)
 
      ! 気圧
      z_Press(k+1) = &
           & z_Press(k)-(Grav*z_Press(k)*DelZ) &
           & /((GasRUniv/ &
           & (MolWtDry + za_MolFr(k,1) * (MolWtWet(1) - MolWtDry))) &
           & *z_Temp(k))

      ! 温位から温度を計算
      z_TakemiTheta(k+1) = Theta_0 + (Theta_tr - Theta_0)*(s_Z(k+1) / TakemiZtr)**1.25
      z_Temp(k+1) = z_TakemiTheta(k+1) * (z_Press(k+1) / 100000.0d0)**(GasRDry / CpDry)

      ! モル比から相対湿度を計算
      z_LLhumidity(k+1) = za_MolFr(k+1,1) * z_Press(k+1) / SvapPress(6,z_Temp(k+1))

      !
      !== Takemi(2007)の基本場の使用上の注意
      !  * 低層で混合比を固定し続けると高度 1.5km までに相対湿度が100%を超える高度がある
      !  * 100% を超えてしまった高度の湿度は, その直下の高度の湿度にすると
      !    Takemi(2007)と同じ湿度の鉛直プロファイルが描ける
      !    * なぞ
      !  * MQ16 は 湿度 99.・・・% という値が出て, 100% を超えなくてもほぼ100%
      !    な値が出てしまう
      !    * Takemi(2007)ではこのような値があるようには見えない
      !      しょうがないのでMQ18 とMQ14 の 100% を超える直前の値の間を取った値を入れておく
      !      * 0.94511507230d0という怪しい数字が上記の値
      if (z_LLhumidity(k+1).gt.1) then
        z_LLhumidity(k+1) = z_LLhumidity(k)
!     if (z_LLhumidity(k+1).gt.0.99) then      !MQ16 用
!       z_LLhumidity(k+1) = 0.94511507230d0    !MQ16 用
     end if

     !
     ! 基本場(II)
     ! 高度 1.5km より上で, 対流圏界面以下(L < k < T)の基本場
     !
     else if (k.le.T.and.k.gt.L) then

      ! k = L での相対湿度
      ! この if 文の初めの k は L+1 から始まり z_TakemiRH(L) が区間外なので値を与える
      z_TakemiRH(L) = z_LLhumidity(L)

      ! 気圧の基本場
      z_Press(k) = &
           & z_Press(k-1)-(Grav*z_Press(k-1)*DelZ) &
           & /((GasRUniv/ &
           & (MolWtDry + za_MolFr(k+1,1) * (MolWtWet(1) - MolWtDry))) &
           & *z_Temp(k-1))

       ! 温位から温度を計算
       z_TakemiTheta(k) = Theta_0 + (Theta_tr - Theta_0)*(s_Z(k) / TakemiZtr)**1.25
       z_Temp(k) = z_TakemiTheta(k) * (z_Press(k) / 100000.0d0)**(GasRDry / CpDry)

       ! モル比と大気の平均分子量
       za_MolFr(k,1) = SvapPress( 6,z_Temp(k))* z_TakemiRH(k)/z_Press(k)
     
       !
       ! 基本場(III)
       ! 対流圏界面より上 (k > T) の場
       !
       else if (k.gt.T) then

       ! 気圧の基本場
       z_Press(k) = &
           & z_Press(k-1)-(Grav*z_Press(k-1)*DelZ) &
           & /((GasRUniv/ &
           & (MolWtDry + za_MolFr(k,1) * (MolWtWet(1) - MolWtDry))) &
           & *z_Temp(k-1))

       ! 温位から温度を計算
       z_TakemiTheta(k) = Theta_tr * exp(Grav * (s_Z(k) - TakemiZtr) / (CpDry * z_Temp(T)))
       z_Temp(k) = z_TakemiTheta(k) * (z_Press(k) / 100000.0d0)**(GasRDry / CpDry)

       ! モル比と大気の平均分子量
       za_MolFr(k,1) = SvapPress( 6,z_Temp(k))* z_TakemiRH(k)/z_Press(k)

     end if
    end do

 end subroutine ECCM_Takemi2007

!!!------------------------------------------------------------------------------!!!
  subroutine HUM_Takemi2007ALL(z_RHall)
    !
    !== 概要
    ! Takemi(2007) での湿度場の与え方の都合上, 高度1.5km以下と高度1.5km以上で
    ! 相対湿度を別々に求めた
    ! 前者は HUM_Takemi2007 or HUM_Takemi2007_TDRY or HUM_Takemi2007_MDRY で求め
    ! 後者は ECCM_Takemi2007 で求めている.
    ! そのため, 全高度での相対湿度 z_RHall を確認のためにも求めるために,
    ! このサブルーチンで両者を結合している
    ! 値の確認が主な目的なため, 基本的に必要ないが, 仮にECCM_MolFrが必要なときに
    ! 相対湿度が必要なこともあり, 一応求めている

    implicit none

    real(8),intent(out) :: z_RHall(DimZMin:DimZMax)
    real(8) :: a_MolFrIn(1:SpcNum)
    real(8) :: Humidit
    real(8) :: z_Tem(DimZMin:DimZMax)
    real(8) :: z_Pre(DimZMin:DimZMax)
    real(8) :: z_MolMea(DimZMin:DimZMax)
    real(8) :: za_MolF(DimZMin:DimZMax)
    real(8) :: ZTR
    real(8) :: z_humidity1(DimZMin:DimZMax) !高度0 - 1.5 km の相対湿度
    real(8) :: z_humidity2(DimZMin:DimZMax) !高度1.5 km 以上の相対湿度
    integer :: k, T, L

    !高度1.5km以下の湿度を取得
    call ECCM_Takemi2007(z_Tem, z_Pre, za_MolF, z_humidity1)

    !高度1.5km以上の湿度を取得
    !三つから選択
    call HUM_Takemi2007(z_humidity2, ZTR, T, L)
!   call HUM_Takemi2007_TDRY(z_humidity2, ZTR, T, L)
!   call HUM_Takemi2007_MDRY(z_humidity2, ZTR, T, L)

    do k = RegZMin+1, DimZMax
      if (k.le.L) then
        z_RHall(k) = z_humidity1(k)

      else if (k.gt.L) then
        z_RHall(k) = z_humidity2(k)

      end if
    end do

  end subroutine HUM_Takemi2007ALL

!!!------------------------------------------------------------------------------!!!
  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_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(inout) :: MolFr(1:SpcNum)    !モル分率
    real(8), intent(out):: DTempDZ
    real(8)             :: MolFrOrig(1:SpcNum) 
    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
    MolFrOrig    = MolFr

    !------------------------------------------------------------
    !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), MolFrOrig(1:SpcNum)) &
!    A = dot_product( Heat(1:SpcNum), MolFr(1:SpcNum)) &
      & / ( GasRUniv * Temp )
    B = dot_product(( Heat(1:SpcNum) ** 2.0d0), MolFrOrig(1:SpcNum)) &
!    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( 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 constants,only: MolWtDry,      &!
      &                 CpDry,         &!
      &                 Grav
    use composition,only: MolWtWet
    use basicset, only: 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
  
end module ECCM
