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

module interface_ChkSpcMol
  interface
     
     subroutine ChkSpcMol(mol, var, epsz, phase, name, method, ssw, sw)
       real(8), intent(in)               :: mol(:) 
       real(8), intent(in)               :: var(:) 
       real(8), intent(in)               :: epsz
       integer, intent(in)               :: phase(:)      
       integer, intent(in)               :: method(:)
       logical, intent(out), allocatable :: ssw(:)
       logical, intent(out)              :: sw
       character(*), intent(in)          :: name(:)
     end subroutine ChkSpcMol
     
  end interface
end module interface_ChkSpcMol



subroutine ChkSpcMol(mol, var, epsz, phase, name, method, ssw, sw)
  use interface_name_cnd2gas
  use interface_name2num
  use interface_display

  implicit none
  real(8), intent(in)               :: mol(:)       !化学種のモル数
  real(8), intent(in)               :: var(:) 
  real(8), intent(in)               :: epsz
  integer, intent(in)               :: phase(:) 
  integer, intent(in)               :: method(:)
  logical, intent(out), allocatable :: ssw(:)
  logical, intent(out)              :: sw
  character(*), intent(in)          :: name(:)

  real(8), allocatable              :: var1(:)
  real(8), allocatable              :: fr(:)
  real(8)                           :: cr
  integer                           :: s, N
  integer                           :: gas(2)
  integer                           :: liq(2)
  integer                           :: sol(2)
  integer                           :: num
  character(20)                     :: gasname

!!!
!!!初期値作成
!!!
  N = size(mol, 1)
  if (allocated(ssw)) deallocate(ssw)
  allocate(ssw(N), var1(N), fr(N))
  ssw = .true.; sw = .false. 
  var1 = mol * var

  gas(1) = 1
  gas(2) = count(phase == 1, 1)
  liq(1) = minloc(phase, 1, phase == 2)
  liq(2) = liq(1) + count(phase == 2, 1) - 1
  sol(1) = minloc(phase, 1, phase == 3)
  sol(2) = sol(1) + count(phase == 3, 1) - 1

!!!
!!!相毎に足し算する
!!!
  do s = 1, N
     if (phase(s) == 1) then 
        fr(s) = sum(var1(gas(1):gas(2))) / sum(var1)         
     elseif (phase(s) == 2) then 
        fr(s) = sum(var1(liq(1):liq(2))) / sum(var1)         
     elseif (phase(s) == 3) then 
        fr(s) = sum(var1(sol(1):sol(2))) / sum(var1)         
     else
        fr(s) = var1(s) / sum(var1) 
     end if
  end do
!  write(*,*) "chkerra, fr", fr

!!!
!!!元素の存在量が EPSZ よりも大きいかどうか判定する
!!!一度に化学種が 1 つしか消去されないようにしてある
!!!元素の消える条件より厳しくするために 1.0d-2 倍してある
!!!
  s = minloc(fr, 1)
  if (fr(s) <= epsz) then

     !この段階では気相が消えないようにする
     if (phase(s) == 2) then 
        ssw(liq(1):liq(2)) = .false.         
     elseif (phase(s) == 3) then 
        ssw(sol(1):sol(2)) = .false.         
     elseif(phase(s) == 4) then
        ssw(s) = .false. 
     end if
     sw = .true. 
!     write(*,*) "chkerra, ssw", ssw

     !消去される凝結相に対応する気相もしきい値より小さいならば気相も消す
     if (method(s) == 3 .OR. method(s) == 4) then 
        call name_cnd2gas(name(s), gasname)
        call name2num(gasname, name, num)

        fr(num) = var1(num) / sum(var1) 
        if (fr(num) <= epsz) ssw(num) = .false.
     end if

!     call display(fr, "fr")
!     write(*,*) "chkerra, ssw", ssw
  end if


!!!
!!!化学ポテンシャルの比較
!!!
  if (sw) then 
     s = minloc(fr, 1)
     if (method(s) == 3 .OR. method(s) == 4) then 
        call name_cnd2gas(name(s), gasname)
        call name2num(gasname, name, num)
        
        cr = abs((var(s) - var(num)) / var(num))
!        if (cr < 1.0d-4) then 
!           write(*,*) "even equiv", name(s), name(num)
!        end if
        
!        write(*,*) var(s), var(num)
     end if
  end if


!!!
!!!念のため
!!!
  if (size(mol,1) /= size(var,1)) then 
     write(*,*) "Dimension of mol and var are not the same"
     stop
  end if

!!!
!!!後片付け
!!!
  deallocate(var1, fr)


end subroutine ChkSpcMol
