!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004, 2005. All rights reserved.
!---------------------------------------------------------------------

!= Module DisturbSet
!
!   * Developer: SUGIYAMA Ko-ichiro
!   * Version: $Id: disturbset.f90,v 1.2.4.1 2006/01/13 06:55:05 sugiyama Exp $ 
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview 
!
!擾乱場の設定を読み込み, 擾乱場の値を計算する. 
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!== Future Plans
!
!  * 設定可能な擾乱のタイプを増やす
!  * エラー処理に gf4f90io を利用
!


module disturbset
  !== 暗黙の型宣言禁止
  implicit none

  !== save 属性
  save

  !== Public Interface
  real(8), allocatable  :: ss_PotTempDisturb(:, :) !温位
  real(8), allocatable  :: ss_ExnerDisturb(:, :)   !エクスナー関数
  real(8), allocatable  :: fs_VelXDisturb(:, :)    !速度 u
  real(8), allocatable  :: sf_VelZDisturb(:, :)    !速度 w
  real(8), allocatable  :: ss_MixRtVDisturb(:, :) !水蒸気混合比
  real(8), allocatable  :: ss_MixRtCDisturb(:, :) !雲水混合比
  real(8), allocatable  :: ss_MixRtRDisturb(:, :)  !雨水混合比


contains
  !== Procedure Interface 
  !
  !=== Initialize module and acquire NAMELIST
  !
  !NAMELIST から擾乱場の初期値生成プログラム disturb の設定を読み取る. 
  !
  subroutine disturbset_calc_init(cfgfile)
    !=== Dependency
    use dc_trace,   only: BeginSub, EndSub 
    use dc_message, only: MessageNotify
    use gridset,    only: DimXMin,         &! 配列の X 方向の下限
      &                   DimXMax,         &! 配列の X 方向の上限
      &                   DimZMin,         &! 配列の Z 方向の下限
      &                   DimZMax,         &! 配列の Z 方向の上限
!      &                  Xmin, XMax,      &! x 座標の始点・終点
!      &                  Zmin, Zmax,      &! z 座標の始点・終点
      &                   s_X,             &! X 座標軸(スカラー格子点)
      &                   s_Z,             &! Z 座標軸(スカラー格子点)
      &                   f_Z               ! Z 座標軸(フラックス格子点)
    use basicset,   only: ss_ExnerBasicZ,  &! エクスナー関数(平均)
      &                   ss_PotTempBasicZ,&! 温位(平均)
      &                   ss_MixRtVBasicZ   ! 水蒸気混合比(平均)
    use physset,    only: PressSfc,        &! 標準圧力 [Pa]
      &                   Pi,              &! 円周率
      &                   MolWtDry,        &! 乾燥成分の分子量 [kg/kg]
      &                   MolWtWet,        &! 湿潤成分の分子量 [kg/kg]
      &                   GasRDry,         &! 乾燥成分の気体定数 [J/K kg]
      &                   CpDry             ! 乾燥成分の平均定圧比熱 [J/K kg]


    !=== Input
    character(*), intent(in) :: cfgfile

    !=== Work
    integer       :: i, k, l
    real(8)       :: RH
    real(8)       :: ss_Temp(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)       :: ss_SvapPress(DimXMin:DimXMax, DimZMin:DimZMax)

    real(8)       :: beta(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)       :: Xr, Xc, Zr, Zc
    real(8)       :: DeltaPotTemp
    real(8)       :: DeltaMixRtV
    real(8)       :: factor


    !=== NAMELIST
    NAMELIST /disturbset/ RH !相対湿度

    call BeginSub("disturbset_calc_init", &
&                 fmt="%c",       &
&                 c1="Initialize values of disturbance. ")

    open (10, FILE=cfgfile)
    read(10, NML=disturbset)
    close(10)

    !==== 配列の割り当て
    allocate( ss_PotTempDisturb(DimXMin:DimXMax, DimZMin:DimZMax), &
         &    ss_ExnerDisturb(DimXMin:DimXMax, DimZMin:DimZMax),   &
         &    fs_VelXDisturb(DimXMin:DimXMax, DimZMin:DimZMax),    &
         &    sf_VelZDisturb(DimXMin:DimXMax, DimZMin:DimZMax),    &  
         &    ss_MixRtVDisturb(DimXMin:DimXMax, DimZMin:DimZMax),  &
         &    ss_MixRtCDisturb(DimXMin:DimXMax, DimZMin:DimZMax),  &
         &    ss_MixRtRDisturb(DimXMin:DimXMax, DimZMin:DimZMax)  )

    !==== 初期化
    ss_ExnerDisturb = 0.0d0
    ss_PotTempDisturb = 0.0d0
    fs_VelXDisturb = 0.0d0
    sf_VelZDisturb = 0.0d0
    ss_MixRtVDisturb = 0.0d0
    ss_MixRtCDisturb = 0.0d0
    ss_MixRtRDisturb  = 0.0d0
    

    !=== Klemp and Wilhelmson (1978) と同様の温位設定
    Xr = 10.0d3
    Xc = 1.20d4
!    Zr = 1.5d3 
!    Zc = 2.0d3
    Zr = 2.0d3
    Zc = 1.5d3 
    DeltaPotTemp = 1.5d0
    
    do k = DimZMin, DimZMax          
      do i = DimXMin, DimXMax
        beta(i,k) =  &
          & ( ( ( s_X(i) - Xc ) / Xr ) ** 2.0d0  &
          &    + ( ( s_Z(k) - Zc ) / Zr ) ** 2.0d0 ) ** 5.0d-1                
      end do
    end do

    where ( beta < 1.0d0 ) 
      ss_PotTempDisturb = &
        & DeltaPotTemp * ( dcos( Pi * 5.0d-1 * beta ) ** 2.0d0 )
    end where
        


    !=== 蒸気圧の設定. 
    !==== 飽和蒸気圧
    ss_Temp = ss_ExnerBasicZ * (ss_PotTempDisturb + ss_PotTempBasicZ)
    ss_SvapPress = ss_svap_ss(ss_Temp)
    
!    DeltaMixRtV = 1.0d-3
!    where ( beta < 1.0d0 ) 
!      ss_MixRtVDisturb = &
!        & DeltaMixRtV * ( dcos( Pi * 5.0d-1 * beta ) ** 2.0d0 )
!    end where
    
    do k = DimZMin, DimZMax          
      factor = RH
      
      ss_MixRtVDisturb(:,k) = &
        & ((MolWtWet/MolWtDry) * ss_SvapPress(:,k) * factor) &
        & / (PressSfc * (ss_ExnerBasicZ(:,k) ** (CpDry / GasRDry)) &
        &                            - ss_SvapPress(:,k) * factor) &
        & - ss_MixRtVBasicZ(:,k)
    end do
    
    where (ss_MixRtVDisturb < 0.0d0) 
      ss_MixRtVDisturb = 0.0d0
    end where

    


    call EndSub("disturbset_calc_init")
    
  end subroutine disturbset_calc_init
  

  function ss_svap_ss(ss_temp)
    use gridset,    only: DimXMin,  &! 配列の X 方向の下限
      &                   DimXMax,  &! 配列の X 方向の上限
      &                   DimZMin,  &! 配列の Z 方向の下限
      &                   DimZMax    ! 配列の Z 方向の上限

    ! 暗黙の型宣言禁止
    implicit none

    ! 入出力変数
    real(8)             :: ss_svap_ss(DimXMin:DimXMax, DimZMin:DimZMax)
                                                          !飽和蒸気圧 [Pa]
    real(8), intent(in) :: ss_temp(DimXMin:DimXMax, DimZMin:DimZMax)
                                                          !温度 [K]
    
    ! Antoine の式を利用. H2O の係数は化学工学便覧参照. 
    real(8), parameter  :: A = 23.1964d0
    real(8), parameter  :: B = 3816.44d0
    real(8), parameter  :: C = -46.13d0
    
    ss_svap_ss = dexp( A - B / (C + ss_temp) ) !Antoine の式
    
  end function ss_svap_ss

end module disturbset
