!c Description:
!c   温度・圧力・初期組成から, 熱力学関数 G と平衡時の組成を導くルーチン
!c   相の存在をチェックするためのルーチン付き
!c
!c Current Code Owner: 
!c   sugiyama@gfd-dennou.org
!c
!c Histry: 
!c   Version    Date          Comment
!c   -------    ----------    --------
!c   1.0        2004-04-25    作成
!c
!c Copyright (C) SUGIYAMA Ko-ichiro, 2003, All rights reserved

module interface_mingibbs_chk
  interface
     
     subroutine mingibbs_chk(temp, pres, emol, struct, property, mol_nxt, mol_ini)
       use prm_type, only: type_property
       real(8), intent(in)               :: temp  
       real(8), intent(in)               :: pres
       real(8), intent(in)               :: emol(:)
       real(8), intent(in)               :: struct(:,:)
       real(8), intent(in), optional     :: mol_ini(:)
       real(8), intent(out)              :: mol_nxt(:)
       type(type_property), intent(inout):: property(:)
     end subroutine mingibbs_chk
     
  end interface
end module interface_mingibbs_chk





subroutine mingibbs_chk(temp, pres, &
     & emol, struct, property, mol_nxt, mol_ini)
  use interface_packing
  use interface_display
  use interface_mingibbs
  use prm_type, only: type_property
  use prm_common, only: zero, debug

  real(8), intent(in)               :: temp  
  real(8), intent(in)               :: pres
  real(8), intent(in)               :: emol(:)
  real(8), intent(in)               :: struct(:,:)
  real(8), intent(in), optional     :: mol_ini(:)
  real(8), intent(out)              :: mol_nxt(:)
  type(type_property), intent(inout):: property(:)

  real(8), allocatable              :: struct_chk(:,:)
  real(8), allocatable              :: emol_chk(:)
  real(8), allocatable              :: mol_chk(:)
  real(8), allocatable              :: mol_store(:)
  type(type_property), allocatable  :: property_chk(:)
  logical, allocatable              :: esw(:)
  logical, allocatable              :: ssw(:)
  logical, allocatable              :: csw(:)
  logical, allocatable              :: ssw_store(:)
  logical, allocatable              :: ssw_tmp(:,:)
  integer                           :: i, l
  integer                           :: smax
  integer                           :: emax
  integer                           :: gas(2)
  integer                           :: liq(2)
  integer                           :: sol(2)
  

!!!
!!!初期化
!!!  
  smax = size(struct, 1)
  emax = size(struct, 2)
  allocate(ssw(smax), ssw_store(smax), ssw_tmp(smax,2), &
       & esw(emax), csw(smax), mol_store(smax))  
  ssw_store = .false. 
  esw = .true. 
  mol_store = 0.0d0
  gas(1) = 1
  gas(2) = count(property%phase == 1, 1)
  liq(1) = minloc(property%phase, 1, property%phase == 2)
  liq(2) = liq(1) + count(property%phase == 2, 1) - 1
  sol(1) = minloc(property%phase, 1, property%phase == 3)
  sol(2) = sol(1) + count(property%phase == 3, 1) - 1


!!!
!!!相を追加しながら, 物質が存在できるかチェック
!!!気相のみの場合を加味して, ループは gas(2) から回す
!!!
  STABLE_PHASE: do i = gas(2), smax
!     write(*,*) "mingibbs_chk: ", i
     ssw = .false.
     csw = .false.
     
     !相番号が 1 は全て .true. 
     do l = gas(1), gas(2)
        ssw(l) = .true.
        if (i == gas(2)) csw(gas(1):gas(2)) = .true.
     end do
     
     !相番号が 2, 3, 4 の扱い
     if (property(i)%phase == 2) then 
        if (i == liq(1)) then 
           do l = liq(1), liq(2) 
              ssw(l) = .true. 
              csw(l) = .true. 
           end do
        else
           cycle STABLE_PHASE
        end if
     elseif (property(i)%phase == 3) then 
        if (i == sol(1)) then 
           do l = sol(1), sol(2) 
              ssw(l) = .true. 
              csw(l) = .true. 
           end do
        else
           cycle STABLE_PHASE
        end if
     elseif (property(i)%phase == 4) then 
        ssw(i) = .true. 
        csw(i) = .true. 
     end if
     
!     write(*,*) "ssw", ssw
!     write(*,*) "csw", csw
     
     !テスト計算する科学種を決める
     call packing(property, ssw, property_chk)
     call packing(struct, ssw, esw, struct_chk)
     call packing(emol, esw, emol_chk)
     if (present(mol_ini)) then 
        call packing(mol_ini, ssw, mol_chk)
        call mingibbs(temp, pres, &
             & emol_chk, struct_chk, property_chk, mol_nxt, mol_chk)
     else
        call mingibbs(temp, pres, &
             & emol_chk, struct_chk, property_chk, mol_nxt)
     end if
     
!     call display(mol_nxt, "計算結果")
     
     ssw_tmp(:,1) = ssw_store
     ssw_tmp(:,2) = mol_nxt > zero
     ssw_store = any(ssw_tmp, 2)
     mol_store = merge(mol_nxt, mol_store, csw) 
     
     !掃除
     deallocate(property_chk, struct_chk, emol_chk)
     if (present(mol_ini)) deallocate(mol_chk)
  end do STABLE_PHASE

  if (debug) call display(mol_store, "mol_store")     
  if (debug) write(*,*) "ssw_store", ssw_store     


!!!
!!!存在できる相を用いて計算
!!!
  call packing(property, ssw_store, property_chk)
  call packing(struct, ssw_store, esw, struct_chk)
  call packing(emol, esw, emol_chk)

  call packing(mol_store, ssw_store, mol_chk)
  call mingibbs(temp, pres, &
       & emol_chk, struct_chk, property_chk, mol_nxt, mol_chk)


!!!
!!!後かたづけ
!!!
  deallocate(struct_chk, emol_chk, property_chk, &
       esw, ssw, ssw_store, ssw_tmp)
  if (present(mol_ini)) deallocate(mol_chk)

end subroutine mingibbs_chk

