!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004, 2005. All rights reserved.
!---------------------------------------------------------------------
!= Subroutine DisturbEnv
!
!   * Developer: SUGIYAMA Ko-ichiro, ODAKA Masatsugu
!   * Version: $Id: disturbenv_mmc.f90,v 1.1 2011-03-30 13:53:07 sugiyama Exp $ 
!   * Tag Name: $Name: arare4-20120911 $
!   * Change History: 
!
!== Overview 
!
!擾乱のデフォルト値を与えるためのルーチン. 
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
! * 火星湿潤対流計算用
!
!== Future Plans
!
!  * 設定可能な擾乱のタイプを増やす
!  * エラー処理に gf4f90io を利用
!

subroutine DisturbEnv_mmc(                                   &
  &    cfgfile, xz_PotTemp, xz_Exner, pz_VelX, xr_VelZ,  &
  &    xz_Km, xz_Kh, xz_DensCloud, xz_SatRatio           &
  &  )
  !
  !擾乱のデフォルト値を与えるためのルーチン. 
  !
  
  !モジュール読み込み
  use dc_message, only: MessageNotify

  use gridset,    only:DimXMin,         &! 配列の X 方向の下限
    &                  DimXMax,         &! 配列の X 方向の上限
    &                  RegXMin,         &! 配列の X 方向の下限
    &                  RegXMax,         &! 配列の X 方向の上限
    &                  DimZMin,         &! 配列の Z 方向の下限
    &                  DimZMax,         &! 配列の Z 方向の上限
    &                  RegZMin,         &! 配列の Z 方向の下限
    &                  RegZMax,         &! 配列の Z 方向の上限
    &                  MarginZ,         &! 計算領域のマージン
    &                  DelZ,            &! Z 方向の刻み幅
    &                  SpcNum,          &! 凝縮性化学種の数
    &                  XMin, XMax,      &! 
    &                  ZMin, ZMax,      &! 
    &                  s_X,             &! X 座標軸(スカラー格子点)
    &                  s_Z               ! Z 座標軸(スカラー格子点)
  use fileset,    only:RandomFile        ! 乱数ファイル
  use basicset,   only:                 &
    &                  SpcWetMolFr,     &!凝縮成分の初期モル比
    &                  MolWtWet,        &!凝縮成分の分子量
    &                  MolWtDry,        &!乾燥成分の分子量
    &                  xz_PotTempBasicZ, &! 基本場の温位
    &                  xz_TempBasicZ,   &! 基本場の温度
    &                  xz_PressBasicZ,  &! 基本場の圧力
    &                  xz_ExnerBasicZ,  &! 無次元圧力
    &                  TempSfc,         &! 地表面温度
    &                  PressSfc,        &! 地表面圧力
    &                  GasRDry,         &! 乾燥成分の気体定数
    &                  CpDry,           &! 乾燥成分の定圧比熱
    &                  Grav              ! 重力加速度
  use ChemData, only:  ChemData_SVapPress_AntoineA, &!Antoine の式の係数
    &                  ChemData_SVapPress_AntoineB   !Antoine の式の係数

  !暗黙の型宣言禁止
  implicit none
  
  !変数定義
  character(*), intent(in) :: cfgfile
  real(8), intent(out)  :: pz_VelX(DimXMin:DimXMax,DimZMin:DimZMax)  
                                    !水平風速の擾乱成分
  real(8), intent(out)  :: xr_VelZ(DimXMin:DimXMax,DimZMin:DimZMax) 
                                    !鉛直風速の擾乱成分 
  real(8), intent(out)  :: xz_Exner(DimXMin:DimXMax,DimZMin:DimZMax)  
                                    !エクスナー関数の擾乱成分 
  real(8), intent(out)  :: xz_PotTemp(DimXMin:DimXMax,DimZMin:DimZMax)  
                                    !温位の擾乱成分 
  real(8), intent(out)  :: xz_Km(DimXMin:DimXMax,DimZMin:DimZMax)
                                    !運動量に対する拡散係数
  real(8), intent(out)  :: xz_Kh(DimXMin:DimXMax,DimZMin:DimZMax)
                                    !熱に対する拡散係数
  real(8), intent(out)  :: xz_DensCloud(DimXMin:DimXMax,DimZMin:DimZMax)
                                    !雲密度
  real(8), intent(out)  :: xz_SatRatio(DimXMin:DimXMax,DimZMin:DimZMax)
                                    !飽和比
  real(8), parameter         :: Pi = 3.1415926535897932385d0 
                                    !円周率
  real(8)       :: XcRate = 0.0d0   !擾乱の中心位置(水平方向)の領域に対する割合
  real(8)       :: XrRate = 0.0d0   !擾乱の半径(水平方向)の領域に対する割合
  real(8)       :: ZcRate = 0.0d0   !擾乱の中心位置(鉛直方向)の領域に対する割合
  real(8)       :: ZrRate = 0.0d0   !擾乱の半径(鉛直方向)の領域に対する割合
  real(8)       :: Xc     = 0.0d0   !擾乱の中心位置(水平方向)
  real(8)       :: Xr     = 0.0d0   !擾乱の半径(水平方向)
  real(8)       :: Zc     = 0.0d0   !擾乱の中心位置(鉛直方向)
  real(8)       :: Zr     = 0.0d0   !擾乱の半径(鉛直方向)
  real(8)       :: Zpos   = 0.0d0   !擾乱の Z 座標 [m] (Therma-Random 用)
  real(8)       :: DelMax = 0.0d0   !温位擾乱の最大値
  real(8)       :: HalfWidth = 0.0d0!温位擾乱の半値幅
  real(8)       :: ShearWidth= 0.0d0!シアー層の幅
  real(8)       :: Random           !ファイルから取得した乱数
  real(8)       :: RandomNum(DimXMin:DimXMax)
  real(8)       :: RandomNum2(DimXMin:DimXMax)
  character(20) :: Type             !温位擾乱のタイプ
  integer       :: i, k, s, n       ! DO ループ変数

  !-------------------------------------------------------------
  ! 初期化
  !-------------------------------------------------------------
  !NAMELIST ファイルの読み込み
  NAMELIST /disturbset/ &
    & Type, DelMax, XrRate, XcRate, ZrRate, ZcRate, &
    & Zpos, HalfWidth, ShearWidth

  open (10, FILE=cfgfile)
  read(10, NML=disturbset)
  close(10)
  
  !初期化
  pz_VelX    = 0.0d0
  xr_VelZ    = 0.0d0
  xz_Exner   = 0.0d0
  xz_PotTemp = 0.0d0
  xz_Km      = 0.0d0
  xz_Kh      = 0.0d0
  xz_DensCloud = 0.0d0
  xz_SatRatio = 0.0d0

  Xr = (XMax - XMin) * XrRate
  Xc = (XMax - XMin) * XcRate
  Zr = (ZMax - ZMin) * ZrRate
  Zc = (ZMax - ZMin) * ZcRate

  !-------------------------------------------------------------
  ! 温位の擾乱
  !-------------------------------------------------------------
  select case(Type)
    
  case ("Kitamori-Gauss")
    ! 北守(2006) の計算設定. 
    ! 指定された座標を中心としたガウス分布の温位擾乱を与える. 
    ! 基本場圧力分布と飽和蒸気圧から飽和比を与える. 

    do k = DimZMin, DimZMax
      do i = DimXMin, DimXMax
        
        xz_PotTemp(i,k) = &
          & DelMax * dexp( - ( (s_X(i) - Xc) / Xr )**2.0d0 * 5.0d-1   &
          &                - ( (s_Z(k) - Zc) / Zr )**2.0d0 * 5.0d-1 ) 
        
        xz_SatRatio(i,k) = &
          & PressSfc * (xz_ExnerBasicZ(i,k) + xz_Exner(i,k))**(CpDry / GasRDry) &
          & * ( exp(  &
          &         - ChemData_SVapPress_AntoineA(12) &
          &         + ChemData_SVapPress_AntoineB(12) &
          &           / ( (xz_ExnerBasicZ(i,k) + xz_Exner(i,k)) &
          &              * (xz_PotTempBasicZ(i,k) + xz_PotTemp(i,k)) ) &
          &         ) &
          &   )
        
      end do
    end do
    
  case ("Thermal-Random")
    ! 下層にランダムな擾乱を与える
    
    open(20,file=RandomFile)
    do i = DimXMin, DimXMax
      read(20,*) random
      RandomNum(i) = random
    end do
    close(20)
    
    do i = DimXMin, DimXMax
      
      !擾乱が全体としてはゼロとなるように調整
      RandomNum2(i) = RandomNum(i)  &
        & - sum( RandomNum(RegXMin+1:RegXMax) ) / real(RegXMax - RegXMin, 8) 
      
      xz_PotTemp(i, maxloc(s_Z, s_Z <= Zpos) - MarginZ - 1) = &
        & DelMax * RandomNum2(i) / xz_ExnerBasicZ(i, maxloc(s_Z, s_Z <= Zpos) - MarginZ - 1)
      
!      write(*,*) "xz_PotTemp()", xz_PotTemp(i, - MarginZ -1)
    end do
    
    
  case ("SK1994")
    ! Skamarock and Klemp (1994) の内部重力波実験
    
    do k = DimZMin, DimZMax
      do i = DimXMin, DimXMax
        xz_PotTemp(i,k) =                             &
          & DelMax * HalfWidth ** 2.0d0               &
          & * sin(pi * s_Z(k) / ZMax)                 &
          & / (HalfWidth ** 2.0d0 +                   &
          & (s_X(i) - XcRate * (XMax - XMin)) ** 2.0d0 )
      end do
    end do
    do k = DimZMin, DimZMax
      do i = DimXMin, DimXMax
        pz_VelX(i,k)    =   20.0d0
      end do
    end do
    
  case ("KH")
     ! KH 不安定実験(温位擾乱を与えるタイプ)
    
    do k = DimZMin, DimZMax
      do i = DimXMin, DimXMax
        pz_VelX(i,k)    =  4.0d0 *                           &
          & (1.0d0 + tanh((s_Z(k) - ZcRate * (ZMax - ZMin)) / ShearWidth))
      end do
    end do
    do k = DimZMin, DimZMax
      do i = DimXMin, DimXMax
        xz_PotTemp(i,k)  =  1.0d0 *                           &
          !          & (1.0d0 + tanh((s_Z(k) - ZcRate * (ZMax - ZMin)) / ShearWidth))
          & sin(s_X(i)*2.0d0*pi/XMax)*sin(s_Z(k)*2.0d0*pi/ZMax)
      end do
    end do
    
  case ("KH-2")
    ! KH 不安定実験(温位擾乱を与えないタイプ)
    
    do k = DimZMin, DimZMax
      do i = DimXMin, DimXMax
        pz_VelX(i,k)    =  3.0d0 *                             &
          & (1.0d0 + tanh((s_Z(k) - ZcRate * (ZMax - ZMin)) / ShearWidth)) &
          & + 0.2d0 * (1.0d0 +                                  &
          & sin(s_X(i)*2.0d0*pi/XMax)*sin(s_Z(k)*2.0d0*pi/ZMax))
      end do
    end do
    
  end select

end subroutine DisturbEnv_mmc
