!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_matrix_const
  interface

     subroutine matrix_const(temp, mol, mu, phase, struct, B)
       real(8), intent(in)     :: temp       
       real(8), intent(in)     :: mol(:)     
       real(8), intent(in)     :: mu(:)      
       real(8), intent(out)    :: B(:)       
       real(8), intent(in)     :: struct(:,:)
       integer, intent(in)     :: phase(:)   
     end subroutine matrix_const

  end interface
end module interface_matrix_const



subroutine matrix_const(temp, mol, mu, phase, struct, B)
  use prm_phys, only: gasr 
  use interface_phase_count
  use interface_phase_pack

  implicit none
  real(8), intent(in)  :: temp         !温度
  real(8), intent(in)  :: mol(:)       !化学種のモル数
  real(8), intent(in)  :: mu(:)        !化学ポテンシャル
  real(8), intent(out) :: B(:)         !定数行列
  real(8), intent(in)  :: struct(:,:)  !化学種に含まれる元素数
  integer, intent(in)  :: phase(:)     !化学種の相番号

  integer, allocatable :: phase_tmp(:)
  integer              :: pmax         !存在する相の数
  integer              :: smax         !存在する化学種の数
  integer              :: emax         !存在する元素の数
  integer              :: enum
  integer              :: p
  integer              :: max
  integer              :: min
  
  
!!!
!!!変数の初期化
!!!
  smax = size(struct, 1)
  emax = size(struct, 2)
  call phase_count(phase, pmax)
  B = 0.0D0

  allocate(phase_tmp(smax))
  phase_tmp = phase
  call phase_pack(phase_tmp)

  
!!!
!!!行列要素の計算
!!!
  do enum = 1, emax
     B(enum) = sum(struct(:,enum) * mol(:) * mu(:))  / (gasr * temp) &
          & + dot_product(struct(:, enum), mol(:))
  end do

  do p = 1, pmax
     min = minloc(phase_tmp, 1, phase_tmp == p)
     max = min + count(phase_tmp == p) - 1
     
     B(emax + p) = sum(mu(min: max) * mol(min: max)) / (gasr * Temp)
  end do


!!!
!!!後片付け
!!!
  deallocate(phase_tmp)
  

end subroutine matrix_const
