!c Description: 
!c   Atreya and Romani (1985) によって温度, 組成の鉛直分布を計算し, 
!c   それを用いて温度減率, 仮温度減率, 静的安定度を計算するプログラム
!c
!c Current Code Owner: 
!c   sugiyama@gfd-dennou.org
!c
!c Copyright (C) SUGIYAMA Ko-ichiro, 2004, All rights reserved

subroutine eccm( MolFrIni, HumVaper, Temp, Pres, MixRat )
  use gridset, only: DimZMin, DimZMax, RegZMin, DelZ
  use physset, only: MolWtDry, MolWtWet, CpDryMol, CpWetMol, &
       &             GasRUniv, TempSfc, PressSfc, Grav

  implicit none

  real(8), intent(in)      :: MolFrIni
  real(8), intent(in)      :: HumVaper
  real(8), intent(out)     :: Temp( DimZMin:DimZMax )
  real(8), intent(out)     :: Pres( DimZMin:DimZMax )
  real(8), intent(out)     :: MixRat( DimZMin:DimZMax )


!!!
!!!設定する必要のある変数
!!!
  ! H2O(g,l) の値
  real(8), parameter       :: Antoine_A = 7.9186968d0
  real(8), parameter       :: Antoine_B = 1636.909d0
  real(8), parameter       :: Antoine_C = 224.92d0


!!!
!!!設定する必要がない変数
!!!
  real(8)                  :: MolFr( DimZMin:DimZMax, 1:2 ) !モル分率
  real(8)                  :: Cp( 1:2 )                     !比熱   (J/K mol)
  real(8)                  :: CpMean( DimZMin:DimZMax )     !比熱   (J/K mol)
  real(8)                  :: MolWtMean( DimZMin:DimZMax )  !分子量 
  real(8)                  :: MolWt( 1:2 )                  !分子量 
  
  real(8)                  :: VapPres           !分圧
  real(8)                  :: SatVapPres        !飽和蒸気圧
  
  real(8)                  :: dtdz
  real(8)                  :: svap
  real(8)                  :: LatentHeat
  real(8)                  :: A, B
  
  integer                  :: i
  integer, parameter       :: d = 1
  integer, parameter       :: w = 2
  

!!!
!!!配列の初期化
!!!
  Temp = TempSfc
  MolFr = 0.0d0
  CpMean = 0.0d0 

  !比熱
  Cp(d) = CpDryMol
  Cp(w) = CpWetMol

  !分子量
  MolWt(d) = MolWtDry
  MolWt(w) = MolWtWet

  !初期モル比
  MolFr(:,w) = MolFrIni 
  MolFr(:,d) = 1.0d0 - MolFr(:,w) 

  !地表面での温度・圧力
  Pres = 0.0d0
  Temp = 0.0d0
  Pres( RegZMin ) = PressSfc
  Temp( RegZMin ) = TempSfc

!!!
!!!断熱減率 dT/dp の計算. 
!!!
  do i = RegZMin, DimZMax - 1 
     
     !前のステップでのモル分率を用いて現在の蒸気圧を計算
     VapPres = MolFr(i-1, w) * Pres(i)
     
     !飽和蒸気圧
     call antoine( Temp(i), svap = SatVapPres )
     
     !飽和蒸気圧から凝結の有無を決める
     if ( VapPres < SatVapPres ) then 
        
        !モル比
        MolFr(i,:) = MolFr(i-1,:)
        
        !平均分子量の計算 (モル数は前のステップと変わらない)
        MolWtMean(i) = dot_product( MolFr(i, d:w), MolWt(d:w) )
        
        !平均比熱の計算 (モル数は前のステップと変わらない)
        CpMean(i) = dot_product( MolFr(i, d:w), Cp(d:w) )
        
        !乾燥断熱減率を計算
        dtdz = - MolWtMean(i) * Grav / CpMean(i)
        
     else
        
        !飽和蒸気圧
        call Antoine( Temp(i), svap ) 

        !潜熱
        call PhaseChangeEntalpy(Temp(i), LatentHeat)
        
        !飽和蒸気圧と圧力から現在のモル数を計算
        MolFr(i,w) = svap / Pres(i)
        MolFr(i,d) = 1.0d0 - MolFr(i,w)
        
        !平均比熱の計算
        CpMean(i)    = dot_product( MolFr(i,d:w), Cp(d:w) )
        
        !平均分子量の計算
        MolWtMean(i) = dot_product( MolFr(i,d:w), MolWt(d:w) )
        
        !係数. 温度 Temp(i) で評価
        A = LatentHeat * MolFr(i,w) / ( GasRUniv * Temp(i) )
        B = ( LatentHeat ** 2.0d0 ) * MolFr(i,w)  &
             & / ( CpMean(i) * GasRUniv * ( Temp(i) ** 2.0d0 ) )

        !断熱温度減率
        dtdz = - MolWtMean(i) * Grav * ( 1.0d0 + A ) / ( CpMean(i) * ( 1.0d0 + B ) )
        
     end if

     !温度を計算
     Temp(i+1) = Temp(i) + dtdz * DelZ

     !圧力を静水圧平衡から計算
     Pres(i+1) = &
          & Pres(i)  - ( Pres(i) * Grav * DelZ * MolWtDry ) / ( GasRUniv * Temp(i) )
     
  end do
  
  do i = RegZMin - 1, DimZMin, - 1 
     
     !モル比
     MolFr(i,:) = MolFr(i+1,:)
     
     !平均分子量の計算 (モル数は前のステップと変わらない)
     MolWtMean(i) = dot_product( MolFr(i+1, d:w), MolWt(d:w) )
     
     !平均比熱の計算 (モル数は前のステップと変わらない)
     CpMean(i) = dot_product( MolFr(i+1, d:w), Cp(d:w) )
     
     !乾燥断熱減率を計算
     dtdz = MolWtMean(i+1) * Grav / CpMean(i)
     
     !温度を計算
     Temp(i) = Temp(i+1) + dtdz * DelZ
     
     !圧力を静水圧平衡から計算
     Pres(i) = &
          & Pres(i+1)  + ( Pres(i+1) * Grav * DelZ * MolWtDry ) / ( GasRUniv * Temp(i+1) )
     
  end do
  

  !--- モル分率から混合比への変換
!  MixRat = &
!       &   ( MolWtWet * MolFr(:,w) * HumVaper ) &
!       & / (   MolWtDry * ( 1.0d0 - MolFr(:,w) * HumVaper ) &
!       &     + MolWtWet * ( MolFr(:,w) * HumVaper ) )
  MixRat = 0.622d0 * MolFr(:,w)
  MixRat( DimZMax ) =   MixRat( DimZMax - 1 )


contains

  
  subroutine Antoine(temp, svap)    
    implicit none
    real(8), intent(in)    :: temp         !温度
    real(8), intent(out)   :: svap         !飽和蒸気圧
    real(8)                :: svap_log
  
    !Antoine の式
    svap_log = (Antoine_A - (Antoine_B / (Antoine_C + Temp - 273.15D0) ) ) &
         & * dlog(10.0D0) + dlog(133.322D0)
    svap = dexp(svap_log)
    
  end subroutine Antoine


  subroutine PhaseChangeEntalpy(temp, LatentHeat)    
    implicit none
    real(8), intent(in)    :: temp         !温度
    real(8), intent(out)   :: LatentHeat   !潜熱
    real(8)                :: svap_log_dt
    
    svap_log_dt =  antoine_B * dlog(10.0D0) &
         & / (antoine_C + temp - 273.15D0)**2.0D0 
    
    LatentHeat = svap_log_dt * GasRUniv * ( Temp ** 2.0d0 ) 
 
  end subroutine PhaseChangeEntalpy
  

end subroutine eccm
