!c Description: !c 各化学種のデータベースファイルを作成するためのプログラム !c 物性データは netCDF ファイルに入力する. !c !c Current Code Owner: !c sugiyama@gfd-dennou.org !c !c Histry: !c Version Date Comment !c ------- ---------- -------- !c 1.0 2003-02-14 作成 !c !c Copyright (C) SUGIYAMA Ko-ichiro, 2003, All rights reserved program obrdata use netcdf use interface_netcdf_def use interface_netcdf_put use interface_netcdf_open use interface_netcdf_close use interface_mol2atm use interface_atm2wt use prm_chem, only: ChemInfo implicit none character(40) :: name !変数名 character(80) :: long_name !long_name 属性値 character(40) :: units !units 属性値 character(40) :: file real(8) :: molwt real(8) :: rvar1 real(8) :: rvar2 real(8) :: rvar3(100, 2) real(8) :: rvar4(5) real(8), allocatable :: rvar5(:,:) real(8), aLlocatable :: enum(:) integer :: ivar1 integer :: ivar2 integer, allocatable :: atmnum(:) character(5) :: temp_unit character(5) :: svap_unit character(40) :: cvar integer :: status integer :: num integer :: len integer :: ncid !ファイルの ID integer :: spcid !次元配列の ID integer :: elmid !次元配列の ID integer :: tempid !次元配列の ID integer :: cef5id !次元配列の ID !!! !!!データ呼び出し(定数参照型モジュールのサブルーチンを呼び出す) !!! call ChemInfo !!! !!!ファイル準備 !!! write(*,*) "化学種の名前を入力して下さい. (例: H2O-g)" read(*,*) file file = trim(file) // '.nc' write(*,*) "creating data file ", trim(file) !ファイル作成 call netcdf_create(file, ncid) !!! !!!化学種に含まれる元素の原子番号, 個数, 分子量を計算 !!! name = file(1:len_trim(file)-5) call mol2atm(name, atmnum, enum) call atm2wt(atmnum, enum, molwt) !!! !!!次元の定義 !!! !化学種に相当する次元 spc を作成 name = "spc" len = 1 call netcdf_def(ncid, name, len, spcid) !元素に相当する次元を作成 name = "atmnum" len = size(atmnum, 1) call netcdf_def(ncid, name, len, elmid) !温度に相当する次元を作成 name = "temp" len = nf90_unlimited call netcdf_def(ncid, name, len, tempid) !飽和蒸気圧の係数を保管するための次元を作成 name = "cef5" len = 5 call netcdf_def(ncid, name, len, cef5id) !!! !!!変数の定義 !!! !モル分率の定義 name = "molwt" long_name = "molecular weight" units = "g/mol" call netcdf_def(ncid, name, long_name, units, nf90_double, spcid) ! !標準温度の定義 ! name = "temp_ref" ! long_name = "reference temperature" ! units = "K" ! call netcdf_def(ncid, name, long_name, units, nf90_double, spcid) ! !標準圧力の定義 ! name = "pres_ref" ! long_name = "reference pressure" ! units = "Pa" ! call netcdf_def(ncid, name, long_name, units, nf90_double, spcid) !標準エントロピーの定義 name = "entr_ref" long_name = "reference entropy" units = "J/ K kg" call netcdf_def(ncid, name, long_name, units, nf90_double, spcid) !標準エンタルピーの定義 name = "enth_ref" long_name = "reference enthalpy" units = "J/ K kg" call netcdf_def(ncid, name, long_name, units, nf90_double, spcid) !水溶液のスイッチ name = "solsw" long_name = "solution switch" units = "1" call netcdf_def(ncid, name, long_name, units, nf90_int, spcid) !計算方法を示すスイッチ name = "method" long_name = "calculation method for chemical potential" units = "1" call netcdf_def(ncid, name, long_name, units, nf90_int, spcid) !原子番号 name = "atmnum" long_name = "atomic number" units = "1" call netcdf_def(ncid, name, long_name, units, nf90_int, elmid) !元素数 name = "enum" long_name = "number of element in the chemical spece" units = "1" call netcdf_def(ncid, name, long_name, units, nf90_double, elmid) !温度 name = "temp" long_name = "temperature" units = "K" call netcdf_def(ncid, name, long_name, units, nf90_double, tempid) !比熱 name = "cp" long_name = "specific heat" units = "J/ K mol" call netcdf_def(ncid, name, long_name, units, nf90_double, tempid) !飽和蒸気圧 name = "svap" long_name = "saturated vaper pressure" units = "Pa" call netcdf_def(ncid, name, long_name, units, nf90_double, tempid) !飽和蒸気圧の係数 name = "amp" long_name = "coefficients for the equation of saturated vaper pressure" units = "1" call netcdf_def(ncid, name, long_name, units, nf90_double, cef5id) !飽和蒸気圧の係数 name = "antoine" long_name = "coefficients for the equation of saturated vaper pressure" units = "1" call netcdf_def(ncid, name, long_name, units, nf90_double, cef5id) !!! !!!グローバル属性 !!! status = nf90_put_att(ncid, nf90_global, "title", "thermodynamical properties for chemical spece") status = nf90_put_att(ncid, nf90_global, "institution", "GFD_Dennou_Club oboro Project") call fdate(name) status = nf90_put_att(ncid, nf90_global, "history", name) !!! !!!定義モード終了 !!! call netcdf_enddef(ncid) !!! !!!情報の追加(その 1). スカラー量 !!! !分子量の入力. 値は化学種名から判定. name = "molwt" call netcdf_put(ncid, name, molwt) !原子番号の入力. name = "atmnum" call netcdf_put(ncid, name, atmnum) !元素数の入力. name = "enum" call netcdf_put(ncid, name, enum) ! !標準温度の入力. ! write(*,*) "標準温度を選択して下さい" ! write(*,*) " 1. 298.15 [K]" ! write(*,*) " 2. 298.15 [K] 以外" ! read(*,*) ivar1 ! if (ivar1 /= 1 .AND. ivar1 /= 2) then ! write(*,*) "選択された番号が間違っています" ! stop ! elseif (ivar1 == 1) then ! rvar1 = 298.15d0 ! elseif (ivar1 == 2) then ! write(*,*) "標準温度を入力して下さい. 単位は K (ケルビン)." ! read(*,*) rvar1 ! end if ! name = "temp_ref" ! call netcdf_put(ncid, name, rvar1) ! !標準圧力の入力. ! write(*,*) "標準圧力を選択して下さい" ! write(*,*) " 1. 100,000 [Pa]" ! write(*,*) " 2. 100.000 [Pa] 以外" ! read(*,*) ivar1 ! if (ivar1 /= 1 .AND. ivar1 /= 2) then ! write(*,*) "選択された番号が間違っています" ! stop ! elseif (ivar1 == 1) then ! rvar1 = 1.0d5 ! elseif (ivar1 == 2) then ! write(*,*) "標準圧力を入力して下さい. 単位は Pa (パスカル)." ! read(*,*) rvar1 ! end if ! name = "pres_ref" ! call netcdf_put(ncid, name, rvar1) !標準状態でのエントロピーの入力. write(*,*) "標準状態でのエントロピーを入力して下さい. 単位は J/ K mol." write(*,*) "化学ポテンシャルを蒸気圧から計算する場合は値が 0 でも OK." read(*,*) rvar1 name = "entr_ref" call netcdf_put(ncid, name, rvar1) !標準状態でのエンタルピーの入力. write(*,*) "標準状態でのエンタルピーを入力して下さい. 単位は J/ K mol." read(*,*) rvar1 name = "enth_ref" call netcdf_put(ncid, name, rvar1) !水溶液のスイッチ write(*,*) "この物質は水溶液を作りますか? [y/n]" read(*,*) cvar if (trim(cvar) /= "y" .AND. trim(cvar) /= "n") then write(*,*) "正しく選択されていません" stop elseif (trim(cvar) == "y") then ivar1 = 1 elseif (trim(cvar) == "n") then ivar1 = 0 end if name = "solsw" call netcdf_put(ncid, name, ivar1) !計算方法を示すスイッチ write(*,*) "化学ポテンシャル計算方法を選択して下さい." write(*,*) " 1. 希ガス" write(*,*) " 2. 比熱をスプライン補間して計算" write(*,*) " 3. 飽和蒸気圧から計算. 補間式は AMP 法" write(*,*) " 4. 飽和蒸気圧から計算. 補間式は Antoine 法(SI 単位系)" write(*,*) " 5. NH4SH(s) の平衡定数から計算" read(*,*) ivar2 if (ivar2 /= 1 .AND. ivar2 /= 2 .AND. ivar2 /= 3 & & .AND. ivar2 /= 4 .AND. ivar2 /= 5 .AND. ivar2 /= 6) then write(*,*) "正しく選択されていません" stop end if name = "method" call netcdf_put(ncid, name, ivar2) !温度と比熱の関係を入力 if (ivar2 == 2) then write(*,*) "温度, 比熱の順に比熱の温度依存性を入力して下さい. " write(*,*) "100K で 200 J/K mol ならば 100 200 と入力して下さい. " write(*,*) "入力が終了したら Ctl-d でストップして下さい. " num = 1 input_cp: do read(*,*, iostat = status) rvar1, rvar2 if (status /= 0) exit input_cp rvar3(num, :) = (/rvar1, rvar2/) num = num + 1 end do input_cp write(*,*) "比熱の入力終了" num = minloc(rvar3(:,1), 1) allocate(rvar5(num-1, 2)) rvar5(:,1) = pack(rvar3(:,1), rvar3(:,1) /= 0.0d0) rvar5(:,2) = pack(rvar3(:,2), rvar3(:,2) /= 0.0d0) name = "temp" call netcdf_put(ncid, name, rvar5(:,1)) name = "cp" call netcdf_put(ncid, name, rvar5(:,2)) elseif (ivar2 == 3 .OR. ivar2 == 4 .OR. ivar2 == 5) then !蒸気圧と比熱の関係を入力 write(*,*) "蒸気圧の入力" write(*,*) " 温度の単位を選択して下さい" write(*,*) " 1. K (ケルビン)" write(*,*) " 2. C (摂氏)" read(*,*) ivar1 if (ivar1 == 1) then temp_unit = "K" elseif(ivar1 == 2) then temp_unit = "C" else write(*,*) "選択された単位は存在しません" end if write(*,*) " 蒸気圧の単位は?" write(*,*) " 1. Pa (パスカル)" write(*,*) " 2. mmHg (水銀柱高度)" read(*,*) ivar1 if (ivar1 == 1) then temp_unit = "Pa" elseif(ivar1 == 2) then svap_unit = "mmHg" else write(*,*) "選択された単位は存在しません" end if write(*,*) "温度, 蒸気圧の順に蒸気圧の温度依存性を入力して下さい. " write(*,*) "100K で 200 J/K mol ならば 100 200 と入力して下さい. " write(*,*) "入力が終了したら Ctl-d でストップして下さい. " num = 1 input_svap: do read(*,*, iostat = status) rvar1, rvar2 if (status /= 0) exit input_svap call convertunits(temp_unit, rvar1) call convertunits(svap_unit, rvar2) rvar3(num, :) = (/rvar1, rvar2/) num = num + 1 end do input_svap write(*,*) "蒸気圧の入力終了" num = minloc(rvar3(:,1), 1) allocate(rvar5(num-1, 2)) rvar5(:,1) = pack(rvar3(:,1), rvar3(:,1) /= 0.0d0) rvar5(:,2) = pack(rvar3(:,2), rvar3(:,2) /= 0.0d0) name = "temp" call netcdf_put(ncid, name, rvar5(:,1)) name = "svap" call netcdf_put(ncid, name, rvar5(:,2)) end if !飽和蒸気圧の係数を入力 if (ivar2 == 3 .OR. ivar2 == 4 .OR. ivar2 == 5) then write(*,*) "AMP 法の係数を入力しますか? [y/n]" read(*,*) cvar if (trim(cvar) == "y") then rvar4 = 0.0d0 write(*,*) "係数を入力して下さい" read(*,*) rvar4(1), rvar4(2), rvar4(3), rvar4(4), rvar4(5) name = "amp" call netcdf_put(ncid, name, rvar4) end if write(*,*) "antoine 法の係数を入力しますか? [y/n]" read(*,*) cvar if (trim(cvar) == "y") then rvar4 = 0.0d0 write(*,*) "係数を入力して下さい" read(*,*) rvar4(1), rvar4(2), rvar4(3) name = "antoine" call netcdf_put(ncid, name, rvar4) end if end if !!! !!!netCDF ファイルのクローズ !!! call netcdf_close(ncid) end program obrdata