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


program main
  use gt4_history    

  implicit none

  real(8), allocatable     :: Temp(:)
  real(8), allocatable     :: Temp_approx(:)
  real(8), allocatable     :: VTemp(:)
  real(8), allocatable     :: VTemp_approx(:)
  real(8), allocatable     :: dtdz(:)
  real(8), allocatable     :: dtdz_approx(:)
  real(8), allocatable     :: Vdtdz(:)
  real(8), allocatable     :: Vdtdz_approx(:)
  real(8), allocatable     :: vstab(:)
  real(8), allocatable     :: vstab_approx(:)
  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 = 2.0d6
  real(8), parameter       :: PresEnd = 1.0d4
  real(8), parameter       :: TempIni  = 4.5d2

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

  integer                  :: i
  character(40)            :: OutputFile


!!!
!!!配列の初期化
!!!
  allocate(Pres( Nstep ), &
       & vstab( Nstep ), vstab_approx( Nstep ), &
       & Temp( Nstep ), Temp_approx( Nstep ), &
       & VTemp( Nstep ), VTemp_approx( Nstep ), &
       & dtdz( NStep ), dtdz_approx( NStep ), &
       & vdtdz( NStep ), vdtdz_approx( NStep ), &
       & Mol( Nstep, s ), MolWt( s ), Cp( s ), &
       & MolFr( Nstep, s ), MolWtMean( NStep), Fr( s ) )
  
  !温度
  Temp = TempIni
  Temp_approx = TempIni
  dtdz = 0.0d0
  dtdz_approx = 0.0d0
  MolFr = 0.0d0
  
  !出力ファイル
  OutputFile = "eccm.nc"
  
  !分子量 (単位は kg/mol)
  MolWt(1) = 2.0d-3   !水素
  MolWt(2) = 4.0d-3   !ヘリウム
  MolWt(3) = 18.0d-3  !水(水蒸気)
  MolWt(4) = 18.0d-3  !水(氷)
  
  !比熱
  Cp = 0.0d0
  Cp(1) = 25.0d0
  Cp(2) = 20.786d0
  Cp(3) = 35.0d0
  
  !初期モル数
  Mol(:,1) = 5.0d-1    !水素
  Mol(:,2) = 9.75d-2   !ヘリウム
!  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) = 8.5d-3    !水(水蒸気) 10 x solar
!  Mol(:,3) = 8.5d-2    !水(水蒸気) 100 x solar
  Mol(:,4) = 0.0d0     !水(氷)
  
  
!!!
!!!圧力刻みの設定
!!!
  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_approx(i) =  - Grav * MolWtMean(i) / CpMean

     else
        !飽和蒸気圧と圧力から現在のモル数を計算
        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) )

        !係数
        A = LatentHeat * Molfr(i,e) / ( GasR * Temp(i) )
        B = ( LatentHeat ** 2.0d0 ) * MolFr(i,e) &
             & / ( CpMean * GasR * ( Temp(i) ** 2.0d0 ) )

        write(*,*) A
        write(*,*) B
        
        !断熱減率
        dtdz_approx(i) =  - Grav * MolWtMean(i) * ( 1 + A - B ) / CpMean
        dtdz(i) = - Grav * MolWtMean(i) * ( 1 + A ) / ( CpMean * ( 1 + B ) )
     end if

     Temp(i+1) = Temp(i) + dtdz(i) * dp * ( - GasR ) * Temp(i) &
          & / ( Pres(i) * Grav * MolWtMean(i) )
     Temp_approx(i+1) = Temp_approx(i) &
          & + dtdz_approx(i) * dp * ( - GasR ) * Temp(i) &
          & / ( Pres(i) * Grav * MolWtMean(i) ) 
  end do

  !決まらない部分を与える
  dtdz( NStep )        = dtdz( NStep - 1 )
  dtdz_approx( NStep ) = dtdz_approx( 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) )
  
  do i = 1, NStep
     !仮温度への変換
!     vTemp(i) = Temp(i) * ( MolWtMean(i) / MolWtDry )
!     vTemp_approx(i) = Temp_approx(i) * ( MolWtMean(i) / MolWtDry )
     vTemp(i) = Temp(i) *  MolWtDry / MolWtMean(i) 
     vTemp_approx(i) = Temp_approx(i) *  MolWtDry / MolWtMean(i) 
  end do
  
  !乾燥気塊の仮温度減率
  VdtdzDry = - Grav * MolWtDry / CpDry
  
  do i = 1, NStep - 1  
     !圧力差
     dp = Pres(i+1) - Pres(i)

     !仮温度減率
     vdtdz(i) = - ( VTemp(i+1) - VTemp(i) ) &
          & * ( Pres(i) * Grav * MolWtMean(i) )/ ( GasR * Temp(i) * dp ) 
     vdtdz_approx(i) = - ( VTemp_approx(i+1) - VTemp_approx(i) ) &
          & * ( Pres(i) * Grav * MolWtMean(i) )/ ( GasR * Temp(i) * dp )
  end do
  vdtdz(NStep) = vdtdz(NStep - 1)
  vdtdz_approx(NStep) = vdtdz_approx(NStep - 1)
  
  !静的安定度
  do i = 1, NStep - 1
     vstab(i) = Grav * ( vdtdz(i) - VdtdzDry ) / VTemp(i)
     vstab_approx(i) = Grav * ( vdtdz_approx(i) - VdtdzDry ) / VTemp(i)
     
     write(*,*) vdtdz(i), VdtdzDry, vstab(i)
  end do
  vstab(NStep) = vstab(NStep - 1)
  vstab_approx(NStep) = vstab_approx(NStep - 1)

  
!!!
!!!単位の換算 (m --> km)
!!!
  dtdz = dtdz * (- 1.0d3 )
  dtdz_approx = dtdz_approx * ( - 1.0d3 )
  vdtdz = vdtdz * (- 1.0d3 )
  vdtdz_approx = vdtdz_approx * ( - 1.0d3 )
  


!!!
!!!物理量の書き出し
!!!
  call HistoryPut('mol', Mol)
  call HistoryPut('temp', temp)
  call HistoryPut('temp_approx', Temp_approx)
  call HistoryPut('vtemp', vtemp)
  call HistoryPut('vtemp_approx', vTemp_approx)
  call HistoryPut('dtdz', dtdz)
  call HistoryPut('dtdz_approx', dtdz_approx)
  call HistoryPut('vdtdz', vdtdz)
  call HistoryPut('vdtdz_approx', vdtdz_approx)
  call HistoryPut('vstab', vstab)
  call HistoryPut('vstab_approx', vstab_approx)


!!!
!!!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)
    real(8)                 :: elm(s)
    real(8)                 :: sizes(20)
    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 ', 'elm ', 'size'/), &
         & dimsizes = (/NStep, s, s, 20/), &
         & longnames = (/"Pressure        ", "chemical species", &
         &               "element         ", "size of name    "/), &
         & units = (/"Pa", "1 ", "1 ", "1 "/), &
         & origin = 0.0,  &
         & interval = 0.0)
    
    spc = (/(i, i = 1, s)/)
    call HistoryPut('spc', real(spc, 4))
    
    elm = (/(i, i = 1, s)/)
    call HistoryPut('elm', real(elm, 4))
    
    call HistoryPut('pres', real(pres, 4))
    
    sizes = (/(i, i = 1, 20)/)
    call HistoryPut('size', real(sizes, 4))    
    
    !変数定義, 温度
    call HistoryAddVariable( &
         & varname="temp", &
         & dims=(/'pres'/), &
         & longname="Temperature", &
         & units='K', &
         & xtype="double" )
    
    !変数定義, 温度
    call HistoryAddVariable( &
         & varname="temp_approx", &
         & dims=(/'pres'/), &
         & longname="Temperature", &
         & units='K', &
         & xtype="double" )
    
    !変数定義, 仮温度
    call HistoryAddVariable( &
         & varname="vtemp", &
         & dims=(/'pres'/), &
         & longname="Temperature", &
         & units='K', &
         & xtype="double" )
    
    !変数定義, 仮温度
    call HistoryAddVariable( &
         & varname="vtemp_approx", &
         & 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_approx", &           
         & 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_approx", &           
         & 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_approx", &           
         & dims=(/'pres'/), &
         & longname="Stability", &
         & units='s^-2', &
         & xtype="double" )
    
  end subroutine gt4f90io_init


end program main
