!c Description:
!c   温度・圧力・初期組成から, 熱力学関数 G と平衡時の組成を導くルーチン
!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_mingibbs
  interface
     
     recursive subroutine mingibbs(temp, pres, emol, struct, property, &
          & mol_nxt, mol_ini, mol_cal)
       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(inout), optional  :: mol_cal(:,:)
       real(8), intent(out)              :: mol_nxt(:)
       type(type_property), intent(inout):: property(:)
     end subroutine mingibbs
     
  end interface
end module interface_mingibbs





recursive subroutine mingibbs(temp, pres, emol, struct, property, &
     & mol_nxt, mol_ini, mol_cal)
  use prm_type, only: type_property
  use prm_common, only: zero, debug, ErrR, ErrA, &
       & obrdebugmode, esw_store, ObrSetEsw
  use ssl2_linear              
  use interface_ChkErr
  use interface_ChkEquiv
  use interface_ChkSpcMol
  use interface_ChkElmMol
  use interface_ChkElmExist  
  use interface_ChkSpcExist
  use interface_matrix_coeff   
  use interface_matrix_const   
  use interface_equivmol       
  use interface_nonzeromol     
  use interface_initmol        
  use interface_gibbs_energy   
  use interface_entropy
  use interface_packing        
  use interface_unpacking      
  use interface_mingibbs       
  use interface_display        
  use interface_def_phase      
  use interface_phase_count    

  implicit none
  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(inout), optional  :: mol_cal(:,:)
  real(8), intent(out)              :: mol_nxt(:)
  type(type_property), intent(inout):: property(:)

  real(8)                           :: gibbs_pre
  real(8)                           :: gibbs_nxt
  real(8), allocatable              :: mol(:)
  real(8), allocatable              :: mol_pre(:)
  real(8), allocatable              :: entr(:)
  real(8), allocatable              :: mu(:)
  real(8), allocatable              :: mu0(:)
  real(8), allocatable              :: A(:)
  real(8), allocatable              :: B(:)
  real(8), allocatable              :: emol_np(:)
  real(8), allocatable              :: struct_np(:,:)
  real(8), allocatable              :: mol_ini_np(:)
  integer                           :: emax
  integer                           :: pmax
  integer                           :: i  
  integer                           :: cal
  logical, allocatable              :: ssw(:)
  logical, allocatable              :: esw(:)
  logical                           :: sw
  logical                           :: cnvsw
  type(type_property), allocatable  :: property_np(:)
  
  !表示
  if (debug) then     
     write(*,*) ""
     write(*,*) "call mingibbs", temp, pres
     write(*,*) ""
     write(*,*) property%name
  end if

  !相番号を付け直す
  call def_phase(property)
  if (debug) write(*,*) "phase", property%phase

  !基準値よりも少なくなった元素を削除
  call ChkElmExist(emol, struct, ssw, esw, esw_store, sw)
  if (sw) then 
     call rcv_mingibbs
     return
  end if
  
  !化学反応によって生成される化学種が存在できるかチェック
  call ChkSpcExist(property%method, property%name, ssw, sw)
  if (sw) then 
     call ChkElmMol(struct, ssw, esw)
     call rcv_mingibbs
     return
  end if
  
  !入力された元素数を元に初期組成を決める
  if (present(mol_ini)) then 
     call initmol(emol, struct, mol, mol_ini)  
  else
     call initmol(emol, struct, mol)  
  end if
  
  !初期組成でのギブス自由エネルギー, エントロピー
  call gibbs_energy(temp, pres, mol, property, &
       & mu = mu, mu0 = mu0, gibbs = gibbs_pre)
  call entropy(temp, pres, mol, property, entr = entr)
  if (debug) call display(mol, "初期組成")

  !凝縮相が単相で存在していないことをチェック
  call ChkEquiv(mu0, property%phase, property%name, ssw, sw)
  if (sw) then
     call ChkElmMol(struct, ssw, esw)
     call rcv_mingibbs
     return
  end if

  !最小化ルーチン実行前の初期化
  emax = size(struct, 2)
  call phase_count(property%phase, pmax)
  i = emax + pmax
  call ssl2_lsx_init(i)
  allocate(A(i*(i+1)/2), B(i))
  
  !最小化ルーチン
  min_gibbs: do 
     
     !係数行列作成
     call matrix_coeff(mol, property%phase, struct, A)	
     
     !定数行列作成
     call matrix_const(temp, mol, mu, property%phase, struct, B)
    
     !連立方程式を解く
     call ssl2_lsx(A, B, sw)
     if (sw) then 
        call rcv_debug_mingibbs  !非正則な場合
        stop
     end if
     
     !連立方程式の解からモル数計算
     call EquivMol(B, struct, property%phase, temp, mu, mol, mol_pre)
     call NonZeroMol(mol, mol_pre)
     if (debug)  call display(mol, "化学種存在量")
     if (present(mol_cal)) then 
        cal = minloc(mol_cal(:,1), 1, mol_cal(:,1) == 0.0d0)
        mol_cal(cal,:) = mol
     end if
     
     !平衡組成でのギブス自由エネルギー, エントロピー
     deallocate(mu, entr)
     call gibbs_energy(temp, pres, mol, property, mu = mu, gibbs = gibbs_nxt)
     call entropy(temp, pres, mol, property, entr = entr)
     
     !化学種のモル数をチェック
     call ChkSpcMol(mol, entr, ErrR, &
          & property%phase, property%name, property%method, ssw, sw)
     if (sw) then 
        if (debug) call display(entr, "entropy")
        call ChkElmMol(struct, ssw, esw)
        call rcv_mingibbs
        return
     end if
     
     !収束条件
     call ChkErr(mol_pre, mol, ErrR, ErrA, cnvsw)
     if (cnvsw) then 
        call unpacking(mol, property%num, zero, mol_nxt)
        return 
     end if
     
     !ループを回す条件
     gibbs_pre = gibbs_nxt
     
  end do min_gibbs
  
contains
  
  subroutine rcv_mingibbs
    if (allocated(A)) deallocate(A, B)    
    call packing(emol, esw, emol_np)
    call packing(struct, ssw, esw, struct_np)
    call packing(property, ssw, property_np)
    
    if (present(mol_ini)) then 
       call packing(mol_ini, ssw, mol_ini_np)
       if (present(mol_cal)) then 
          call mingibbs(temp, pres, emol_np, &
               & struct_np, property_np, mol_nxt, mol_ini_np, &
               & mol_cal=mol_cal)
       else
          call mingibbs(temp, pres, emol_np, &
               & struct_np, property_np, mol_nxt, mol_ini_np)
       end if
    else    
       call mingibbs(temp, pres, emol_np, &
            & struct_np, property_np, mol_nxt)
    end if
  end subroutine rcv_mingibbs
  
  

  subroutine rcv_debug_mingibbs

    call ObrDebugMode

    if (present(mol_ini)) then 
       call mingibbs(temp, pres, emol, struct, property, mol, mol_ini=mol_ini)
    else
       call mingibbs(temp, pres, emol, struct, property, mol)
    end if

  end subroutine rcv_debug_mingibbs
  
  
end subroutine mingibbs
