!c Description: !c 断熱曲線の計算プログラム !c !c Current Code Owner: !c sugiyama@gfd-dennou.org !c !c Histry: !c Version Date Comment !c ------- ---------- -------- !c 1.0 2003-03-31 新規作成 !c !c Copyright (C) SUGIYAMA Ko-ichiro, 2003, All rights reserved program obradbt use gtool_history use prm_common, only: obrinit, struct, property, emol_ini, & & smax, emax, step, pres_ini, pres_end, temp_ini, cnv, pseudosw use interface_mingibbs use interface_mingibbs_chk use interface_entropy use interface_mol2emol use interface_rmcnd use interface_molfraction use interface_exec_time implicit none real(8), parameter :: temp_chg = -0.5d0 real(8), allocatable :: pres_map(:) real(8), allocatable :: entr_map(:) real(8), allocatable :: temp_map(:) real(8), allocatable :: emol_map(:,:) real(8), allocatable :: mol_map(:,:) real(8), allocatable :: molfr_map(:,:) real(8), allocatable :: mol(:) real(8), allocatable :: molfr(:) real(8), allocatable :: emol(:) real(8) :: entr1 real(8) :: entr2 real(8) :: entr3 real(8) :: entr_gas real(8) :: temp1 real(8) :: temp2 real(8) :: temp3 integer :: i, l integer :: time1(8) integer :: time2(8) integer :: time(2) logical :: esw !!! !!!実行開始時間を計測 !!! call date_and_time(values = time1) !!! !!!データ呼び出し !!! call obrinit !!! !!!配列の割り当て !!! allocate(pres_map(step), & & entr_map(step), & & temp_map(step), & & mol_map(step, smax), & & molfr_map(step, smax), & & emol_map(step, emax), & & mol(smax)) !!! !!!NetCDF ファイルの用意 !!! call output_gtool4_init !!! !!!計算に用いる圧力の設定 !!! esw = .false. pres_map = pres_ini do i = 2, step ! pres_map(i) = pres_ini & ! & * 10.0d0**((i - 1) * dlog10(pres_end / pres_ini) / (step - 1)) pres_map(i) = pres_ini + (i - 1) * (pres_end - pres_ini) / (step - 1) end do call HistoryPut('p', real(pres_map, 4)) !!! !!!初期温度, 初期圧力を別途計算 !!! l = 1 temp2 = temp_ini call mingibbs_chk(temp2, pres_map(l), emol_ini, struct, property, mol) call entropy(temp2, pres_map(l), mol, & & property, entr_all=entr2, entr_gas=entr_gas) call store_data(l, esw) !!! !!!温度を変化させながら計算 !!! CHANGE_PRES: do l = 2, step ! write(*,*) "START CHANGE_PRES" !前の圧力での収束温度でのエントロピーだけ別途計算 temp1 = temp_map(l - 1) + 2.0d0 call mingibbs_chk(temp1, pres_map(l), & & emol_map(l-1,:), struct, property, mol, mol_map(l-1,:)) call entropy(temp1, pres_map(l), mol, property, entr_all = entr1) !念のためのチェック if (entr1 < entr_map(l-1)) then write(*,*) "the entropy is lower than the pervious entropy" write(*,*) "temp ", temp1 write(*,*) "entr ", entr1, entr_map(l-1) stop end if !ステップ1: エントロピーの値を挟む温度を求める i = 1 STEP1: do !温度定義 temp3 = temp1 + 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) !ループ終了条件 if (entr3 < entr_map(l-1)) exit STEP1 !ループを返すための操作 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, entr_gas=entr_gas) !ループ終了条件. if (abs((entr2 - entr_map(l-1))/entr_map(l-1)) <= cnv) then temp_map(l) = temp2 call store_data(l, esw) 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 ! ! temp3 で計算しなおし ! 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=entr2, entr_gas=entr_gas) esw = .true. call store_data(l, esw) esw = .false. 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 end do CHANGE_PRES !!! !!!NetCDF ファイルを閉じる !!! call HistoryClose !!! !!!実行開始時間を計測 !!! call date_and_time(values = time2) call exec_time(time1, time2, time) write(*,*) "計算実行時間: ", time(1), "分", time(2), "秒" contains subroutine output_gtool4_init use prm_common, only: smax, emax, planet, AdbtFile use gtool_history implicit none real(8) :: spc(smax) real(8) :: elm(emax) integer :: num integer :: i character(1000) :: dataset !ファイルに書き出すべき初期値を変数に入れる write(dataset, *) "pseudosw= ", pseudosw, " ", & & planet, property%name num = len_trim(dataset) !大域属性, 座標軸の設定 call HistoryCreate(& & file=AdbtFile, & & title="Equilibrium Calculation", & & source=dataset(:num+1), & & institution="GFD_Dennou_Club OBORO Project", & & dims=(/'spc ', 'elm ', 'pres'/), & & dimsizes=(/int(smax, 4), int(emax, 4), 0/), & & longnames=(/"chemical species", "element ", "Pressure "/), & & units=(/"1 ", "1 ", "Pa"/), & & origin=0.0, & & interval=0.0) !次元の定義 spc = (/(i, i = 1, smax)/) call HistoryPut('spc', real(spc, 4)) elm = (/(i, i = 1, emax)/) call HistoryPut('elm', real(elm, 4)) !変数定義, エントロピー call HistoryAddVariable( & & varname="entr", & & dims=(/'pres'/), & & longname="Entropy", & & units='J', & & xtype="double" ) !変数定義, 温度 call HistoryAddVariable( & & varname="temp", & & dims=(/'pres'/), & & longname="Temperature", & & units='K', & & xtype="double" ) !変数定義, モル call HistoryAddVariable( & & varname="mol", & & dims=(/'spc ', 'pres'/), & & longname="Chemical Species Abundance", & & units='mol', & & xtype="double" ) !変数定義, モル比 call HistoryAddVariable( & & varname="molfr", & & dims=(/'spc ', 'pres'/), & & longname="Mol Fraction", & & units='1', & & xtype="double" ) !変数定義, 元素量 call HistoryAddVariable( & & varname="emol", & & dims=(/'elm ', 'pres'/), & & longname="Element Abundance", & & units='mol', & & xtype="double" ) end subroutine output_gtool4_init subroutine store_data(l, sw) use prm_common, only: ErrCode implicit none integer :: l logical :: sw ! エントロピーが収束しない場合の処理. 欠損値として扱う if (sw) then temp_map(l) = temp2 mol_map(l, :) = ErrCode molfr_map(l,:) = ErrCode entr_map(l) = entr_map(l-1) emol_map(l,:) = emol_map(l-1,:) else temp_map(l) = temp2 mol_map(l, :) = mol call molfraction(mol, property%phase, molfr) molfr_map(l,:) = molfr(:) if (pseudosw) then entr_map(l) = entr_gas write(*,*) entr_gas call rmcnd(mol, property) call mol2emol(mol, struct, emol) emol_map(l, :) = emol(:) else entr_map(l) = entr2 emol_map(l, :) = emol_ini !元素量は保存 end if end if call HistoryPut('pres', pres_map(l)) call HistoryPut('temp', temp_map(l)) call HistoryPut('mol', mol_map(l,:)) call HistoryPut('molfr', molfr_map(l,:)) call HistoryPut('entr', entr_map(l)) call HistoryPut('emol', emol_map(l,:)) end subroutine store_data end program obradbt