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

subroutine plot_cp
  use prm_common, only: property
  use interface_cp_noblegas
  use interface_cp_spl
  use interface_list_output
  
  implicit none
  integer               :: iws, i, s, M2
  integer, parameter    :: M1 = 10000
  real                  :: dt
  real(8)               :: temp, cp
  real(8), allocatable  :: temp_tmp(:)
  real(8), allocatable  :: cp_tmp(:)
  real, allocatable     :: temp1_map(:)
  real, allocatable     :: cp1_map(:)
  real, allocatable     :: temp2_map(:)
  real, allocatable     :: cp2_map(:)
  real, parameter       :: temp_high = 1.1d3
  real, parameter       :: temp_low = 2.0d1


  !グラフのオープン
  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, size(property, 1)
     if (property(s)%method /= 2) cycle plot
     
     allocate(temp1_map(M1), cp1_map(M1))
     temp1_map = 0.0D0; cp1_map = 0.0D0
     dt = (temp_high - temp_low) / (M1 - 1)
     
     do i = 1, M1
        !比熱の計算
        temp = temp_low + (i-1) * dt       
        
        if (property(s)%method == 1) then 
           call cp_noblegas(temp, cp=cp)
        elseif (property(s)%method == 2) then 
           call cp_spl(temp, property(s), cp=cp)
        else
           cycle plot
        end if
        
        temp1_map(i) = real(temp, 4)
        cp1_map(i) = real(cp, 4)
     end do
     
     
     !比熱のデータを呼び出す
     call list_output(property(s)%cp, temp=temp_tmp, data=cp_tmp)  

     M2 = size(cp_tmp, 1)
     allocate(temp2_map(M2), cp2_map(M2))
     temp2_map = real(temp_tmp, 4)
     cp2_map = real(cp_tmp, 4)
     
     ! グラフ作成
     call grfrm
     call usspnt(M1, temp1_map, cp1_map)
     call usspnt(M2, temp2_map, cp2_map)
     
     call grstrn(1)
     call uspfit
     call grstrf
     call ussttl('Temperature', 'K', 'Specific Heat', 'J/ K mol')
     call usdaxs
     
     call sgspli(25)
     call sgplu(M1, temp1_map, cp1_map)
     call uusmkt(5)
     call uusmki(3)
     call uumrk(M2, temp2_map, cp2_map)
     
     call uxsttl('T', trim(property(s)%name), 0.0)

     deallocate(temp1_map, cp1_map, temp2_map, cp2_map, temp_tmp, cp_tmp)

  end do plot


  !グラフのクローズ
  call grcls

end subroutine plot_cp