!c Description: 
!c   Atreya and Romani (1985) によって温度, 組成の鉛直分布を計算し, 
!c   それを用いて温度減率, 仮温度減率, 静的安定度を計算するプログラム
!c
!c Current Code Owner: 
!c   sugiyama@gfd-dennou.org
!c
!c Copyright (C) SUGIYAMA Ko-ichiro, 2004, All rights reserved

program main
  use gt4_history    

  implicit none

!!!
!!!設定する必要のある変数
!!!
  integer, parameter       :: NStep      = 2000 
  real(8), parameter       :: PresIni    = 2.0d6
  real(8), parameter       :: PresEnd    = 1.0d4
  real(8), parameter       :: TempIni    = 441.0d0
  real(8), parameter       :: Grav       = 23.2d0      !木星重力  

  real(8), parameter       :: factor     = 1.0d0 !太陽系元素存在度に対する比
  real(8), parameter       :: humidity   = 1.0d0 !湿度 < 1 
  character(40), parameter :: OutputFile = "eccm_jupiter_x1_441.0K-2.0d6Pa-test3.nc"  
  
  !真なら飽和蒸気圧から, 偽なら潜熱から, 物性値を計算
  logical, parameter       :: psw = .true. 
  
  ! H2O(g,l) の値
  real(8), parameter       :: solar     = 1.5d-3 
  real(8), parameter       :: W_MolWt   = 18.0d-3
  real(8), parameter       :: W_Cp      = 33.5d0
  real(8), parameter       :: Antoine_A = 7.9186968d0
  real(8), parameter       :: Antoine_B = 1636.909d0
  real(8), parameter       :: Antoine_C = 224.92d0
  real(8), parameter       :: CC_C      = 50161184813.21968d0       
  real(8), parameter       :: CC_LatentHeat = 40.66d3     !潜熱 (J/mol)
  
  ! NH3(g,s)の値
!  real(8), parameter       :: solar     = 1.8d-4
!  real(8), parameter       :: W_MolWt   = 17.0d-3 
!  real(8), parameter       :: W_Cp      = 33.76d0   !200 K
!  real(8), parameter       :: Antoine_A = 9.96382d0
!  real(8), parameter       :: Antoine_B = 1617.907d0
!  real(8), parameter       :: Antoine_C = 272.55d0
!  real(8), parameter       :: CC_C = 11440514731.94403d0    
!  real(8), parameter       :: CC_LatentHeat = 29.01d3     !潜熱 (J/mol)
  
!  ! CH4(g,s)の値
!  real(8), parameter       :: solar     = 6.0d-4
!  real(8), parameter       :: W_MolWt   = 16.0d-3
!  real(8), parameter       :: W_Cp      = 33.25d0   !100 K 
!  real(8), parameter       :: Antoine_A = 7.6954d0
!  real(8), parameter       :: Antoine_B = 532.2d0
!  real(8), parameter       :: Antoine_C = 275.0d0
!  real(8), parameter       :: CC_C = 570230814.6256214d0       
!  real(8), parameter       :: CC_LatentHeat = 9.118d3     !潜熱 (J/mol)o

  !非凝縮成分 (H/He = 0.095, 1 x solar)
  real(8), parameter        :: D_MolWt = 2.32d-3  
  real(8), parameter        :: D_Cp    = 27.66d0


!!!
!!!設定する必要がない変数
!!!
  real(8), allocatable     :: Temp(:)
  real(8), allocatable     :: dtdz(:)
  real(8), allocatable     :: dtdz_low(:)
  real(8), allocatable     :: dtdz_high(:)
  real(8), allocatable     :: stab(:)
  real(8), allocatable     :: stab_low(:)
  real(8), allocatable     :: stab_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), allocatable     :: LatentHeat(:)     !潜熱   (J/mol)
  
  real(8)                  :: VapPres           !分圧
  real(8)                  :: SatVapPres        !飽和蒸気圧
  
  real(8)                  :: dp
  real(8)                  :: dtdzDry
  
  real(8)                  :: temp0, temp1, temp2, temp3
  real(8)                  :: pres0, pres1, pres2, pres3
  real(8)                  :: dt1, dt2, dt3, dt4, dt
  real(8)                  :: dtdp
  
  integer                  :: i
  logical, allocatable     :: lcl(:)
  
  integer, parameter       :: d = 1
  integer, parameter       :: w = 2
  
  real(8), parameter       :: GasR = 8.31d0        !気体定数 (J/k mol)    

  integer                  :: p
  real(8)                  :: temp_sat
  real(8)                  :: MolWtMean_sat
  real(8)                  :: CpMean_sat
  real(8)                  :: LatentHeat_sat 
  real(8)                  :: dtdz_high_sat
  real(8), allocatable     :: dtdz_sat(:)
  real(8), allocatable     :: dtdz_sat2(:)
  real(8), allocatable     :: dtdz_app1(:)
  real(8), allocatable     :: dtdz_app11(:)
  real(8), allocatable     :: dtdz_app2(:)
  real(8), allocatable     :: dtdz_app22(:)
  real(8), allocatable     :: dtdz_ratio1(:)
  real(8), allocatable     :: dtdz_ratio2(:)
  real(8), allocatable     :: stab0(:)
  real(8), allocatable     :: stab1(:)
  real(8), allocatable     :: stab2(:)
  real(8), allocatable     :: MolFr_sat(:)


!!!
!!!配列の初期化
!!!
  allocate(Pres( Nstep ), &
       & stab( Nstep ), stab_low( Nstep ), stab_high( Nstep ), &
       & Temp( Nstep ),     &
       & dtdz( NStep ), dtdz_low( NStep ), dtdz_high( NStep ), &
       & MolWt( d:w ), MolWtMean( NStep ), &
       & Cp( d:w ), CpMean( NStep ), LatentHeat( NStep ), &
       & Mol( 1:3 ), MolFr( 0:Nstep, d:w ), lcl(NStep), &
       & MolFr_sat( NStep ) )
  
  Temp = TempIni
  dtdz = 0.0d0; dtdz_low = 0.0d0; dtdz_high = 0.0d0
  stab = 0.0d0; stab_low = 0.0d0; stab_high = 0.0d0
  MolFr = 0.0d0; Mol = 0.0d0
  lcl = .false. 
  LatentHeat = 0.0d0
  
  !分子量 (単位は kg/mol)
  MolWt(d) = D_MolWt
  MolWt(w) = W_MolWt

  !比熱
  Cp(d) = D_Cp
  Cp(w) = W_Cp

  !初期モル比
  MolFr(:,w) = solar * factor
  MolFr(:,d) = 1 - MolFr(:,w) 

  !解析用
  do p = 1, Nstep
     MolFr_sat(p) = 1.0d-3 &
          & + ( p - 1.0d0 ) * (1.0d0 - 1.0d-3 ) / ( Nstep - 1.0d0 ) 
!          & + ( p - 1.0d0 ) * (1.0d-1 - 1.0d-3 ) / ( Nstep - 1.0d0 ) 
  end do
  
  
!!!
!!!圧力刻みの設定
!!!
  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-1, w) * Pres(i)
     
     !飽和蒸気圧
     if (psw) then 
        call antoine( Temp(i), svap = SatVapPres )
     else
        call CC( Temp(i), svap = SatVapPres )
     end if
     
     write(*,*) i, VapPres < SatVapPres, VapPres, SatVapPres
     
     
     !飽和蒸気圧から凝結の有無を決める
     if ( VapPres < SatVapPres ) then 
        
        !モル比
        MolFr(i,:) = MolFr(i-1,:)
        
        !平均比熱の計算 (モル数は前のステップと変わらない)
        CpMean(i) = dot_product( MolFr(i, d:w), Cp(d:w) )
        
        !乾燥断熱減率を計算
        dtdp =  GasR * Temp(i) / ( CpMean(i) * pres(i) )
        
        !温度プロファイルを計算
        Temp(i+1) = Temp(i)  +  dtdp * dp
        
     else
        
        !凝結高度ではスイッチを入れる. 
        lcl(i) = .true. 

        !湿潤断熱減率をルンゲクッタ法を用いて計算
        temp0 = temp(i)
        pres0 = pres(i)
        call equiv_dtdp(Temp0, Pres0, dtdp)        
        dt1 = dtdp * dp
        
        Temp1 = Temp0 + dt1 / 2.0d0  
        Pres1 = Pres0 + dp / 2.0d0
        call equiv_dtdp(Temp1, Pres1, dtdp)
        dt2 = dtdp * dp
        
        Temp2 = Temp0 + dt2  / 2.0d0  
        Pres2 = Pres0 + dp / 2.0d0
        call equiv_dtdp(Temp2, Pres2, dtdp)
        dt3 = dtdp * dp

        Temp3 = Temp0 + dt3
        Pres3 = Pres0 + dp
        call equiv_dtdp(Temp3, Pres3, dtdp)
        dt4 = dtdp * dp
        
        dt = (dt1 + 2.0d0 * dt2 + 2.0d0 * dt3 + dt4) / 6.0d0
        Temp(i+1) = Temp(i) + dt

     end if
  end do
  lcl(NStep) = lcl(NStep - 1) !決まらない部分
  

!!!
!!!温度減率の計算
!!!
  !温度減率
  do i = 1, NStep
     if ( lcl(i) ) then      
        call equiv_dtdz( Temp(i), Pres(i), &
             & MolFr(i,:), MolWtMean(i), CpMean(i), LatentHeat(i), dtdz(i) )
     else        
        !モル比
        MolFr(i,:) = MolFr(i-1,:)

        !平均分子量の計算 (モル数は前のステップと変わらない)
        MolWtMean(i) = dot_product( MolFr(i, d:w), MolWt(d:w) )
        
        !平均比熱の計算 (モル数は前のステップと変わらない)
        CpMean(i) = dot_product( MolFr(i, d:w), Cp(d:w) )
        
        !乾燥断熱減率を計算
        dtdz(i) = - MolWtMean(i) * Grav / CpMean(i)
     end if
     
  end do
  
  !温度減率 (上限, 下限)
  do i = 1, NStep

     if ( lcl(i) ) then

        dtdz_high(i) = - MolWtMean(i) * Grav * Temp(i) / LatentHeat(i)    
        
        dtdz_low(i)  = - MolWtMean(i) * Grav *  &
             &          ( 1  &
             &            + LatentHeat(i) * Molfr(i,w) / ( GasR * Temp(i) )&
             &            - ( LatentHeat(i) ** 2.0d0 ) * MolFr(i,w)  &
             &                / ( CpMean(i) * GasR * ( Temp(i) ** 2.0d0 ) )&
             &           ) &
             & / CpMean(i)
     else
        dtdz_high(i) = dtdz(i)
        dtdz_low(i)  = dtdz(i)
     end if
  end do


!!!
!!!静的安定度. 
!!!
  do i = 1, NStep
     !乾燥断熱減率
     dtdzDry   =  - MolWtMean(i) * Grav / CpMean(i)

     !現在は凝結高度以外で LatentHeat = 0 なので.... 
     if ( lcl(i) ) then 
        stab(i) =   &
             &   Grav * ( dtdz(i) - dtdzDry ) / Temp(i) &
             & - ( Grav * MolFr(i,w) * ( MolWt(w) - MolWt(d) ) / MolWtMean(i) ) &
             &   * ( &
             &       ( LatentHeat(i) * dtdz(i) / ( GasR * Temp(i) ** 2.0d0 ) )&
             &     + ( MolWtMean(i) * Grav / ( GasR * Temp(i) )  ) &
             &     )
        
!        stab_high(i) =   &
!             &   Grav * ( dtdz_high(i) - dtdzDry ) / Temp(i) &
!             & - ( Grav * MolFr(i,w) * ( MolWt(w) - MolWt(d) ) / MolWtMean(i) ) &
!             &   * ( &
!             &       ( LatentHeat(i) * dtdz_high(i) / ( GasR * Temp(i) ** 2.0d0 ) ) &
!             &     + ( MolWtMean(i) * Grav / ( GasR * Temp(i) )  ) &
!             &     )

        stab_high(i) =   &
             &   Molwt(w) * ( Grav ** (2.0d0) ) &
             &   * ( 1.0d0 - Cp(w) * Temp(i) / LatentHeat(i) ) &
             &   / ( Cp(w) * Temp(i) )

        stab_low(i) &
             & = ( & 
             &       MolWtMean(i) * ( Grav ** 2.0d0 ) * MolFr(i,w) &
             &       / ( Cp(d) * Temp(i) ) &
             &     ) &
             &   * ( &
             &        LatentHeat(i) / ( GasR * Temp(i) ) &
             &      - CpMean(i) / GasR &
             &     ) &
             &   * ( &
             &        LatentHeat(i) / ( CpMean(i) * Temp(i) ) &
             &      + ( MolWt(w) - MolWt(d) ) / MolWtMean(i) &
             &     ) 
     else
        stab(i)      = 0.0d0
        stab_low(i)  = 0.0d0
        stab_high(i) = 0.0d0
     end if
  end do


!!!
!!!単位の換算 (m --> km), 符号の付け換え
!!!
  dtdz       = dtdz * (- 1.0d3 )
  dtdz_low   = dtdz_low * ( - 1.0d3 )
  dtdz_high  = dtdz_high * ( - 1.0d3 )
  

!!!
!!!物理量の書き出し
!!!
  call HistoryPut('molfr', MolFr(1:NStep,:))
  call HistoryPut('molwt', MolWtMean)
  call HistoryPut('cp', CpMean)
  call HistoryPut('temp',  temp)
  call HistoryPut('dtdz',  dtdz)
  call HistoryPut('dtdz_low',  dtdz_low)
  call HistoryPut('dtdz_high',  dtdz_high)
  call HistoryPut('stab', stab)
  call HistoryPut('stab_low', stab_low)
  call HistoryPut('stab_high', stab_high)


!!!
!!!モル比をパラメタとした計算
!!!
  !凝結温度
  Caltemp: do i = 1, Nstep
     if ( lcl(i) ) then 
        temp_sat       = temp(i)
        MolWtMean_sat  = MolWtMean(i)  
        CpMean_sat     = CpMean(i)     
        LatentHeat_sat = LatentHeat(i)    
        exit Caltemp
     end if
  end do Caltemp

  write(*,*) "temp       ", temp_sat 
  write(*,*) "MolWtMean  ", MolWtMean_sat 
  write(*,*) "CpMean     ", CpMean_sat
  write(*,*) "LatentHeat ", LatentHeat_sat 

  allocate( dtdz_sat(Nstep),    dtdz_sat2(Nstep), &
       &    dtdz_app1(Nstep),   dtdz_app11(Nstep), &
       &    dtdz_app2(Nstep),   dtdz_app22(Nstep), &
       &    dtdz_ratio1(Nstep), dtdz_ratio2(Nstep), &
       &    stab0(Nstep), stab1(Nstep), stab2(Nstep) )

  do p = 1, Nstep

     dtdzDry   = - MolWtMean_sat * Grav / CpMean_sat
     dtdz_sat( p )  = &
          & - ( MolWtMean_sat * Grav / CpMean_sat ) &
          &   * ( 1.0d0 + LatentHeat_sat * MolFr_sat(p) / ( GasR * Temp_sat ))&
          &   / ( 1.0d0 + ( LatentHeat_sat ** 2.0d0 * MolFr_sat(p) &
          &             / ( CpMean_sat * GasR * Temp_sat ** 2.0d0 )) )
     dtdz_sat2( p ) = - MolWtMean_sat * Grav * Temp_sat / LatentHeat_sat     

     stab0( p ) =  &
          &   Grav * ( dtdz_sat(p) - dtdzDry ) / Temp_sat &
          & - ( Grav * MolFr_sat(p) * ( MolWt(w) - MolWt(d) ) &
          &        / MolWtMean_sat ) &
          &   * ( &
          &       ( LatentHeat_sat * dtdz_sat(p) &
          &          / ( GasR * Temp_sat ** 2.0d0 ) ) &
          &     + ( MolWtMean_sat * Grav / ( GasR * Temp_sat )  ) &
          &     )
    
     stab1( p ) = &
          &   ( & 
          &       MolWt(d) * ( Grav ** 2.0d0 ) * MolFr_sat(p) &
          &       / ( Cp(d) * Temp_sat ) &
          &     ) &
          &   * ( &
          &        LatentHeat_sat / ( GasR * Temp_sat ) &
          &      - Cp(d) / GasR &
          &     ) &
          &   * ( &
          &        LatentHeat_sat / ( Cp(d) * Temp_sat ) &
          &      + ( MolWt(w) - MolWt(d) ) / MolWt(d)    &
          &     ) 

!     stab1( p ) = &
!          &   ( & 
!          &       MolWtMean_sat * ( Grav ** 2.0d0 ) * MolFr_sat(p) &
!          &       / ( Cp(d) * Temp_sat ) &
!          &     ) &
!          &   * ( &
!          &        LatentHeat_sat / ( GasR * Temp_sat ) &
!          &      - CpMean_sat / GasR &
!          &     ) &
!          &   * ( &
!          &        LatentHeat_sat / ( CpMean_sat * Temp_sat ) &
!          &      + ( MolWt(w) - MolWt(d) ) / MolWtMean_sat &
!          &     ) 

     stab2( p ) =   &
          &   Molwt(w) * ( Grav ** (2.0d0) ) &
          &   * ( 1.0d0 - Cp(w) * Temp_sat / LatentHeat_sat ) &
          &   / ( Cp(w) * Temp_sat )

!     stab2( p ) =  &
!          &   Grav * ( dtdz_high_sat - dtdzDry ) / Temp_sat &
!          & - ( Grav * MolFr_sat(p) * ( MolWt(w) - MolWt(d) ) &
!          &        / MolWtMean_sat ) &
!          &   * ( &
!          &       ( LatentHeat_sat * dtdz_high_sat &
!          &          / ( GasR * Temp_sat ** 2.0d0 ) ) &
!          &     + ( MolWtMean_sat * Grav / ( GasR * Temp_sat )  ) &
!          &     )
  end do


  write(*,*) "OK! ",  CpMean_sat, MolWtMean_sat

  do p = 1, Nstep

     CpMean_sat = Cp(d) * ( 1.0d0 - MolFr_sat(p) ) + Cp(w) * MolFr_sat(p)
     MolWtMean_sat = MolWt(d) * ( 1.0d0 - MolFr_sat(p) ) &
          &           + MolWt(w) * MolFr_sat(p) 
     
     write(*,*) MolFr_sat(p), CpMean_sat, MolWtMean_sat
     
     dtdzDry   = - MolWtMean_sat * Grav / CpMean_sat
     dtdz_sat( p )  = &
          & - ( MolWtMean_sat * Grav / CpMean_sat ) &
          &   * ( 1.0d0 + LatentHeat_sat * MolFr_sat(p) / ( GasR * Temp_sat ))&
          &   / ( 1.0d0 + ( LatentHeat_sat ** 2.0d0 * MolFr_sat(p) &
          &             / ( CpMean_sat * GasR * Temp_sat ** 2.0d0 )) )
     dtdz_sat2( p ) = - MolWtMean_sat * Grav * Temp_sat / LatentHeat_sat     


     dtdz_app1( p ) = &
          &  LatentHeat_sat * MolFr_sat(p) / ( GasR * Temp_sat )
     
     dtdz_app11( p ) = 1.0d0 + dtdz_app1(p)
     
     dtdz_app2( p ) = &
          & LatentHeat_sat ** 2.0d0 * MolFr_sat(p) &
          &  / ( CpMean_sat * GasR * Temp_sat ** 2.0d0 )

     dtdz_app22( p ) = 1.0d0 + dtdz_app2(p)

!     dtdz_ratio1( p ) = dtdz_app1(p) / dtdz_app2(p) 

     dtdz_ratio1( p ) =  MolWtMean_sat * Grav &
          & * (   LatentHeat_sat ** 2.0d0 &
          &        / ( CpMean_sat * GasR * Temp_sat **2.0d0 ) &
          &     - LatentHeat_sat  / ( GasR * Temp_sat ) ) &
          & / ( CpMean_sat * ( dtdz_app2( p ) ** 2.0d0 ) )

     dtdz_ratio2( p ) =  dtdz_app11(p) / dtdz_app22(p) 

  end do
  


!!!
!!!物理量の書き出し
!!!
  call HistoryPut('stab0', stab0)
  call HistoryPut('stab1', stab1)
  call HistoryPut('stab2', stab2)
  call HistoryPut('sdtdz', dtdz_sat)
  call HistoryPut('sdtdz2', dtdz_sat2)
  call HistoryPut('dtdz_app1',   dtdz_app1)
  call HistoryPut('dtdz_app11',  dtdz_app11)
  call HistoryPut('dtdz_app2',   dtdz_app2)
  call HistoryPut('dtdz_app22',  dtdz_app22)
  call HistoryPut('dtdz_ratio1', dtdz_ratio1)
  call HistoryPut('dtdz_ratio2', dtdz_ratio2)


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


!!!
!!!積分定数を決める. 普段はコメントアウトしておく
!!!
  
!  call Antoine( 88.09d0, svap = SatVapPres )
!  write(*,*) "Ant", SatVapPres
  
!  call CC( 88.09d0, svap = SatVapPres )
!  write(*,*) "CC", SatVapPres

  
contains

  subroutine equiv_dtdp(TempTmp, PresTmp, dtdpTmp)
    real(8), intent(in)            :: TempTmp
    real(8), intent(in)            :: PresTmp
    real(8), intent(out)           :: dtdpTmp
    real(8)              :: svap
    real(8)              :: LatentHeat
    real(8)              :: A
    real(8)              :: B
    real(8)              :: MolFr0(d:w)
    real(8)              :: Cp0    
    
    !飽和蒸気圧
    if (psw) then 
       call Antoine(TempTmp, svap)
       call PhaseChangeEntalpy(TempTmp, LatentHeat)
    else
       call CC(TempTmp, svap, LatentHeat)
    end if
    
    !飽和蒸気圧と圧力から現在のモル数を計算
    MolFr0(w) = svap / PresTmp
    MolFr0(d) = 1.0d0 - MolFr0(w)
    
    !平均比熱の計算
    Cp0    = dot_product( MolFr0(d:w), Cp(d:w) )
    
    !係数. 温度 Temp で評価
    A = LatentHeat * MolFr0(w) / ( GasR * TempTmp )
    B = ( LatentHeat ** 2.0d0 ) * MolFr0(w)  &
         & / ( Cp0 * GasR * ( TempTmp ** 2.0d0 ) )
    
    !断熱減率
    dtdpTmp = GasR * TempTmp * ( 1 + A ) / ( Cp0 * PresTmp * ( 1 + B ) )
    
  end subroutine equiv_dtdp


  subroutine equiv_dtdz(TempTmp, PresTmp, &
       & MolFrTmp, MolWtTmp, CpTmp, LatentHeat, dtdzTmp)
    real(8), intent(in)            :: TempTmp
    real(8), intent(in)            :: PresTmp
    real(8), intent(out)           :: MolWtTmp
    real(8), intent(out)           :: MolFrTmp(d:w)
    real(8), intent(out)           :: CpTmp
    real(8), intent(out)           :: LatentHeat
    real(8), intent(out)           :: dtdzTmp
    real(8)                        :: svap
    real(8)                        :: A
    real(8)                        :: B
    
    
    !飽和蒸気圧
    if (psw) then 
       call Antoine(TempTmp, svap)
       call PhaseChangeEntalpy(TempTmp, LatentHeat)
    else
       call CC(TempTmp, svap, LatentHeat)
    end if

    !飽和蒸気圧と圧力から現在のモル数を計算
    MolFrTmp(w) = svap / PresTmp
    MolFrTmp(d) = 1.0d0 - MolFrTmp(w)
    
    !平均分子量
    MolWtTmp = dot_product( MolFrTmp(d:w), MolWt(d:w) )
    
    !平均比熱の計算
    CpTmp    = dot_product( MolFrTmp(d:w), Cp(d:w) )
    
    !係数. 温度 Temp で評価
    A = LatentHeat * MolFrTmp(w) / ( GasR * TempTmp )
    B = ( LatentHeat ** 2.0d0 ) * MolFrTmp(w)  &
         & / ( CpTmp * GasR * ( TempTmp ** 2.0d0 ) )

    !断熱減率
    dtdzTmp = - MolWtTmp * Grav * ( 1 + A ) / ( CpTmp * ( 1 + B ) )

  end subroutine equiv_dtdz


  subroutine CC(temp, svap, LatentHeat)    
    implicit none
    real(8), intent(in)            :: temp         !温度
    real(8), intent(out)           :: svap         !飽和蒸気圧
    real(8), intent(out),optional  :: LatentHeat   !潜熱
  
    !クラウジウスクラペイロンを積分した式
    svap =  CC_C * dexp( - CC_LatentHeat / ( GasR * temp) ) 
    
    !潜熱は定数
    if ( present(LatentHeat) ) LatentHeat = CC_LatentHeat

  end subroutine CC


  subroutine Antoine(temp, svap)    
    implicit none
    real(8), intent(in)    :: temp         !温度
    real(8), intent(out)   :: svap         !飽和蒸気圧
    real(8)                :: svap_log
  
    !Antoine の式
    svap_log = (Antoine_A - (Antoine_B / (Antoine_C + Temp - 273.15D0) ) ) &
         & * dlog(10.0D0) + dlog(133.322D0)
    svap = dexp(svap_log)
    
  end subroutine Antoine


  subroutine PhaseChangeEntalpy(temp, LatentHeat)    
    implicit none
    real(8), intent(in)    :: temp         !温度
    real(8), intent(out)   :: LatentHeat   !潜熱
    real(8)                :: svap_log_dt
  
    svap_log_dt =  antoine_B * dlog(10.0D0) &
         & / (antoine_C + temp - 273.15D0)**2.0D0 

    LatentHeat = svap_log_dt * GasR * ( Temp ** 2.0d0 ) 
 
  end subroutine PhaseChangeEntalpy


  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  ', 'solar' /), &
         & dimsizes = (/ NStep, w, Nstep /), &
         & longnames = (/"Pressure        ", &
         &               "chemical species", &
         &               "Mol ratio       "/), &
         & units = (/"Pa", "1 ", "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 HistoryAddAttr( 'pres', 'positive', 'down' )
    
    call HistoryPut('solar', real(molfr_sat, 4))
    
    
    !変数定義, 温度
    call HistoryAddVariable( &
         & varname="temp", &
         & dims=(/'pres'/), &
         & longname="Temperature", &
         & units='K', &
         & 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="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="stab0", &           
         & dims=(/'solar'/), &
         & longname="Stability", &
         & units='s^-2', &
         & xtype="double" )

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

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

    !変数定義, 温度減率
    call HistoryAddVariable( &
         & varname="sdtdz", &           
         & dims=(/'solar'/), &
         & longname="dT/dz (non-approx)", &
         & units='Pa', &
         & xtype="double" )

    !変数定義, 温度減率
    call HistoryAddVariable( &
         & varname="sdtdz2", &           
         & dims=(/'solar'/), &
         & longname="dT/dz (max)", &
         & units='Pa', &
         & xtype="double" )

    !変数定義, 温度減率
    call HistoryAddVariable( &
         & varname="dtdz_app1", &           
         & dims=(/'solar'/), &
         & longname="Contribution (1)", &
         & units='1', &
         & xtype="double" )

    !変数定義, 温度減率
    call HistoryAddVariable( &
         & varname="dtdz_app11", &           
         & dims=(/'solar'/), &
         & longname="Contribution (2)", &
         & units='1', &
         & xtype="double" )

    !変数定義, 温度減率
    call HistoryAddVariable( &
         & varname="dtdz_app2", &           
         & dims=(/'solar'/), &
         & longname="Contribution (3)", &
         & units='1', &
         & xtype="double" )

    !変数定義, 温度減率
    call HistoryAddVariable( &
         & varname="dtdz_app22", &           
         & dims=(/'solar'/), &
         & longname="Contribution (3)", &
         & units='1', &
         & xtype="double" )

    !変数定義, 温度減率
    call HistoryAddVariable( &
         & varname="dtdz_ratio1", &           
         & dims=(/'solar'/), &
         & longname="Contribution (3)", &
         & units='1', &
         & xtype="double" )

    !変数定義, 温度減率
    call HistoryAddVariable( &
         & varname="dtdz_ratio2", &           
         & dims=(/'solar'/), &
         & longname="Contribution (3)", &
         & units='1', &
         & xtype="double" )

 end subroutine gt4f90io_init


end program main
