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

!= Module DisturbSet
!
!   * Developer: SUGIYAMA Ko-ichiro
!   * Version: $Id: disturbset.f90,v 1.2 2005/04/22 15:01:52 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 方向の上限
      &                   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]
    use ChemEquiv, only:  ss_SvapPress_ss  ! 飽和蒸気圧 [Pa]
    use ConvertUnit, only:ss_Vap2MixRt_ss  ! 蒸気圧から混合比への変換

    !=== 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
    real(8)       :: Xr, Xc, Zr, Zc
    real(8)       :: DeltaPotTemp
    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 = 2.0d3
    Zc = 1.5d3 
    DeltaPotTemp = 1.5d0
    
    do i = DimXMin, DimXMax
      do k = DimZMin, DimZMax          
        beta =  &
          & ( ( ( s_X(i) - Xc ) / Xr ) ** 2.0d0  &
          &    + ( ( s_Z(k) - Zc ) / Zr ) ** 2.0d0 ) ** 5.0d-1 
        
        if ( beta < 1.0d0 ) then 
          ss_PotTempDisturb(i,k) = &
            & DeltaPotTemp * ( dcos( Pi * 5.0d-1 * beta ) ** 2.0d0 )
        end if
        
      end do
    end do

    !=== 風速場 (Ogura and Takahashi, 1971)
!    ss_PotTempDisturb = 0.0d0 
!    do k = DimZMin, DimZMax          
!      if ((2.0d3 - f_Z(k)) < 0.0d0) cycle
!      if (f_Z(k) < 0.0d0)           cycle
!      ss_PotTempDisturb(:,k) = 1.0d1 * f_Z(k)/1.0d3 * (2.0d0 - f_Z(k)/1.0d3)
!    end do

    !=== 蒸気圧の設定. 
    !==== 飽和蒸気圧
    ss_Temp = ss_ExnerBasicZ * (ss_PotTempDisturb + ss_PotTempBasicZ)
    ss_SvapPress = ss_SvapPress_ss(ss_Temp) * RH

    !==== 混合比 (Ogura and Takahashi, 1971)
    ss_Vap2MixRt_ss_ss(ss_SVapPress, ss_ExnerBasicZ)
    
    where (ss_MixRtVDisturb < 0.0d0) 
      ss_MixRtVDisturb = 0.0d0
    end where

    call EndSub("disturbset_calc_init")
    
  end subroutine disturbset_calc_init
  
end module disturbset
