! ECCM の再計算
!
! 内部で利用する式は 中島 (1998) を参照


program main
  use gt4_history    

  implicit none

  real(8), allocatable     :: Temp(:)
  real(8), allocatable     :: VTemp(:)
  real(8), allocatable     :: dtdz(:)
  real(8), allocatable     :: dtdz_low(:)
  real(8), allocatable     :: dtdz_high(:)
  real(8), allocatable     :: Vdtdz(:)
  real(8), allocatable     :: Vdtdz_low(:)
  real(8), allocatable     :: Vdtdz_high(:)
  real(8), allocatable     :: stab(:)
  real(8), allocatable     :: stab_low(:)
  real(8), allocatable     :: stab_high(:)
  real(8), allocatable     :: vstab(:)
  real(8), allocatable     :: vstab_low(:)
  real(8), allocatable     :: vstab_high(:)
  real(8), allocatable     :: Pres(:)
  real(8), allocatable     :: Mol(:)
  real(8), allocatable     :: MolWt(:)          !分子量 (kg/mol)
  real(8), allocatable     :: MolWtMean(:)      !分子量 (kg/mol)
  real(8), allocatable     :: MolFr(:,:)        !モル分率
  real(8), allocatable     :: Cp(:)             !比熱   (J/K mol)
  real(8), allocatable     :: CpMean(:)         !比熱   (J/K mol)

  real(8)                  :: VapPres           !分圧
  real(8)                  :: SatVapPres        !飽和蒸気圧
  real(8)                  :: SatVapPresDT      !飽和蒸気圧の温度勾配

  real(8)                  :: A
  real(8)                  :: B
  real(8)                  :: C
  real(8)                  :: dp
  real(8)                  :: VdtdzDry

  integer                  :: i
  logical, allocatable     :: lcl(:)


  real(8), parameter       :: GasR = 8.31d0        !気体定数 (J/k mol)
  real(8), parameter       :: LatentHeat = 40.66d3 !潜熱 (J/mol)
  real(8), parameter       :: Grav = 23.2d0        !木星重力
  real(8), parameter       :: solar  = 8.51d-4   !H/O

  real(8), parameter       :: PresIni = 1.0d6
  real(8), parameter       :: PresEnd = 1.0d4

  integer, parameter       :: NStep = 2000 
  integer, parameter       :: d = 1
  integer, parameter       :: w = 2

  character(40), parameter :: OutputFile = "eccm_x20.nc"  
  real(8), parameter       :: TempIni    = 4.5d2
  real(8), parameter       :: factor     = 20.0d0



!!!
!!!配列の初期化
!!!
  allocate(Pres( Nstep ), &
       & vstab( Nstep ), vstab_low( Nstep ), vstab_high( Nstep ), &
       & stab( Nstep ), stab_low( Nstep ), stab_high( Nstep ), &
       & Temp( Nstep ), vTemp( Nstep ),    &
       & dtdz( NStep ), dtdz_low( NStep ), dtdz_high( NStep ), &
       & vdtdz( NStep ), vdtdz_low( NStep ), vdtdz_high( NStep ), &
       & MolWt( d:w ), MolWtMean( NStep ), &
       & Cp( d:w ), CpMean( NStep ), &
       & Mol( 1:3 ), MolFr( Nstep, d:w ), lcl(NStep) )
  
  Temp = TempIni; vTemp = TempIni
  dtdz = 0.0d0; dtdz_low = 0.0d0; dtdz_high = 0.0d0
  vdtdz = 0.0d0; vdtdz_low = 0.0d0; vdtdz_high = 0.0d0
  vstab = 0.0d0; vstab_low = 0.0d0; vstab_high = 0.0d0
  stab = 0.0d0;  stab_low = 0.0d0;  stab_high = 0.0d0
  MolFr = 0.0d0; Mol = 0.0d0
  lcl = .false. 
  
  !分子量 (単位は kg/mol)
  MolWt(d) = 2.0d-3  * ( 0.5d0 / ( 0.5d0 + 0.0975d0) ) &
       &     + 4.0d-3 * ( 0.0975d0 / ( 0.5d0 + 0.0975d0) ) 
  MolWt(w)= 18.0d-3  !水(水蒸気)
  
  !比熱
  Cp(d) = 29.0d0  * ( 0.5d0 / ( 0.5d0 + 0.0975d0) ) &
       &  + 20.786d0 * ( 0.0975d0 / ( 0.5d0 + 0.0975d0) ) 
  Cp(w) = 33.5d0
  
!  !初期モル
!  Mol(d)   = 1.0d0
!  Mol(w)   = 2 * factor * solar / (1 - 2.0d0 * factor * solar)
!  Mol(w+1) = 0.04875d0
  
  !初期モル比
  MolFr(:,w) = 1.5d-3 * factor
  MolFr(:,d) = 1 - MolFr(:,w) 

  
  
!!!
!!!圧力刻みの設定
!!!
  Pres = PresIni
  do i = 2, NStep     
     Pres(i) = PresIni &
          & * 10.0d0**( ( i - 1 ) * dlog10( PresEnd / PresIni ) &
          & / ( NStep - 1 ) )
  end do
  

!!!
!!!NetCDF ファイルの用意
!!!
  call gt4f90io_init

  
!!!
!!!断熱減率の計算
!!!
  do i = 1, NStep - 1
     
     !圧力変化
     dp = Pres(i+1) - Pres(i)
     
     !前のステップでのモル分率を用いて現在の蒸気圧を計算
     VapPres = MolFr(i, w) * Pres(i)
     
     !飽和蒸気圧
     call Amp( Temp(i), svap = SatVapPres )
     
     if ( VapPres < SatVapPres ) then 
        !平均分子量 (モル数は前のステップと変わらない)
        MolWtMean(i) = dot_product( MolFr(i, d:w), MolWt(d:w) )
        
        !平均比熱の計算 (モル数は前のステップと変わらない)
        CpMean(i) = dot_product( MolFr(i, d:w), Cp(d:w) )
        
        !乾燥断熱減率を計算
        dtdz(i)      = - Grav * MolWtMean(i) / CpMean(i)
        dtdz_low(i)  =  dtdz(i)
        dtdz_high(i) =  dtdz(i)

     else
        lcl(i) = .true. 
        
        !飽和蒸気圧と圧力から現在のモル数を計算
        MolFr(i,w) = SatVapPres / Pres(i)
        MolFr(i,d) = 1.0d0 - MolFr(i,w)
        
        !平均分子量
        MolWtMean(i) = dot_product( MolFr(i, d:w), MolWt(d:w) )
        
        !平均比熱の計算
        CpMean(i) = dot_product( MolFr(i, d:w), Cp(d:w) )

        !係数. 温度 Temp(i) で評価
        A = LatentHeat * Molfr(i,w) / ( GasR * Temp(i) )
        B = ( LatentHeat ** 2.0d0 ) * MolFr(i,w)  &
             & / ( CpMean(i) * GasR * ( Temp(i) ** 2.0d0 ) )

!        write(*,*) A, B

        !断熱減率
        dtdz(i)      = - MolWtMean(i) * Grav * ( 1 + A ) &
             &           / ( CpMean(i) * ( 1 + B ) )
        dtdz_low(i)  = - MolWtMean(i) * Grav * ( 1 + A - B ) / CpMean(i)
        dtdz_high(i) = - MolWtMean(i) * Grav * Temp(i) / LatentHeat
        
!        write(*,*) dtdz(i), dtdz_low(i), dtdz_high(i)
     end if
     
     Temp(i+1) = Temp(i) &
          & - dtdz(i) * dp * GasR * Temp(i) &
          &   / ( Pres(i) * Grav * MolWtMean(i) )
  end do
  
  !決まらない部分を与える
  dtdz( NStep )       = dtdz( NStep - 1 )
  dtdz_low( NStep )   = dtdz_low( NStep - 1 )
  dtdz_high( NStep )  = dtdz_high( NStep - 1 )
  Temp( NStep )       = Temp( NStep - 1 )
  MolFr( NStep, : )   = MolFr( NStep - 1, : )
  MolWtMean( NStep )  = MolWtMean( NStep - 1 )
  CpMean( NStep )     = CpMean( NStep - 1 )
  
  
!!!
!!!仮温度への変換
!!!
!!!平均分子量, モル分率は既に求まっている  
  
  !乾燥気塊の仮温度減率
  VdtdzDry = - Grav * MolWt(d) / Cp(d)

  do i = 1, NStep - 1
     call Amp( Temp(i), svap = SatVapPres, svap_dt = SatVapPresDT )
     
     vdtdz(i)      = dtdz(i)      * MolWt(d) / MolWtMean(i)      
     vdtdz_low(i)  = dtdz_low(i)  * MolWt(d) / MolWtMean(i)      
     vdtdz_high(i) = dtdz_high(i) * MolWt(d) / MolWtMean(i)      
     
     if ( lcl(i) ) then 

!        vdtdz(i) = vdtdz(i) &
!             & + MolWt(d) * Pres(i) * Grav &
!             &  * ( MolWtMean(i+1) - MolWtMean(i) ) &
!             &  / ( MolWtMean(i) * GasR * ( Pres(i+1) - Pres(i) ) )

        vdtdz_low(i) = vdtdz_low(i) &
             & + MolWt(d) * Pres(i) * Grav &
             &  * ( MolWtMean(i+1) - MolWtMean(i) ) &
             &  / ( MolWtMean(i) * GasR * ( Pres(i+1) - Pres(i) ) )

        vdtdz_high(i) = vdtdz_high(i) &
             & + MolWt(d) * Pres(i) * Grav &
             &  * ( MolWtMean(i+1) - MolWtMean(i) ) &
             &  / ( MolWtMean(i) * GasR * ( Pres(i+1) - Pres(i) ) )


        
        vdtdz(i) = vdtdz(i) &
             & - MolWt(d) * ( MolWt(w) - MolWt(d) ) &
             &   * LatentHeat * MolFr(i,w) * dtdz(i) &
             &   / ( (MolWtMean(i) ** 2.0d0) * GasR * Temp(i) ) &
             & - MolWt(d) * ( MolWt(w) - MolWt(d) ) &
             &   * Grav * MolFr(i,w) &
             &   / ( MolWtMean(i) * GasR ) 

!        vdtdz_low(i) = vdtdz_low(i) &
!             & - MolWt(d) * ( MolWt(w) - MolWt(d) ) &
!             &   * LatentHeat * MolFr(i,w) * dtdz_low(i) &
!             &   / ( (MolWtMean(i) ** 2.0d0) * GasR * Temp(i) ) &
!             & - MolWt(d) * ( MolWt(w) - MolWt(d) ) &
!             &   * Grav * MolFr(i,w) &
!             &   / ( MolWtMean(i) * GasR ) 

!        vdtdz_high(i) = vdtdz_high(i) &
!             & - MolWt(d) * ( MolWt(w) - MolWt(d) ) &
!             &   * LatentHeat * MolFr(i,w) * dtdz_high(i) &
!             &   / ( (MolWtMean(i) ** 2.0d0) * GasR * Temp(i) ) &
!             & - MolWt(d) * ( MolWt(w) - MolWt(d) ) &
!             &   * Grav * MolFr(i,w) &
!             &   / ( MolWtMean(i) * GasR ) 

!        write(*,*) vdtdz(i)

!        vdtdz(i) = vdtdz(i) &
!             & - MolWt(d) * ( MolWt(w) - MolWt(d) ) * Temp(i) &
!             &   * dtdz(i) * SatVapPresDT &
!             &   / ( (MolWtMean(i) ** 2.0d0) * Pres(i) ) &
!             & - MolWt(d) * ( MolWt(w) - MolWt(d) ) &
!             &   * Grav * MolFr(i,w) &
!             &   / ( MolWtMean(i) * GasR ) 

!        write(*,*) vdtdz(i)

!        vdtdz_low(i) = vdtdz_low(i) &
!             & + MolWt(d) * Pres(i) * Grav &
!             &  * ( MolWtMean(i+1) - MolWtMean(i) ) &
!             &  / ( MolWtMean(i) * GasR * ( Pres(i+1) - Pres(i) ) )
        
!        vdtdz_high(i) = vdtdz_high(i) &
!             & + MolWt(d) * Pres(i) * Grav &
!             &  * ( MolWtMean(i+1) - MolWtMean(i) ) &
!             &  / ( MolWtMean(i) * GasR * ( Pres(i+1) - Pres(i) ) )

!        vdtdz(i) =  vdtdz(i) &
!             & - MolWt(d) * ( MolWt(w) - MolWt(d) ) &
!             &   * SatVapPres * Grav      &
!             &   / ( MolWtMean(i) * Pres(i) * GasR) &
!             & - MolWt(d) * ( MolWt(w) - MolWt(d) ) * Temp(i)  &
!             &   * dtdz(i) * SatVapPresDT  &
!             &   / ( ( MolWtMean(i) ** 2.0d0 ) * Pres(i) ) 
        
!        vdtdz_low(i) = vdtdz_low(i)  &
!             & - MolWt(d) * ( MolWt(w) - MolWt(d) ) &
!             &   * SatVapPres * Grav      &
!             &   / ( MolWtMean(i) * Pres(i) * GasR) &
!             & - MolWt(d) * ( MolWt(w) - MolWt(d) ) * Temp(i)  &
!             &   * dtdz_low(i) * SatVapPresDT  &
!             &   / ( ( MolWtMean(i) ** 2.0d0 ) * Pres(i) ) 

!        vdtdz_high(i) = vdtdz_high(i) &
!             & - MolWt(d) * ( MolWt(w) - MolWt(d) ) &
!             &   * SatVapPres * Grav      &
!             &   / ( MolWtMean(i) * Pres(i) * GasR) &
!             & - MolWt(d) * ( MolWt(w) - MolWt(d) ) * Temp(i)  &
!             &   * dtdz_high(i) * SatVapPresDT  &
!             &   / ( ( MolWtMean(i) ** 2.0d0 ) * Pres(i) ) 
     end if
  end do
  vdtdz(NStep)      = vdtdz(NStep - 1)
  vdtdz_low(NStep)  = vdtdz_low(NStep - 1)
  vdtdz_high(NStep) = vdtdz_high(NStep - 1)
  

!!!
!!!仮温度
!!!
  VTemp = Temp(1) * MolWt(d) / MolWtMean(1)  !初期化

  do i = 1, NStep - 1
     vTemp(i+1) = vTemp(i) &
          & - vdtdz(i) * &
          &   (  GasR * Temp(i) * ( pres(i+1) - pres(i) ) ) &
          &   / ( pres(i) * Grav * MolWtMean(i)  )
  end do
  vtemp(NStep)      = vtemp(NStep-1)


!!!
!!!静的安定度
!!!
  do i = 1, NStep
     stab(i) = Grav * ( dtdz(i) - vdtdzDry ) / Temp(i)
     stab_low(i) &
          & =   Grav * ( dtdz_low(i) - vdtdzDry ) / Temp(i)
     stab_high(i) &
          & =   Grav * ( dtdz_high(i) - vdtdzDry ) / Temp(i)
  end do

  do i = 1, NStep
     vstab(i) = Grav * ( vdtdz(i) - VdtdzDry ) / VTemp(i)
     vstab_low(i) &
          & =   Grav * ( vdtdz_low(i) - VdtdzDry ) / VTemp(i)
     vstab_high(i) &
          & =   Grav * ( vdtdz_high(i) - VdtdzDry ) / VTemp(i)
  end do
  
  
!!!
!!!単位の換算 (m --> km), 符号の付け換え
!!!
  dtdz       = dtdz * (- 1.0d3 )
  dtdz_low   = dtdz_low * ( - 1.0d3 )
  dtdz_high  = dtdz_high * ( - 1.0d3 )
  vdtdz      = vdtdz * (- 1.0d3 )
  vdtdz_low  = vdtdz_low * ( - 1.0d3 )
  vdtdz_high = vdtdz_high * ( - 1.0d3 )
  


!!!
!!!物理量の書き出し
!!!
  call HistoryPut('molfr', MolFr)
  call HistoryPut('molwt', MolWtMean)
  call HistoryPut('cp', CpMean)
  call HistoryPut('temp',  temp)
  call HistoryPut('vtemp', vtemp)
  call HistoryPut('dtdz',  dtdz)
  call HistoryPut('vdtdz', vdtdz)
  call HistoryPut('stab', vstab)
  call HistoryPut('vstab', vstab)

!  call HistoryPut('temp_low',  Temp_low)
!  call HistoryPut('vtemp_low', vTemp_low)
  call HistoryPut('dtdz_low',  dtdz_low)
  call HistoryPut('vdtdz_low', vdtdz_low)
  call HistoryPut('stab_low', vstab_low)
  call HistoryPut('vstab_low', vstab_low)

!  call HistoryPut('temp_high',  Temp_high)
!  call HistoryPut('vtemp_high', vTemp_high)
  call HistoryPut('dtdz_high',  dtdz_high)
  call HistoryPut('vdtdz_high', vdtdz_high)
  call HistoryPut('stab_high', vstab_high)
  call HistoryPut('vstab_high', vstab_high)


!!!
!!!NetCDF ファイルを閉じる
!!!
  call HistoryClose


contains

  subroutine Antoine(temp, svap)    
    implicit none
    real(8), intent(in)    :: temp         !温度
    real(8), intent(out)   :: svap         !飽和蒸気圧
    real(8)                :: svap_log
  
    !以下は水での値
    real(8), parameter     :: AA = 8.184254d0
    real(8), parameter     :: AB = 1791.3d0
    real(8), parameter     :: AC = 238.1d0
    
    !Antoine の式
    svap_log = (AA - (AB / (AC + Temp - 273.15D0) ) ) &
         & * dlog(10.0D0) + dlog(133.322D0)
    svap = dexp(svap_log)

  end subroutine Antoine


  subroutine Amp(temp, svap, svap_dt)    
    implicit none
    real(8), intent(in)              :: temp         !温度
    real(8), intent(out)             :: svap         !飽和蒸気圧
    real(8), intent(out), optional   :: svap_dt      !飽和蒸気圧
    real(8)                          :: svap_log
    real(8)                          :: svap_log_dt
  
    !以下は水での値
    real(8), parameter     :: A = -5631.1206d0
    real(8), parameter     :: B =  -8.363602d0
    real(8), parameter     :: C = 8.2312d0
    real(8), parameter     :: D = -0.03861449d0
    real(8), parameter     :: E = 2.77494d-5

    svap_log = &
         & A / temp &
         & + B &
         & + (C * dlog(temp)) &
         & + (D * temp) &
         & + (E * (temp**2)) &
         & + dlog(1.0D-1)

    svap_log_dt = &
         & - (A / temp**2.0D0) &
         & + (C / temp) &
         & + (D) &
         & + (2.0D0 * E * temp) 
    
    svap = dexp(svap_log) 
    if ( present(svap_dt) ) svap_dt = dexp( svap_log_dt )
    
  end subroutine Amp


  
  subroutine gt4f90io_init

    implicit none
    real(8)                 :: spc(d:w)
    integer                 :: i 
    integer                 :: num
    character(1000)         :: dataset
    
    write(dataset, *) "; PresIni=", pres(1), &
         & "; PresEnd=", pres(NStep), &
         & "; TempIni=", temp(1), &
         & "; step=", NStep
    num = len_trim(dataset)
    
    call HistoryCreate(&
         & file = OutputFile, &
         & title = "Equilibrium Calculation by ECCM", &
         & source = dataset(:num+1), &
         & institution = "GFD_Dennou_Club oboro Project", &
         & dims = (/ 'pres', 'spc ' /), &
         & dimsizes = (/ NStep, w /), &
         & longnames = (/"Pressure        ", "chemical species"/), &
         & units = (/"Pa", "1 "/), &
         & origin = 0.0,  &
         & interval = 0.0)
    
    spc = (/(i, i = d, w)/)
    call HistoryPut('spc', real(spc, 4))
    
    call HistoryPut('pres', real(pres, 4))
    
    !変数定義, 温度
    call HistoryAddVariable( &
         & varname="temp", &
         & dims=(/'pres'/), &
         & longname="Temperature", &
         & units='K', &
         & xtype="double" )
    
    !変数定義, 温度
!    call HistoryAddVariable( &
!         & varname="temp_low", &
!         & dims=(/'pres'/), &
!         & longname="Temperature", &
!         & units='K', &
!         & xtype="double" )
    
    !変数定義, 温度
!    call HistoryAddVariable( &
!         & varname="temp_high", &
!         & dims=(/'pres'/), &
!         & longname="Temperature", &
!         & units='K', &
!         & xtype="double" )
    
    !変数定義, 仮温度
    call HistoryAddVariable( &
         & varname="vtemp", &
         & dims=(/'pres'/), &
         & longname="Temperature", &
         & units='K', &
         & xtype="double" )
    
    !変数定義, 仮温度
!    call HistoryAddVariable( &
!         & varname="vtemp_low", &
!         & dims=(/'pres'/), &
!         & longname="Temperature", &
!         & units='K', &
!         & xtype="double" )
    
    
!    !変数定義, 仮温度
!    call HistoryAddVariable( &
!         & varname="vtemp_high", &
!         & dims=(/'pres'/), &
!         & longname="Temperature", &
!         & units='K', &
!         & xtype="double" )
    
!    !変数定義, モル
!    call HistoryAddVariable( &
!         & varname="mol", &           
!         & dims=(/'pres', 'spc '/), &
!         & longname="Chemical Species Abundance", &
!         & units='mol', &
!         & xtype="double" )

    !変数定義, モル比
    call HistoryAddVariable( &
         & varname="molfr", &           
         & dims=(/'pres', 'spc '/), &
         & longname="Mol fraction", &
         & units='mol', &
         & xtype="double" )

    !変数定義, 平均分子量
    call HistoryAddVariable( &
         & varname="molwt", &           
         & dims=(/'pres'/), &
         & longname="Mean Mol Weight", &
         & units='mol', &
         & xtype="double" )

    !変数定義, モル
    call HistoryAddVariable( &
         & varname="cp", &           
         & dims=(/'pres'/), &
         & longname="Mean Specific Heat", &
         & units='mol', &
         & xtype="double" )

    !変数定義, 温度減率
    call HistoryAddVariable( &
         & varname="dtdz", &           
         & dims=(/'pres'/), &
         & longname="Lapse Rate", &
         & units='K/km', &
         & xtype="double" )
    
    !変数定義, 温度減率
    call HistoryAddVariable( &
         & varname="dtdz_low", &           
         & dims=(/'pres'/), &
         & longname="Lapse Rate", &
         & units='K/km', &
         & xtype="double" )
    
    !変数定義, 温度減率
    call HistoryAddVariable( &
         & varname="dtdz_high", &           
         & dims=(/'pres'/), &
         & longname="Lapse Rate", &
         & units='K/km', &
         & xtype="double" )
    
    !変数定義, 温度減率
    call HistoryAddVariable( &
         & varname="vdtdz", &           
         & dims=(/'pres'/), &
         & longname="Lapse Rate", &
         & units='K/km', &
         & xtype="double" )
    
    !変数定義, 温度減率
    call HistoryAddVariable( &
         & varname="vdtdz_low", &           
         & dims=(/'pres'/), &
         & longname="Lapse Rate", &
         & units='K/km', &
         & xtype="double" )
    
    !変数定義, 温度減率
    call HistoryAddVariable( &
         & varname="vdtdz_high", &           
         & dims=(/'pres'/), &
         & longname="Lapse Rate", &
         & units='K/km', &
         & xtype="double" )
     
    !変数定義, 温度減率
    call HistoryAddVariable( &
         & varname="stab", &           
         & dims=(/'pres'/), &
         & longname="Stability", &
         & units='s^-2', &
         & xtype="double" )
    
    !変数定義, 温度減率
    call HistoryAddVariable( &
         & varname="stab_low", &           
         & dims=(/'pres'/), &
         & longname="Stability", &
         & units='s^-2', &
         & xtype="double" )

    !変数定義, 温度減率
    call HistoryAddVariable( &
         & varname="stab_high", &           
         & dims=(/'pres'/), &
         & longname="Stability", &
         & units='s^-2', &
         & xtype="double" )
     
    !変数定義, 温度減率
    call HistoryAddVariable( &
         & varname="vstab", &           
         & dims=(/'pres'/), &
         & longname="Stability", &
         & units='s^-2', &
         & xtype="double" )
    
    !変数定義, 温度減率
    call HistoryAddVariable( &
         & varname="vstab_low", &           
         & dims=(/'pres'/), &
         & longname="Stability", &
         & units='s^-2', &
         & xtype="double" )

    !変数定義, 温度減率
    call HistoryAddVariable( &
         & varname="vstab_high", &           
         & dims=(/'pres'/), &
         & longname="Stability", &
         & units='s^-2', &
         & xtype="double" )
    
 end subroutine gt4f90io_init


end program main