!c Description:
!c   アンモニア水溶液に特化させた
!c   ギブス自由エネルギーの極小点を探るプログラム
!c
!c Current Code Owner: 
!c   sugiyama@gfd-dennou.org
!c
!c Histry: 
!c   Version    Date          Comment
!c   -------    ----------    --------
!c   1.0        2004-07-03    新規作成
!c
!c Copyright (C) SUGIYAMA Ko-ichiro, 2004, All rights reserved

program OboroPhase
  use dcl
  use prm_common, only: obrinit, property, smax, emol_ini, struct
  use interface_gibbs_energy
  use interface_mingibbs
  
  implicit none
  real(8), parameter          :: mol1 =  4.1322d-1
  real(8), parameter          :: mol2 =  9.75d-2
  real(8), parameter          :: mol3 =  8.51d-2
  real(8), parameter          :: mol4 =  1.12d-3
  real(8), parameter          :: mol_min1 =  1.0d-9
  real(8), parameter          :: mol_max1 =  1.0d-2
  real(8), parameter          :: mol_min2 =  1.0d-11
  real(8), parameter          :: mol_max2 =  1.0d-3
  real(8), parameter          :: temp1 = 396.45d0
  real(8), parameter          :: temp2 = 396.5d0
  integer, parameter          :: step = 100
  integer, parameter          :: step_t = 51
  integer, allocatable        :: spc(:)
  real(8), allocatable        :: mol(:)
  real(8), allocatable        :: mol_store(:,:)
  real(8), allocatable        :: mol_map(:,:,:)
  real(8), allocatable        :: mol_cal(:,:,:)
  real(8), allocatable        :: gibbs_eqv_store(:)
  real(8), allocatable        :: gibbs_gas_store(:)
  real(8), allocatable        :: gibbs_map(:,:)
  real(8), allocatable        :: temp_map(:)
  real(8)                     :: temp
  real(8)                     :: pres
  real(8)                     :: gibbs
  real, allocatable           :: gibbs_plot(:,:,:)
  real, allocatable           :: xaxis(:)
  real, allocatable           :: yaxis(:)
  real                        :: xmin, xmax
  real                        :: ymin, ymax
  integer                     :: i, j
  integer                     :: t
  integer                     :: n1(step_t), n2(step_t)
  character(40)               :: tpinfo
  character(40)               :: molinfo

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

!!!
!!!データの初期化
!!!
  pres = 1551515.125d0

  allocate(gibbs_plot(step, step, step_t), &
       & gibbs_map(step, step), &
       & gibbs_eqv_store(step_t), &       
       & gibbs_gas_store(step_t), &
       & temp_map(step_t), &
       & spc(smax), &
       & mol_map(step, step, smax), &
       & mol_store(step_t, smax), &
       & mol(smax), &
       & mol_cal(1000, smax, step_t), &
       & xaxis(step), yaxis(step))
  gibbs_plot = 0.0d0
  mol_cal = 0.0d0
  mol_map = 0.0d0
  spc = (/(i, i = 1, smax)/)
  
  mol_map(:,:,1) = mol1
  mol_map(:,:,2) = mol2
  do i = 1, step
     mol_map(i,:,5) = mol_min1 &
          & * 10.0d0**((i-1) * dlog10(mol_max1 / mol_min1) / step)
     mol_map(i,:,3) = mol3 - mol_map(i,:,5)
     xaxis(i) = real(mol_map(i,1,5), 4)
  end do
  do j = 1, step
     mol_map(:,j,6) = mol_min2 &
          & * 10.0d0**((j-1) * dlog10(mol_max2 / mol_min2) / step)
     mol_map(:,j,4) = mol4 - mol_map(:,j,6)
     yaxis(j) = real(mol_map(1,j,6), 4)
  end do
  
  
!!!
!!!ギブス自由エネルギーの計算
!!!
  do t = 1, step_t
     temp = temp2 - (t - 1) *  (temp2 - temp1) / (step_t - 1)
     temp_map(t) = temp

     !組成を変化させてギブス自由エネルギーを計算
     do i = 1, step
        write(*,*) temp, i
        do j = 1, step
           call gibbs_energy(temp, pres, mol_map(i,j,:), property, gibbs=gibbs)
           gibbs_map(i,j) = gibbs
        end do
     end do
     
     !平衡組成を計算
     mol_map(1,1,3) = mol3 
     mol_map(1,1,4) = mol4 
     mol_map(1,1,5) = 0.0d0
     mol_map(1,1,6) = 0.0d0
     call mingibbs(temp, pres, emol_ini, struct, property, &
          & mol, mol_ini=mol_map(1,1,:), mol_cal=mol_cal(:,:,t))
     call gibbs_energy(temp, pres, mol, property, gibbs=gibbs)
     gibbs_eqv_store(t) = gibbs
     mol_store(t,:) = mol
     SEARCH1: do i = 1, 100
        if (mol_cal(i,5,t) < mol_max1) then 
           n1(t) = i 
           exit SEARCH1
        end if
     end do SEARCH1
     SEARCH2: do i = 1, 100
        if (mol_cal(i,5,t) < mol_min1) then 
           n2(t) = i - 1
           exit SEARCH2
        end if
     end do SEARCH2
!     n1(t) = minloc(mol_cal(:,5,t), 1, mol_cal(:,5,t) < mol_min1) - 1 
     write(*,*) gibbs

     !平衡組成を計算     
     mol_map(1,1,3) = mol3
     mol_map(1,1,4) = mol4
     mol_map(1,1,5) = 1.0d-40
     mol_map(1,1,6) = 1.0d-40
     call gibbs_energy(temp, pres, mol_map(1,1,:), property, gibbs=gibbs)
     gibbs_gas_store(t) = gibbs
     write(*,*) gibbs

     !気相のみのギブス自由エネルギーの差として保存
     gibbs_map = gibbs_map - gibbs
     gibbs_plot(:,:,t) = real(gibbs_map, 4)
  end do

!!!
!!!プロットするための配列を作る
!!!
  xmin = real(mol_min1, 4)
  xmax = real(mol_max1, 4)
  ymin = real(mol_min2, 4)
  ymax = real(mol_max2, 4)
  write(*,*) n1
  write(*,*) n2


!!!
!!!値の表示
!!!
  !グラフのオープン
  call DclOpenGraphics()
!  call DclGetParm( 'GLOBAL:rmiss', rmiss )
  call DclSetParm( 'GRAPH:LSOFTF', .false. )
!  call DclSetFrameMargin( 0.05, 0.05, 0.05, 0.05 )
  call DclSetFrameMargin( 0.1, 0.1, 0.1, 0.1 )

  do t = 1, step_t
     write(tpinfo, *) "temp= ", real(temp_map(t), 4), "  ", &
          & "pres= ", real(pres, 4)
     
     if (mol_store(t,5) <= 1.0d-40) then
        molinfo = "OBORO: ALL GAS"
     else        
        molinfo = "OBORO: H2O-NH3(liq)"
     end if

     call DclSetFrameTitle( 'checking for gibbs free energy plane', 't',  0.,  0., 0.03, 1 )
     call DclSetFrameTitle( 'H2, He, H2O, NH3', 'b', -1.,  1., 0.02, 2 )
     call DclSetFrameTitle( tpinfo,  'b',  0.,  0., 0.02, 3 )
     call DclSetFrameTitle( molinfo,   'b',  1., -1., 0.02, 4 )     

     call DclNewFrame  
     call DclSetWindow( xmin, xmax, ymin, ymax )
     call DclSetViewPort( 0.2, 0.8, 0.2, 0.8 )
     call DclSetTransNumber( DCL_LOG_LOG )
     call DclSetTransFunction
     
     call DclDrawAxisLog( 'bt', 1, 9 )
     call DclDrawTitle( 'b', 'H2O(l)', 0.0 )
     
     call DclDrawAxisLog( 'lr', 1, 9 )
     call DclDrawTitle( 'l', 'NH3(l)', 0.0 )
     
     call DclSetXGrid( xaxis )
     call DclSetYGrid( yaxis )
     
     if (minval(gibbs_plot(:,:,t)) < 0.0d0) &
          & call DclSetContourLevel( minval(gibbs_plot(:,:,t)), 0., -10.)

     call DclSetShadeLevel( minval(gibbs_plot), maxval(gibbs_plot), -20. )
     call DclShadeContour( gibbs_plot(:,:,t) )
     
     call DclDrawContour( gibbs_plot(:,:,t) )

     call DclDrawMarker(&
          & real(mol_cal(n1(t):n2(t),5,t), 4), &
          & real(mol_cal(n1(t):n2(t),6,t), 4), index=31, type=3)
!     call DclDrawLine(&
!          & real(mol_cal(1:n1(t),5,t), 4), &
!          & real(mol_cal(1:n1(t),6,t), 4), index=21)
     
     call DclDrawTitle( 'T', 'gibbs - gibbs(gas)', sw=2 )    

  end do

  call DclCloseGraphics

end program OboroPhase
