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

program OboroPlot  
  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        :: entr_tmp(:)       !エントロピー
  real, allocatable           :: entr_plot(:,:)       !エントロピー
  real(8), allocatable        :: temp_tmp(:)       !温度
  real, allocatable           :: temp_plot(:,:)       !温度
  real(8), allocatable        :: mol_tmp(:,:)  
  real(8), allocatable        :: mol(:,:)
  real, allocatable           :: mol_plot(:,:,:)
  real(8), allocatable        :: molfr_tmp(:,:)
  real, allocatable           :: molfr_plot(:,:,:)
  real(8), allocatable        :: emol_tmp(:,:)
  real, allocatable           :: emol_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        :: cldns_tmp(:,:)
  real, allocatable           :: cldns_plot(:,:,:)
  real, allocatable           :: cldns_all_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 = 'entr'
     call netcdf_get(ncid, name, entr_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 = 'cldns'
     call netcdf_get(ncid, name, cldns_tmp) 
     name = 'vstab'
     call netcdf_get(ncid, name, vstab_tmp) 

     name = 'mol'
     call netcdf_get(ncid, name, mol_tmp)
     name = 'molfr'
     call netcdf_get(ncid, name, molfr_tmp)
     name = 'emol'
     call netcdf_get(ncid, name, emol_tmp) 

     name = 'sname'
     call netcdf_get(ncid, name, sname_tmp) 

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

     !配列の割り当て
!     x_num = size(mol_tmp,1)
     x_nums(i) = size(mol_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), &
             & sname_plot(y_num,fnum), &
             & entr_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), &
             & cldns_all_plot(x_num,fnum), &
             & cldns_plot(x_num,y_num,fnum), &
             & mol(x_num,y_num), &
             & mol_plot(x_num,y_num,fnum), &
             & molfr_plot(x_num,y_num,fnum), &       
             & emol_plot(x_num,z_num,fnum), &
             & false(x_num))
        pres_plot  = 0.0d0
        sname_plot = "" 
        entr_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 
        cldns_all_plot  = 0.0d0 
        cldns_plot = 0.0d0 
        mol        = 0.0d0 
        mol_plot   = 0.0d0 
        molfr_plot = 0.0d0        
        emol_plot  = 0.0d0
     end if

!     if (size(pres_plot, 1) /= x_num) then 
!        write(*,*) "data point is not same"
!        stop
!     end if

     pres_plot(1:x_nums(i),i)    = real(pres_tmp,4)
     sname_plot(:,i)             = sname_tmp
     entr_plot(1:x_nums(i),i)    = real(entr_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)
!     mol_plot(1:x_nums(i),:,i)   = real(mol_tmp,4)
     molfr_plot(1:x_nums(i),:,i) = real(molfr_tmp,4)
     emol_plot(1:x_nums(i),:,i)  = real(emol_tmp,4)
     cldns_plot(1:x_nums(i),:,i) = real(cldns_tmp,4)
     cldns_all_plot(1:x_nums(i),i) = real(sum(cldns_tmp, 2),4)

     where (cldns_all_plot <= 1.0d-6)
        cldns_all_plot = 1.0d-6
     end where
     
     where (vstab_plot <= 1.0d-8)
        vstab_plot = 1.0d-8
     end where    
     
     do s = 1, y_num
        if (maxval(mol_tmp(:,s)) == minval(mol_tmp(:,s))) then
           mol_plot(1:x_nums(i),s,i) =  real(mol_tmp(:,s), 4)
           cycle
        end if
        zero = minval(mol_tmp(:,s))
        min  = minval(mol_tmp(:,s), mol_tmp(:,s) > zero)
        write(*,*) zero        
        write(*,*) min

        false = real(min, 4)
        mol(1:x_nums(i),s) =  merge(mol_tmp(:,s), &
             & false(1:x_nums(i)), mol_tmp(:,s) /= zero)
        mol_plot(1:x_nums(i),s,i) =  real(mol(1:x_nums(i),s), 4)
     end do

     deallocate(pres_tmp, sname_tmp, temp_tmp, vtemp_tmp, entr_tmp, &
          & dtdz_tmp, vdtdz_tmp, cldns_tmp, vstab_tmp, &
          & mol_tmp, molfr_tmp, emol_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 DclNewFrame
  call DclSetWindow(&
       & minval(cldns_all_plot, cldns_all_plot > 0.0d0), &
       & maxval(cldns_all_plot, cldns_all_plot > 0.0d0), &
       & maxval(pres_plot, pres_plot > 0.0d0), &
       & minval(pres_plot, pres_plot > 0.0d0))
!  call DclSetWindow(minval(cldns_all_plot)*0.999, &
!       & maxval(cldns_all_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('cldns all', 'pres', 'g/l', 'Pa')
  call DclDrawScaledAxis

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

  !エントロピー
  call DclNewFrame
  call DclSetWindow(&
       & minval(entr_plot, entr_plot > 0.0d0), &
       & maxval(entr_plot, entr_plot > 0.0d0), &
       & maxval(pres_plot, pres_plot > 0.0d0), &
       & minval(pres_plot, pres_plot > 0.0d0))
!  call DclSetWindow(minval(entr_plot)*0.999, maxval(entr_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('entr', 'pres', 'J/K mol', 'Pa')
  call DclDrawScaledAxis

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

  !モル数変化
  do s = 1, y_num
     write(*,*) sname_plot(s,1)
     write(*,*) minval(mol_plot(:,s,:), mol_plot(:,s,:) > 0.0d0)
     write(*,*) maxval(mol_plot(:,s,:), mol_plot(:,s,:) > 0.0d0)
     if (minval(mol_plot(:,s,:), mol_plot(:,s,:) > 0.0d0) &
          & == maxval(mol_plot(:,s,:), mol_plot(:,s,:) > 0.0d0)) cycle

     if (minval(mol_plot(:,s,:), mol_plot(:,s,:) > 0.0d0) * 1.5d0 &
          & >= maxval(mol_plot(:,s,:), mol_plot(:,s,:) > 0.0d0)) cycle
     
     call DclNewFrame
     call DclSetWindow(&
          & minval(mol_plot(:,s,:), mol_plot(:,s,:) > 0.0d0), &
          & maxval(mol_plot(:,s,:), mol_plot(:,s,:) > 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('mol', 'pres', 'mol', 'Pa')
     call DclDrawScaledAxis
     
     do i = 1, fnum
        index = i * 10 + 2
        call DclDrawLine(mol_plot(1:x_nums(i),s,i), pres_plot(1:x_nums(i),i), index=index) 
     end do
     call DclDrawTitle( 'T', trim(sname_plot(s,1)), sw=2 )
  end do


  !グラフのクローズ
  call DCLCloseGraphics
  

end program OboroPlot

