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

subroutine plot_svap
  use prm_common, only: property
  use interface_svap_amp
  use interface_svap_antoine 
  use interface_list_output

  implicit none
  real                  :: dt
  real(8)               :: temp
  real(8)               :: svap  
  real(8), allocatable  :: temp_tmp(:)
  real(8), allocatable  :: svap_tmp(:)
  real, parameter       :: temp_high = 1.0d3
  real, parameter       :: temp_low = 5.0d1
  real, allocatable     :: temp1_map(:)  
  real, allocatable     :: temp2_map(:)  
  real, allocatable     :: temp3_map(:)  
  real, allocatable     :: svap1_map(:)  
  real, allocatable     :: svap2_map(:)  
  real, allocatable     :: svap3_map(:)  
  real, allocatable     :: svap4_map(:)  
  integer               :: i, iws, s
  integer               :: N1, N2  
  integer               :: M2, M3
  integer               :: type
  integer, parameter    :: M1 = 10000


  !グラフのオープン
  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)
  CALL GLLSET( 'LMISS', .TRUE. )    !欠損値処理
  
  plot: do s = 1, size(property, 1)

     if (property(s)%method /= 3 .AND. property(s)%method /= 4) cycle plot
     
     allocate(temp1_map(M1), svap1_map(M1))
     temp1_map = 0.0D0; svap1_map = 0.0D0
     dt = (temp_high - temp_low) / (M1 - 1)
     
     !飽和蒸気圧の計算
     do i = 1, M1
        temp = temp_low + (i-1) * dt       
        
        if (property(s)%method == 3) then 
           call svap_amp(temp, property(s)%svap_cef, svap, type=type)
        elseif (property(s)%method == 4) then 
           call svap_antoine(temp, property(s)%svap_cef, svap, type=type)
        else
           cycle plot
        end if
        
        temp1_map(i) = real(temp, 4)
        svap1_map(i) = real(svap, 4)
     end do
     
     
     !データの取得
     call list_output(property(s)%svap, temp=temp_tmp, data=svap_tmp)  
     M2 = size(temp_tmp, 1)
     allocate(temp2_map(M2), svap2_map(M2))
     temp2_map = real(temp_tmp, 4)
     svap2_map = real(svap_tmp, 4)
     
     
     !データ点が存在する部分だけ取り出す
     N1 = minloc(temp1_map, 1, temp1_map >= temp2_map(1))
     N2 = maxloc(temp1_map, 1, temp1_map <= temp2_map(M2))
     M3 = N2 - N1 + 1
     allocate(temp3_map(M3), svap3_map(M3))
     temp3_map = temp1_map(N1:N2) 
     svap3_map = svap1_map(N1:N2) 
     
     
     !上限温度以下, 上限温度以上はカット
     N1 = minloc(temp1_map, 1, temp1_map >= property(s)%svap_cef%low%temp)
     N2 = maxloc(temp1_map, 1, temp1_map <= property(s)%svap_cef%high%temp)
     if (N2 == 0) N2 = M1
     allocate(svap4_map(M1))
     svap4_map = svap1_map
     svap4_map(1:N1)  = svap1_map(N1)
     svap4_map(N2:M1) = svap1_map(N2) 
     
     
     !グラフ作成
     call grfrm
     call usspnt(M1, temp1_map, svap1_map)
     call usspnt(M2, temp2_map, svap2_map)
     
     call grstrn(2)
     call uspfit
     call grstrf
     call ussttl('Temperature', 'K', 'Saturated pressure', 'Pa')
     call usdaxs
     
     call sgspli(25)
     call sgplu(M1, temp1_map, svap1_map)
     call uusmkt(5)
     call uusmki(3)
     call uumrk(M2, temp2_map, svap2_map)
     call uxsttl('T', trim(property(s)%name), 0.0)

     
     call grfrm
     call usspnt(M3, temp3_map, svap3_map)
     call usspnt(M2, temp2_map, svap2_map)

     call grstrn(2)
     call uspfit
     call grstrf
     call ussttl('Temperature', 'K', 'Saturated pressure', 'Pa')
     call usdaxs
     
     call sgspli(25)
     call sgplu(M3, temp3_map, svap3_map)
     call uusmkt(5)
     call uusmki(3)
     call uumrk(M2, temp2_map, svap2_map)
     call uxsttl('T', trim(property(s)%name), 0.0)     


!     call grfrm
!     call usspnt(M1, temp1_map, svap4_map)
!     call usspnt(M2, temp2_map, svap2_map)
!
!     call grstrn(2)
!     call uspfit
!     call grstrf
!     call ussttl('Temperature', 'K', 'Saturated pressure', 'Pa')
!     call usdaxs
     
!     call sgspli(25)
!     call sgplu(M1, temp1_map, svap4_map)
!     call uusmkt(5)
!     call uusmki(3)
!     call uumrk(M2, temp2_map, svap2_map)
!     call uxsttl('T', trim(property(s)%name), 0.0)     


     deallocate(temp1_map, temp2_map, temp3_map, temp_tmp, &
          & svap1_map, svap2_map, svap3_map, svap4_map, svap_tmp)

  end do plot


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

end subroutine plot_svap

