subroutine disturbset_calc_init(cfgfile)
!=== Dependency
!=== Input
character(*), intent(in) :: cfgfile
!=end
!=== Work
character(80) :: ExnerDisturbType !無次元圧力の擾乱 (DisturbType.f90 参照)
real(8) :: ExnerSigmaX !擾乱の幅 (X 方向)
real(8) :: ExnerSigmaZ !擾乱の幅 (Z 方向)
real(8) :: ExnerMax !擾乱の最大値
real(8) :: ExnerCenterX !擾乱の中心 (X 方向)
real(8) :: ExnerCenterZ !擾乱の中心 (X 方向)
character(80) :: PotTempDisturbType !温位の擾乱 (DisturbType.f90 参照)
real(8) :: PotTempSigmaX !擾乱の幅 (X 方向)
real(8) :: PotTempSigmaZ !擾乱の幅 (X 方向)
real(8) :: PotTempMax !擾乱の最大値
real(8) :: PotTempCenterX !擾乱の中心 (X 方向)
real(8) :: PotTempCenterZ !擾乱の中心 (X 方向)
character(80) :: VelDisturbType !速度の擾乱 (DisturbType.f90 参照)
real(8) :: VelX !速度 u
real(8) :: VelZ !速度 w
integer :: i, k
real(8) :: mu
real(8) :: sigma
!=end
!=begin
!=== NAMELIST
NAMELIST /disturbset/ ExnerDisturbType, ExnerMax, ExnerSigmaX, ExnerSigmaZ, ExnerCenterX, ExnerCenterZ, PotTempDisturbType, PotTempMax, PotTempSigmaX, PotTempSigmaZ, PotTempCenterX, PotTempCenterZ, VelDisturbType, VelX, VelZ
!=end
call BeginSub("disturbset_calc_init", fmt="%c", c1="Initialize values of disturbance. ")
open (10, FILE=cfgfile)
read(10, NML=disturbset)
close(10)
!==== 変数の値のチェック
if ( ExnerSigmaX > 1 ) then
call MessageNotify("Message", "disturbset_calc_init", "SigmaX is > 1")
stop
end if
if ( ExnerSigmaZ > 1 ) then
call MessageNotify("Message", "disturbset_calc_init", "SigmaZ is > 1")
stop
end if
if ( PotTempSigmaX > 1 ) then
call MessageNotify("Message", "disturbset_calc_init", "SigmaX is > 1")
stop
end if
if ( PotTempSigmaZ > 1 ) then
call MessageNotify("Message", "disturbset_calc_init", "SigmaZ is > 1")
stop
end if
if ( ExnerCenterX > 1 ) then
call MessageNotify("Message", "disturbset_calc_init", "CenterX is > 1")
stop
end if
if ( ExnerCenterZ > 1 ) then
call MessageNotify("Message", "disturbset_calc_init", "CenterZ is > 1")
stop
end if
if ( PotTempCenterX > 1 ) then
call MessageNotify("Message", "disturbset_calc_init", "CenterX is > 1")
stop
end if
if ( PotTempCenterZ > 1 ) then
call MessageNotify("Message", "disturbset_calc_init", "CenterZ is > 1")
stop
end if
!==== 配列の割り当て, 初期化
allocate( ss_PotTempDisturb(DimXMin:DimXMax, DimZMin:DimZMax), ss_ExnerDisturb(DimXMin:DimXMax, DimZMin:DimZMax), fs_VelXDisturb(DimXMin:DimXMax, DimZMin:DimZMax), sf_VelZDisturb(DimXMin:DimXMax, DimZMin:DimZMax) )
ss_ExnerDisturb = 0.0d0
ss_PotTempDisturb = 0.0d0
fs_VelXDisturb = 0.0d0
sf_VelZDisturb = 0.0d0
!-----------------------------------------------------------
!--- エクスナー関数の擾乱
!-----------------------------------------------------------
select case (ExnerDisturbType)
case ("Exner_Gauss_X")
!--- X 方向のガウシアン分布
mu = ( Xmax - Xmin ) * ExnerCenterX
sigma = ( Xmax - Xmin ) * ExnerSigmaX
do i = DimXMin, DimXMax
ss_ExnerDisturb(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_ExnerDisturb(:,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_ExnerDisturb(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_ExnerDisturb(:,k) = ss_ExnerDisturb(:,k) * dexp(-((s_Z(k) - mu) / sigma)**2.0d0 / 2.0d0)
end do
end select
!-----------------------------------------------------------
!--- 仮温位の擾乱
!-----------------------------------------------------------
select case (PotTempDisturbType)
case ("PotTemp_Gauss_X")
!--- X 方向のガウシアン分布
mu = ( Xmax - Xmin ) * PotTempCenterX
sigma = ( Xmax - Xmin ) * PotTempSigmaX
do i = DimXMin, DimXMax
ss_PotTempDisturb(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_PotTempDisturb(:,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_PotTempDisturb(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_PotTempDisturb(:,k) = ss_PotTempDisturb(:,k) * dexp(-((s_Z(k) - mu) / sigma)**2.0d0 / 2.0d0)
end do
! where (ss_PotTempDisturb < PotTempMax * 1.0d-4)
! ss_PotTempDisturb = 0.0d0
! end where
case default
ss_PotTempDisturb = 0.0d0
end select
!----------------------------------------------------------
!--- 速度
!----------------------------------------------------------
select case (VelDisturbType)
case ("Zero")
fs_VelXDisturb = 0.0d0
sf_VelZDisturb = 0.0d0
case ("VelX_test")
fs_VelXDisturb = VelX
fs_VelXDisturb(-5:10, :) = 0.0d0
fs_VelXDisturb(90:105, :) = 0.0d0
sf_VelZDisturb = 0.0d0
case ("VelX_Uniform")
!--- 速度 u 一定
fs_VelXDisturb = VelX
sf_VelZDisturb = 0.0d0
case ("VelZ_Uniform")
!--- 速度 w 一定
fs_VelXDisturb = 0.0d0
sf_VelZDisturb = VelZ
case ("Solid_Rotation")
!--- 領域中心の剛体回転
do i = DimXMin, DimXMax
fs_VelXDisturb(i,:) = 2.0d0*pi / 1000.0d0 * ( f_Z(:) - (Zmax + Zmin) / 2.0D0 )
end do
do k = DimZMin, DimZMax
sf_VelZDisturb(:,k) = - 2.0d0 * pi / 1000.0d0 * ( s_X(:) - (Xmax + Xmin) / 2.0D0 )
end do
end select
call EndSub("disturbset_calc_init")
end subroutine disturbset_calc_init