!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004, 2005. All rights reserved.
!---------------------------------------------------------------------
!= Subroutine DisturbEnv_3d
!
!   * Developer: SUGIYAMA Ko-ichiro, ODAKA Masatsugu
!   * Version: $Id: disturbenv_3d.f90,v 1.4 2007/10/01 08:56:45 odakker Exp $ 
!   * Tag Name: $Name: arare4-20071012 $
!   * Change History: 
!
!== Overview 
!
!擾乱のデフォルト値を与えるためのルーチン. 
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!== Future Plans
!
!  * 設定可能な擾乱のタイプを増やす
!  * エラー処理に gf4f90io を利用
!

subroutine DisturbEnv_3d(                           &
  &    cfgfile, xyz_PotTemp, xyz_Exner,          &
  &    pyz_VelX, xqz_VelY, xyr_VelZ,             &
  &    xyz_Km, xyz_Kh                            &
  &  )
  !
  !擾乱のデフォルト値を与えるためのルーチン. 
  !
  
  !モジュール読み込み
  use dc_types, only : DP
  use dc_iounit,  only: FileOpen
  use dc_message, only: MessageNotify

  use gridset_3d, only:DimXMin,         &! 配列の X 方向の下限
    &                  DimXMax,         &! 配列の X 方向の上限
    &                  DimYMin,         &! 配列の Y 方向の下限
    &                  DimYMax,         &! 配列の Y 方向の上限
    &                  DimZMin,         &! 配列の Z 方向の下限
    &                  DimZMax,         &! 配列の Z 方向の上限
    &                  Zmargin,         &!  
    &                  SpcNum,          &! 凝縮性化学種の数
    &                  XMin, XMax,      &! 
    &                  YMin, YMax,      &! 
    &                  ZMin, ZMax,      &! 
    &                  x_X,             &! X 座標軸(スカラー格子点)
    &                  y_Y,             &! X 座標軸(スカラー格子点)
    &                  z_Z,             &! Z 座標軸(スカラー格子点)
    &                  xyz_X, xyz_Y, xyz_Z, &
    &                  x_dx, y_dy, z_dz
  use fileset_3d, only:RandomFile        ! 乱数ファイル
  use basicset_3d,only:                 &
    &                  SpcWetMolFr,     &!凝縮成分の初期モル比
    &                  MolWtWet,        &!凝縮成分の分子量
    &                  MolWtDry,        &!乾燥成分の分子量
    &                  xyz_TempBasicZ,  &! 基本場の温度
    &                  xyz_PressBasicZ, &! 基本場の圧力
    &                  xyz_ExnerBasicZ   ! 無次元圧力
!    &                  xyza_MixRtBasicZ  ! 基本場の混合比

  use xyz_module, only: BoundaryXCyc_xyz, &
    &                   BoundaryYCyc_xyz, &
    &                   BoundaryZSym_xyz 
!  use ECCM_3d,     only:  ECCM_MolFr


  !暗黙の型宣言禁止
  implicit none
  
  !変数定義
  character(*), intent(in) :: cfgfile
  real(DP), intent(out)  :: pyz_VelX(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)  
                                    !水平風速の擾乱成分
  real(DP), intent(out)  :: xqz_VelY(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)  
                                    !水平風速の擾乱成分
  real(DP), intent(out)  :: xyr_VelZ(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) 
                                    !鉛直風速の擾乱成分 
  real(DP), intent(out)  :: xyz_Exner(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)  
                                    !エクスナー関数の擾乱成分 
  real(DP), intent(out)  :: xyz_PotTemp(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)  
                                    !温位の擾乱成分 
!  real(DP), intent(out)  :: xyza_MixRt(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax,SpcNum)  
                                    !凝縮成分の混合比(擾乱成分)
  real(DP), intent(out)  :: xyz_Km(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                    !運動量に対する拡散係数
  real(DP), intent(out)  :: xyz_Kh(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                    !熱に対する拡散係数
  real(DP), parameter         :: Pi = 3.1415926535897932385d0 
                                    !円周率
  real(DP)       :: Humidity        !相対湿度
  real(DP)       :: XcRate          !擾乱の中心位置(水平方向)の領域に対する割合
  real(DP)       :: XrRate          !擾乱の半径(水平方向)の領域に対する割合
  real(DP)       :: YcRate          !擾乱の中心位置(水平方向)の領域に対する割合
  real(DP)       :: YrRate          !擾乱の半径(水平方向)の領域に対する割合
  real(DP)       :: ZcRate          !擾乱の中心位置(鉛直方向)の領域に対する割合
  real(DP)       :: ZrRate          !擾乱の半径(鉛直方向)の領域に対する割合
  real(DP)       :: Xc              !擾乱の中心位置(水平方向)
  real(DP)       :: Xr              !擾乱の半径(水平方向)
  real(DP)       :: Yc              !擾乱の中心位置(水平方向)
  real(DP)       :: Yr              !擾乱の半径(水平方向)
  real(DP)       :: Zc              !擾乱の中心位置(鉛直方向)
  real(DP)       :: Zr              !擾乱の半径(鉛直方向)
  real(DP)       :: Xpos            !擾乱の X 座標 [m] (Therma-Random 用)
  real(DP)       :: Ypos            !擾乱の Y 座標 [m] (Therma-Random 用)
  real(DP)       :: Zpos            !擾乱の Z 座標 [m] (Therma-Random 用)
  real(DP)       :: beta(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax)
                                    !擾乱の最大値に対する割合
  real(DP)       :: DelMax          !温位擾乱の最大値
  real(DP)       :: Random          !ファイルから取得した乱数
  real(DP)       :: RandomNum(DimXMin:DimXMax,DimYMin:DimYMax)
  real(DP)       :: RandomNum2(DimXMin:DimXMax,DimYMin:DimYMax)
  real(DP)       :: RandomMean       
  character(20)  :: Type            !温位擾乱のタイプ
!  real(DP)       :: xyza_MolFr(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum)
                                    !モル比

  integer       :: i, j, k, s, n    ! DO ループ変数
  integer       :: unit             ! 装置番号

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

  call FileOpen(unit, file=cfgfile, mode='r')
  read(unit, NML=disturbset)
  close(unit)

  !初期化
  pyz_VelX    = 0.0d0
  xqz_VelY    = 0.0d0
  xyr_VelZ    = 0.0d0
  xyz_Exner   = 0.0d0
  xyz_PotTemp = 0.0d0
!  xyza_MixRt  = 0.0d0
  xyz_Km      = 0.0d0
  xyz_Kh      = 0.0d0

  Xr = minval( x_X, 1, x_X > (XMax - XMin) * XrRate )
  Xc = minval( x_X, 1, x_X > (XMax - XMin) * XcRate )
  Yr = minval( y_Y, 1, y_Y > (YMax - YMin) * YrRate )
  Yc = minval( y_Y, 1, y_Y > (YMax - YMin) * YcRate )
  Zr = minval( z_Z, 1, z_Z > (ZMax - ZMin) * ZrRate )
  Zc = minval( z_Z, 1, z_Z > (ZMax - ZMin) * ZcRate )

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

  case ("Thermal-KW1978")
    ! 指定された領域内に温位擾乱を与える (Klemp and Wilhelmson, 1978) 
    
    beta =                               &
      & (                                       &
      &      ( ( xyz_X - Xc ) / Xr ) ** 2.0d0  &
      &    + ( ( xyz_Y - Yc ) / Yr ) ** 2.0d0  &
      &    + ( ( xyz_Z - Zc ) / Zr ) ** 2.0d0  &
      &  ) ** 5.0d-1
    
    where ( beta < 1.0d0 )
      xyz_PotTemp = DelMax * ( dcos( Pi * 5.0d-1 * beta ) ** 2.0d0 ) &
        &          / xyz_ExnerBasicZ
    end where
    

  case ("Thermal-Gauss")
    ! 指定された座標を中心としたガウス分布の温位擾乱を与える. 

    xyz_PotTemp = &
      & DelMax    &
      & * dexp( - ( (xyz_X - Xc) / Xr )**2.0d0 * 5.0d-1   &
      &         - ( (xyz_Y - Yc) / Yr )**2.0d0 * 5.0d-1   &
      &         - ( (xyz_Z - Zc) / Yr )**2.0d0 * 5.0d-1 ) &
      & / xyz_ExnerBasicZ

    where ( sign(1.0d0, DelMax) * xyz_PotTemp < DelMax * 1.0d-4 )
      xyz_PotTemp = 0.0d0 
    end where

  case ("Thermal-GaussXZ")
    ! 指定された座標を中心としたガウス分布の温位擾乱を与える. 

    xyz_PotTemp = &
      & DelMax    &
      & * dexp( - ( (xyz_X - Xc) / Xr )**2.0d0 * 5.0d-1   &
      &         - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) &
      & / xyz_ExnerBasicZ

    where ( sign(1.0d0, DelMax) * xyz_PotTemp < DelMax * 1.0d-4 )
      xyz_PotTemp = 0.0d0 
    end where
  
  case ("Thermal-GaussYZ")
    ! 指定された座標を中心としたガウス分布の温位擾乱を与える. 

    xyz_PotTemp = &
      & DelMax    &
      & * dexp( - ( (xyz_Y - Yc) / Yr )**2.0d0 * 5.0d-1   &
      &         - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) &
      & / xyz_ExnerBasicZ

    where ( sign(1.0d0, DelMax) * xyz_PotTemp < DelMax * 1.0d-4 )
      xyz_PotTemp = 0.0d0 
    end where
 
  case ("Adv-XZ-X")
    ! X 方向移流のテスト
    ! 指定された座標を中心としたガウス分布の温位擾乱を与える. 

    xyz_PotTemp = &
      & DelMax    &
      & * dexp( - ( (xyz_X - Xc) / Xr )**2.0d0 * 5.0d-1   & 
      &         - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) 

    where ( sign(1.0d0, DelMax) * xyz_PotTemp < DelMax * 1.0d-4 )
      xyz_PotTemp = 0.0d0 
    end where

    where ( xyz_PotTemp > DelMax  )
      xyz_PotTemp = DelMax
    end where

    pyz_VelX = 20.0d0

  case ("Adv-XZ-Z")
    ! Z 方向移流のテスト
    ! 指定された座標を中心としたガウス分布の温位擾乱を与える. 

    xyz_PotTemp = &
      & DelMax    &
      & * dexp( - ( (xyz_X - Xc) / Xr )**2.0d0 * 5.0d-1   & 
      &         - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) 

    where ( sign(1.0d0, DelMax) * xyz_PotTemp < DelMax * 1.0d-4 )
      xyz_PotTemp = 0.0d0 
    end where

    where ( xyz_PotTemp > DelMax  )
      xyz_PotTemp = DelMax
    end where

    xyr_VelZ = 20.0d0

  case ("Adv-YZ-Y")
    ! Y 方向移流のテスト
    ! 指定された座標を中心としたガウス分布の温位擾乱を与える. 

    xyz_PotTemp = &
      & DelMax    &
      & * dexp( - ( (xyz_Y - Yc) / Yr )**2.0d0 * 5.0d-1   & 
      &         - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) 

    where ( sign(1.0d0, DelMax) * xyz_PotTemp < DelMax * 1.0d-4 )
      xyz_PotTemp = 0.0d0 
    end where

    where ( xyz_PotTemp > DelMax  )
      xyz_PotTemp = DelMax
    end where

    xqz_VelY = 20.0d0

  case ("Adv-YZ-Z")
    ! Z 方向移流のテスト
    ! 指定された座標を中心としたガウス分布の温位擾乱を与える. 

    xyz_PotTemp = &
      & DelMax    &
      & * dexp( - ( (xyz_Y - Yc) / Yr )**2.0d0 * 5.0d-1   & 
      &         - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) 

    where ( sign(1.0d0, DelMax) * xyz_PotTemp < DelMax * 1.0d-4 )
      xyz_PotTemp = 0.0d0 
    end where

    where ( xyz_PotTemp > DelMax  )
      xyz_PotTemp = DelMax
    end where

    xyr_VelZ = 20.0d0

  case ("Adv-XZ-XZ")
    ! XZ 方向移流のテスト
    ! 指定された座標を中心としたガウス分布の温位擾乱を与える. 

    xyz_PotTemp = &
      & DelMax    &
      & * dexp( - ( (xyz_X - Xc) / Xr )**2.0d0 * 5.0d-1   & 
      &         - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) 

    where ( sign(1.0d0, DelMax) * xyz_PotTemp < DelMax * 1.0d-4 )
      xyz_PotTemp = 0.0d0 
    end where

    where ( xyz_PotTemp > DelMax  )
      xyz_PotTemp = DelMax
    end where

    pyz_VelX = 20.0d0
    xyr_VelZ = 20.0d0

  case ("Adv-YZ-YZ")
    ! YZ 方向移流のテスト
    ! 指定された座標を中心としたガウス分布の温位擾乱を与える. 

    xyz_PotTemp = &
      & DelMax    &
      & * dexp( - ( (xyz_Y - Yc) / Yr )**2.0d0 * 5.0d-1   & 
      &         - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) 

    where ( sign(1.0d0, DelMax) * xyz_PotTemp < DelMax * 1.0d-4 )
      xyz_PotTemp = 0.0d0 
    end where

    where ( xyz_PotTemp > DelMax  )
      xyz_PotTemp = DelMax
    end where

    xqz_VelY = 20.0d0
    xyr_VelZ = 20.0d0

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

    xyz_Exner = &
      & DelMax  &
      & * dexp( - ( (xyz_X - Xc) / Xr )**2.0d0 * 5.0d-1 ) &
      & * dexp( - ( (xyz_Y - Yc) / Yr )**2.0d0 * 5.0d-1 ) &
      & * dexp( - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) 

    where ( xyz_Exner < DelMax * 1.0d-4)
      xyz_Exner = 0.0d0 
    end where
    
  case ("Exner-GaussX")
    ! 指定された場所に, ガウシアンな圧力の擾乱を与える. 
    ! (X 方向一様)

    xyz_Exner = &
      & DelMax  &
      & * dexp( - ( (xyz_Y - Yc) / Yr )**2.0d0 * 5.0d-1 ) &
      & * dexp( - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) 

    where ( xyz_Exner < DelMax * 1.0d-4)
      xyz_Exner = 0.0d0 
    end where

  case ("Exner-GaussY")
    ! 指定された場所に, ガウシアンな圧力の擾乱を与える. 
    ! (Y 方向一様)

    xyz_Exner = &
      & DelMax  &
      & * dexp( - ( (xyz_X - Xc) / Xr )**2.0d0 * 5.0d-1 ) &
      & * dexp( - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) 

    where ( xyz_Exner < DelMax * 1.0d-4)
      xyz_Exner = 0.0d0 
    end where

  case ("Exner-GaussZ")
    ! 指定された場所に, ガウシアンな圧力の擾乱を与える. 
    ! (Z 方向一様)

    xyz_Exner = &
      & DelMax  &
      & * dexp( - ( (xyz_X - Xc) / Xr )**2.0d0 * 5.0d-1 ) &
      & * dexp( - ( (xyz_Y - Yc) / Yr )**2.0d0 * 5.0d-1 ) 

    where ( xyz_Exner < DelMax * 1.0d-4)
      xyz_Exner = 0.0d0 
    end where
    
  case ("Thermal-Random")
    ! 下層にランダムな擾乱を与える

    call FileOpen(unit, file=RandomFile, mode='r')
    do j = DimYMin, DimYMax
      do i = DimXMin, DimXMax
        read(unit,*) random
        RandomNum(i,j) = random
    end do
    end do
    close(unit)

    ! 乱数の平均値を計算
    RandomMean = 0.0d0

    do j = DimYMin, DimYMax
      do i = DimXMin, DimXMax
        RandomMean = RandomMean + &
          & (RandomNum(i,j)*x_dx(i)*y_dy(j))/((XMax-Xmin)*(YMax-YMin))
      end do
    end do

    do j = DimYMin, DimYMax
      do i = DimXMin, DimXMax
        !擾乱が全体としてはゼロとなるように調整
        RandomNum2(i,j) = RandomNum(i,j)  - RandomMean

        xyz_PotTemp(i,j,maxloc(z_Z, z_Z <= Zpos) - Zmargin - 1) = &
          & DelMax * RandomNum2(i,j) / xyz_ExnerBasicZ(i,j,maxloc(z_Z, z_Z <= Zpos) - Zmargin - 1)
        
      end do
    end do

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

    i = ( DimXMax - DimXMin - 10) / 2 
    j = ( DimYMax - DimYMin - 10) / 2 
    k = minloc( z_Z, 1, z_Z > 2.5d4 ) - Zmargin
    n = int( 5.0d3 / z_dz(1) )

    xyz_PotTemp(i-n:i,j-n:j,k-n:k) = DelMax


  case ("SK1989")
    ! Skamarock and Klemp (1989) の Cold-bubble 実験
    
    xyz_PotTemp = 0.0d0

    beta =                                 &
      & sqrt(                                   &
      &      ( ( xyz_X - Xc ) / Xr ) ** 2.0d0  &
      &    + ( ( xyz_Y - Yc ) / Yr ) ** 2.0d0  &
      &    + ( ( xyz_Z - Zc ) / Zr ) ** 2.0d0  &
      &  )

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

  end select

  !-------------------------------------------------------------
  ! 蒸気圧の設定. 
  !-------------------------------------------------------------
!  if (Humidity /= 0.0d0) then 
!    do i = DimXMin, DimXMax      
!      do j = DimYMin, DimYMax      
!        call ECCM_MolFr( SpcWetMolFr(1:SpcNum), Humidity, &
!          &              xyz_TempBasicZ(i,j,:),           &
!          &              xyz_PressBasicZ(i,j,:), xyza_MolFr(i,j,:,:) )
!      end do
!    end do
    
    !気相のモル比を混合比に変換
!    do s = 1, SpcNum
!      xyza_MixRt(:,:,:,s) =  &
!        & xyza_MolFr(:,:,:,s) &
!        & * MolWtWet(s) / MolWtDry - xyza_MixRtBasicZ(:,:,:,s)
!    end do
!  end if
  
  
  !値が小さくなりすぎないように最低値を与える
!  where (xza_MixRt <= 1.0d-20 )
!    xza_MixRt = 1.0d-20
!  end where
  
  !境界条件
!  do s =1, SpcNum
!    call BoundaryXCyc_xyz( xyza_MixRt(:,:,:,s) )
!    call BoundaryYCyc_xyz( xyza_MixRt(:,:,:,s) )
!    call BoundaryZSym_xyz( xyza_MixRt(:,:,:,s) )
!  end do

end subroutine DisturbEnv_3d