!c Description: 
!c   環境場の計算
!c
!c Current Code Owner: 
!c   sugiyama@gfd-dennou.org
!c
!c Copyright (C) SUGIYAMA Ko-ichiro, 2004, All rights reserved


subroutine DistrbType(ss_ExnerDistrb, ss_ThetaDistrb, &
     &               fs_VelXDistrb,  sf_VelZDistrb  )

  use gridset, only: s_X, f_X, s_Z, f_Z, Xmin, Xmax, Zmin, Zmax, &
       &             DimXMin, DimXMax, DimZMin, DimZMax,      &
       &             RegXMin, RegXMax, RegZMin, RegZMax
  use gacchanprm, only: MolWtDry, CpDry, &
       & ExnerDistrbType, ExnerMax, ExnerSigmaX, ExnerSigmaZ, &
       &                  ExnerCenterX, ExnerCenterZ, &
       & ThetaDistrbType, ThetaMax, ThetaSigmaX, ThetaSigmaZ, &
       &                  ThetaCenterX, ThetaCenterZ, &
       & VelDistrbType, VelX, VelZ, AnglVel
  use if_avr

  implicit none

  real(8), intent(out)  :: ss_ThetaDistrb(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(out)  :: ss_ExnerDistrb(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(out)  :: fs_VelXDistrb(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(out)  :: sf_VelZDistrb(DimXMin:DimXMax, DimZMin:DimZMax)

  integer :: i, k
  real(8) :: mu
  real(8) :: sigma

  !--- 初期化
  ss_ExnerDistrb = 0.0d0
  ss_ThetaDistrb = 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 (ThetaDistrbType)

  case ("Theta_Gauss_X")
     !--- X 方向のガウシアン分布
     mu = ( Xmax - Xmin ) * ThetaCenterX
     sigma = ( Xmax - Xmin ) * ThetaSigmaX

     do i = DimXMin, DimXMax
        ss_ThetaDistrb(i,:) = &
             & ThetaMax * dexp(-((s_X(i) - mu) / sigma)**2.0d0 / 2.0d0) 
     end do
     
  case ("Theta_Gauss_Z")
     !--- Z 方向のガウシアン分布
     mu = ( Zmax - Zmin ) * ThetaCenterX
     sigma = ( Zmax - Zmin ) * ThetaSigmaZ
     
     do k = DimZMin, DimZMax
        ss_ThetaDistrb(:,k) = &
             & ThetaMax * dexp(-((s_Z(k) - mu) / sigma)**2.0d0 / 2.0d0) 
     end do
     
  case ("Theta_Gauss_XZ")     
     !--- 同心円状のガウシアン分布
     mu = ( Xmax - Xmin ) * ThetaCenterX
     sigma = ( Xmax - Xmin ) * ThetaSigmaX
     
     do i = DimXMin, DimXMax
        ss_ThetaDistrb(i,:) = &
             & ThetaMax * dexp(-((s_X(i) - mu) / sigma)**2.0d0 / 2.0d0) 
     end do
     
     mu = ( Zmax - Zmin ) * ThetaCenterZ
     sigma = ( Zmax - Zmin ) * ThetaSigmaZ

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

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

  case ("Zero")
     fs_VelXDistrb = 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 ("Rotate")
     !--- 領域の中心を軸とする回転流
     do k = DimZMin, DimZMax
        fs_VelXDistrb(:,k) = - AnglVel * (f_Z(k) - (ZMax - ZMin) * 5.0d-1) 
     end do
     do i = DimXMin, DimXMax
        sf_VelZDistrb(i,:) = AnglVel * (f_X(i) - (XMax - XMin) * 5.0d-1)
     end do

  case ("RotateQuarter")
     !--- 領域中に 4 個の回転流を置く(周期境界を用いても境界でシアができない)
     fs_VelXDistrb = 0.0d0
     sf_VelZDistrb = 0.0d0
     do k = RegZMin, RegZMax/2
        fs_VelXDistrb(DimXMin:RegXMax/2,k)  =   &
             & - AnglVel * (f_Z(k) - (ZMax - ZMin) * 2.5d-1) 
     end do
     do i = RegXMin, RegXMax/2
        sf_VelZDistrb(i,RegZMin:RegZMax/2) =    &
             &   AnglVel * (f_X(i) - (XMax - XMin) * 2.5d-1)
     end do

     fs_VelXDistrb(RegXMax/2+1:RegXMax,RegZMax/2+1:RegZMax) =    &
          & fs_VelXDistrb(RegXMin:RegXMax/2,RegZMin:RegZMax/2)
     sf_VelZDistrb(RegXMax/2+1:RegXMax,RegZMax/2+1:RegZMax) =    &
          & sf_VelZDistrb(RegXMin:RegXMax/2,RegZMin:RegZMax/2)

     fs_VelXDistrb(RegXMax/2+1:RegXMax,RegZMin:RegZMax/2) =    &
          & - fs_VelXDistrb(RegXMin:RegXMax/2,RegZMin:RegZMax/2)
     sf_VelZDistrb(RegXMax/2+1:RegXMax,RegZMin:RegZMax/2) =    &
          & - sf_VelZDistrb(RegXMin:RegXMax/2,RegZMin:RegZMax/2)

     fs_VelXDistrb(RegXMin:RegXMax/2,RegZMax/2+1:RegZMax) =    &
          & - fs_VelXDistrb(RegXMin:RegXMax/2,RegZMin:RegZMax/2)
     sf_VelZDistrb(RegXMin:RegXMax/2,RegZMax/2+1:RegZMax) =    &
          & - sf_VelZDistrb(RegXMin:RegXMax/2,RegZMin:RegZMax/2)

  end select
        
   
end subroutine DistrbType

