subroutine DisturbEnv( cfgfile, xz_PotTemp, xz_Exner, pz_VelX, xr_VelZ, xza_MixRt, xz_Km, xz_Kh, xz_DensCloud, xz_SatRatio )
!
!擾乱のデフォルト値を与えるためのルーチン.
!
!モジュール読み込み
use dc_message, only: MessageNotify
use gridset, only:DimXMin, DimXMax, RegXMin, RegXMax, DimZMin, DimZMax, RegZMin, RegZMax, MarginZ, DelZ, SpcNum, XMin, XMax, ZMin, ZMax, s_X, s_Z ! Z 座標軸(スカラー格子点)
use fileset, only:RandomFile ! 乱数ファイル
use basicset, only: SpcWetMolFr, MolWtWet, MolWtDry, xz_PotTempBasicZ, xz_TempBasicZ, xz_PressBasicZ, xz_ExnerBasicZ, xza_MixRtBasicZ, TempSfc, PressSfc, GasRDry, CpDry, Grav ! 重力加速度
use ChemData, only: ChemData_SVapPress_AntoineA, ChemData_SVapPress_AntoineB !Antoine の式の係数
use Boundary, only: BoundaryXCyc_xza , BoundaryZSym_xza
use ECCM, only: ECCM_MolFr
!暗黙の型宣言禁止
implicit none
!変数定義
character(*), intent(in) :: cfgfile
real(8), intent(out) :: pz_VelX(DimXMin:DimXMax,DimZMin:DimZMax)
!水平風速の擾乱成分
real(8), intent(out) :: xr_VelZ(DimXMin:DimXMax,DimZMin:DimZMax)
!鉛直風速の擾乱成分
real(8), intent(out) :: xz_Exner(DimXMin:DimXMax,DimZMin:DimZMax)
!エクスナー関数の擾乱成分
real(8), intent(out) :: xz_PotTemp(DimXMin:DimXMax,DimZMin:DimZMax)
!温位の擾乱成分
real(8), intent(out) :: xza_MixRt(DimXMin:DimXMax,DimZMin:DimZMax, SpcNum)
!凝縮成分の混合比(擾乱成分)
real(8), intent(out) :: xz_Km(DimXMin:DimXMax,DimZMin:DimZMax)
!運動量に対する拡散係数
real(8), intent(out) :: xz_Kh(DimXMin:DimXMax,DimZMin:DimZMax)
!熱に対する拡散係数
real(8), intent(out) :: xz_DensCloud(DimXMin:DimXMax,DimZMin:DimZMax)
!雲密度
real(8), intent(out) :: xz_SatRatio(DimXMin:DimXMax,DimZMin:DimZMax)
!飽和比
real(8), parameter :: Pi = 3.1415926535897932385d0
!円周率
real(8) :: Humidity !相対湿度
real(8) :: XcRate !擾乱の中心位置(水平方向)の領域に対する割合
real(8) :: XrRate !擾乱の半径(水平方向)の領域に対する割合
real(8) :: ZcRate !擾乱の中心位置(鉛直方向)の領域に対する割合
real(8) :: ZrRate !擾乱の半径(鉛直方向)の領域に対する割合
real(8) :: Xc !擾乱の中心位置(水平方向)
real(8) :: Xr !擾乱の半径(水平方向)
real(8) :: Zc !擾乱の中心位置(鉛直方向)
real(8) :: Zr !擾乱の半径(鉛直方向)
real(8) :: Xpos !擾乱の X 座標 [m] (Therma-Random 用)
real(8) :: Zpos !擾乱の Z 座標 [m] (Therma-Random 用)
real(8) :: beta(DimXMin:DimXMax, DimZMin:DimZMax)
!擾乱の最大値に対する割合
real(8) :: DelMax !温位擾乱の最大値
real(8) :: HalfWidth !温位擾乱の半値幅
real(8) :: ShearWidth !シアー層の幅
real(8) :: Random !ファイルから取得した乱数
real(8) :: RandomNum(DimXMin:DimXMax)
real(8) :: RandomNum2(DimXMin:DimXMax)
character(20) :: Type !温位擾乱のタイプ
real(8) :: xz_MolFr(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
!モル比
integer :: i, k, s, n ! DO ループ変数
!-------------------------------------------------------------
! 初期化
!-------------------------------------------------------------
!NAMELIST ファイルの読み込み
NAMELIST /disturbset/ Type, DelMax, XrRate, XcRate, ZrRate, ZcRate, Humidity, Xpos, Zpos, HalfWidth, ShearWidth
open (10, FILE=cfgfile)
read(10, NML=disturbset)
close(10)
!初期化
pz_VelX = 0.0d0
! pz_VelX = 20.0d0
! xr_VelZ = 0.0d0
xr_VelZ = 0.0d0
xz_Exner = 0.0d0
xz_PotTemp = 0.0d0
xza_MixRt = 0.0d0
xz_Km = 0.0d0
xz_Kh = 0.0d0
xz_DensCloud = 0.0d0
xz_SatRatio = 0.0d0
! Xr = minval( s_X, 1, s_X > (XMax - XMin) * XrRate )
Xr = (XMax - XMin) * XrRate
! Xc = minval( s_X, 1, s_X > (XMax - XMin) * XcRate )
Xc = (XMax - XMin) * XcRate
! Zr = minval( s_Z, 1, s_Z > (ZMax - ZMin) * ZrRate )
Zr = (ZMax - ZMin) * ZrRate
! Zc = minval( s_Z, 1, s_Z > (ZMax - ZMin) * ZcRate )
Zc = (ZMax - ZMin) * ZcRate
!-------------------------------------------------------------
! 温位の擾乱
!-------------------------------------------------------------
select case(Type)
case ("Thermal-KW1978")
! 指定された領域内に温位擾乱を与える (Klemp and Wilhelmson, 1978)
do k = DimZMin, DimZMax
do i = DimXMin, DimXMax
beta(i,k) = ( ( ( s_X(i) - Xc ) / Xr ) ** 2.0d0 + ( ( s_Z(k) - Zc ) / Zr ) ** 2.0d0 ) ** 5.0d-1
end do
end do
where ( beta < 1.0d0 )
xz_PotTemp = DelMax * ( dcos( Pi * 5.0d-1 * beta ) ** 2.0d0 ) / xz_ExnerBasicZ
end where
case ("Thermal-Gauss")
! 指定された座標を中心としたガウス分布の温位擾乱を与える.
do k = DimZMin, DimZMax
do i = DimXMin, DimXMax
xz_PotTemp(i,k) = DelMax * dexp( - ( (s_X(i) - Xc) / Xr )**2.0d0 * 5.0d-1 - ( (s_Z(k) - Zc) / Zr )**2.0d0 * 5.0d-1 ) / xz_ExnerBasicZ(i,k)
end do
end do
where ( sign(1.0d0, DelMax) * xz_PotTemp < DelMax * 1.0d-4 )
xz_PotTemp = 0.0d0
end where
case ("Kitamori-Gauss")
! 北守(2006) の計算設定.
! 指定された座標を中心としたガウス分布の温位擾乱を与える.
! 基本場圧力分布と飽和蒸気圧から飽和比を与える.
do k = DimZMin, DimZMax
do i = DimXMin, DimXMax
! do k = RegZMin, RegZMax
! do i = RegXMin, RegXMax
! xz_PotTemp(i,k) = &
! & DelMax * dexp( - ( (s_X(i) - Xc) / Xr )**2.0d0 * 5.0d-1 &
! & - ( (s_Z(k) - Zc) / Zr )**2.0d0 * 5.0d-1 ) &
! & / xz_ExnerBasicZ(i,k)
xz_PotTemp(i,k) = DelMax * dexp( - ( (s_X(i) - Xc) / Xr )**2.0d0 * 5.0d-1 - ( (s_Z(k) - Zc) / Zr )**2.0d0 * 5.0d-1 )
! xz_SatRatio(i,k) = &
!! & PressSfc * (TempSfc / xz_TempBasicZ(i,k) )**(Grav / CpDry ) &
! & xz_PressBasicZ(i,k) &
! & / exp( ChemData_SVapPress_AntoineA(12) - &
! & ChemData_SVapPress_AntoineB(12) / ( xz_TempBasicZ(i,k) &
! & * (1.0d0 + xz_PotTemp(i,k) / TempSfc) ) )
xz_SatRatio(i,k) = PressSfc * (xz_ExnerBasicZ(i,k) + xz_Exner(i,k))**(CpDry / GasRDry) * ( exp( - ChemData_SVapPress_AntoineA(12) + ChemData_SVapPress_AntoineB(12) / ( (xz_ExnerBasicZ(i,k) + xz_Exner(i,k)) * (xz_PotTempBasicZ(i,k) + xz_PotTemp(i,k)) ) ) )
end do
end do
case ("Thermal-MixRt-Gauss")
! 指定された座標を中心としたガウス分布の温位擾乱と混合比擾乱を与える.
! ただし, 混合比擾乱の最大値は, 基本場の混合比の 0.01 倍とする.
do k = DimZMin, DimZMax
do i = DimXMin, DimXMax
xz_PotTemp(i,k) = DelMax * dexp( - ( (s_X(i) - Xc) / Xr )**2.0d0 * 5.0d-1 - ( (s_Z(k) - Zc) / Zr )**2.0d0 * 5.0d-1 ) / xz_ExnerBasicZ(i,k)
do s = 1, SpcNum
xza_MixRt(i,:,s) = maxval( xza_MixRtBasicZ(:,:,s) ) * 1.0d-1 * dexp( - ( (s_X(i) - Xc) / Xr )**2.0d0 * 5.0d-1 - ( (s_Z(k) - Zc) / Zr )**2.0d0 * 5.0d-1 ) / xz_ExnerBasicZ(i,:)
end do
end do
end do
! where ( sign(1.0d0, DelMax) * xz_PotTemp < DelMax * 1.0d-4 )
! xz_PotTemp = 0.0d0
! end where
case ("Exner-Gauss")
! 指定された場所に, ガウシアンな温位の擾乱を与える.
do k = DimZMin, DimZMax
do i = DimXMin, DimXMax
xz_Exner(i,k) = DelMax * dexp( - ( (s_X(i) - Xc) / Xr )**2.0d0 * 5.0d-1 - ( (s_Z(k) - Zc) / Zr )**2.0d0 * 5.0d-1 )
end do
end do
! where ( xz_Exner < DelMax * 1.0d-4)
! xz_Exner = 0.0d0
! end where
case ("Thermal-Random")
! 下層にランダムな擾乱を与える
open(20,file=RandomFile)
do i = DimXMin, DimXMax
read(20,*) random
RandomNum(i) = random
! write(*,*) RandomNum(i)
end do
close(20)
do i = DimXMin, DimXMax
!擾乱が全体としてはゼロとなるように調整
RandomNum2(i) = RandomNum(i) - sum( RandomNum(RegXMin+1:RegXMax) ) / real(RegXMax - RegXMin, 8)
! write(*,*) "RandomNum2", RandomNum2(i)
! xz_PotTemp(i, maxloc(s_Z, s_Z <= 90.0d3) - MarginZ) = &
! & DelMax * RandomNum2(i) / xz_ExnerBasicZ(i, maxloc(s_Z, s_Z <= 100.0d3) - MarginZ)
xz_PotTemp(i, maxloc(s_Z, s_Z <= Zpos) - MarginZ - 1) = DelMax * RandomNum2(i) / xz_ExnerBasicZ(i, maxloc(s_Z, s_Z <= Zpos) - MarginZ - 1)
write(*,*) "xz_PotTemp()", xz_PotTemp(i, - MarginZ -1)
end do
case ("HS2001")
! Hueso and Sanchez-Lavega を模した設定
i = ( DimXMax - DimXMin - 10) / 2
k = minloc( s_Z, 1, s_Z > 2.5d4 ) - MarginZ
n = int( 5.0d3 / DelZ )
xz_PotTemp(i-n:i,k-n:k) = DelMax
! case ("SK1989")
! ! Skamarock and Klemp (1989) の Cold-bubble 実験
!
! xz_PotTemp = 0.0d0
!
! do k = DimZMin, DimZMax
! do i = DimXMin, DimXMax
!
! beta(i,k) = &
! & sqrt( &
! & ( ( s_X(i) - Xc ) / Xr ) ** 2.0d0 &
! & + ( ( s_Z(k) - Zc ) / Zr ) ** 2.0d0 &
! & )
! end do
! end do
!
! where ( beta < 1.0d0 )
! xz_PotTemp = 0.5d0*DelMax*(cos(pi*beta) + 1.0d0)
! end where
case ("SK1994")
! Skamarock and Klemp (1994) の内部重力波実験
do k = DimZMin, DimZMax
do i = DimXMin, DimXMax
xz_PotTemp(i,k) = DelMax * HalfWidth ** 2.0d0 * sin(pi * s_Z(k) / ZMax) / (HalfWidth ** 2.0d0 + (s_X(i) - XcRate * (XMax - XMin)) ** 2.0d0 )
end do
end do
do k = DimZMin, DimZMax
do i = DimXMin, DimXMax
pz_VelX(i,k) = 20.0d0
end do
end do
case ("KH")
! KH 不安定実験(温位擾乱を与えるタイプ)
do k = DimZMin, DimZMax
do i = DimXMin, DimXMax
pz_VelX(i,k) = 4.0d0 * (1.0d0 + tanh((s_Z(k) - ZcRate * (ZMax - ZMin)) / ShearWidth))
end do
end do
do k = DimZMin, DimZMax
do i = DimXMin, DimXMax
xz_PotTemp(i,k) = 1.0d0 * sin(s_X(i)*2.0d0*pi/XMax)*sin(s_Z(k)*2.0d0*pi/ZMax)
end do
end do
case ("KH-2")
! KH 不安定実験(温位擾乱を与えないタイプ)
do k = DimZMin, DimZMax
do i = DimXMin, DimXMax
pz_VelX(i,k) = 3.0d0 * (1.0d0 + tanh((s_Z(k) - ZcRate * (ZMax - ZMin)) / ShearWidth)) + 0.2d0 * (1.0d0 + sin(s_X(i)*2.0d0*pi/XMax)*sin(s_Z(k)*2.0d0*pi/ZMax))
end do
end do
end select
!-------------------------------------------------------------
! 蒸気圧の設定.
!-------------------------------------------------------------
if (Humidity /= 0.0d0) then
do i = DimXMin, DimXMax
call ECCM_MolFr( SpcWetMolFr(1:SpcNum), Humidity, xz_TempBasicZ(i,:), xz_PressBasicZ(i,:), xz_MolFr(i,:,:) )
end do
!気相のモル比を混合比に変換
do s = 1, SpcNum
xza_MixRt(:,:,s) = xz_MolFr(:,:,s) * MolWtWet(s) / MolWtDry - xza_MixRtBasicZ(:,:,s)
end do
end if
!値が小さくなりすぎないように最低値を与える
! where (xza_MixRt <= 1.0d-20 )
! xza_MixRt = 1.0d-20
! end where
!境界条件
call BoundaryXCyc_xza( xza_MixRt )
call BoundaryZSym_xza( xza_MixRt )
! xza_MixRt = xza_BoundaryXCyc_xza( xza_MixRt )
! xza_MixRt = xza_BoundaryZSym_xza( xza_MixRt )
end subroutine DisturbEnv