! 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     :: Mol(:,:)          !モル数 (mol)
  real(8), allocatable     :: Mol_low(:,:)     
  real(8), allocatable     :: Mol_high(:,:)    
  real(8), allocatable     :: MolFr(:,:)        !モル分率
  real(8), allocatable     :: MolFr_low(:,:)   
  real(8), allocatable     :: MolFr_high(:,:)  
  real(8), allocatable     :: MolWtMean(:)      !平均分子量
  real(8), allocatable     :: MolWtMean_low(:) 
  real(8), allocatable     :: MolWtMean_high(:)
  real(8), allocatable     :: Pres(:)
  real(8), allocatable     :: MolWt(:)          !分子量 (kg/mol)
  real(8), allocatable     :: Cp(:)             !比熱   (J/K mol)
  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


!!!
!!!配列の初期化
!!!
  !出力ファイル
  OutputFile = "eccm.nc"

  !配列の割り当て
  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 ), Mol_low( Nstep, s ), Mol_high( Nstep, s ), &
       & MolFr( Nstep, s ), MolFr_low( Nstep, s ), MolFr_high( Nstep, s ), &
       & MolWtMean( NStep), MolWtMean_low( NStep), MolWtMean_high( NStep), &
       & MolWt( s ), Cp( s ), Fr( s ) )

  !値を初期化する
  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
  Mol = 0.0d0; Mol_low = 0.0d0; Mol_high = 0.0d0; 
  MolFr = 0.0d0; MolFr_low = 0.0d0; MolFr_high = 0.0d0; 
  MolWtMean = 0.0d0; MolWtMean_low = 0.0d0; MolWtMean_high = 0.0d0; 
  
  !分子量 (単位は 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     !水(氷)
  Mol_low = Mol; Mol_high = Mol
  
  
!!!
!!!圧力刻みの設定
!!!
  Pres = PresIni
  do i = 2, NStep     
     Pres(i) = PresIni &
          & * 10.0d0**( ( i - 1 ) * dlog10( PresEnd / PresIni ) &
          & / ( NStep - 1 ) )
  end do
  

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


!!!
!!!断熱減率の計算
!!!  
  call lapserate(Temp, dtdz, Mol, MolFr, MolWtMean, 1)
  call lapserate(Temp_low, dtdz_low, Mol_low, MolFr_low, MolWtMean_low, 2)
  call lapserate(Temp_high, dtdz_high, Mol_high, MolFr_high, MolWtMean_high, 3)

!!!
!!!安定度の計算
!!!  
  call stability(Temp, Mol, MolWtMean, vTemp, vdtdz, vstab, 1)
  call stability(Temp_low, Mol_low, MolWtMean_low, vTemp_low, vdtdz_low, vstab_low, 2)
  call stability(Temp_high, Mol_high, MolWtMean_high, vTemp_high, vdtdz_high, vstab_high, 3)
  

!!!
!!!単位の換算 (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)
    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_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



!!!
!!!断熱減率の計算
!!!
  subroutine lapserate(Temp_anal, dtdz_anal, Mol_anal, &
       & MolFr_anal, MolWtMean_anal, case)
    implicit none
    
    real(8), intent(inout) :: Temp_anal(:)
    real(8), intent(out)   :: dtdz_anal(:)
    real(8), intent(inout) :: Mol_anal(:,:)
    real(8), intent(out)   :: MolFr_anal(:,:)
    real(8), intent(out)   :: MolWtMean_anal(:)
    integer, intent(in)    :: case 

    integer              :: i


    do i = 1, NStep - 1       

       !圧力変化
       dp = Pres(i+1) - Pres(i)

       !前のステップでの気相の全モル数
       MolGas = sum(Mol_anal(i, 1:v), 1)
       
       !前のステップでのモル分率
       MolFr_anal(i, 1:v) = Mol_anal(i, 1:v) / MolGas
       
       !前のステップでのモル分率を用いて現在の蒸気圧を計算
       VapPres = MolFr_anal(i, e) * Pres(i)
       
       !飽和蒸気圧
       call Amp( Temp_anal(i), SatVapPres )

       if ( VapPres < SatVapPres ) then 
          !平均分子量 (モル数は前のステップと変わらない)
          MolWtMean_anal(i) = dot_product( MolFr_anal(i, 1:v), MolWt(1:v) )
          
          !平均比熱の計算 (モル数は前のステップと変わらない)
          CpMean = dot_product( MolFr_anal(i, 1:v), Cp(1:v) )
          
          !乾燥断熱減率を計算
          dtdz_anal(i) = - Grav * MolWtMean_anal(i) / CpMean
          
       else
          !飽和蒸気圧と圧力から現在のモル数を計算
          Mol_anal(i,e) = SatVapPres * ( MolGas - Mol_anal(i,e) ) &
               & / ( Pres(i) - SatVapPres) 
          
          !現在の気相の全モル数
          MolGas = sum(Mol_anal(i, 1:v), 1)
          
          !現在のモル分率
          MolFr_anal(i, 1:v) = Mol_anal(i,1:v) / MolGas
          
          !平均分子量
          MolWtMean_anal(i) = dot_product( MolFr_anal(i, 1:v), MolWt(1:v) )
        
          !平均比熱の計算
          CpMean = dot_product( MolFr_anal(i, 1:v), Cp(1:v) )
          
          !係数
          A = LatentHeat * Molfr(i,e) / ( GasR * Temp_anal(i) )
          B = ( LatentHeat ** 2.0d0 ) * MolFr_anal(i,e) &
               & / ( CpMean * GasR * ( Temp_anal(i) ** 2.0d0 ) )
          
          if (case == 1) then 
             !断熱減率 (一般形)
             dtdz_anal(i) &
                  & = - Grav * MolWtMean_anal(i) * ( 1 + A ) &
                  & / ( CpMean * ( 1 + B ) )             
          elseif (case == 2) then 
             !断熱減率 (凝結成分が少ない近似)
             dtdz_anal(i)  &
                  & = - Grav * MolWtMean_anal(i) * ( 1 + A - B ) / CpMean
          elseif (case == 3) then 
             !断熱減率 (凝結成分が多い近似)
             dtdz_anal(i) &
                  & = - Grav * MolWtMean_anal(i) * A / ( CpMean * B )
          end if
       end if
       
       Temp_anal(i+1) = Temp_anal(i) &
            & + dtdz_anal(i) * dp * ( - GasR ) * Temp_anal(i) &
            & / ( Pres(i) * Grav * MolWtMean_anal(i) )       
    end do
    
    !決まらない部分を与える
    dtdz_anal( NStep )      = dtdz_anal( NStep - 1 )
    Mol_anal( NStep, : )    = Mol_anal( NStep - 1, : )
    MolFr_anal( NStep, : )  = MolFr_anal( NStep - 1, : )
    MolWtMean_anal( NStep ) = MolWtMean_anal( NStep - 1 )
    
  end subroutine lapserate
     

  subroutine stability(Temp_anal, Mol_anal, MolWtMean_anal, vTemp_anal, vdtdz_anal, vstab_anal, case)
    real(8), intent(in)  :: Temp_anal(:)
    real(8), intent(in)  :: Mol_anal(:,:)
    real(8), intent(in)  :: MolWtMean_anal(:)
    real(8), intent(out) :: vTemp_anal(:)
    real(8), intent(out) :: vdtdz_anal(:)
    real(8), intent(out) :: vstab_anal(:)
    integer, intent(in)  :: case
    
    
    !乾燥気塊の分子量. 最上層でのモル数を利用
    MolWtDry = MolWtMean_anal( NStep )
    
    !乾燥気塊の比熱. 非凝縮性成分で考える
    Fr(1:v-1) = Mol_anal(NStep, 1:v-1) / sum( Mol_anal(NStep, 1:v-1) )
    CpDry = dot_product( Fr(1:v-1), Cp(1:v-1) )
    
    !仮温度への変換
    vTemp_anal = Temp_anal *  MolWtDry / MolWtMean_anal
    
    !乾燥気塊の仮温度減率
    VdtdzDry = - Grav * MolWtDry / CpDry
    
    do i = 1, NStep - 1  
       !圧力差
       dp = Pres(i+1) - Pres(i)

       !仮温度減率
       if (case == 1) then 
          vdtdz_anal(i) = - ( VTemp_anal(i+1) - VTemp_anal(i) ) &
               & * ( Pres(i) * Grav * MolWtMean_anal(i) ) &
               & / ( GasR * Temp_anal(i) * dp ) 
       elseif (case == 2) then 
          vdtdz_anal(i) = - ( VTemp_anal(i+1) - VTemp_anal(i) ) &
               & * ( Pres(i) * Grav * MolWtMean_anal(i) ) &
               & / ( GasR * Temp_anal(i) * dp )
       elseif (case == 3) then 
          vdtdz_anal(i) = - ( VTemp_anal(i+1) - VTemp_anal(i)) &
               & * ( Pres(i) * Grav * MolWtMean_anal(i) ) &
               & / ( GasR * Temp_anal(i) * dp )
       end if
    end do
    vdtdz_anal(NStep) = vdtdz_anal(NStep - 1)
    
    !静的安定度
    vstab_anal = Grav * ( vdtdz_anal - VdtdzDry ) / VTemp_anal
    
  end subroutine stability


end program main
