!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-02-07    作成
!c
!c Copyright (C) SUGIYAMA Ko-ichiro, 2003, All rights reserved

module prm_common
  use prm_type, only: type_property
  use prm_file, only: InputFile

  implicit none
  save
  real(8)                           :: temp_ini    !初期温度
  real(8)                           :: pres_ini    !初期圧力
  real(8)                           :: pres_end    !計算終了圧力
  real(8)                           :: grav        !重力
  real(8)                           :: zero        !最小値
  real(8)                           :: ErrR        !許容相対誤差 
  real(8)                           :: ErrElmR     !許容相対誤差(エントロピー)
  real(8)                           :: ErrA        !許容絶対誤差
  real(8)                           :: SVect       !探索ベクトル変化
  real(8)                           :: CSeed       !凝縮物の種
  real(8), allocatable              :: struct(:,:) !構造式の係数
  real(8), allocatable              :: emol_ini(:) !元素のモル数
  integer                           :: step        !計算のステップ数
  integer                           :: smax        !化学種の数
  integer                           :: emax        !化学種の数
  integer, allocatable              :: atmnum(:)   !原子番号
  character(10)                     :: planet      !惑星名
  character(80)                     :: PropertyDir !物性データファイル置場
  character(80)                     :: AdbtFile    !アウトプットファイル
  character(80)                     :: AtmosFile   !アウトプットファイル
  logical                           :: pseudosw    !擬湿潤のスイッチ
  logical                           :: debug
  logical                           :: solarsw
  logical, allocatable              :: esw_store(:)
  type(type_property), allocatable  :: property(:) !物性値

contains

  subroutine ObrSetErrA(mol)

    implicit none
    real(8), intent(in) :: mol(:)

    !本来は ErrA = ErrR * sum(entr_i * mol_i) / entr_i となるのだが
    !配列にするのは面倒なのでスカラーとした. entr_i の値はおよそ
    !ファクター倍くらいしか大きさが変わらないので, 以下で OK だろう. 
    ! ErrA = ErrR * sum(mol_i) * bar(entr) / entr_i 
    !      = ErrR * sum(mol_i) 
    !念のため 0.1 倍してある. 

    ErrA = sum(mol) * ErrR * 1.0d-1   
    
  end subroutine ObrSetErrA


  subroutine ObrSetEsw(mol, entr)    
    use interface_mol2emol
    
    implicit none
    real(8), intent(inout) :: mol(smax)    
    real(8), intent(inout) :: entr(smax)
    real(8), allocatable   :: emol(:)
    integer                :: i

    !気相以外のエントロピーの値をゼロに
    !凝縮相のモル数をゼロにしてからエントロピーを計算すると
    !対数項が小さくなりすぎるので駄目
    do i = 1, smax
       if (mol(i) <= zero) entr(i) = 0.0d0
    end do

    !エントロピーの比が許容相対誤差に満たないもの, 
    !系から除去された化学種のモル数をゼロに設定する.
    do i = 1, smax
       if (mol(i) * entr(i) <= ErrElmR * dot_product(mol, entr)) then
          mol(i) = 0.0d0
       end if
    end do
    
    !元素数を計算
    call mol2emol(mol, struct, emol)
    
    !元素が 0 なもののスイッチを入れる
    esw_store = emol /= 0.0d0
    write(*,*) "prm_common, esw_store", esw_store
    
  end subroutine ObrSetEsw
  

!!!
!!!デバッグモードのスイッチを入れる
!!!
  subroutine ObrDebugMode
    debug = .true.
    
  end subroutine ObrDebugMode
  

  subroutine obrinit
    use interface_def_gravity
    use interface_def_property
    use interface_def_atmnum
    use interface_def_phase
    use interface_def_struct
    use interface_def_emol
    use interface_def_cp_cef
    use interface_def_svap_cef
    
    implicit none
    character(10)            :: sname(100)     !物質名
    real(8)                  :: srate(20)
    real(8)                  :: emol(20)
    integer                  :: l
    
    NAMELIST /ControlData/ &
         & ErrR, SVect, CSeed, zero, debug, &
         & solarsw, PropertyDir
    NAMELIST /AdbtData/ &
         & temp_ini, pres_ini, pres_end, &
         & step, pseudosw, AdbtFile
    NAMELIST /AtmosData/ &
         & planet, AtmosFile
    NAMELIST /Names/ sname
    NAMELIST /SolarRate/ srate
    NAMELIST /ElementMoles/ emol
    

!!!
!!!Namelist ファイルから変数を取り出す
!!!    
    open (10, FILE=InputFile)
    read(10, NML=ControlData)
    close(10)
    
    open (10, FILE=InputFile)
    read(10, NML=AtmosData)
    close(10)
    
    open (10, FILE=InputFile)
    read(10, NML=AdbtData)
    close(10)
    
    sname = "NotKnown"
    open (10, FILE=InputFile)
    read(10, NML=Names)
    close(10)   

    srate = 0.0d0
    open (10, FILE=InputFile)
    read(10, NML=SolarRate)
    close(10)   

    emol = 0.0d0
    open (10, FILE=InputFile)
    read(10, NML=ElementMoles)
    close(10)   
    
!!!
!!!重力の値を得る
!!!
    call def_gravity(planet, grav)

!!!
!!!物性値を求める
!!!
    smax = count(sname /= "NotKnown")
    call def_property(sname(1:smax), propertydir, property)
    
!!!
!!!各物質に通し番号を振る
!!!
    property%num = (/(l, l = 1, size(property, 1))/)
    
!!!
!!!元素番号及び元素数
!!!
    call def_atmnum(property, atmnum)
    emax = size(atmnum, 1)

!!!
!!! esw_store の初期化
!!!
    allocate(esw_store(emax))
    esw_store = .true.

!!!
!!! ErrA の初期化
!!!
    ErrA = ErrR                !とりあえず一緒にしておく
    ErrElmR = ErrR * 1.0d3    !元素のスイッチが先に入るように

!!!
!!!相毎に番号を割り当てる. 
!!!  気体は 1, 溶液は 2, 固溶体は 3, それ以外は 4
!!!
    call def_phase(property)

!!!
!!!構造式の係数
!!!
    call def_struct(property, struct)

!!!
!!!元素のモル数
!!!
    if (solarsw) then 
       call def_emol(atmnum, srate, emol_ini)
    else
       call def_emol(emol, atmnum, emol_ini)
    end if

!!!
!!!比熱のスプライン補間係数
!!!
    call def_cp_cef(property)

!!!
!!!飽和蒸気圧の係数
!!!
    call def_svap_cef(property, propertydir)

  end subroutine obrinit

end module prm_common
