!= Module ECCM
!
! Authors::   SUGIYAMA Koichiro, ODAKA Masatsugu
! Version::   $Id: eccm.f90,v 1.21 2007/04/19 10:46:03 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
!

subroutine ECCM2( a_MolFrIni, Humidity, z_Temp, z_Press, z_MolWtMean, za_MolFr )

  !モジュール読み込み
  use dc_message, only: MessageNotify

  use gridset,  only: DimZMin,       &! 配列の Z 方向の下限
    &                 DimZMax,       &! 配列の Z 方向の上限 
    &                 RegZMin,       &!
    &                 SpcNum,        &!
    &                 DelZ            !
  use basicset,only:  MolWtDry,      &!
    &                 MolWtWet,      &!
    &                 CpDryMol,      &!
    &                 SpcWetID,      &!
    &                 TempSfc,       &!
    &                 PressSfc,      &!
    &                 Grav            !
  use chemcalc, only: SvapPress       !,     &!
!    &                 LatentHeatPerMol, &!
!    &                 ReactHeatNH4SHPerMol
  use moistset, only: CondNum,          &!凝結過程の数
    &                 IdxCG,            &!凝結過程(蒸気)の配列添え字
    &                 IdxCC,            &!凝結過程(雲)の配列添え字
    &                 GasNum,           &!気体の数
    &                 IdxNH3,           &!NH3(蒸気)の配列添え字
    &                 IdxH2S             !H2S(蒸気)の配列添え字
  use ChemData, only: GasRUniv        
  use MoistFunc,only: DelMolFrNH4SH 

  !暗黙の型宣言禁止
  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), intent(out):: za_MolFr(DimZMin:DimZMax, 1:SpcNum) 
  !モル分率
!  real(8)             :: a_MolWt(0:SpcNum)       !分子量の作業配列
!  real(8)             :: a_MolFr(0:SpcNum)       !モル比の作業配列
  real(8)             :: SatPress                !飽和蒸気圧
  real(8)             :: VapPress                !蒸気圧
  real(8)             :: DelMolFr
  integer             :: k, s
    
  !-------------------------------------------------------------
  ! 配列の初期化
  !-------------------------------------------------------------
  !分子量の配列を用意
!  a_MolWt(0)        = MolWtDry
!  a_MolWt(1:SpcNum) = MolWtWet(1:SpcNum)
  
  write(*,*) a_MolFrIni(1:SpcNum) 
  write(*,*) Humidity 

  !初期モル比
  za_MolFr = 0.0d0
  za_MolFr(RegZMin, 1:SpcNum)   = a_MolFrIni(1:SpcNum) 
  
  !初期の平均分子量
!  a_MolFr(0)           = 1.0d0 - sum(a_MolFrIni)
!  a_MolFr(1:SpcNum)    = a_MolFrIni
  z_MolWtMean(RegZMin) = dot_product(MolWtWet, a_MolFrIni) + (1.0d0 - sum(a_MolFrIni)) * MolWtDry
    
  !地表面での温度(RegZMin は, 高度 DelZ / 2 に相当)
  z_Temp          = 1.0d-60
  z_Temp(RegZMin) = TempSfc - Grav * z_MolWtMean(RegZMin) &
    &               / CpDryMol * ( DelZ * 5.0d-1 )
  
  !地表面での圧力(RegZMin は, 高度 DelZ / 2 に相当)
  z_Press           = 1.0d-60
  z_Press(RegZMin)  = &
    & PressSfc *((TempSfc / z_Temp(RegZMin)) ** (- CpDryMol /  GasRUniv))

  !-----------------------------------------------------------
  ! (1) 乾燥断熱線に沿った温度を決める
  ! (2) 静水圧平衡から圧力を求める
  ! (3) (1),(2) の温度圧力に対して, とある相対湿度となるモル比を決める
  !-----------------------------------------------------------    
  DtDz: do k = RegZMin, DimZMax-1
    
    !(1)乾燥断熱線に沿って k+1 での温度を計算
    z_Temp(k+1) = z_Temp(k) - Grav * z_MolWtMean(k) / 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
    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
    
    !------------------------------------------------------------
    !温度勾配を計算
    !------------------------------------------------------------
    !モル比
!    a_MolFr(0) = 1.0d0 - sum( za_MolFr(k+1, 1:SpcNum) )
!    a_MolFr(1:SpcNum) = za_MolFr(k+1, 1:SpcNum)
    
    !平均分子量
!    z_MolWtMean(k+1) = dot_product(a_MolWt(0:SpcNum), a_MolFr(0:SpcNum))
    z_MolWtMean(k+1) = dot_product(MolWtWet, za_MolFr(k+1,:)) + (1.0d0 - sum(za_MolFr(k+1,:))) * MolWtDry    
    
  end do DtDz
  
end subroutine ECCM2
