!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2005. All rights reserved.
!---------------------------------------------------------------------

!= Subroutine ECCM
!
!   * Developer: SUGIYAMA Ko-ichiro 
!   * Version: $Id: eccm.f90,v 1.1.2.1 2006/01/13 06:10:20 sugiyama Exp $
!   * Tag Name: $Name: arare3j $
!   * Change History: 
!
!== Overview 
!
!断熱的に上昇する気塊の温度減率を計算し, 静水圧平衡から圧力を決める
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!== Future Plans
!

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

  implicit none

  real(8), intent(in)      :: MolFrIni
  real(8), intent(out)     :: Temp( DimZMin:DimZMax )
  real(8), intent(out)     :: Pres( DimZMin:DimZMax )
!  real(8), intent(out)     :: MixRat( DimZMin:DimZMax )
  real(8), intent(out)     :: MolFrV( DimZMin:DimZMax )
  real(8), intent(out)     :: LatentHeat( DimZMin:DimZMax )

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


!!!
!!!設定する必要がない変数
!!!
  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)                  :: A, B
  
  integer                  :: i
  integer, parameter       :: d = 1
  integer, parameter       :: w = 2
  

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

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

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

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

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

!!!
!!!断熱減率 dT/dz の計算. 
!!!

  do i = RegZMin, RegZMin
     
     !前のステップでのモル分率を用いて現在の蒸気圧を計算
!     VapPres = MolFr(i-1, w) * PressSfc
     VapPres = PressSfc
     
     !飽和蒸気圧
     call antoine( TempSfc, 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( TempSfc, svap ) 

        !潜熱
        call PhaseChangeEntalpy(TempSfc, LatentHeat(i))
        
        !飽和蒸気圧と圧力から現在のモル数を計算
!        MolFr(i,w) = svap / PressSfc
        MolFr(i,w) = 1.0d0   ! 凝結成分が主成分の場合
        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) )
        
        !係数. 温度 TempSfc で評価
        A = LatentHeat(i) * MolFr(i,w) / ( GasRUniv * TempSfc )
        B = ( LatentHeat(i) ** 2.0d0 ) * MolFr(i,w)  &
             & / ( CpMean(i) * GasRUniv * ( TempSfc ** 2.0d0 ) )

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

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

     !圧力を静水圧平衡から計算
     Pres(i+1) = &
       & PressSfc  - ( PressSfc * Grav * DelZ * 0.5d0 * MolWtDry ) & 
       & / ( GasRUniv * 0.5d0 * (TempSfc + Temp(i+1)))

     !モル比 (DimZMax でも値を持つようにするため)
     MolFr(i+1,:) = MolFr(i,:)
     
  end do

  do i = RegZMin+1, DimZMax-1
     
     !前のステップでのモル分率を用いて現在の蒸気圧を計算
!     VapPres = MolFr(i-1, w) * Pres(i)
     VapPres = 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(i))
        
        !飽和蒸気圧と圧力から現在のモル数を計算
!        MolFr(i,w) = svap / Pres(i)
        MolFr(i,w) = 1.0d0   ! 凝結成分が主成分の場合
        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(i) * MolFr(i,w) / ( GasRUniv * Temp(i) )
        B = ( LatentHeat(i) ** 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 * 0.5d0 * (Temp(i) + Temp(i+1)))

     !モル比 (DimZMax でも値を持つようにするため)
     MolFr(i+1,:) = MolFr(i,:)
     
  end do
  
!  do i = RegZMin - 1, DimZMin, - 1 
  do i = RegZMin, 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

  !=== 出力用
  MolFrV = MolFr(:,w)
  LatentHeat(:) = LatentHeat(:) / MolWt(w)
  

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_log = (Antoine_A - (Antoine_B / (Antoine_C + Temp) ) )
    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 
    svap_log_dt =  antoine_B / (antoine_C + temp)**2.0D0 
    
    LatentHeat = svap_log_dt * GasRUniv * ( Temp ** 2.0d0 ) 
 
  end subroutine PhaseChangeEntalpy
  

end subroutine eccm
