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

!= Subroutine ECCM2
!
!   * Developer: SUGIYAMA Ko-ichiro 
!   * Version: $Id: eccm.f90,v 1.1.2.1 2006/01/13 06:10:20 sugiyama Exp $
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview 
!
!断熱的に上昇する気塊の温度減率を計算し, 静水圧平衡から圧力を決める
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!  * 比熱は定数, 平均分子量は配列
!  * 積分は 1 次精度であるが, 十分にステップ数を細かくすれば OK. 
!    * 水のみの場合でルンゲクッタを用いた場合と値は同じ. 
!
!== Future Plans
!

subroutine eccm2( MolFrIni)
  use gridset, only: RegZMin, RegZMax, SpcNum, DelZ, f_Z
  use dataset, only: MolWtDry, MolWtWet, &
    &                CpDryMol, CpWetMol,SpcWetID, &
    &                TempSfc, PressSfc, Grav
  use chemdata, only: GasRUniv
  use chemcalc, only: SvapPress, LatentHeat, ChemPotOneSpc, ChemPotNH4SH, &
    &                 EntropyOneSpc, EntropyNH4SH  
  
  implicit none

  real(8), intent(in) :: MolFrIni(1:SpcNum)
  real(8)             :: Temp( RegZMin:RegZMax+1 )
  real(8)             :: Press( RegZMin:RegZMax+1 )
  real(8)             :: DTempDz( RegZMin:RegZMax )
  real(8)             :: DTempDzDry( RegZMin:RegZMax )
  real(8)             :: MolFr( RegZMin-1:RegZMax+1, 0:SpcNum ) !モル分率
  real(8)             :: Cp(0:SpcNum) 
  real(8)             :: CpMean( RegZMin-1:RegZMax+1) 
  real(8)             :: MolWt( 0:SpcNum) 
  real(8)             :: MolWtMean( RegZMin-1:RegZMax+1) 
  real(8)             :: dM( RegZMin:RegZMax) 
  real(8)             :: dp( RegZMin:RegZMax) 
  real(8)             :: dMdz( RegZMin:RegZMax) 
  real(8)             :: Stab( RegZMin:RegZMax) 
  real(8)             :: LatentHeatAll(SpcNum )  
  real(8)             :: SvapPressAll(SpcNum)        !飽和蒸気圧
  real(8)             :: VapPress                   !分圧
  real(8)             :: A, B
  real(8)             :: VapPressNH3
  real(8)             :: VapPressH2S
  real(8)             :: MolFrNH3
  real(8)             :: MolFrH2S
  real(8)             :: num
  real(8)             :: del, Q
  integer             :: i, s
  integer, parameter  :: ID_NH3   = 7
  integer, parameter  :: ID_H2S   = 9
  integer, parameter  :: ID_NH4SH = 9
  integer, parameter  :: Num_NH3   = 2
  integer, parameter  :: Num_H2S   = 3

  
!  write(*,*) 'Sub ECCM2, MolFrIni: ', MolFrIni

!!!
!!!配列の初期化
!!!
  !初期モル比
  MolFr = 0.0d0
  MolFr(RegZMin, 1:SpcNum)   = MolFrIni(1:SpcNum) 
  MolFr(RegZMin, 0)          = 1.0d0 - sum(MolFrIni)
  MolFr(RegZMin-1, 1:SpcNum) = MolFrIni(1:SpcNum) 
  MolFr(RegZMin-1, 0)        = 1.0d0 - sum(MolFrIni)

  !地表面での温度・圧力
  Press = 0.0d0
  Temp = 0.0d0
  Press(RegZMin) = PressSfc
  Temp(RegZMin)  = TempSfc
  Cp(0)        = CpDryMol
  Cp(1:SpcNum) = CpWetMol
!  CpMean       = 0.0d0
  MolWt(0)        = MolWtDry
  MolWt(1:SpcNum) = MolWtWet
  MolWtMean       = 0.0d0
  num = 0.0d0


!!!
!!!断熱減率 dT/dz の計算. 
!!!
  do i = RegZMin, RegZMax+1
    write(*,*) i, Temp(i), Press(i)

    !初期化
    Q = 0.0d0

    !飽和蒸気圧
    SvapPressAll(1:SpcNum) = SvapPress( Temp(i) )
    
    !潜熱    
    LatentHeatAll(1:SpcNum)  = LatentHeat( Temp(i) )

    !------------------------------------------------------------
    !NH4SH 以外の化学種の平衡条件
    !------------------------------------------------------------
    do s = 1, SpcNum      
      if (s /= Num_H2S ) then 
        !前のステップでのモル分率を用いて現在の蒸気圧を計算
        VapPress = MolFr(i-1, s) * Press(i)     

        !飽和蒸気圧から凝結の有無を決める
        if ( VapPress < SvapPressAll(s) ) then         
          !モル比
          MolFr(i,s) = MolFr(i-1,s)
          
          !凝結していないので潜熱なし.
          LatentHeatAll(s) = 0.0d0
          
        else      

          !飽和蒸気圧と圧力から現在のモル数を計算
          MolFr(i,s) = max(SvapPressAll(s) / Press(i), 1.0d-16)
          
        end if
      end if
    end do
    
    !------------------------------------------------------------
    !NH4SH の平衡条件
    !------------------------------------------------------------
    num = 0.0d0
    Adjast: do
      del = 1.0d-9 * num
      MolFrNH3 = MolFr(i, Num_NH3) - del
      MolFrH2S = MolFr(i, Num_H2S) - del
      VapPressNH3 = MolFrNH3 * Press(i)     
      VapPressH2S = MolFrH2S * Press(i)     

      if (MolFrH2S <= 1.0d-16) then 
        MolFr(i,Num_H2S) = max(MolFrH2S, 1.0d-16)
        LatentHeatAll(Num_H2S) = 0.0d0
        exit Adjast
      end if

      if (MolFrNH3 <= 1.0d-16) then 
        MolFr(i,Num_NH3) = max(MolFrNH3, 1.0d-16)
        LatentHeatAll(Num_NH3) = 0.0d0
        exit Adjast
      end if

      if ( ChemPotNH4SH(ID_NH3, ID_H2S, Temp(i))  &
        &  >=   ChemPotOneSpc(ID_NH3, Temp(i), VapPressNH3) &
        &     + ChemPotOneSpc(ID_H2S, Temp(i), VapPressH2S) ) then 
        
        MolFr(i,Num_NH3) = max(MolFrNH3, 1.0d-16)
        MolFr(i,Num_H2S) = max(MolFrH2S, 1.0d-16)

        !反応熱解放 (これで OK だと思うが.... )
        del = MolFr(i-1,Num_H2S) - MolFr(i,Num_H2S)
        Q = 10834.0d0 * GasRUniv * del / DelZ / CpDryMol
        
        LatentHeatAll(Num_NH3) = 0.0d0
        LatentHeatAll(Num_H2S) = 0.0d0
        
        if (num == 0.0d0) then 
          LatentHeatAll(Num_NH3) = 0.0d0
        end if
        exit Adjast
        
      end if
      
      num = num + 1.0d0
    end do Adjast

    !------------------------------------------------------------
    !次のステップの温度, 圧力を求める
    !------------------------------------------------------------
    !モル比
    MolFr(i,0) = 1.0d0 - sum(MolFr(i,1:SpcNum))
    
    !平均分子量
    MolWtMean(i) = dot_product(MolWt, MolFr(i,:))

    !係数. 温度 Temp(i) で評価
    A = dot_product(LatentHeatAll(1:SpcNum), MolFr(i,1:SpcNum)) &
      &  / ( GasRUniv * Temp(i) )
    B = dot_product((LatentHeatAll(1:SpcNum) ** 2.0d0), MolFr(i,1:SpcNum)) &
      &  / ( CpDryMol * GasRUniv * ( Temp(i) ** 2.0d0 ) )
    
    !断熱温度減率
    DTempDz(i) = - Grav * MolWtMean(i) * ( 1.0d0 + A ) &
      &          / ( CpDryMol * ( 1.0d0 + B ) )        &
      &          + Q
    DTempDzDry(i) = - Grav * MolWtMean(i) / CpDryMol
    
    !温度を計算
    Temp(i+1) = Temp(i) + DTempDz(i) * DelZ

    !圧力を静水圧平衡から計算
    Press(i+1) = &
      & Press(i)  - ( Press(i) * Grav * DelZ * MolWtMean(i)) &
      & / ( GasRUniv * Temp(i) )
    
    !モル比 (RegZMax でも値を持つようにするため)
    MolFr(i+1,:) = MolFr(i,:)
    MolFr(i+1,0) = 1.0d0 - sum(MolFr(i,1:SpcNum))

  end do
  

!!!
!!!安定度の計算
!!!
  !分子量, 圧力差
  do i = RegZMin+1, RegZMax
    dM(i) = ( MolWtMean(i + 1) - MolWtMean(i - 1) ) * 5.0d-1
    dp(i) = ( Press(i + 1) - Press(i - 1) ) * 5.0d-1
  end do
  dM(RegZMin) = dM(RegZMin+1)
  dp(RegZMin) = dp(RegZMin+1)
  
  do i = RegZMin, RegZMax
    dMdz(i) = - dM(i) * Press(i) * Grav * MolWtMean(i) &
      &        / ( dp(i) * GasRUniv * Temp(i) )
  
    Stab(i) = Grav * ( DTempDz(i) - DTempDzDry(i) ) / Temp(i)  &
      &    - Grav *  dMdz(i) / MolWtMean(i) 

    write(*,*) "dM, dp   ", dM(i), dp(i)
    write(*,*) "dMdz     ", dMdz(i)
    write(*,*) "Deldtdzd ", DTempDz(i), DTempDzDry(i),DTempDz(i)-DTempDzDry(i) 
    write(*,*) "Stab     ", Stab(i)    
  end do

  where (Stab < 1.0d-7) 
    Stab = 1.0d-7
  end where
  

  call HistoryFile()  


contains

  subroutine HistoryFile()
    use gt4_history
    use fileset,  only:              &
      &                 exptitle,    &
      &                 expsrc,      &
      &                 expinst

    implicit none
    
    real :: ID(0:SpcNum)
    
    ID(0) = 1.0
    ID(1:SpcNum) = real(SpcWetID, 4)

    !-----------------------------------------------------------
    ! ヒストリー作成
    !-----------------------------------------------------------
    call HistoryCreate(                              &
      & file = 'arare-eccm.nc',                      &
      & title = exptitle,                            &
      & source = expsrc,                             &
      & institution = expinst,                       &
      & dims=(/'p','s'/),                            &
      & dimsizes=(/RegZMax+1, SpcNum+1/),              &
      & longnames=(/'Pressure    ',                  &
      &             'Species Num '/),                &
      & units=(/'Pa','1 '/), origin=0.0,               &
      & interval=0.0 )

    !-----------------------------------------------------------  
    ! 変数出力
    !-----------------------------------------------------------
    call HistoryPut('p', real(Press(RegZMin:RegZMax),4) )
    call HistoryPut('s', ID)
    call HistoryAddAttr( 'p', 'positive', 'down' )  

!    !圧力
!    call HistoryAddVariable(                         &
!      & varname='Press', dims=(/'p'/),               &
!      & longname='Pressure',                         &
!      & units='Pa',                                  &
!      & xtype='double' )

    !温度
    call HistoryAddVariable(                         &
      & varname='temp', dims=(/'p'/),                &
      & longname='Temperature',                      &
      & units='K',                                   &
      & xtype='double' )

    !モル比
    call HistoryAddVariable(                         &
      & varname='molfr', dims=(/'p','s'/),           &
      & longname='Mole Fraction',                    &
      & units='1',                                   &
      & xtype='double' )

    !温度減率
    call HistoryAddVariable(                         &
      & varname='dtdz', dims=(/'p'/),             &
      & longname='Adiabatic Lapse Rate',             &
      & units='K/km',                                &
      & xtype='double' )

    !安定度
    call HistoryAddVariable(                         &
      & varname='stab', dims=(/'p'/),                &
      & longname='Static Stability',                 &
      & units='s^-2',                                &
      & xtype='double' )

!    call HistoryPut( 'Press',   Press(RegZMin:RegZMax))
    call HistoryPut( 'temp',    Temp(RegZMin:RegZMax))
    call HistoryPut( 'molfr',   MolFr(RegZMin:RegZMax,:))
    call HistoryPut( 'dtdz',  - DTempDz(RegZMin:RegZMax) * 1.0d3)
    call HistoryPut( 'stab',    Stab(RegZMin:RegZMax))

    call HistoryClose

  end subroutine HistoryFile

end subroutine eccm2
