!c Description:
!c   デバッグモード. 
!c
!c Current Code Owner: 
!c   sugiyama@gfd-dennou.org
!c
!c Histry: 
!c   Version    Date          Comment
!c   -------    ----------    --------
!c   1.0        2004-04-26    新規作成
!c
!c Copyright (C) SUGIYAMA Ko-ichiro, 2004, All rights reserved

program OboroDebug
  use gtool_history
  use dcl
  use prm_type, only: type_property
  use prm_common, only: obrinit, struct, property, emol_ini, &
       & smax, emax, pres_ini, pres_end, temp_ini, cnv, min, &
       & pseudosw, AdbtFile
  use interface_netcdf_open
  use interface_netcdf_get
  use interface_netcdf_close
  use interface_mingibbs_chk
  use interface_entropy
  use interface_gibbs_energy
  use interface_molfraction
  use interface_name_cnd2gas 
  use interface_svap_antoine
  use interface_svap_amp
  
  implicit none
  real(8), allocatable        :: pres_map(:)      
  real(8), allocatable        :: entr_map(:)
  real(8), allocatable        :: temp_map(:)
  real(8), allocatable        :: mol(:)
  real(8), allocatable        :: mol_map(:,:)
  real(8), allocatable        :: mol_tmp(:,:)
  real(8), allocatable        :: mol_store(:,:)
  real(8), allocatable        :: emol_map(:,:)
  real(8), allocatable        :: emol_tmp(:,:)
  real(8), allocatable        :: molfr(:)
  real(8), allocatable        :: molfr_store(:,:)
  real(8), allocatable        :: entr(:)
  real(8), allocatable        :: entr_store(:,:)
  real(8), parameter          :: temp_chg = -0.5d0
  real(8)                     :: temp_store(100)     
  real(8)                     :: entr_all_store(100) 
  real(8)                     :: entr_err_store(100) 
  real(8)                     :: entr_gas_store(100)
  real(8)                     :: err
  real(8)                     :: gibbs
  real(8)                     :: gibbs_store(100)
  real(8)                     :: svap1
  real(8)                     :: temp1, temp2, temp3  
  real(8)                     :: entr1, entr2, entr3  
  real(8)                     :: entr_gas
  real, allocatable           :: temp_plot(:)      
  real, allocatable           :: entr_plot(:,:)
  real, allocatable           :: entr_all_plot(:)
  real, allocatable           :: entr_err_plot(:)
  real, allocatable           :: entr_gas_plot(:)
  real, allocatable           :: gibbs_plot(:)
  real, allocatable           :: mol_plot(:,:)
  real, allocatable           :: presfr_plot(:,:)
  real, allocatable           :: svap_plot(:,:)
  real                        :: temp_bs(2), entr_bs(2)  
  integer                     :: i, l, s
  integer                     :: step
  integer                     :: ncid
  integer                     :: iws
  character(10)               :: name

  
!!!
!!!データ呼び出し
!!!
  call obrinit


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

  !変数の取得. 取り出す変数は pres, spc, elm, temp, mol, molfr 
  name = 'pres'
  call netcdf_get(ncid, name, pres_map)
  name = 'temp'
  call netcdf_get(ncid, name, temp_map)
  name = 'entr'
  call netcdf_get(ncid, name, entr_map)
  name = 'mol'
  call netcdf_get(ncid, name, mol_tmp)
  name = 'emol'
  call netcdf_get(ncid, name, emol_tmp) 
  
  !配列の割り当て
  step = size(pres_map)
  allocate(mol_map(step, smax), emol_map(step, emax))
  mol_map   = transpose(mol_tmp)
  emol_map  = transpose(emol_tmp)

  !ファイルを閉じる
  call netcdf_close(ncid)
  
!!!
!!!ギブス最小化法を行うための初期値を作成する
!!!
  !変数初期化
  allocate(mol(smax), mol_store(100,smax), &
       & molfr(smax), molfr_store(100, smax), &
       & entr(smax),  entr_store(100, smax))
  temp_store = 0.0d0
  mol = 0.0d0; mol_store = 0.0d0
  molfr = 0.0d0; molfr_store = 0.0d0
  entr = 0.0d0; entr_store = 0.0d0
  
  !計算圧力の入力
  write(*,*) "■圧力を入力して下さい (ステップ数)"
  read(*,*) l

  !入力された圧力に近い格子点を決める
  write(*,*) "この圧力に関して計算します: ", pres_map(l)
  write(*,*) "前のステップでの収束温度: ", temp_map(l-1)
  write(*,*) "前のステップでのエントロピー: ", entr_map(l-1)
  write(*,*) "前のステップでの元素組成: ", emol_map(l-1,:)
  write(*,*) ""  

!!!
!!!ギブス自由エネルギー最小化
!!!
  !ステップ1: エントロピーの値を挟む温度を求める
  i = 0
  STEP1: do 
     !温度定義
     temp3 = temp_map(l - 1) + i * temp_chg
        
     !エントロピーの計算
     call mingibbs_chk(temp3, pres_map(l), &
          & emol_map(l-1,:), struct, property, mol, mol_map(l-1,:))
     call entropy(temp3, pres_map(l), mol, property, &
          & entr=entr, entr_all = entr3, entr_gas = entr_gas)
     call gibbs_energy(temp3, pres_map(l), mol, property, gibbs=gibbs)
     call molfraction(mol, property%phase, molfr)
        
     !ループ終了条件
     if (entr3 < entr_map(l-1)) then 
        exit STEP1
     end if
     
     !ループを返すための操作
     i = i + 1
     entr1 = entr3
     temp1 = temp3

     !値の保存
     err = abs(entr3 - entr_map(l-1))/entr_map(l-1)
     temp_store(i) = temp3
     entr_all_store(i) = entr3
     entr_err_store(i) = err
     entr_gas_store(i) = entr_gas
     gibbs_store(i) = gibbs
     entr_store(i,:) = entr(:)
     mol_store(i,:) = mol(:)
     molfr_store(i,:) = molfr(:)
  end do STEP1
     
  !ステップ 2: 温度の中点でのエントロピーを求め, 1 ステップ前での
  !            エントロピーと同じになるまで収束させる.
  STEP2: do 
     !中間点での値を計算する
     temp2 = (temp1 + temp3) / 2.0d0
     
     call mingibbs_chk(temp2, pres_map(l), &
          & emol_map(l-1,:), struct, property, mol, mol_map(l-1,:))
     call entropy(temp2, pres_map(l), mol, property, &
          & entr=entr, entr_all=entr2, entr_gas = entr_gas)
     call gibbs_energy(temp2, pres_map(l), mol, property, gibbs=gibbs)
     call molfraction(mol, property%phase, molfr)

     !値の保存
     err = abs(entr2 - entr_map(l-1))/entr_map(l-1)
     i = i + 1
     temp_store(i) = temp2
     entr_all_store(i) = entr2
     entr_err_store(i) = err
     entr_gas_store(i) = entr_gas
     entr_store(i,:) = entr(:)
     gibbs_store(i) = gibbs
     mol_store(i,:) = mol(:)
     molfr_store(i,:) = molfr(:)
     
     !ループ終了条件. 
     if (abs((entr2 - entr_map(l-1))/entr_map(l-1)) <= cnv) then
        write(*,*) "(", l, ") ", "Entropy is now converged" 
        write(*,*) "temp = ", temp2
        exit STEP2
     end if
     
     !数値的に温度が等しい場合にはメッセージを出力. 
     if (temp1 == temp2 .OR. temp2 == temp3) then
        temp2 = temp3
        write(*,*) "Conserved Entropy is slightly changed", &
             & entr_map(l-1), " -> ", entr2
        write(*,*) "(", l, ") ", "Entropy is now converged" 
        exit STEP2
     end if
     
     !ステップを進めるための処理
     if (entr3 < entr_map(l-1) .and. entr_map(l-1) < entr2) then
        temp1 = temp2
        entr1 = entr2
     elseif (entr2 <= entr_map(l-1) .and. entr_map(l-1) < entr1) then
        temp3 = temp2
        entr3 = entr2
     end if
  end do STEP2

  
!!!
!!!プロットするための配列を作る
!!!
  step = minloc(temp_store, 1, temp_store == 0.0d0) -1

  allocate(temp_plot(step), &
       & entr_plot(step,smax), entr_err_plot(step), entr_gas_plot(step), &
       & entr_all_plot(step), gibbs_plot(step), mol_plot(step,smax), &
       & presfr_plot(step,smax), svap_plot(step,smax))
  temp_plot = 0.0d0; entr_plot = 0.0d0
  entr_all_plot = 0.0d0; entr_err_plot = 0.0d0
  molfr = 0.0d0
  presfr_plot = 0.0d0
  svap_plot = 0.0d0
  
  temp_plot = real(temp_store(1:step), 4)
  entr_plot(:,:) = real(entr_store(1:step,:), 4)
  entr_all_plot = real(entr_all_store(1:step), 4)
  entr_err_plot = real(entr_err_store(1:step), 4)
  entr_gas_plot = real(entr_gas_store(1:step), 4)
  gibbs_plot = real(gibbs_store, 4)
  mol_plot(:,:) = real(mol_store(1:step,:), 4)
  presfr_plot(:,:) = real(molfr_store(1:step,:) * pres_map(l), 4)
  
  entr_bs = entr_map(l-1)
  temp_bs(1) = minval(temp_plot)
  temp_bs(2) = maxval(temp_plot)


!!!
!!!飽和蒸気圧
!!!
  do i = 1, step
     PLOT: do s = 1, smax
        if (property(s)%method == 3) then 
           call svap_amp(temp_store(i), property(s)%svap_cef, svap1)
        elseif (property(s)%method == 4) then 
           call svap_antoine(temp_store(i), property(s)%svap_cef, svap1)
        else
           cycle plot
        end if
        
        svap_plot(i,s) = real(svap1, 4) 
     end do PLOT
  end do


  
!!!
!!!値の表示
!!!
  !グラフのオープン
  call DclOpenGraphics()
!  call DclDivideFrame( 'y', 3, 2 )
  call DclSetFrameMargin( 0.05, 0.05, 0.05, 0.05 )

  !グラフ作成
  call DclNewFrame
  call DclScalingPoint(temp_plot, entr_all_plot)
!  call DclScalingPoint(temp_bs, entr_bs) 
!  call DclScalingPoint(temp_bs, entr_gas_plot) 
  call DclFitScalingParm
  call DclSetTransFunction

  call DclSetTitle('Temperature', 'entropy', 'K', 'J/K')
  call DclDrawScaledAxis
  
  call DclDrawMarker(temp_plot, entr_all_plot, index=5, type=2)
!  call DclDrawMarker(temp_plot, entr_gas_plot, index=3, type=1)
!  call DclDrawLine(temp_bs, entr_bs, index=21) 

  call DclDrawTitle( 'T', 'entropy', sw=2 )
  

  !グラフ作成
  call DclNewFrame
  call DclSetTransNumber(2)
  call DclScalingPoint(temp_plot, entr_err_plot)
  call DclFitScalingParm
  call DclSetTransFunction

  call DclSetTitle('Temperature', 'entropy error', 'K', '')
  call DclDrawScaledAxis
  call DclDrawMarker(temp_plot, entr_err_plot, index=5)
  call DclDrawTitle( 'T', 'entropy error', sw=2 )

  !化学ポテンシャルのプロット
  call DclNewFrame
  call DclScalingPoint(temp_plot, gibbs_plot)
  call DclFitScalingParm
  call DclSetTransFunction
  
  call DclSetTitle('Temperature', 'K', 'entropy', 'J/K')
  call DclDrawScaledAxis
  call DclDrawMarker(temp_plot, gibbs_plot, index=5)
  call DclDrawTitle( 'T', 'gibbs', sw=2 )

  !モル数
  do i = 1, smax
     if (maxval(mol_plot(:,i)) <= minval(mol_plot(:,i)) * 1.001d0 ) cycle

     call DclNewFrame
     call DclScalingPoint(temp_plot, mol_plot(:,i))
     call DclFitScalingParm
     call DclSetTransFunction
     
     call DclSetTitle('Temperature', 'mol', 'K', 'mol')
     call DclDrawScaledAxis
     call DclDrawMarker(temp_plot, mol_plot(:,i), index=5)
!     call DclDrawLine(temp_plot, mol_plot(:,i), index=22) 
     call DclDrawTitle( 'T', property(i)%name, sw=2 )
  end do
  

  !気相の分圧. 飽和蒸気圧も一緒にプロット
  do i = 1, count(property%phase == 1)
     call DclNewFrame
     call DclSetTransNumber(1)
     call DclSetParm('INTERPRET_MISSING_VALUE', .true. )
     call DclScalingPoint(temp_plot, presfr_plot(:,i))

     do s = count(property%phase == 1) + 1, smax
        call name_cnd2gas(property(s)%name, name)
        if (trim(name) /= trim(property(i)%name)) cycle
        call DclScalingPoint(temp_plot, svap_plot(:,s))
     end do

     call DclFitScalingParm
     call DclSetTransFunction

     call DclSetTitle('partial pressure', 'pressure', 'Pa', 'Pa')
     call DclDrawScaledAxis

     call DclDrawMarker(temp_plot, presfr_plot(:,i), index=5)
     do s = count(property%phase == 1) + 1, smax
        call name_cnd2gas(property(s)%name, name)
        if (trim(name) /= trim(property(i)%name)) cycle
        l = s * 10 + 2
        call DclDrawLine(temp_plot, svap_plot(:,s), index=l) 
     end do
     call DclDrawTitle( 'T', property(i)%name, sw=2 )
  end do


  !グラフのクローズ
  call DCLCloseGraphics


end program OboroDebug
