!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_entr0
  use prm_common, only: property
  use prm_phys, only: pres_ref, temp_ref
  use interface_entropy

  implicit none
  integer                 :: i, iws, s, N
  integer, parameter      :: M = 1000
  real                    :: low(2)
  real                    :: high(2)
  real                    :: limit(2)
  real, allocatable       :: entr0_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 = 400.0d0
  real(8), allocatable    :: entr0(:)
  real(8), allocatable    :: mol(:)
  logical                 :: sw1 
  logical                 :: sw2
  
  
!!!
!!!変数の初期化. 標準圧力での 1 mol での値をプロット. 
!!!
  N = size(property, 1)
  allocate(entr0_map(N, M), temp_map(M), mol(N))
  pres = pres_ref  
  mol = 1.0D0   !モル当たりの量を計算

  
!!!
!!!比熱の計算
!!!
  entr0_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 entropy(temp, pres, mol, property, entr0=entr0)

     temp_map(i) = real(temp, 4)
     entr0_map(:,i) = real(entr0, 4)

     deallocate(entr0)
  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
     sw1 = property(s)%ref%entr /= 0.0d0    
     sw2 = 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(entr0_map(s,:)) 
     limit(2) = maxval(entr0_map(s,:)) 

     call grfrm
     call usspnt(M, temp_map, entr0_map(s, :))
     if(sw1) call usspnt(1, real(temp_ref,4), real(property(s)%ref%entr, 4))
     if(sw2) call usspnt(2, low, limit)
     if(sw2) call usspnt(2, high, limit)
     
     call grstrn(1)
     call uspfit
     call grstrf
     call ussttl('Temperature', 'K', 'Entropy', 'J/ K mol')
     call usdaxs
     
     call sgspli(25)
     call sgplu(M, temp_map, entr0_map(s, :))
     if(sw1) call uusmkt(5)
     if(sw1) call uusmki(3)
     if(sw1) call uumrk(1, real(temp_ref,4), real(property(s)%ref%entr, 4))
     if(sw2) call sgspli(41)
     if(sw2) call sgplu(2, low, limit)
     if(sw2) call sgspli(41)
     if(sw2) call sgplu(2, high, limit)
     
     call uxsttl('T', trim(property(s)%name), 0.0)
  end do plot


!!!
!!!グラフのクローズ
!!!
    
  call grcls



end subroutine Plot_entr0


