!
!Description:
!  出力結果を表示するためのプログラム
!
! Current Code Owner: sugiyama@gfd-dennou.org
!
! Code Description: 
!   Language:          Fortran 90
!   Software Standaeds: "気象庁 標準コーディング"
!
! Modules used: 
!
! Copyright (C) SUGIYAMA Ko-ichiro, 2004, All rights reserved

program OboroPlotMinimum2
  use dcl
  use interface_netcdf_open
  use interface_netcdf_close
  use interface_netcdf_get
  use interface_display

  implicit none
  real(8), allocatable        :: pres_tmp(:)       !圧力
  real, allocatable           :: pres_plot(:)      !圧力
  real(8), allocatable        :: vstab_tmp(:)
  real(8), allocatable        :: vstab_low_tmp(:)
  real(8), allocatable        :: vstab_high_tmp(:)
  real, allocatable           :: vstab_plot(:,:)  
  real(8), allocatable        :: dtdz_tmp(:)
  real(8), allocatable        :: dtdz_low_tmp(:)
  real(8), allocatable        :: dtdz_high_tmp(:)
  real, allocatable           :: dtdz_plot(:,:)  
  real(8), allocatable        :: vdtdz_tmp(:)
  real(8), allocatable        :: vdtdz_low_tmp(:)
  real(8), allocatable        :: vdtdz_high_tmp(:)
  real, allocatable           :: vdtdz_plot(:,:)  
  integer                     :: i, s
  integer                     :: x_num
  integer                     :: y_num
  integer                     :: z_num
  integer                     :: ncid
  integer                     :: fnum
  integer                     :: index
  integer, allocatable        :: x_nums(:)
  character(80)               :: fname
  character(20), allocatable  :: sname_tmp(:)
  character(20), allocatable  :: sname_plot(:,:)
  character(80)               :: name
  character(1)                :: replay
  logical                     :: disp


  write(*,*) "============================================"
  write(*,*) "======NetCDF ファイルを指定して下さい======="
  write(*,*) "============================================"
  read(*,*) fname


!!!
!!!NetCDF ファイルのオープン
!!!次元は pres, spc, elm であることを知っている
!!!
  !ファイルのオープン
  call netcdf_open(trim(fname), ncid)
  
  !変数の取得. 取り出す変数は pres, spc, elm, temp, mol, molfr 
  name = 'pres'
  call netcdf_get(ncid, name, pres_tmp)
  
  name = 'vstab'
  call netcdf_get(ncid, name, vstab_tmp) 
     
  name = 'vstab_low'
  call netcdf_get(ncid, name, vstab_low_tmp) 
     
  name = 'vstab_high'
  call netcdf_get(ncid, name, vstab_high_tmp) 
  
  name = 'dtdz'
  call netcdf_get(ncid, name, dtdz_tmp) 
     
  name = 'dtdz_low'
  call netcdf_get(ncid, name, dtdz_low_tmp) 
     
  name = 'dtdz_high'
  call netcdf_get(ncid, name, dtdz_high_tmp) 
  
  name = 'vdtdz'
  call netcdf_get(ncid, name, vdtdz_tmp) 
     
  name = 'vdtdz_low'
  call netcdf_get(ncid, name, vdtdz_low_tmp) 
     
  name = 'vdtdz_high'
  call netcdf_get(ncid, name, vdtdz_high_tmp) 
  
  !ファイルを閉じる
  call netcdf_close(ncid)
  
  !配列の割り当て
  x_num = size(pres_tmp,1)

  allocate( pres_plot(x_num), vstab_plot(x_num,3), &
       dtdz_plot(x_num,3), vdtdz_plot(x_num,3) ) 
  pres_plot  = 0.0d0
  vstab_plot = 0.0d0 
  dtdz_plot = 0.0d0 
  vdtdz_plot = 0.0d0 

  
  pres_plot(:)    = real(pres_tmp,4)
  vstab_plot(:,1) = real(vstab_tmp,4)
  vstab_plot(:,2) = real(vstab_low_tmp,4)
  vstab_plot(:,3) = real(vstab_high_tmp,4)
  dtdz_plot(:,1)  = real(dtdz_tmp,4)
  dtdz_plot(:,2)  = real(dtdz_low_tmp,4)
  dtdz_plot(:,3)  = real(dtdz_high_tmp,4)
  vdtdz_plot(:,1) = real(vdtdz_tmp,4)
  vdtdz_plot(:,2) = real(vdtdz_low_tmp,4)
  vdtdz_plot(:,3) = real(vdtdz_high_tmp,4)
  
  where (vstab_plot <= 1.0d-8)
     vstab_plot = 1.0d-8
  end where
  
  deallocate(pres_tmp, &
       & vstab_tmp, vstab_low_tmp, vstab_high_tmp, &
       & dtdz_tmp,  dtdz_low_tmp,  dtdz_high_tmp, &
       & vdtdz_tmp, vdtdz_low_tmp, vdtdz_high_tmp )


!!!
!!!プロット
!!!  
  !グラフのオープン
  call DclOpenGraphics()
  call DclSetParm('WINDOW_HEIGHT', 600)
  call DclSetParm('WINDOW_WIDTH', 600)
  call DCLSetParm('DRAW_CORNERMARK', .false.)
  call DclSetFrameMargin( 0.05, 0.05, 0.05, 0.05 )
  

  call DclNewFrame
  call DclSetWindow(&
       & minval(dtdz_plot, dtdz_plot > 0.0d0), &
       & maxval(dtdz_plot, dtdz_plot > 0.0d0), &
       & maxval(pres_plot, pres_plot > 0.0d0), &
       & minval(pres_plot, pres_plot > 0.0d0))
  call DclSetViewPort( 0.2, 0.8, 0.2, 0.8 )
  call DclSetTransNumber(2)
  call DclSetTransFunction
  
  call DclSetTitle('dtdz', 'pres', 'K/km', 'Pa')
  call DclDrawScaledAxis

  do i = 1, 3
     index = i * 10 + 2
     call DclDrawLine(dtdz_plot(:,i), pres_plot, index=index) 
  end do
  call DclDrawTitle( 'T', 'lapse rate', sw=2 )
  

  call DclNewFrame
  call DclSetWindow(&
       & minval(vdtdz_plot, vdtdz_plot > 0.0d0), &
       & maxval(vdtdz_plot, vdtdz_plot > 0.0d0), &
       & maxval(pres_plot, pres_plot > 0.0d0), &
       & minval(pres_plot, pres_plot > 0.0d0))
  call DclSetViewPort( 0.2, 0.8, 0.2, 0.8 )
  call DclSetTransNumber(2)
  call DclSetTransFunction
  
  call DclSetTitle('vdtdz', 'pres', 'K/km', 'Pa')
  call DclDrawScaledAxis

  do i = 1, 3
     index = i * 10 + 2
     call DclDrawLine(vdtdz_plot(:,i), pres_plot, index=index) 
  end do
  call DclDrawTitle( 'T', 'virtual lapse rate', sw=2 )
  

  !安定成層 (log-scaler plot)
  call DclNewFrame
  call DclSetWindow(&
       & minval(vstab_plot, vstab_plot > 0.0d0), &
       & maxval(vstab_plot, vstab_plot > 0.0d0), &
       & maxval(pres_plot, pres_plot > 0.0d0), &
       & minval(pres_plot, pres_plot > 0.0d0))
  call DclSetViewPort( 0.2, 0.8, 0.2, 0.8 )
  call DclSetTransNumber(2)
  call DclSetTransFunction
  
  call DclSetTitle('vstab', 'pres', '', 'Pa')
  call DclDrawScaledAxis

  do i = 1, 3
     index = i * 10 + 2
     call DclDrawLine(vstab_plot(:,i), pres_plot, index=index) 
  end do
  call DclDrawTitle( 'T', 'stability', sw=2 )
  

  !安定成層 (log-log plot)
  call DclNewFrame
  call DclSetWindow(&
       & minval(vstab_plot, vstab_plot > 0.0d0), &
       & maxval(vstab_plot, vstab_plot > 0.0d0), &
       & maxval(pres_plot, pres_plot > 0.0d0), &
       & minval(pres_plot, pres_plot > 0.0d0))
  call DclSetViewPort( 0.2, 0.8, 0.2, 0.8 )
  call DclSetTransNumber(4)
  call DclSetTransFunction
  
  call DclSetTitle('vstab', 'pres', '', 'Pa')
  call DclDrawScaledAxis

  do i = 1, 3
     index = i * 10 + 2
     call DclDrawLine(vstab_plot(:,i), pres_plot, index=index) 
  end do
  call DclDrawTitle( 'T', 'stability', sw=2 )

  
  !グラフのクローズ
  call DCLCloseGraphics
  

end program OboroPlotMinimum2
