!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004, 2005. All rights reserved.
!---------------------------------------------------------------------
!= Subroutine DisturbEnv
!
!   * Developer: SUGIYAMA Ko-ichiro, ODAKA Masatsugu
!   * Version: $Id: disturbenv.f90,v 1.14 2006/11/16 15:54:28 sugiyama Exp $ 
!   * Tag Name: $Name: arare4-20061224 $
!   * Change History: 
!
!== Overview 
!
!擾乱のデフォルト値を与えるためのルーチン. 
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!== Future Plans
!
!  * 設定可能な擾乱のタイプを増やす
!  * エラー処理に gf4f90io を利用
!

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

  use gridset,    only:DimXMin,         &! 配列の X 方向の下限
    &                  DimXMax,         &! 配列の X 方向の上限
    &                  RegXMin,         &! 配列の X 方向の下限
    &                  RegXMax,         &! 配列の X 方向の上限
    &                  DimZMin,         &! 配列の Z 方向の下限
    &                  DimZMax,         &! 配列の 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_TempBasicZ,   &! 基本場の温度
    &                  xz_PressBasicZ,  &! 基本場の圧力
    &                  xz_ExnerBasicZ,  &! 無次元圧力
    &                  xza_MixRtBasicZ   ! 基本場の混合比
  use Boundary, only:  xza_BoundaryXCyc_xza , &!
    &                  xza_BoundaryZSym_xza 
  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)  :: xza_MixRt(DimXMin:DimXMax,DimZMin:DimZMax, SpcNum)  
                                    !凝縮成分の混合比(擾乱成分)
  real(8), intent(out)  :: xz_Km(DimXMin:DimXMax,DimZMin:DimZMax)
                                    !運動量に対する拡散係数
  real(8), intent(out)  :: xz_Kh(DimXMin:DimXMax,DimZMin:DimZMax)
                                    !熱に対する拡散係数
  real(8), parameter         :: Pi = 3.1415926535897932385d0 
                                    !円周率
  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)       :: Xpos             !擾乱の X 座標 [m] (Therma-Random 用)
  real(8)       :: Zpos             !擾乱の Z 座標 [m] (Therma-Random 用)
  real(8)       :: beta(DimXMin:DimXMax, DimZMin:DimZMax)
                                    !擾乱の最大値に対する割合
  real(8)       :: DelMax           !温位擾乱の最大値
  real(8)       :: Random           !ファイルから取得した乱数
  real(8)       :: RandomNum(DimXMin:DimXMax)
  real(8)       :: RandomNum2(DimXMin:DimXMax)
  character(20) :: Type             !温位擾乱のタイプ
  real(8)       :: xz_MolFr(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                    !モル比

  integer       :: i, k, s, n       ! DO ループ変数



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

  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
  xza_MixRt  = 0.0d0
  xz_Km      = 0.0d0
  xz_Kh      = 0.0d0

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

  case ("Thermal-KW1978")
    ! 指定された領域内に温位擾乱を与える (Klemp and Wilhelmson, 1978) 
    
    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 ) &
        &          / xz_ExnerBasicZ
    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 ) &
        & / xz_ExnerBasicZ(i,:)
    end do
    
    do k = DimZMin, DimZMax
      xz_PotTemp(:,k) =  &
        & xz_PotTemp(:,k) * dexp( - ( (s_Z(k) - Zc) / Zr )**2.0d0 * 5.0d-1 ) &
        & / xz_ExnerBasicZ(:,k)
    end do
    
    where ( sign(1.0d0, DelMax) * xz_PotTemp < DelMax * 1.0d-4 )
      xz_PotTemp = 0.0d0 
    end where
    

  case ("Thermal-MixRt-Gauss")
    ! 指定された座標を中心としたガウス分布の温位擾乱と混合比擾乱を与える. 
    ! ただし, 混合比擾乱の最大値は, 基本場の混合比の 0.01 倍とする.
    
    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 ) &
        & / xz_ExnerBasicZ(i,:)
      do s = 1, SpcNum
        xza_MixRt(i,:,s) =                                     &
          & maxval( xza_MixRtBasicZ(:,:,s) ) * 1.0d-1          &
          & * dexp( - ( (s_X(i) - Xc) / Xr )**2.0d0 * 5.0d-1 ) &
          & / xz_ExnerBasicZ(i,:)
      end do
    end do
    
    do k = DimZMin, DimZMax
      xz_PotTemp(:,k) =  &
        & xz_PotTemp(:,k) * dexp( - ( (s_Z(k) - Zc) / Zr )**2.0d0 * 5.0d-1 ) &
        & / xz_ExnerBasicZ(:,k)
      do s = 1, SpcNum
        xza_MixRt(:,k,s) =                                     &
          & xza_MixRt(:,k,s)                                   &
          & * dexp( - ( (s_Z(k) - Zc) / Zr )**2.0d0 * 5.0d-1 ) &
          & / xz_ExnerBasicZ(:,k)
      end do
    end do
    
!    where ( sign(1.0d0, DelMax) * 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
      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 <= 90.0d3) - MarginZ) = &
!        & DelMax * RandomNum2(i) / xz_ExnerBasicZ(i, maxloc(s_Z, s_Z <= 100.0d3) - MarginZ)

      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)

    end do
 

  case ("HS2001")
    ! Hueso and Sanchez-Lavega を模した設定

    i = ( DimXMax - DimXMin - 10) / 2 
    k = minloc( s_Z, 1, s_Z > 2.5d4 ) - MarginZ
    n = int( 5.0d3 / DelZ )

    xz_PotTemp(i-n:i,k-n:k) = DelMax


  case ("SK1989")
    ! Skamarock and Klemp (1989) の Cold-bubble 実験
    
    Xr = (XMax - XMin) * XrRate
    Xc = (XMax - XMin) * XcRate
    Zr = (ZMax - ZMin) * ZrRate
    Zc = (ZMax - ZMin) * ZcRate
    
    xz_PotTemp = 0.0d0

    do k = DimZMin, DimZMax
      do i = DimXMin, DimXMax

        beta(i,k) =                                 &
          & sqrt(                                   &
          &      ( ( s_X(i) - Xc ) / Xr ) ** 2.0d0  &
          &    + ( ( s_Z(k) - Zc ) / Zr ) ** 2.0d0  &
          &  )
      end do
    end do

    where ( beta < 1.0d0 )
      xz_PotTemp = 0.5d0*DelMax*(cos(pi*beta) + 1.0d0)
    end where

  end select

  !-------------------------------------------------------------
  ! 蒸気圧の設定. 
  !-------------------------------------------------------------
  if (Humidity /= 0.0d0) then 
    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
      xza_MixRt(:,:,s) =                                                    &
        & xz_MolFr(:,:,s) * MolWtWet(s) / MolWtDry - xza_MixRtBasicZ(:,:,s)
    end do
  end if
  
  
  !値が小さくなりすぎないように最低値を与える
!  where (xza_MixRt <= 1.0d-20 )
!    xza_MixRt = 1.0d-20
!  end where
  
  !境界条件
  xza_MixRt = xza_BoundaryXCyc_xza( xza_MixRt )
  xza_MixRt = xza_BoundaryZSym_xza( xza_MixRt )
  

end subroutine DisturbEnv
