!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
