subroutine DisturbEnv_3d( cfgfile, xyz_PotTemp, xyz_Exner, pyz_VelX, xqz_VelY, xyr_VelZ, xyza_MixRt, xyz_Km, xyz_Kh )
!
!擾乱のデフォルト値を与えるためのルーチン.
!
!モジュール読み込み
use dc_types, only : DP
use dc_iounit, only: FileOpen
use dc_message, only: MessageNotify
use gridset_3d, only:DimXMin, DimXMax, DimYMin, DimYMax, DimZMin, DimZMax, Zmargin, SpcNum, XMin, XMax, YMin, YMax, ZMin, ZMax, x_X, y_Y, 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 ) = 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