!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