!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 gt4_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, ErrR, zero, & & 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_name_cnd2gas use interface_display implicit none real(8), allocatable :: pres_map(:) real(8), allocatable :: temp_map(:) real(8), allocatable :: entr_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 :: mu(:) real(8), allocatable :: mu_store(:,:) real(8), parameter :: temp_chg = -0.1d0 real(8) :: temp_store(100) real(8) :: entr_all_store(100) real(8) :: entr_err_store(100) real(8) :: err real(8) :: gibbs real(8) :: gibbs_store(100) real(8) :: temp1, temp2, temp3 real(8) :: entr1, entr2, entr3 real(8) :: pres_region real(8) :: min real, allocatable :: temp_plot(:) real, allocatable :: entr_all_plot(:) real, allocatable :: entr_err_plot(:) real, allocatable :: gibbs_plot(:) real, allocatable :: mol_plot(:,:) real, allocatable :: mu_plot(:,:) real, allocatable :: false(:) 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 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), & & mu_store(100,smax)) temp_store = 0.0d0 mol = 0.0d0; mol_store = 0.0d0 !計算圧力の入力 write(*,*) "■どのくらいの圧力領域を計算しますか (Pa)" ! read(*,*) pres_region !計算圧力の入力 write(*,*) "■どの圧力で計算しますか (l)" ! l = minloc(pres_map, 1, pres_map > pres_region) ! do i = l - 5, l + 5 ! write(*,*) i, pres_map(i) ! end do read(*,*) l !入力された圧力に近い格子点を決める write(*,*) "この圧力に関して計算します: ", pres_map(l) !!! !!!ギブス自由エネルギー最小化 !!! !ステップ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_all = entr3) call gibbs_energy(temp3, pres_map(l), mol, property, mu=mu, gibbs=gibbs) !値の保存 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 gibbs_store(i) = gibbs mol_store(i,:) = mol(:) mu_store(i,:) = mu(:) deallocate(mu) !ループ終了条件 if (entr3 < entr_map(l-1)) then exit STEP1 end if !ループを返すための操作 i = i + 1 entr1 = entr3 temp1 = temp3 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_all=entr2) call gibbs_energy(temp2, pres_map(l), mol, property, mu=mu, gibbs=gibbs) !値の保存 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 gibbs_store(i) = gibbs mol_store(i,:) = mol(:) mu_store(i,:) = mu(:) deallocate(mu) !ループ終了条件. if (abs((entr2 - entr_map(l-1))) <= ErrR * entr_map(l-1) ) then write(*,*) "(", l, ") ", "Entropy is now converged" exit STEP2 end if !数値的に温度が等しい場合にはメッセージを出力. if (temp1 == temp2 .OR. temp2 == temp3) then 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), false(step), & & entr_err_plot(step), gibbs_plot(step), & & entr_all_plot(step), mol_plot(step,smax), mu_plot(step,smax)) temp_plot = 0.0d0 entr_all_plot = 0.0d0; entr_err_plot = 0.0d0 mu_plot = 0.0d0 temp_plot = real(temp_store(1:step), 4) entr_all_plot = real(entr_all_store(1:step), 4) entr_err_plot = real(entr_err_store(1:step), 4) gibbs_plot = real(gibbs_store(1:step), 4) mu_plot(:,:) = real(mu_store(1:step,:), 4) mol_plot(:,:) = real(mol_store(1:step,:), 4) temp_bs(1) = minval(temp_plot) temp_bs(2) = maxval(temp_plot) entr_bs(1) = real(entr_map(l-1), 4) entr_bs(2) = real(entr_map(l-1), 4) !!! !!!結果を標準出力へ表示 !!! write(*,*) "temp" write(*,*) " ", temp_map(l) write(*,*) "pres" write(*,*) " ", pres_map(l) call display(emol_map(l,:), "emol") call display(mol_map(l,:), "mol") ! write(tpinfo, *) "temp= ", real(temp_map(t), 4), " ", & ! & "pres= ", real(pres, 4) !!! !!!値の表示 !!! !グラフのオープン call DclOpenGraphics() ! call DclDivideFrame( 'y', 3, 2 ) call DclSetFrameMargin( 0.05, 0.05, 0.05, 0.05 ) ! call DclSetFrameMargin( 0.1, 0.1, 0.1, 0.1 ) !グラフ作成 ! call DclSetFrameTitle( 'OBORO Debugging Mode', & ! & '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 DclScalingPoint(temp_plot, entr_all_plot) call DclScalingPoint(temp_bs, entr_bs) call DclFitScalingParm call DclSetTransFunction call DclSetTitle('Temperature', 'entropy', 'K', 'J/K') call DclDrawScaledAxis call DclDrawMarker(temp_plot, entr_all_plot, index=5) call DclDrawLine(temp_bs, entr_bs, index=21) call DclDrawTitle( 'T', 'entropy', sw=2 ) !グラフ作成 ! call DclSetFrameTitle( 'OBORO Debugging Mode', & ! & '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 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 DclSetFrameTitle( 'OBORO Debugging Mode', & ! & '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 DclScalingPoint(temp_plot, gibbs_plot) call DclFitScalingParm call DclSetTransFunction call DclSetTitle('Temperature', 'gibbs', 'K', '') call DclDrawScaledAxis call DclDrawMarker(temp_plot, gibbs_plot, index=5) call DclDrawTitle( 'T', 'gibbs', sw=2 ) !モル数 do i = 1, smax ! call DclSetFrameTitle( 'OBORO Debugging Mode', & ! & '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 ) if (maxval(mol_plot(:,i)) == minval(mol_plot(:,i))) 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 DclDrawTitle( 'T', property(i)%name, sw=2 ) end do !気相の化学ポテンシャル. 凝縮相の化学ポテンシャルも同時にプロット do i = 1, count(property%phase == 1) ! call DclSetFrameTitle( 'OBORO Debugging Mode', & ! & '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 DclSetTransNumber(1) call DclSetParm('INTERPRET_MISSING_VALUE', .true. ) call DclScalingPoint(temp_plot, mu_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, mu_plot(:,s)) end do call DclFitScalingParm call DclSetTransFunction call DclSetTitle('Temperature', 'chemical potential', 'K', 'J/ K mol') call DclDrawScaledAxis call DclDrawMarker(temp_plot, mu_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 DclDrawMarker(temp_plot, mu_plot(:,s), type=2, index=l) end do call DclDrawTitle( 'T', property(i)%name, sw=2 ) end do !グラフのクローズ call DCLCloseGraphics contains subroutine obrsort(var1, var2) real(8), intent(inout) :: var1 real(8), intent(inout) :: var2 end subroutine obrsort end program OboroDebug