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

subroutine plot_mu0
  use prm_common, only: property
  use prm_phys, only: pres_ref, temp_ref
  use interface_gibbs_energy
  use interface_name2num

  implicit none
  integer                 :: i, iws, s, N, l, k1, k2
  integer, parameter      :: M = 10000
  real                    :: low(2)
  real                    :: high(2)
  real                    :: limit(2)
  real, allocatable       :: mu0_map(:,:)
  real, allocatable       :: temp_map(:)
  real(8)                 :: temp
  real(8)                 :: pres
  real(8)                 :: dt
  real(8), parameter      :: temp_low = 10.0d0
  real(8), parameter      :: temp_high = 300.0d0
  real(8), allocatable    :: mu0(:)
  real(8), allocatable    :: mol(:)
  character(20)           :: name
  logical                 :: sw
  
  
!!!
!!!変数の初期化. 標準圧力での 1 mol での値をプロット. 
!!!
  N = size(property, 1)
  allocate(mu0_map(N, M), temp_map(M), mol(N))
  pres = pres_ref  
  mol = 1.0D0   !モル当たりの量を計算

  
!!!
!!!比熱の計算
!!!
  mu0_map = 0.0D0; temp_map = 0.0D0
  dt = (temp_high - temp_low) / (M - 1)

  do i = 1, M
     temp = temp_low + (i - 1) * dt
     
     call gibbs_energy(temp, pres, mol, property, mu0=mu0)

     temp_map(i) = real(temp, 4)
     mu0_map(:,i) = real(mu0, 4)

     deallocate(mu0)
  end do


!!!
!!!グラフのオープン
!!!
  write(*,*) "workstation ID (i) ?"
  call sgpwsn
  read(*,*) iws
  call sglset('LCORNER', .FALSE.)
  call swiset('IWIDTH', 600)
  call swiset('IHEIGHT', 600)
  call gropn(iws)
  call slmgn(0.05, 0.05, 0.01, 0.01)
  
  
!!!  
!!! グラフ作成
!!!
  plot: do s = 1, N
     sw = property(s)%method == 3 &
          & .OR. property(s)%method == 4 .OR. property(s)%method == 5
     
     low(1) = real(property(s)%svap_cef%low%temp,4)
     low(2) = real(property(s)%svap_cef%low%temp,4) + dt
     high(1) = real(property(s)%svap_cef%high%temp,4)
     high(2) = real(property(s)%svap_cef%high%temp,4) + dt
     limit(1) = minval(mu0_map(s,:)) 
     limit(2) = maxval(mu0_map(s,:)) 
     
     call grfrm
     call usspnt(M, temp_map, mu0_map(s, :))
     if(sw) call usspnt(2, low, limit)
     if(sw) call usspnt(2, high, limit)
     
     call grstrn(1)
     call uspfit
     call grstrf
     call ussttl('Temperature', 'K', 'Mu0', 'J/ K mol')
     call usdaxs
     
     call sgspli(25)
     call sgplu(M, temp_map, mu0_map(s, :))
     if(sw) call sgspli(41)
     if(sw) call sgplu(2, low, limit)
     if(sw) call sgspli(41)
     if(sw) call sgplu(2, high, limit)
     
     call uxsttl('T', trim(property(s)%name), 0.0)
  end do plot


  do s = 1, N
     name = property(s)%name
     l = len_trim(name)
     write(*,*) name
     if (name(l:l) /= "g") cycle

     name = name(1:l-1) // 'l'
     call name2num(name, property%name, k1)
     
     name = name(1:l-1) // 's'
     call name2num(name, property%name, k2)

     write(*,*) k1, k2
     if (k1 == 0 .AND. k2 == 0) cycle

     call grfrm
     if(k1 /= 0) call usspnt(M, temp_map, mu0_map(k1, :))
     if(k2 /= 0) call usspnt(M, temp_map, mu0_map(k2, :))
     
     call grstrn(1)
     call uspfit
     call grstrf
     call ussttl('Temperature', 'K', 'Mu0', 'J/ K mol')
     call usdaxs
     
     call sgspli(25)
     if(k1 /= 0) call sgplu(M, temp_map, mu0_map(k1, :))
     if(k1 /= 0) call sgplu(M, temp_map, mu0_map(k2, :))
     
     call uxsttl('T', trim(property(s)%name), 0.0)
  end do
  
  
!!!
!!!グラフのクローズ
!!!
  call grcls

end subroutine Plot_mu0