!= Module Constants
!
! Authors::   SUGIYAMA Ko-ichiro, ODAKA Masatsugu
! Version::   $Id: constants.f90,v 1.1 2011-04-26 04:08:10 sugiyama Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2007. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview 
!
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!== Future Plans
!

module constants
  !
  !

  !モジュール読み込み
  use dc_types,   only: DP
  use dc_iounit,  only: FileOpen
  use dc_message, only: MessageNotify   !メッセージ出力
  use timeset,    only: DelTimeShort, DelTimeLong
  use ChemData,   only: GasRUniv,             &!普遍気体定数
    &                   ChemData_OneSpcID,    &!化学種の ID
    &                   ChemData_CpPerMolRef, &!標準状態での比熱
    &                   ChemData_MolWt         !分子量
  use axesset, only:    x_dx, y_dy, z_dz

  !暗黙の型宣言禁止
  implicit none

  private

  !Public Interface
  real(DP), save, public :: Nu  = 0.0d0           !音波減衰項の減衰係数
  real(DP), save, public :: NuH   = 0.0d0         !数値粘性の係数 (水平方向)
  real(DP), save, public :: NuV   = 0.0d0         !数値粘性の係数 (鉛直方向)
  real(DP), save, public :: NuH2  = 0.0d0         !数値粘性の係数 (水平方向)
  real(DP), save, public :: NuV2  = 0.0d0         !数値粘性の係数 (鉛直方向)
  real(DP), save, public :: Grav = 9.8d0          !重力 [m/s^2]
  real(DP), save, public :: PressBasis = 965.0d0  !温位の基準圧力 [Pa]
  real(DP), save, public :: TempSfc = 300.0d0     !地表面温度 [K]
  real(DP), save, public :: PressSfc = 965.0d0    !地表面圧力 [Pa]
  real(DP), save, public :: CpDry  = 0.0d0        !乾燥成分の定圧比熱 [J/K kg]
  real(DP), save, public :: CpDryMol = 0.0d0      !乾燥成分の定圧比熱 [J/K kg]
  real(DP), save, public :: CvDry = 0.0d0         !乾燥成分の定積比熱 [J/K kg]
  real(DP), save, public :: MolWtDry = 0.0d0      !乾燥成分の分子量   [kg/mol]
  real(DP), save, public :: GasRDry  = 0.0d0      !乾燥成分の気体定数 [J/K kg]
  real(DP), save, public :: DayTime = 86400.0d0   ! 1 日の長さ [s]

  public constants_init

contains

!!!-----------------------------------------------------------------!!!
  subroutine constants_Init(cfgfile)

    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    character(*), intent(in) :: cfgfile         !NAMELIST ファイル
    integer                  :: SpcDryNum       !乾燥成分の化学種の数
    character(20)            :: SpcDrySymbol(5) !乾燥成分の化学種名
    real(DP)                 :: SpcDryMolFr(5)  !乾燥成分の化学種の存在度
    integer, allocatable     :: SpcDryID(:)     !乾燥成分の化学種のID
    real(DP), allocatable    :: PropertyDry(:)  !作業配列
    integer                  :: s, n1, n2, n3   !作業変数
    integer                  :: unit            !装置番号
    logical                  :: flag
    real(DP)                 :: AlphaSound = 1.0e-7  !音波減衰項の係数
    real(DP)                 :: DelXMin, DelYMin, DelZMin
    real(DP)                 :: AlphaNumDiff = 1.0d-4
    real(DP)                 :: AlphaNumDiffFact = 1.0d0  !速度に対する粘性を上げる場合は数字を 1 以上にする. 
    
    !NAMELIST の定義
    NAMELIST /constants_nml/ &
      & Grav, PressBasis, TempSfc, PressSfc, & 
      & SpcDrySymbol, SpcDryMolFr, DayTime, CpDry, MolWtDry, AlphaSound, AlphaNumDiff, AlphaNumDiffFact

    SpcDrySymbol = '' 
    SpcDryMolFr  = 0.0d0
    
    !ファイルオープン. 情報取得. 
    call FileOpen(unit, file=cfgfile, mode='r')
    read(unit, NML=constants_nml)
    close(unit)

    !-------------------------------------------------------------------
    ! 音波減衰項の減衰率   Min(DelX, DelZ) ** 2.0 に比例
    !
    DelXMin = minval(x_dx)
    DelYMin = minval(y_dy)
    DelZMin = minval(z_dz)
    Nu = AlphaSound * ( Min(DelXMin, DelZMin) ** 2.0d0 ) / DelTimeShort
    
    ! CReSS マニュアルでは, 2 次精度中心差分の場合, Alpha < 1/8 くらいが適当と述べている. 
    ! deepconv では NuH として 500 以上の値が入っていたそうだ,
    NuH = AlphaNumDiff * ( SQRT(DelXMin*DelYMin) ** 2.0d0 ) / DelTimeLong
    NuV = AlphaNumDiff * ( DelZMin ** 2.0d0 ) / DelTimeLong
    NuH2 = NuH * AlphaNumDiffFact 
    NuV2 = NuV * AlphaNumDiffFact 

    if (CpDry /= 0.0d0 .AND. MolWtDry /= 0.0d0) then 
    ! namelist から CpDry の値が入力された場合

      !定圧モル比熱
      CpDryMol = CpDry *  MolWtDry 
      
      !気体定数
      GasRDry = GasRUniv / MolWtDry
      
      !定積比熱
      CvDry    = CpDry - GasRDry

    elseif (SpcDryMolFr(1) /= 0.0d0) then 
    ! namelist から モル比が入力された場合

      flag = .true.

      !----------------------------------------------------------
      ! 乾燥成分の物性値の初期化
      !
      !乾燥成分の個数を数える
      SpcDryNum = count(SpcDrySymbol /= "")
      
      !化学種の ID を取得    
      allocate(SpcDryID(SpcDryNum))    
      do s = 1, SpcDryNum
        SpcDryID(s) = ChemData_OneSpcID( SpcDrySymbol(s) )
      end do
      
      !作業配列の準備
      allocate(PropertyDry(SpcDryNum))
      
      !分子量
      do s = 1, SpcDryNum
        PropertyDry(s) = ChemData_MolWt(SpcDryID(s))
      end do
      MolWtDry = dot_product(PropertyDry, SpcDryMolFr(1:SpcDryNum)) 
      
      !定圧比熱(モル当量)
      do s = 1, SpcDryNum    
        PropertyDry(s) = ChemData_CpPerMolRef(SpcDryID(s))
      end do
      CpDryMol = dot_product(PropertyDry, SpcDryMolFr(1:SpcDryNum)) 
      
      !定圧比熱
      CpDry    = CpDryMol / MolWtDry
      
      !気体定数
      GasRDry = GasRUniv / MolWtDry
      
      !定積比熱
      CvDry    = CpDry - GasRDry
    end if
    
    !----------------------------------------------------------
    ! 確認
    !----------------------------------------------------------
    call MessageNotify( "M", "constants_init", "Nu = %f", d=(/Nu/) )
    call MessageNotify( "M", "constants_init", "NuH = %f", d=(/NuH/) )
    call MessageNotify( "M", "constants_init", "NuV = %f", d=(/NuV/) )
    call MessageNotify( "M", "constants_init", "NuH2 = %f", d=(/NuH2/) )
    call MessageNotify( "M", "constants_init", "NuV2 = %f", d=(/NuV2/) )
    call MessageNotify( "M", &
      & "constants_init", "Grav = %f",     d=(/Grav/) )
    call MessageNotify( "M", &
      & "constants_init", "PressBasis = %f", d=(/PressBasis/))
    call MessageNotify( "M", &
      & "constants_init", "TempSfc = %f",  d=(/TempSfc/) )
    call MessageNotify( "M", &
      & "constants_init", "PressSfc = %f", d=(/PressSfc/) )
    if (flag) then 
      do s = 1, SpcDryNum
        call MessageNotify( "M", &
          &  "constants_init", "SpcDryID = %d",      i=(/SpcDryID(s)/))
        call MessageNotify( "M", &
          &  "constants_init", "SpcDrySymbol = %c", c1=trim(SpcDrySymbol(s)))
        call MessageNotify( "M", &
          &  "constants_init", "SpcDryMolFr = %f",   d=(/SpcDryMolFr(s)/))
      end do
    end if

    call MessageNotify( "M", "constants_init", "CpDry    = %f",    d=(/CpDry/) )
    call MessageNotify( "M", "constants_init", "CpDryMol = %f", d=(/CpDryMol/) )
    call MessageNotify( "M", "constants_init", "CvDry    = %f",    d=(/CvDry/) )
    call MessageNotify( "M", "constants_init", "GasRDry  = %f",  d=(/GasRDry/) )
    call MessageNotify( "M", "constants_init", "MolWtDry = %f", d=(/MolWtDry/) )
    call MessageNotify( "M", "constants_init", "DayTime  = %f",  d=(/DayTime/)  )

  end subroutine Constants_Init

end module Constants
