!
!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 OboroPlotMinimum
  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        :: temp_tmp(:)       !温度
  real, allocatable           :: temp_plot(:,:)       !温度
  real(8), allocatable        :: dtdz_tmp(:)
  real, allocatable           :: dtdz_plot(:,:)
  real(8), allocatable        :: vdtdz_tmp(:)
  real, allocatable           :: vdtdz_plot(:,:)
  real(8), allocatable        :: vtemp_tmp(:)
  real, allocatable           :: vtemp_plot(:,:)
  real(8), allocatable        :: vstab_tmp(:)
  real, allocatable           :: vstab_plot(:,:)  
  real, allocatable           :: n_plot(:,:)  
  real(8), allocatable        :: false(:)
  real(8)                     :: zero
  real(8)                     :: min
  integer                     :: i, s
  integer                     :: x_num
  integer                     :: y_num
  integer                     :: z_num
  integer                     :: ncid
  integer                     :: fnum
  integer                     :: index
  integer, allocatable        :: x_nums(:)
  character(80), allocatable  :: fname(:)
  character(20), allocatable  :: sname_tmp(:)
  character(20), allocatable  :: sname_plot(:,:)
  character(80)               :: name
  character(1)                :: replay
  logical                     :: disp


  write(*,*) "============================================"
  write(*,*) "======ファイル数を入力して下さい============"
  write(*,*) "============================================"
  read(*,*) fnum
  allocate(fname(fnum), x_nums(fnum))

  write(*,*) "============================================"
  write(*,*) "======NetCDF ファイルを指定して下さい======="
  write(*,*) "============================================"
!  fname(1) = "obradbt.nc"
  do i = 1, fnum
     read(*,*) fname(i) 
  end do

  write(*,*) "============================================"
  write(*,*) "======図を 1 枚毎に出力しますか?============"
  write(*,*) "====== say Y or N =========================="
  write(*,*) "============================================"
  read(*,*) replay
  if (replay == "Y" .OR. replay == "y") then 
     disp = .false.
  elseif (replay == "N" .OR. replay == "n") then 
     disp = .true.
  else
     write(*,*) "Unknown Answer: ", replay
     stop
  end if
  

!!!
!!!NetCDF ファイルのオープン
!!!次元は pres, spc, elm であることを知っている
!!!
  do i = 1, fnum
     !ファイルのオープン
     call netcdf_open(trim(fname(i)), ncid)

     !変数の取得. 取り出す変数は pres, spc, elm, temp, mol, molfr 
     name = 'pres'
     call netcdf_get(ncid, name, pres_tmp)

     name = 'temp'
     call netcdf_get(ncid, name, temp_tmp)
     name = 'vtemp'
     call netcdf_get(ncid, name, vtemp_tmp)

     name = 'dtdz'
     call netcdf_get(ncid, name, dtdz_tmp) 
     name = 'vdtdz'
     call netcdf_get(ncid, name, vdtdz_tmp) 

     name = 'vstab'
     call netcdf_get(ncid, name, vstab_tmp) 

     !ファイルを閉じる
     call netcdf_close(ncid)

     !配列の割り当て
!     x_num = size(mol_tmp,1)
     x_nums(i) = size(temp_tmp,1)
     x_num = 100000
!     y_num = size(mol_tmp,2)
!     z_num = size(emol_tmp,2)

     if (.not. allocated(pres_plot)) then 
        allocate(pres_plot(x_num,fnum), &
             & temp_plot(x_num,fnum), &
             & vtemp_plot(x_num,fnum), &
             & dtdz_plot(x_num,fnum), &
             & vdtdz_plot(x_num,fnum), &
             & vstab_plot(x_num,fnum), &
             & n_plot(x_num,fnum), &
             & false(x_num))
        pres_plot  = 0.0d0
        temp_plot  = 0.0d0 
        vtemp_plot = 0.0d0 
        dtdz_plot  = 0.0d0
        vdtdz_plot = 0.0d0 
        vstab_plot = 0.0d0 
        n_plot     = 0.0d0 
     end if

     pres_plot(1:x_nums(i),i)    = real(pres_tmp,4)
     temp_plot(1:x_nums(i),i)    = real(temp_tmp,4)
     vtemp_plot(1:x_nums(i),i)   = real(vtemp_tmp,4)
     dtdz_plot(1:x_nums(i),i)    = real(dtdz_tmp,4)
     vdtdz_plot(1:x_nums(i),i)   = real(vdtdz_tmp,4)
     vstab_plot(1:x_nums(i),i)   = real(vstab_tmp,4)
     
     where (vstab_plot <= 1.0d-8)
        vstab_plot = 1.0d-8
     end where    
     
     deallocate(pres_tmp, temp_tmp, vtemp_tmp,  &
          & dtdz_tmp, vdtdz_tmp, vstab_tmp )
  end do

  !N2 --> N へ変換
  n_plot = sqrt(vstab_plot)



!!!
!!!プロット
!!!  
  !グラフのオープン
  call DclOpenGraphics()
  if (disp) call DclDivideFrame( 'y', 3, 2 )
  if (.NOT. disp) call DclSetParm('WINDOW_HEIGHT', 600)
  if (.NOT. disp) call DclSetParm('WINDOW_WIDTH', 600)
  if (.NOT. disp) call DCLSetParm('DRAW_CORNERMARK', .false.)
  call DclSetFrameMargin( 0.05, 0.05, 0.05, 0.05 )
  
  !温度
  call DclNewFrame
  call DclSetWindow(&
       & minval(temp_plot, temp_plot > 0.0d0), &
       & maxval(temp_plot, temp_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('temp', 'pres', 'K', 'Pa')
  call DclDrawScaledAxis

  do i = 1, fnum
     index = i * 10 + 2
     call DclDrawLine(temp_plot(1:x_nums(i),i), pres_plot(1:x_nums(i),i), index=index) 
  end do
  call DclDrawTitle( 'T', 'temperature', sw=2 )
  

  !仮温度
  call DclNewFrame
  call DclSetWindow(&
       & minval(vtemp_plot, vtemp_plot > 0.0d0), &
       & maxval(vtemp_plot, vtemp_plot > 0.0d0), &
       & maxval(pres_plot, pres_plot > 0.0d0), &
       & minval(pres_plot, pres_plot > 0.0d0))
!  call DclSetWindow(minval(vtemp_plot), maxval(vtemp_plot), &
!       & maxval(pres_plot), minval(pres_plot))
  call DclSetViewPort( 0.2, 0.8, 0.2, 0.8 )
  call DclSetTransNumber(2)
  call DclSetTransFunction
     
  call DclSetTitle('vtemp', 'pres', 'K', 'Pa')
  call DclDrawScaledAxis

  do i = 1, fnum
     index = i * 10 + 2
     call DclDrawLine(vtemp_plot(1:x_nums(i),i), pres_plot(1:x_nums(i),i), index=index) 
  end do
  call DclDrawTitle( 'T', 'virtual temperature', sw=2 )
  

  !断熱減率
  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 DclSetWindow(minval(dtdz_plot)*0.999, maxval(dtdz_plot)*1.001, &
!       & maxval(pres_plot)*1.001, minval(pres_plot)*0.999)
  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, fnum
     index = i * 10 + 2
     call DclDrawLine(dtdz_plot(1:x_nums(i),i), pres_plot(1:x_nums(i),i), 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 DclSetWindow(minval(vdtdz_plot)*0.999, maxval(vdtdz_plot)*1.001, &
!       & maxval(pres_plot)*1.001, minval(pres_plot)*0.999)
  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, fnum
     index = i * 10 + 2
     call DclDrawLine(vdtdz_plot(1:x_nums(i),i), pres_plot(1:x_nums(i),i), 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 DclSetWindow(minval(vstab_plot)*0.999, maxval(vstab_plot)*1.001, &
!       & maxval(pres_plot)*1.001, minval(pres_plot)*0.999)
  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, fnum
     index = i * 10 + 2
     call DclDrawLine(vstab_plot(1:x_nums(i),i), pres_plot(1:x_nums(i),i), 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 DclSetWindow(minval(vstab_plot)*0.999, maxval(vstab_plot)*1.001, &
!       & maxval(pres_plot)*1.001, minval(pres_plot)*0.999)
  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, fnum
     index = i * 10 + 2
     call DclDrawLine(vstab_plot(1:x_nums(i),i), pres_plot(1:x_nums(i),i), index=index) 
  end do
  call DclDrawTitle( 'T', 'stability', sw=2 )
  
 
  !ブラントバイサラ振動数
  call DclNewFrame
  call DclSetWindow(&
       & minval(n_plot, n_plot > 0.0d0), &
       & maxval(n_plot, n_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('n', 'pres', '', 'Pa')
  call DclDrawScaledAxis

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

  !ブラントバイサラ振動数 (log-log plot)
  call DclNewFrame
  call DclSetWindow(&
       & minval(n_plot, n_plot > 0.0d0), &
       & maxval(n_plot, n_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('n', 'pres', '', 'Pa')
  call DclDrawScaledAxis

  do i = 1, fnum
     index = i * 10 + 2
     call DclDrawLine(n_plot(1:x_nums(i),i), pres_plot(1:x_nums(i),i), index=index) 
  end do
  call DclDrawTitle( 'T', 'stability', sw=2 )
  
  !グラフのクローズ
  call DCLCloseGraphics
  

end program OboroPlotMinimum

