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


program main
  use gt4_history    

  implicit none

  real(8), allocatable     :: Temp(:)
  real(8), allocatable     :: Temp_low(:)
  real(8), allocatable     :: Temp_high(:)
  real(8), allocatable     :: VTemp(:)
  real(8), allocatable     :: VTemp_low(:)
  real(8), allocatable     :: VTemp_high(:)
  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     :: vstab(:)
  real(8), allocatable     :: vstab_low(:)
  real(8), allocatable     :: vstab_high(:)
  real(8), allocatable     :: Pres(:)
  real(8), allocatable     :: Mol(:,:)          !モル数 (mol)
  real(8), allocatable     :: MolWt(:)          !分子量 (kg/mol)
  real(8), allocatable     :: MolFr(:,:)        !モル分率
  real(8), allocatable     :: Cp(:)             !比熱   (J/K mol)
  real(8), allocatable     :: MolWtMean(:)      !平均分子量
  real(8), allocatable     :: Fr(:)

  real(8)                  :: VapPres           !分圧
  real(8)                  :: SatVapPres        !飽和蒸気圧

  real(8)                  :: MolGas            !気相の全モル数
  real(8)                  :: MolWtDry          !非凝結成分の平均分子量
  real(8)                  :: MolWtWet          !凝結成分の平均分子量

  real(8)                  :: A
  real(8)                  :: B
  real(8)                  :: dp
  real(8)                  :: CpMean
  real(8)                  :: CpDry
  real(8)                  :: VdtdzDry

  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       :: PresIni = 1.0d6
  real(8), parameter       :: PresEnd = 1.0d4
  real(8), parameter       :: TempIni  = 3.45d2

  integer, parameter       :: NStep = 2000 
  integer, parameter       :: e = 3
  integer, parameter       :: v = 3
  integer, parameter       :: s = 4

  integer                  :: i
  logical, allocatable     :: lcl(:)
  character(40)            :: OutputFile


!!!
!!!配列の初期化
!!!
  allocate(Pres( Nstep ), &
       & vstab( Nstep ), vstab_low( Nstep ), vstab_high( Nstep ), &
       & Temp( Nstep ), Temp_low( Nstep ), Temp_high( Nstep ), &
       & VTemp( Nstep ), VTemp_low( Nstep ), VTemp_high( Nstep ), &
       & dtdz( NStep ), dtdz_low( NStep ), dtdz_high( NStep ), &
       & vdtdz( NStep ), vdtdz_low( NStep ), vdtdz_high( NStep ), &
       & Mol( Nstep, s ), MolWt( s ), Cp( s ), &
       & MolFr( Nstep, s ), MolWtMean( NStep), Fr( s ), lcl(NStep) )
  
  Temp = TempIni; Temp_low = TempIni; Temp_high = TempIni
  vTemp = TempIni; vTemp_low = TempIni; vTemp_high = 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
  MolFr = 0.0d0
  lcl = .false. 
  
  !分子量 (単位は kg/mol)
  MolWt(1) = 2.0d-3   !水素
  MolWt(2) = 4.0d-3   !ヘリウム
  MolWt(3) = 18.0d-3  !水(水蒸気)
  MolWt(4) = 18.0d-3  !水(氷)
  MolWtWet = 18.0d-3  !水
  
  !比熱
  Cp = 0.0d0
  Cp(1) = 29.0d0
  Cp(2) = 20.786d0
  Cp(3) = 33.5d0
  
  !初期モル数
!  Mol(:,1) = 5.0d-1    !水素
  Mol(:,1) = 0.491203999999999d0
!  Mol(:,2) = 9.75d-2   !ヘリウム
  Mol(:,2) = 0.0974999999999999d0
!  Mol(:,3) = 8.5d-6    !水(水蒸気) 0.01 x solar
!  Mol(:,3) = 8.5d-5    !水(水蒸気) 0.1 x solar
  Mol(:,3) = 8.5d-4    !水(水蒸気) 1 x solar
!  Mol(:,3) = 0.00425500d0    !水(水蒸気) 5 x solar
!  Mol(:,3) = 8.5d-3    !水(水蒸気) 10 x solar
!  Mol(:,3) = 8.5d-3 * 2.0d0   !水(水蒸気) 20 x solar
!  Mol(:,3) = 8.5d-3 * 5.0d0   !水(水蒸気) 20 x solar
!  Mol(:,3) = 8.5d-2    !水(水蒸気) 100 x solar
  Mol(:,4) = 0.0d0     !水 

  !出力ファイル
  OutputFile = "eccm_x1.nc"  
  
  
!!!
!!!圧力刻みの設定
!!!
  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)
     
     !前のステップでの気相の全モル数
     MolGas = sum(Mol(i, 1:v), 1)
     
     !前のステップでのモル分率
     MolFr(i, 1:v) = Mol(i, 1:v) / MolGas
     
     !前のステップでのモル分率を用いて現在の蒸気圧を計算
     VapPres = MolFr(i, e) * Pres(i)
     
     !飽和蒸気圧
     call Amp( Temp(i), SatVapPres )
     
     if ( VapPres < SatVapPres ) then 
        !平均分子量 (モル数は前のステップと変わらない)
        MolWtMean(i) = dot_product( MolFr(i, 1:v), MolWt(1:v) )
        
        !平均比熱の計算 (モル数は前のステップと変わらない)
        CpMean = dot_product( MolFr(i, 1:v), Cp(1:v) )

        !乾燥断熱減率を計算
        dtdz(i)      = - Grav * MolWtMean(i) / CpMean
        dtdz_low(i)  =  dtdz(i)
        dtdz_high(i) =  dtdz(i)

     else
        lcl(i) = .true. 

        !飽和蒸気圧と圧力から現在のモル数を計算
        Mol(i,e) = SatVapPres * ( MolGas - Mol(i,e) ) &
             & / ( Pres(i) - SatVapPres) 
        
        !現在の気相の全モル数
        MolGas = sum(Mol(i, 1:v), 1)
        
        !現在のモル分率
        MolFr(i, 1:v) = Mol(i,1:v) / MolGas
        
        !平均分子量
        MolWtMean(i) = dot_product( MolFr(i, 1:v), MolWt(1:v) )
        
        !平均比熱の計算
        CpMean = dot_product( MolFr(i, 1:v), Cp(1:v) )

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

        !断熱減率
        dtdz(i) = - Grav * MolWtMean(i) * ( 1 + A ) / ( CpMean * ( 1 + B ) )
        dtdz_low(i) =  - Grav * MolWtMean(i) * ( 1 + A - B ) / CpMean
        dtdz_high(i) =  - Grav * MolWtMean(i) * A / ( CpMean * B )

     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 )
  
  
!!!
!!!仮温度への変換
!!!
!!!平均分子量, モル分率は既に求まっている  
  
  !乾燥気塊の分子量. 最上層でのモル数を利用
  MolWtDry = MolWtMean( NStep )
  
  !乾燥気塊の比熱. 非凝縮性成分で考える
  Fr(1:v-1) = Mol(NStep, 1:v-1) / sum( Mol(NStep, 1:v-1) )
  CpDry = dot_product( Fr(1:v-1), Cp(1:v-1) )
  
  !乾燥気塊の仮温度減率
  VdtdzDry = - Grav * MolWtDry / CpDry
  
!  do i = 1, NStep
!     vtemp(i) =  temp(i) * MolWtDry / MolWtMean(i) 
!  end do
  
  
  !仮温度減率
  do i = 1, NStep
     
     vdtdz(i) = dtdz(i) *  MolWtDry / MolWtMean(i) 
     vdtdz_low(i) = dtdz_low(i) *  MolWtDry / MolWtMean(i) 
     vdtdz_high(i) = dtdz_high(i) *  MolWtDry / MolWtMean(i) 
     
     if ( lcl(i) )  then 
        call Amp( Temp(i), SatVapPres )
        
        vdtdz(i) = vdtdz(i)  &
             &     - MolWtDry * pres(i) * Grav &
             &        * ( MolWtDry - MolWtWet ) * SatVapPres &
             &        / ( MolWtMean(i) * GasR * ( pres(i) ** 2.0d0 ) )

        vdtdz_low(i) = vdtdz_low(i)  &
             &     - MolWtDry * pres(i) * Grav &
             &        * ( MolWtDry - MolWtWet ) * SatVapPres &
             &        / ( MolWtMean(i) * GasR * ( pres(i) ** 2.0d0 ) )
 
        vdtdz_high(i) = vdtdz_high(i)  &
             &     - MolWtDry * pres(i) * Grav &
             &        * ( MolWtDry - MolWtWet ) * SatVapPres &
             &        / ( MolWtMean(i) * GasR * ( pres(i) ** 2.0d0 ) )
     end if
  end do
  
!  do i = 1, NStep-1
!     vdtdz(i) = &
!          &   dtdz(i) * MolWtDry / MolWtMean(i) &
!          & + pres(i) * Grav * MolWtDry &
!          &     * ( MolWtMean(i + 1) - MolWtMean(i) ) &
!          & / ( GasR * MolWtMean(i) * ( pres(i+1) - pres(i) ) )

!     vdtdz_low(i) = &
!          &  dtdz_low(i) * MolWtDry / MolWtMean(i) &
!          & + pres(i) * Grav * MolWtDry &
!          &     * ( MolWtMean(i + 1) - MolWtMean(i) ) &
!          & / ( GasR * MolWtMean(i) * ( pres(i+1) - pres(i) ) )

!     vdtdz_high(i) = &
!          &   dtdz_high(i) * MolWtDry / MolWtMean(i) &
!          & + pres(i) * Grav * MolWtDry &
!          &     * ( MolWtMean(i + 1) - MolWtMean(i) ) &
!          & / ( GasR * MolWtMean(i) * ( pres(i+1) - pres(i) ) )

!     vdtdz_high(i) = &
!          & ( ( vTemp(i+1) - vTemp(i) )  &
!          &      *  Pres(i) * Grav * MolWtMean(i) )  &
!          & / ( ( pres(i+1) - pres(i) ) * ( - GasR ) * Temp(i) ) 

!  end do
!  vdtdz(NStep)      = vdtdz(NStep-1)
!  vdtdz_low(NStep)  = vdtdz_low(NStep-1)
!  vdtdz_high(NStep) = vdtdz_high(NStep-1)

  
  !仮温度
  VTemp = TempIni * MolWtDry / MolWtMean(1)  !初期化
  VTemp_low = TempIni * MolWtDry / MolWtMean(1)  !初期化
  VTemp_high = TempIni * MolWtDry / 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)  )
     VTemp_low(i+1) = VTemp(i) &
          & + vdtdz_low(i) * &
          &   ( ( - GasR ) * Temp(i) * ( pres(i+1) - pres(i) ) ) &
          &   / ( pres(i) * Grav * MolWtMean(i)  )
     VTemp_high(i+1) = VTemp(i) &
          & + vdtdz_high(i) * &
          &   ( ( - GasR ) * Temp(i) * ( pres(i+1) - pres(i) ) ) &
          &   / ( pres(i) * Grav * MolWtMean(i)  )

  end do
  vtemp(NStep)      = vtemp(NStep-1)
  vtemp_low(NStep)  = vtemp_low(NStep-1)
  vtemp_high(NStep) = vtemp_high(NStep-1)


  
  !静的安定度
  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('mol', Mol)
  call HistoryPut('temp',  temp)
  call HistoryPut('vtemp', vtemp)
  call HistoryPut('dtdz',  dtdz)
  call HistoryPut('vdtdz', vdtdz)
  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('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('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)    
    implicit none
    real(8), intent(in)    :: temp         !温度
    real(8), intent(out)   :: svap         !飽和蒸気圧
    real(8)                :: svap_log
  
    !以下は水での値
    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 = dexp(svap_log)
    
  end subroutine Amp


  
  subroutine gt4f90io_init

    implicit none
    real(8)                 :: spc(s)
    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, s /), &
         & longnames = (/"Pressure        ", "chemical species"/), &
         & units = (/"Pa", "1 "/), &
         & origin = 0.0,  &
         & interval = 0.0)
    
    spc = (/(i, i = 1, s)/)
    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="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="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