!
!Description:
!  化学ポテンシャルをプロットして平衡状態を調べるためのプログラム
!
! Current Code Owner: sugiyama@gfd-dennou.org
!
! Histry: 
!
! Version    Date             Comment
! -------    ----------       --------
! 1.0        2004-07-26       新規作成
!
! Code Description: 
!   Language:          Fortran 90
!   Software Standaeds: "気象庁 標準コーディング"
!
! Modules used: 
!
! Copyright (C) SUGIYAMA Ko-ichiro, 2000, All rights reserved

program ObrChkMu  
  use dcl
  use interface_netcdf_open
  use interface_netcdf_close
  use interface_netcdf_get
  use interface_display  
  use interface_name_cnd2gas
  use interface_name2num

  implicit none
  real(8), allocatable        :: mu_tmp(:,:) 
  real(8), allocatable        :: mu(:,:)
  real, allocatable           :: mu_plot(:,:)
  real, allocatable           :: mu_err_plot(:,:)
  real(8), allocatable        :: pres(:) 
  real, allocatable           :: pres_plot(:)
  real(8), allocatable        :: vap(:)
  real(8), allocatable        :: svap(:)
  real(8), allocatable        :: molfr(:,:) 
  real(8), allocatable        :: molfr_tmp(:,:) 
  real, allocatable           :: molfr_plot(:,:) 
  real(8), allocatable        :: mol(:,:) 
  real(8), allocatable        :: mol_tmp(:,:) 
  real(8), allocatable        :: temp(:)
  real(8), allocatable        :: false8(:)
  real, allocatable           :: false(:)
  real(8), parameter          :: zero = 1.0d-40
  integer                     :: x
  integer                     :: y
  integer                     :: l
  integer                     :: ncid  
  integer                     :: k, k1, k2
  integer                     :: s, s1, s2
  integer                     :: index
  character(80)               :: fname
  character(1)                :: replay  
  character(20)               :: name
  character(20), allocatable  :: sname(:)
  logical                     :: disp


  write(*,*) "============================================"
  write(*,*) "======NetCDF ファイルを指定して下さい======="
  write(*,*) "============================================"
  read(*,*) fname  
!  fname = "obradbt-1000.nc"

  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 であることを知っている
!!!
  !ファイルのオープン
  call netcdf_open(trim(fname), ncid)
  
  !変数の取得. 
  name = 'mu'
  call netcdf_get(ncid, name, mu_tmp)
  write(*,*) "mu"

  name = 'pres'
  call netcdf_get(ncid, name, pres)
  write(*,*) "pres"

  name = 'temp'
  call netcdf_get(ncid, name, temp)
  write(*,*) "temp"

  name = 'molfr'
  call netcdf_get(ncid, name, molfr_tmp)
  write(*,*) "molfr"

  name = 'mol'
  call netcdf_get(ncid, name, mol_tmp)
  write(*,*) "mol"

  name = 'sname'
  call netcdf_get(ncid, name, sname)
  write(*,*) "sname"

  !ファイルのクローズ
  call netcdf_close(ncid)


!!!
!!!プロットする配列を作る
!!!
  y = size(mu_tmp, 1)  
  x = size(mu_tmp, 2)
  allocate(mu(x,y), molfr_plot(x,y), &
       & mu_plot(x,y), mu_err_plot(x,y), pres_plot(x), &
       & vap(x), svap(x), molfr(x,y), mol(x,y), false(x), false8(x))
  mu = transpose(mu_tmp)
  molfr = transpose(molfr_tmp)
  mol = transpose(mol_tmp)
  mu_plot = real(mu, 4)  
  pres_plot = real(pres, 4)
  molfr_plot = real(molfr, 4)

!!!
!!!欠損値処理
!!!
  false = 999.0d0
  false8 = 999.0d0
  do s = 1, y
     l = len_trim(sname(s))
     name = sname(s)
     mu_plot(:,s) = merge(mu_plot(:,s), false, mol(:,s) /= zero)
  end do
  

!!!
!!!プロット
!!!  
  !グラフのオープン
  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 )

  
!!!
!!!各化学ポテンシャルをプロット
!!!  
  do s = 1, y
     l = len_trim(sname(s))
     name = sname(s)

     if (name(l:l) /= "g") cycle
!     if (minval(mu_plot(:,s), mu_plot(:,s) > 0.0d0) &
!          & == maxval(mu_plot(:,s), mu_plot(:,s) > 0.0d0)) cycle
     if (trim(sname(s)) == "NH4SH-s") cycle
     
     name = name(1:l-1) // 'l'
     call name2num(name, sname, k1)

     name = name(1:l-1) // 's'
     call name2num(name, sname, k2)

     call grfrm
     call uslset('lyinv', .true.)  
     call gllset( 'LMISS', .TRUE. )

     call usspnt(x, mu_plot(:,s), pres_plot)
     if (k1 /= 0) call usspnt(x, mu_plot(:,k1), pres_plot)
     if (k2 /= 0) call usspnt(x, mu_plot(:,k2), pres_plot)

     call grstrn(2)
     call uspfit
     call grstrf
     call ussttl('equilibrium constant', '1', 'Pressure', 'Pa')
     call usdaxs
     
     call sgspli(25)
     call sgplu(x, mu_plot(:,s), pres_plot)
    
     if (k2 /= 0) call sgspli(45)
     if (k2 /= 0) call sgplu(x, mu_plot(:,k2), pres_plot)

     if (k1 /= 0) call sgspli(35)
     if (k1 /= 0) call sgplu(x, mu_plot(:,k1), pres_plot)

     call uxsttl('T', trim(sname(s)), 0.0)


     !化学ポテンシャルの相対誤差
     if (k1 == 0 .AND. k2 == 0) cycle

     call grfrm
     call uslset('lyinv', .true.)  
     call gllset( 'LMISS', .TRUE. )
     
     if (k1 /= 0) then 
        mu(:,k1) = abs((mu(:,k1) - mu(:,s)) / mu(:,s))
        mu_err_plot(:,k1) = real(merge(mu(:,k1), false8, mol(:,k1) /= zero),4)
        where (mu_err_plot(:,k1) <= 0.0d0)
           mu_err_plot(:,k1) = 1.0d-20
        end where
     end if

     if (k2 /= 0) then 
        mu(:,k2) = abs((mu(:,k2) - mu(:,s)) / mu(:,s))
        mu_err_plot(:,k2) = real(merge(mu(:,k2), false8, mol(:,k2) /= zero),4)
        where (mu_err_plot(:,k2) <= 0.0d0)
           mu_err_plot(:,k2) = 1.0d-20
        end where
     end if

     if (k1 /= 0) call usspnt(x, mu_err_plot(:,k1), pres_plot)
     if (k2 /= 0) call usspnt(x, mu_err_plot(:,k2), pres_plot)
     
     call grstrn(4)
     call uspfit
     call grstrf
     call ussttl('equilibrium constant', '1', 'Pressure', 'Pa')
     call usdaxs
     
     if (k1 /= 0) call sgspli(35)
     if (k1 /= 0) call sgplu(x, mu_err_plot(:,k1), pres_plot)

     if (k2 /= 0) call sgspli(45)
     if (k2 /= 0) call sgplu(x, mu_err_plot(:,k2), pres_plot)
     
     call uxsttl('T', trim(sname(s)), 0.0)

  end do


  GRAPH2: do s = 1, y
     if (trim(sname(s)) /= "NH4SH-s") cycle GRAPH2
     name = 'NH3-g'
     call name2num(name, sname, s1)
     name = 'H2S-g'
     call name2num(name, sname, s2)
     
     vap = log(molfr(:,s1) * molfr(:,s2) * pres **2.0d0) - log(1.0d-2)
     svap = 61.781d0 - (10834.0d0 / temp)
     
     call grfrm
     call uslset('lyinv', .true.)
     call usspnt(x, real(vap, 4), pres_plot)
     call usspnt(x, real(svap, 4), pres_plot)
     
     call grstrn(2)
     call uspfit
     call grstrf
     call ussttl('equilibrium constant', '1', 'Pressure', 'Pa')
     call usdaxs
     
     call sgspli(25)
     call sgplu(x, real(vap, 4), pres_plot)
     call sgspli(35)
     call sgplu(x, real(svap, 4), pres_plot)
     
     call uxsttl('T', trim(sname(s)), 0.0)
  end do GRAPH2


  call DclCloseGraphics()
  
  
end program ObrChkMu

