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

module interface_initmol
  interface
     
     subroutine initmol(emol, struct, mol, mol_ini)
       real(8), intent(in)                :: emol(:)   
       real(8), intent(in)                :: struct(:,:) 
       real(8), intent(in), optional      :: mol_ini(:)  
       real(8), intent(out), allocatable  :: mol(:)  
     end subroutine initmol
     
  end interface
end module interface_initmol




subroutine initmol(emol, struct, mol, mol_ini)
  use ssl2_linear
  use interface_display
  use prm_common, only: zero, CSeed, debug, ErrA
  
  implicit none
  real(8), intent(in)                :: emol(:)   
  real(8), intent(in)                :: struct(:,:) 
  real(8), intent(in), optional      :: mol_ini(:)  
  real(8), intent(out), allocatable  :: mol(:)  
  
  real(8), allocatable   :: A(:,:)
  real(8), allocatable   :: B(:)
  real(8), allocatable   :: FA(:,:)
  real(8), allocatable   :: FB(:)
  real(8), allocatable   :: emol_tmp(:)
  real(8), parameter     :: fct = 9.9d-1
  integer                :: i, e
  integer                :: N, M, num
  integer, allocatable   :: order(:)
  

  N = size(struct, 1)
  M = size(struct, 2)
  num = 1
  allocate(A(M, M), B(M), FA(M, M), FB(M), order(M), mol(N), emol_tmp(M))
  order = 0
  mol = 0.0d0

  
!!!
!!!係数行列の作成
!!!
  do e = 1, M
     
     use: do i = 1, N
        if (struct(i, e) >= 1 .AND. .not. any(order == i)) then 
           A(num, :) = struct(i, :)
           order(e) = i
           exit use
        end if
     end do use
     
     num = num + 1
  end do

  !転置行列に変換
  A = transpose(A)
  FA = A
  
  
!!!
!!!係数行列に使われなかった化学種の存在量を適当に与える
!!!
  not_use: do i = 1, N
     if (any(order == i)) cycle not_use
     
     emol_tmp =  emol / struct(i, :) 
     num = minloc(emol_tmp, 1, emol_tmp > 0.0d0)
     
     mol(i) = (emol(num) * CSeed) / struct(i, num)
     
!     if (present(mol_ini) .AND. mol_ini(i) > ErrA) then 
     if (present(mol_ini) .AND. mol_ini(i) > zero) then 
        mol(i) = mol_ini(i)
     end if
  end do not_use
  

!!!
!!!定数行列の作成. 元素モルから上記で与えた分の元素量を差し引く
!!!
  LOOP: do
     do e = 1, M
        B(e) = emol(e) - dot_product(struct(:, e), mol)
     end do
     
     !もしも B が負になったならば, mol に適当なファクターをかける
     !溶液のモル比を保存するため, mol 全部をファクター倍. 
     if (minval(B) > 0.0d0) then 
        exit LOOP
     else
        mol = mol * fct
     end if
  end do LOOP
  FB = B
  
  
!!!
!!!行列式を解き, 残りの化学種存在量を決める
!!!
  call ssl2_lax_init(M)
  call ssl2_lax(A, B)

  do e = 1, M
     mol(order(e)) = B(e)
  end do


!!!
!!!値が負になっていないかどうかのチェック
!!!  
  if (minval(mol) < 0.0d0) then 
     write(*,*) "ERROR, initmol.f90: mol < 0.0"
     call display(emol, "emol")
     call display(mol, "mol")
     call display(A, "A")
     call display(B, "B")     
     if (present(mol_ini)) display(mol_ini, "mol_ini")
     stop
  end if


!!!
!!!後かたづけ
!!!
  deallocate(A, B, emol_tmp, order)

end subroutine Initmol


