!c Description:
!c  ギブス自由エネルギーの計算
!c
!c Current Code Owner: 
!c   sugiyama@gfd-dennou.org
!c
!c Histry: 
!c   Version    Date          Comment
!c   -------    ----------    --------
!c   1.0        2003-02-06    作成
!c
!c Copyright (C) SUGIYAMA Ko-ichiro, 2003, All rights reserved

module interface_gibbs_energy
  interface

     subroutine gibbs_energy(temp, pres, mol, property, mu0, mu, gibbs)
       use prm_type, only: type_property
       real(8), intent(in)                         :: temp      
       real(8), intent(in)                         :: pres      
       real(8), intent(in)                         :: mol(:)    
       real(8), intent(out), optional, allocatable :: mu0(:)    
       real(8), intent(out), optional, allocatable :: mu(:)     
       real(8), intent(out), optional              :: gibbs     
       type(type_property), intent(in)             :: property(:)
     end subroutine gibbs_energy

  end interface
end module interface_gibbs_energy




subroutine gibbs_energy(temp, pres, mol, property, mu0, mu, gibbs)
  use prm_phys, only: gasr, pres_ref
  use prm_type, only: type_property
  use interface_molfraction
  use interface_cp_spl
  use interface_cp_noblegas
  use interface_svap_amp
  use interface_svap_antoine
  use interface_phase2alpha
  use interface_name_cnd2gas
  use interface_name2num
  use interface_stop
  
  implicit none
  real(8), intent(in)                         :: temp    
  real(8), intent(in)                         :: pres    
  real(8), intent(in)                         :: mol(:)  
  real(8), intent(out), optional, allocatable :: mu0(:)
  real(8), intent(out), optional, allocatable :: mu(:) 
  real(8), intent(out), optional              :: gibbs 
  type(type_property), intent(in)             :: property(:) 
  
  real(8)                        :: mu0_gas       !気相の標準化学ポテンシャル
  real(8)                        :: cp_int        !比熱の温度積分
  real(8)                        :: cp_dvt_int    !比熱/温度の温度積分
  real(8)                        :: svap          !飽和蒸気圧
  real(8), allocatable           :: mu0_tmp(:)    !標準化学ポテンシャル
  real(8), allocatable           :: mu_tmp(:)     !化学ポテンシャル
  real(8), allocatable           :: molfr(:)      !モル比
  real(8), allocatable           :: alpha(:)      !係数
  integer                        :: s
  integer                        :: num
  integer                        :: type
  integer                        :: smax
  character(40)                  :: gasname  


!!!
!!!変数の初期化
!!!  
  smax = size(property, 1)
  allocate(mu0_tmp(smax), mu_tmp(smax))
  mu0_tmp = 0.0D0
  mu_tmp = 0.0D0
  call phase2alpha(property%phase, alpha)
  call molfraction(mol, property%phase, molfr)

!!!
!!!標準化学ポテンシャル, 化学ポテンシャルの計算
!!!
  CalcMu0: do s = 1, smax
!!!
!!!変数の初期化
!!!
     mu0_gas = 0.0D0
     cp_int = 0.0D0
     cp_dvt_int = 0.0D0
     svap = 0.0D0
     
!!!
!!!ラベルが "NobleGas" の場合
!!!
     if (property(s)%method == 1) then 
        !比熱, 比熱/温度の積分
        call cp_nobleGas(temp, cp_int=cp_int, cp_dvt_int=cp_dvt_int)        

        !標準化学ポテンシャルの計算        
        mu0_tmp(s) =  property(s)%ref%enth - temp * property(s)%ref%entr  &
             & + cp_int - temp * cp_dvt_int                     
        
!!!
!!!ラベルが "Gas" の場合
!!!
     elseif (property(s)%method == 2) then
        !比熱, 比熱/温度の積分
        call cp_spl(temp, property(s), cp_int=cp_int, cp_dvt_int=cp_dvt_int)
        
        !標準化学ポテンシャルの計算
        mu0_tmp(s) = property(s)%ref%enth - temp * property(s)%ref%entr &
             & + cp_int - temp * cp_dvt_int
              
!!!
!!!ラベルが飽和蒸気圧(AMP)の場合
!!!
     elseif (property(s)%method == 3) then 
        !飽和蒸気圧の計算
        call svap_amp(temp, property(s)%svap_cef, svap, type=type)
        
        if (type == 1) then
           mu0_tmp(s) = - (property(s)%svap_cef%low%A * temp ** 2.0d0)&
                & / (2.0d0 * property(s)%svap_cef%low%temp) &
                & + property(s)%svap_cef%low%B
           
        elseif (type == 2) then 
           mu0_tmp(s) = - (property(s)%svap_cef%high%A * (temp ** 2.0d0)) &
                & / 2.0d0 - property(s)%svap_cef%high%B * temp &
                & + property(s)%svap_cef%high%C
           
        elseif (type == 0) then 
           !対応する気相のエントロピーを得る
           call name_cnd2gas(property(s)%name, gasname)
           call name2num(gasname, property%name, num)
           
           mu0_gas = mu0_tmp(num)

           !標準化学ポテンシャルの計算
           mu0_tmp(s)  = mu0_gas + gasr * temp * dlog(svap/pres_ref)

        else
           write(*,*) "gibbs_entropy.f90 at AMP: unknown type ", type
        end if

!!!
!!!ラベルが飽和蒸気圧(Antoine)の場合
!!!
     elseif (property(s)%method == 4) then 
        !飽和蒸気圧の計算
        call svap_antoine(temp, property(s)%svap_cef, svap, type=type)
           
        if (type == 1) then 
           mu0_tmp(s) = - (property(s)%svap_cef%low%A * temp ** 2.0d0)&
                & / (2.0d0 * property(s)%svap_cef%low%temp) &
                & + property(s)%svap_cef%low%B
           
        elseif (type == 2) then 
           mu0_tmp(s) = - (property(s)%svap_cef%high%A * (temp ** 2.0d0)) &
                & / 2.0d0 - property(s)%svap_cef%high%B * temp &
                & + property(s)%svap_cef%high%C

        elseif (type == 0) then 
           !対応する気相のエントロピーを得る
           call name_cnd2gas(property(s)%name, gasname)
           call name2num(gasname, property%name, num)
           mu0_gas = mu0_tmp(num)
           
           !標準化学ポテンシャルの計算
           mu0_tmp(s)  = mu0_gas + gasr * temp * dlog(svap/pres_ref)
           
        else
           write(*,*) "gibbs_entropy.f90 at AMP: unknown type ", type

        end if


!!!
!!!ラベルが NH4SH の場合
!!!
     elseif (property(s)%method == 5) then

        !対応する気相のエントロピーを得る
        gasname = 'H2S-g'
        call name2num(gasname, property%name, num)
        mu0_gas = mu0_tmp(num)
        
        !対応する気相のエントロピーを得る
        gasname = 'NH3-g'
        call name2num(gasname, property%name, num)
        mu0_gas = mu0_gas + mu0_tmp(num)        
        
        !標準化学ポテンシャルの計算
        if (temp < property(s)%svap_cef%low%temp) then 
           mu0_tmp(s) = - (property(s)%svap_cef%low%A * temp ** 2.0d0)&
                & / (2.0d0 * property(s)%svap_cef%low%temp) &
                & + property(s)%svap_cef%low%B

        elseif (temp > property(s)%svap_cef%high%temp) then 
           mu0_tmp(s) = - (property(s)%svap_cef%high%A * (temp ** 2.0d0)) &
                & / 2.0d0 - property(s)%svap_cef%high%B * temp &
                & + property(s)%svap_cef%high%C
           
        else
           mu0_tmp(s)  = mu0_gas &
                & + gasr * temp * (61.781d0 - (10834.0d0 / temp)) &
                & - gasr * temp * dlog(1.0d12)
        end if



!!!
!!!ラベルが unknown な場合
!!!
     else
        call stop(property(s)%method, 'gibbs_energy.f90', 'unkown method')
     end if

  end do CalcMu0


!!!
!!!化学ポテンシャルの計算
!!!
  mu_tmp = mu0_tmp + gasr * temp * dlog(molfr) &
       & + alpha * gasr * temp * dlog(pres / pres_ref)


!!!
!!!出力
!!!
  !標準化学ポテンシャル
  if (present(mu0)) then 
     allocate(mu0(smax))
     mu0 = mu0_tmp
  end if

  !化学ポテンシャル
  if (present(mu)) then 
     allocate(mu(smax))
     mu = mu_tmp
  end if

  !ギブス自由エネルギー
  if (present(gibbs)) gibbs = dot_product(mol, mu_tmp)
  

!!!
!!!後始末
!!!
  deallocate(mu0_tmp, mu_tmp, molfr, alpha)

end subroutine gibbs_energy


