!Description:
!  化学種名から含まれる元素名とその個数を判定する
!
! Current Code Owner: sugiyama@gfd-dennou.org
!
! History: 
!   Version    Date          Comment
!   -------    ----------    -------
!   1.0        2003-01-19    新規作成
!   1.1        2003-01-23    原子記号を原子番号に変換するように
!   1.2        2003-01-24    インターフェイス文も含める
!
! Copyright (C) SUGIYAMA Ko-ichiro, 2003, All rights reserved

!!!
!!!インターフェイス文
!!!
module interface_mol2atm
  interface
     subroutine mol2atm(sname, anum, enum)
       character(40), intent(in)            :: sname
       integer, allocatable, intent(out)    :: anum(:)
       real(8), allocatable, intent(out)    :: enum(:)
     end subroutine mol2atm
  end interface
end module interface_mol2atm


!!!
!!!サブルーチン本文
!!!
subroutine mol2atm(sname, anum, enum)
  use prm_chem, only: atmsymbol, amax

  implicit none
  character(40), intent(in)            :: sname
  integer, allocatable, intent(out)    :: anum(:)
  real(8), allocatable, intent(out)    :: enum(:)

  integer, allocatable       :: anum_tmp(:)
  real(8), allocatable       :: enum_tmp(:)
  integer, parameter         :: M = 20
  integer                    :: num1, num2, num3
  character(2), allocatable  :: symbol(:)
  logical, allocatable       :: sw(:)


!!!
!!!変数の初期化
!!!
  allocate(anum_tmp(M), enum_tmp(M), symbol(M), sw(M))
  anum_tmp = 0
  enum_tmp = 0.0d0
  sw = .false. 
  num1 = 1
  num2 = 1


!!!
!!!元素名とその数を化学種名から判定する. 
!!!
  do num1 = 1, len_trim(sname)

     if (ichar('A') <= ichar(sname(num1:num1)) &
          & .AND. ichar(sname(num1:num1)) <= ichar('Z')) then 
        symbol(num2) = sname(num1:num1)
        enum_tmp(num2) = 1.0d0
        num2 = num2 + 1

     elseif (ichar('a') <= ichar(sname(num1:num1)) &
          & .AND. ichar(sname(num1:num1)) <= ichar('z')) then 
        symbol(num2-1) = sname(num1-1:num1)
        
     else
        symbol(num2-1) = sname(num1-1:num1-1)
        read(sname(num1:num1), *) enum_tmp(num2-1)  

     end if
  end do


!!!
!!!元素名に重複がないかの判定
!!!
  num3 = minloc(enum_tmp, 1) - 1
  sw(1: num3) = .true.

  do num1 = 1, num3 - 1
     do num2 = num1 + 1, num3
        
        if (symbol(num1) == symbol(num2)) then 
           enum_tmp(num1) = enum_tmp(num1) + enum_tmp(num2) 
           sw(num2) = .false. 
        end if

     end do
  end do


!!!
!!!元素名から原子番号を求める
!!!
  do num1 = 1, M 
     atm: do num2 = 1, amax
        if(symbol(num1) == atmsymbol(num2)) then 
           anum_tmp(num1) = num2
           exit atm
        end if
     end do atm
  end do


!!!
!!!出力用の変数を作成
!!!
  allocate(anum(count(sw)), enum(count(sw)))
  anum = pack(anum_tmp, sw)
  enum = pack(enum_tmp, sw)


!!!
!!!後かたづけ
!!!
  deallocate(anum_tmp, enum_tmp, symbol, sw)

end subroutine mol2atm
