!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Subroutine DistrbType
!
!   * Developer: SUGIYAMA Ko-ichiro (sugiyama@gfd-dennou.org)
!   * Version: $Id: distrbtype.f90,v 1.3 2005/02/02 03:51:31 odakker Exp $ 
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview 
!
!典型的な擾乱場を計算するためのサブルーチン. タイプとしては以下が設定可能. 
!   * 無次元圧力: 擾乱なし or X 方向ガウシアン分布 or Z 方向ガウシアン分布 or XZ 方向ガウシアン分布 
!   * 温位: 擾乱なし or X 方向ガウシアン分布 or Z 方向ガウシアン分布 or XZ 方向ガウシアン分布 
!   * 速度: 擾乱なし or X 方向速度一定 or Z 方向速度一定 
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!== Future Plans
!
!設定可能な擾乱のタイプを増やす
!
                                                                  !=end

subroutine DistrbType(ss_ExnerDistrb, ss_PotTempDistrb, &
     &               fs_VelXDistrb,  sf_VelZDistrb  )
                                                                  !=begin
  !== Dependency
  use physprm 
  use gridset, only: &
       & s_X, s_Z, f_X, f_Z, Xmin, Xmax, Zmin, Zmax, &
       & DimXMin, DimXMax, DimZMin, DimZMax
  use distrbprm, only: &
       & MolWtDry, CpDry, &
       & ExnerDistrbType, ExnerMax, ExnerSigmaX, ExnerSigmaZ, &
       &                  ExnerCenterX, ExnerCenterZ, &
       & PotTempDistrbType, PotTempMax, PotTempSigmaX, PotTempSigmaZ, &
       &                  PotTempCenterX, PotTempCenterZ, &
       & VelDistrbType, VelX, VelZ
                                                                  !=end
  implicit none
                                                                  !=begin
  !== Output
  real(8), intent(out)  :: ss_PotTempDistrb(DimXMin:DimXMax, DimZMin:DimZMax)   !温位
  real(8), intent(out)  :: ss_ExnerDistrb(DimXMin:DimXMax, DimZMin:DimZMax)   !エクスナー関数
  real(8), intent(out)  :: fs_VelXDistrb(DimXMin:DimXMax, DimZMin:DimZMax)    !速度 u
  real(8), intent(out)  :: sf_VelZDistrb(DimXMin:DimXMax, DimZMin:DimZMax)    !速度 w
                                                                  !=end
  integer :: i, k
  real(8) :: mu
  real(8) :: sigma

  !--- 初期化
  ss_ExnerDistrb = 0.0d0
  ss_PotTempDistrb = 0.0d0
  fs_VelXDistrb = 0.0d0
  sf_VelZDistrb = 0.0d0
  

  !-----------------------------------------------------------
  !--- エクスナー関数の擾乱
  !-----------------------------------------------------------
  select case (ExnerDistrbType)

  case ("Exner_Gauss_X")
     !--- X 方向のガウシアン分布
     mu = ( Xmax - Xmin ) * ExnerCenterX
     sigma = ( Xmax - Xmin ) * ExnerSigmaX

     do i = DimXMin, DimXMax
        ss_ExnerDistrb(i,:) = &
             & ExnerMax * dexp(-((s_X(i) - mu) / sigma )**2.0d0 / 2.0d0) 
     end do

  case ("Exner_Gauss_Z")    
     !--- Z 方向のガウシアン分布
     mu = ( Zmax - Zmin ) * ExnerCenterZ
     sigma = ( Zmax - Zmin ) * ExnerSigmaZ

     do k = DimZMin, DimZMax
        ss_ExnerDistrb(:,k) = &
             & ExnerMax * dexp(-((s_Z(k) - mu) / sigma)**2.0d0 / 2.0d0) 
     end do

  case ("Exner_Gauss_XZ")
     !--- 同心円状のガウシアン分布
     mu = ( Xmax - Xmin ) * ExnerCenterX
     sigma = ( Xmax - Xmin ) * ExnerSigmaX

     do i = DimXMin, DimXMax
        ss_ExnerDistrb(i,:) = &
             & ExnerMax * dexp(-((s_X(i) - mu) / sigma)**2.0d0 / 2.0d0) 
     end do

     mu = ( Zmax - Zmin ) * ExnerCenterZ
     sigma = ( Zmax - Zmin ) * ExnerSigmaZ
     
     do k = DimZMin, DimZMax
        ss_ExnerDistrb(:,k) = ss_ExnerDistrb(:,k) &
             & * dexp(-((s_Z(k) - mu) / sigma)**2.0d0 / 2.0d0) 
     end do
  end select

  !-----------------------------------------------------------
  !--- 仮温位の擾乱
  !-----------------------------------------------------------
  select case (PotTempDistrbType)

  case ("PotTemp_Gauss_X")
     !--- X 方向のガウシアン分布
     mu = ( Xmax - Xmin ) * PotTempCenterX
     sigma = ( Xmax - Xmin ) * PotTempSigmaX

     do i = DimXMin, DimXMax
        ss_PotTempDistrb(i,:) = &
             & PotTempMax * dexp(-((s_X(i) - mu) / sigma)**2.0d0 / 2.0d0) 
     end do
     
  case ("PotTemp_Gauss_Z")
     !--- Z 方向のガウシアン分布
     mu = ( Zmax - Zmin ) * PotTempCenterX
     sigma = ( Zmax - Zmin ) * PotTempSigmaZ
     
     do k = DimZMin, DimZMax
        ss_PotTempDistrb(:,k) = &
             & PotTempMax * dexp(-((s_Z(k) - mu) / sigma)**2.0d0 / 2.0d0) 
     end do
     
  case ("PotTemp_Gauss_XZ")     
     !--- 同心円状のガウシアン分布
     mu = ( Xmax - Xmin ) * PotTempCenterX
     sigma = ( Xmax - Xmin ) * PotTempSigmaX
     
     do i = DimXMin, DimXMax
        ss_PotTempDistrb(i,:) = &
             & PotTempMax * dexp(-((s_X(i) - mu) / sigma)**2.0d0 / 2.0d0) 
     end do
     
     mu = ( Zmax - Zmin ) * PotTempCenterZ
     sigma = ( Zmax - Zmin ) * PotTempSigmaZ

     do k = DimZMin, DimZMax
        ss_PotTempDistrb(:,k) =  &
             & ss_PotTempDistrb(:,k) &
             &   * dexp(-((s_Z(k) - mu) / sigma)**2.0d0 / 2.0d0) 
     end do

!     where (ss_PotTempDistrb < PotTempMax * 1.0d-4)
!        ss_PotTempDistrb = 0.0d0 
!     end where
 
   case default

     ss_PotTempDistrb = 0.0d0

  end select

     
  !----------------------------------------------------------
  !--- 速度
  !----------------------------------------------------------
  select case (VelDistrbType)

  case ("Zero")
     fs_VelXDistrb = 0.0d0
     sf_VelZDistrb = 0.0d0

  case ("VelX_test")
     fs_VelXDistrb = VelX
     fs_VelXDistrb(-5:10, :) = 0.0d0
     fs_VelXDistrb(90:105, :) = 0.0d0
     sf_VelZDistrb = 0.0d0

  case ("VelX_Uniform")
     !--- 速度 u 一定
     fs_VelXDistrb = VelX
     sf_VelZDistrb = 0.0d0

  case ("VelZ_Uniform")
     !--- 速度 w 一定
     fs_VelXDistrb = 0.0d0
     sf_VelZDistrb = VelZ

  case ("Solid_Rotation")
     !--- 領域中心の剛体回転

     do i = DimXMin, DimXMax   
       fs_VelXDistrb(i,:) =   2.0d0*pi/1000.0d0 * ( f_Z(:) - (Zmax + Zmin)/2.0D0 )
     end do 

     do k = DimZMin, DimZMax
       sf_VelZDistrb(:,k) = - 2.0d0*pi/1000.0d0 * ( s_X(:) - (Xmax + Xmin)/2.0D0 )
     end do     

  end select
   
   
end subroutine DistrbType

