!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-27    作成
!c
!c Copyright (C) SUGIYAMA Ko-ichiro, 2003, All rights reserved

module interface_entropy
  interface

     subroutine entropy(temp, pres, mol, property, &
          & entr0, entr, entr_all, entr_gas)
       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 :: entr0(:)    
       real(8), intent(out), optional, allocatable :: entr(:)      
       real(8), intent(out), optional              :: entr_all     
       real(8), intent(out), optional              :: entr_gas     
       type(type_property), intent(in)             :: property(:)  
     end subroutine entropy

  end interface
end module interface_entropy


subroutine entropy(temp, pres, mol, property, &
     & entr0, entr, entr_all, entr_gas)
  use prm_type, only: type_property
  use prm_phys, only: gasr, pres_ref
  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 :: entr0(:)
  real(8), intent(out), optional, allocatable :: entr(:)
  real(8), intent(out), optional              :: entr_all   
  real(8), intent(out), optional              :: entr_gas   
  type(type_property), intent(in)             :: property(:)

  real(8)                     :: entr0_gas
  real(8)                     :: cp_dvt_int
  real(8)                     :: svap
  real(8)                     :: svap_dt
  real(8), allocatable        :: entr0_tmp(:)
  real(8), allocatable        :: entr_tmp(:)
  real(8), allocatable        :: molfr(:)
  real(8), allocatable        :: alpha(:)
  integer                     :: s
  integer                     :: num
  integer                     :: type
  integer                     :: smax
  integer                     :: smax_gas
  character(40)               :: gasname

  
!!!
!!!変数の初期化
!!!
  smax = size(mol, 1)
  smax_gas = count(property%phase == 1)
  allocate(entr0_tmp(smax), entr_tmp(smax))
  entr0_tmp = 0.0d0
  entr_tmp = 0.0d0

  call phase2alpha(property%phase, alpha)
  call molfraction(mol, property%phase, molfr)

!!!
!!!単位モル当たりのエントロピーの計算
!!!
  Calc: do s = 1, smax

!!!
!!!変数の初期化
!!!
     entr0_gas = 0.0d0
     cp_dvt_int = 0.0d0
     svap       = 0.0d0 
     svap_dt    = 0.0d0 

!!!
!!!ラベルが "NobleGas" の場合
!!!
     if (property(s)%method == 1) then 
        call cp_noblegas(temp, cp_dvt_int=cp_dvt_int)
        entr0_tmp(s) = property(s)%ref%entr + cp_dvt_int

!!!
!!!ラベルが "Gas" の場合
!!!
     elseif (property(s)%method == 2) then
        call cp_spl(temp, property(s), cp_dvt_int=cp_dvt_int)
        entr0_tmp(s) = property(s)%ref%entr + cp_dvt_int
        
!!!
!!!ラベルが飽和蒸気圧 (AMP) の場合
!!!
     elseif (property(s)%method == 3) then 
        !飽和蒸気圧の 1 階微分を得る
        call svap_amp(temp, property(s)%svap_cef, &
             & svap, svap_dt=svap_dt, type=type)

        if (type == 1) then 
           entr0_tmp(s) =  (property(s)%svap_cef%low%A * temp) &
                & / property(s)%svap_cef%low%temp
            
        elseif (type == 2) then 
           entr0_tmp(s) =  property(s)%svap_cef%high%A * temp &
                & + property(s)%svap_cef%high%B
           
        elseif (type == 0) then            
           !対応する気相のエントロピーを得る
           call name_cnd2gas(property(s)%name, gasname)
           call name2num(gasname, property%name, num)
           entr0_gas = entr0_tmp(num)
           
           !エントロピー計算
           entr0_tmp(s) = entr0_gas &
                & - gasr * dlog(svap) - gasr * temp * dlog(svap_dt) &
                & + gasr * dlog(pres_ref)

        else
           write(*,*) "entropy.f90 at AMP: unknown type: ", type

        end if
        
           
!!!
!!!ラベルが飽和蒸気圧 (Antoine) の場合
!!!
     elseif (property(s)%method == 4) then 
        !飽和蒸気圧の 1 階微分を得る
        call svap_antoine(temp, property(s)%svap_cef, &
             & svap, svap_dt=svap_dt, type=type)
        
        if (type == 1) then 
           entr0_tmp(s) =  (property(s)%svap_cef%low%A * temp) &
                & / property(s)%svap_cef%low%temp
           
        elseif (type == 2) then 
           entr0_tmp(s) =  property(s)%svap_cef%high%A * temp &
                & + property(s)%svap_cef%high%B

        elseif (type == 0) then           
           !対応する気相のエントロピーを得る
           call name_cnd2gas(property(s)%name, gasname)
           call name2num(gasname, property%name, num)
           entr0_gas = entr0_tmp(num)
           
           !エントロピー計算
           entr0_tmp(s) = entr0_gas &
                & - gasr * dlog(svap) - gasr * temp * dlog(svap_dt) &
                & + gasr * dlog(pres_ref)

        else
           write(*,*) "entropy.f90 at Antoine: unknown type: ", type
        end if

           
!!!
!!!ラベルが NH4SH の場合
!!!
     elseif (property(s)%method == 5) then
        !NH3-g の値を求める
        gasname = 'NH3-g'
        call name2num(gasname, property%name, num)
        entr0_gas = entr0_tmp(num)
        
        !H2S-g の値を求める
        gasname = 'H2S-g'
        call name2num(gasname, property%name, num)
        entr0_gas = entr0_gas + entr0_tmp(num)

        !エントロピー計算
        if (temp < property(s)%svap_cef%low%temp) then 
           entr0_tmp(s) =  (property(s)%svap_cef%low%A * temp) &
                & / property(s)%svap_cef%low%temp

        elseif (temp > property(s)%svap_cef%high%temp) then 
           entr0_tmp(s) =  property(s)%svap_cef%high%A * temp &
                & + property(s)%svap_cef%high%B
           
        else
           entr0_tmp(s) = entr0_gas - gasr * 61.781d0  + gasr * dlog(1.0d12)
        end if
        
!!!
!!!ラベルが unknown の場合
!!!
     else
        call stop(property(s)%method, 'entropy.f90', 'unknown method')
     end if
        
  end do Calc


!!!
!!!各化学種のモルエントロピーを計算
!!!
  entr_tmp = entr0_tmp - gasr * dlog(molfr) &
       & - alpha * gasr * dlog(pres / pres_ref)


!!!  
!!!出力
!!!
  !モル当たりのエントロピー
  if (present(entr0)) then 
     allocate(entr0(smax))
     entr0 = entr0_tmp
  end if

  !モルエントロピー
  if (present(entr)) then 
     if (allocated(entr)) deallocate(entr)
     allocate(entr(smax))
     entr = entr_tmp
  end if

  !混合大気のエントロピーの計算
  if (present(entr_all)) entr_all = dot_product(mol, entr_tmp)  

  !気相だけのエントロピーの計算
  if (present(entr_gas)) &
       & entr_gas = dot_product(mol(1:smax_gas), entr_tmp(1:smax_gas))

!!!
!!!後始末
!!!
  deallocate(entr0_tmp, entr_tmp, molfr, alpha)


end subroutine Entropy

