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

module interface_equivmol
  interface

     subroutine EquivMol(B, struct, phase, temp, mu, mol, mol_pre)
       real(8), intent(in)     :: temp            !温度
       real(8), intent(inout)  :: B(:)            !解行列
       real(8), intent(in)     :: mu(:)           !化学ポテンシャル
       real(8), intent(in)     :: struct(:,:)     !化学種中に含まれる元素数
       real(8), intent(inout)  :: mol(:)          !化学種のモル数
       real(8), intent(out), allocatable :: mol_pre(:)   
       integer, intent(in)  :: phase(:)        !相番号
     end subroutine EquivMol
     
  end interface
end module interface_equivmol





subroutine equivmol(B, struct, phase, temp, mu, mol, mol_pre)
  use prm_phys, only: gasr
  use interface_phase_pack
  
  implicit none
  real(8), intent(in)     :: temp            !温度
  real(8), intent(inout)  :: B(:)            !解行列
  real(8), intent(in)     :: mu(:)           !化学ポテンシャル
  real(8), intent(in)     :: struct(:,:)     !化学種中に含まれる元素数
  real(8), intent(inout)  :: mol(:)          !化学種のモル数
  real(8), intent(out), allocatable :: mol_pre(:)   
  integer, intent(in)     :: phase(:)        !相番号

  real(8), allocatable    :: B_phi(:)        !作業用
  real(8), allocatable    :: struct_dot_B(:) !struct と B の内積
  integer                 :: i 
  integer                 :: smax, emax
  integer, allocatable    :: phase_tmp(:)     


!!!
!!!初期化
!!!  
  !配列の大きさを得る
  smax = size(struct, 1)
  emax = size(struct, 2)

  !配列の割り当て
  if (allocated(mol_pre)) deallocate(mol_pre)
  allocate(mol_pre(smax), &
       & B_phi(smax), &
       & struct_dot_B(smax))
  
  !初期値を代入
  mol_pre = mol

  allocate(phase_tmp(smax))
  phase_tmp = phase
  call phase_pack(phase_tmp)
 
  
!!!
!!!平衡組成計算の前準備
!!!
  !B_phi には相番号に相当する行の値を入れる
  do i = 1, smax
     B_phi(i) = B(emax+ phase_tmp(i))
  end do
  
  !struct と B の内積をとる
  do i = 1, smax
     struct_dot_B(i) = sum(B(1: emax) * struct(i, 1: emax)) 
  end do
  

!!!
!!!平衡計算
!!!
  !平衡組成の計算(仮)
  mol = &
       & - mol * mu / (gasr * temp) &
       & + mol * B_phi & 
       & + mol * struct_dot_B
  
  
!!!  
!!!後始末
!!!
  deallocate(B_phi, struct_dot_B, phase_tmp)
  
end subroutine equivmol





