!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_dvt_int
  use prm_phys, only: temp_ref
  use prm_common, only: property
  use interface_cp_noblegas
  use interface_cp_spl
  use interface_list_output
  
  implicit none
  integer               :: iws, i, s, M2, L
  integer, parameter    :: M1 = 10000
  real                  :: dt
  real(8)               :: temp, cp_dvt_int
  real(8), allocatable  :: temp_tmp(:)
  real(8), allocatable  :: cp_dvt_int_tmp(:)
  real, allocatable     :: temp1_map(:)
  real, allocatable     :: cp_dvt_int1_map(:)
  real, allocatable     :: temp2_map(:)
  real, allocatable     :: cp_dvt_int2_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), cp_dvt_int1_map(M1))
     temp1_map = 0.0D0; cp_dvt_int1_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_dvt_int=cp_dvt_int)
        elseif (property(s)%method == 2) then 
           call cp_spl(temp, property(s), cp_dvt_int=cp_dvt_int)
        else
           cycle plot
        end if
        
        temp1_map(i) = real(temp, 4)
        cp_dvt_int1_map(i) = real(cp_dvt_int, 4)
     end do
     
     
     !比熱のデータを呼び出す
     call list_output(property(s)%spl, temp=temp_tmp, dat_dvt_int=cp_dvt_int_tmp)
     M2 = size(cp_dvt_int_tmp, 1)
     allocate(temp2_map(M2), cp_dvt_int2_map(M2))

     do i = 1, M2
        L = maxloc(temp_tmp, 1, temp_tmp <= temp_ref)
        temp2_map = real(temp_tmp, 4)
        if (i <= L) then 
           cp_dvt_int2_map(i) = - real(sum(cp_dvt_int_tmp(i:L-1)), 4)
        else
           cp_dvt_int2_map(i) = real(sum(cp_dvt_int_tmp(L:i-1)), 4)
        end if
     end do
     
     
     ! グラフ作成
     call grfrm
     call usspnt(M1, temp1_map, cp_dvt_int1_map)
     call usspnt(M2, temp2_map, cp_dvt_int2_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, cp_dvt_int1_map)
     call uusmkt(5)
     call uusmki(3)
     call uumrk(M2, temp2_map, cp_dvt_int2_map)
     
     call uxsttl('T', trim(property(s)%name), 0.0)

     deallocate(temp1_map, cp_dvt_int1_map, temp2_map, cp_dvt_int2_map, temp_tmp, cp_dvt_int_tmp)

  end do plot


  !グラフのクローズ
  call grcls

end subroutine plot_cp_dvt_int