!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004, 2005. All rights reserved.
!---------------------------------------------------------------------
!= Subroutine DisturbEnv
!
!   * 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 を利用
!

subroutine DisturbEnv( cfgfile,         &
  & pz_VelX,     xr_VelZ,               &
  & xz_Exner,    xz_PotTemp,   xz_MixRt )
  !
  !デフォルトの擾乱を与えるためのルーチン. 
  !
  
  !モジュール読み込み
  use gridset,    only:DimXMin,         &! 配列の X 方向の下限
    &                  DimXMax,         &! 配列の X 方向の上限
    &                  DimZMin,         &! 配列の Z 方向の下限
    &                  DimZMax,         &! 配列の Z 方向の上限
    &                  RegZMin,         &! 物理領域の Z 方向の下限
    &                  SpcNum,          &! 凝縮性化学種の数
    &                  XMin, XMax,      &! 
    &                  ZMin, ZMax,      &! 
    &                  s_X,             &! X 座標軸(スカラー格子点)
    &                  s_Z               ! Z 座標軸(スカラー格子点)
  use fileset,    only:RandomFile        ! 乱数ファイル
  use basicset,   only:                  &
    &                  SpcWetMolFr,     &!凝縮成分の初期モル比
    &                  MolWtWet,        &!凝縮成分の分子量
    &                  xz_TempBasicZ,   &! 基本場の温度
    &                  xz_PressBasicZ,  &! 基本場の圧力
    &                  xz_MixRtBasicZ    ! 基本場の混合比
  use Boundary, only:  xz_Boundary_Cyc_Sym 
  use MoistFunc,only:  xz_MolFr2MixRt
  use ECCM,     only:  ECCM_MolFr


  !暗黙の型宣言禁止
  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_MixRt(DimXMin:DimXMax,DimZMin:DimZMax, SpcNum)  
                                    !凝縮成分の混合比(擾乱成分)
  real(8), parameter         :: Pi = 3.1415926535897932385d0 
                                    !円周率
  integer       :: i, k, s
  real(8)       :: Humidity         !相対湿度
  real(8)       :: XcRate           !擾乱の中心位置(水平方向)の領域に対する割合
  real(8)       :: XrRate           !擾乱の半径(水平方向)の領域に対する割合
  real(8)       :: ZcRate           !擾乱の中心位置(鉛直方向)の領域に対する割合
  real(8)       :: ZrRate           !擾乱の半径(鉛直方向)の領域に対する割合
  real(8)       :: Xc               !擾乱の中心位置(水平方向)
  real(8)       :: Xr               !擾乱の半径(水平方向)
  real(8)       :: Zc               !擾乱の中心位置(鉛直方向)
  real(8)       :: Zr               !擾乱の半径(鉛直方向)
  real(8)       :: beta(DimXMin:DimXMax, DimZMin:DimZMax)
                                    !擾乱の最大値に対する割合
  real(8)       :: DelMax           !温位擾乱の最大値
  real(8)       :: Random           !ファイルから取得した乱数
  character(20) :: Type             !温位擾乱のタイプ
  real(8)       :: xz_MolFr(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                    !モル比

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

  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_MixRt   = 0.0d0

  !-------------------------------------------------------------
  ! 温位の擾乱
  !-------------------------------------------------------------
  select case(Type)

  case ("Thermal")
    ! 指定された場所に温位の 1 つの熱的サーマルを与える. 
    
    Xr = (XMax - XMin) * XrRate
    Xc = (XMax - XMin) * XcRate
    Zr = (ZMax - ZMin) * ZrRate
    Zc = (ZMax - ZMin) * ZcRate
    
    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 )
      xz_PotTemp = DelMax * ( dcos( Pi * 5.0d-1 * beta ) ** 2.0d0 )
    end where


  case ("Thermal-Gauss")
    ! 指定された場所に, ガウシアンな温位の擾乱を与える. 

    Xr = (XMax - XMin) * XrRate
    Xc = (XMax - XMin) * XcRate
    Zr = (ZMax - ZMin) * ZrRate
    Zc = (ZMax - ZMin) * ZcRate
    
    do i = DimXMin, DimXMax
      xz_PotTemp(i,:) = &
        & DelMax * dexp( - ( (s_X(i) - Xc) / Xr )**2.0d0 * 5.0d-1 ) 
    end do
    
    do k = DimZMin, DimZMax
      xz_PotTemp(:,k) =  &
        & xz_PotTemp(:,k) * dexp( - ( (s_Z(k) - Zc) / Zr )**2.0d0 * 5.0d-1 ) 
    end do
    
    where ( xz_PotTemp < DelMax * 1.0d-4)
      xz_PotTemp = 0.0d0 
    end where
    
  case ("Exner-Gauss")
    ! 指定された場所に, ガウシアンな温位の擾乱を与える. 

    Xr = (XMax - XMin) * XrRate
    Xc = (XMax - XMin) * XcRate
    Zr = (ZMax - ZMin) * ZrRate
    Zc = (ZMax - ZMin) * ZcRate
    
    do i = DimXMin, DimXMax
      xz_Exner(i,:) = &
        & DelMax * dexp( - ( (s_X(i) - Xc) / Xr )**2.0d0 * 5.0d-1 ) 
    end do
    
    do k = DimZMin, DimZMax
      xz_Exner(:,k) =  &
        & xz_Exner(:,k) * dexp( - ( (s_Z(k) - Zc) / Zr )**2.0d0 * 5.0d-1 ) 
    end do
    
    where ( xz_Exner < DelMax * 1.0d-4)
      xz_Exner = 0.0d0 
    end where
    
  case ("Thermal-Random")
    ! 最下層にランダムな擾乱を与える

    open(20,file=RandomFile)
    do i = DimXMin, DimXMax
      read(20,*) random
      xz_Exner(i, RegZMin+1) = DelMax * random
    end do
    close(20)
    
  end select

  
  !-------------------------------------------------------------
  ! 蒸気圧の設定. 
  !-------------------------------------------------------------
  do i = DimXMin, DimXMax      
    call ECCM_MolFr( SpcWetMolFr(1:SpcNum), Humidity, xz_TempBasicZ(i,:), &
      &              xz_PressBasicZ(i,:), xz_MolFr(i,:,:) )
  end do
    
  !気相のモル比を混合比に変換
  do s = 1, SpcNum
    xz_MixRt(:,:,s) = xz_MolFr2MixRt(MolWtWet(s), xz_MolFr(:,:,s)) - xz_MixRtBasicZ(:,:,s)
  end do

  !値が小さくなりすぎないように最低値を与える
  where (xz_MixRt <= 1.0d-20 )
    xz_MixRt = 1.0d-20
  end where
  
  !境界条件
  do s = 1, SpcNum
    xz_MixRt(:,:,s) = xz_Boundary_Cyc_Sym(xz_MixRt(:,:,s))
  end do

end subroutine DisturbEnv