subroutine DisturbEnvMPI( cfgfile, xz_PotTemp, xz_Exner, pz_VelX, xr_VelZ, xza_MixRt, xz_Km, xz_Kh )
!
!擾乱のデフォルト値を与えるためのルーチン.
!
!モジュール読み込み
use dc_message, only: MessageNotify
use gridset, only:DimXMin, DimXMax, RegXMin, RegXMax, DimZMin, DimZMax, MarginX, MarginZ, DelZ, SpcNum, XMin, XMax, ZMin, ZMax, s_X, s_Z ! Z 座標軸(スカラー格子点)
use filesetMPI, only:RandomFile ! 乱数ファイル
use basicset, only: SpcWetMolFr, MolWtWet, MolWtDry, xz_TempBasicZ, xz_PressBasicZ, xz_ExnerBasicZ, xza_MixRtBasicZ ! 基本場の混合比
use BoundaryMPI, only: BoundaryXCyc_xza , BoundaryZSym_xza
use ECCM, only: ECCM_MolFr
use mpiset, only: nprocs, myrank
!暗黙の型宣言禁止
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), 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) :: Random !ファイルから取得した乱数
! real(8) :: RandomNum(DimXMin:DimXMax)
real(8), allocatable :: RandomNum(:)
real(8), allocatable :: RandomNum2(:)
character(20) :: Type !温位擾乱のタイプ
real(8) :: xz_MolFr(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
!モル比
real(8) :: XMaxMPI, XMinMPI
integer :: RegXMaxMPI, RegXMinMPI
integer :: DimXMaxMPI, DimXMinMPI
real(8) :: s_Xmpi(DimXMin:DimXMax)
integer :: i, k, s, n ! DO ループ変数
!-------------------------------------------------------------
! 初期化
!-------------------------------------------------------------
!NAMELIST ファイルの読み込み
NAMELIST /disturbset/ Type, DelMax, XrRate, XcRate, ZrRate, ZcRate, Humidity, Xpos, Zpos
open (10, FILE=cfgfile)
read(10, NML=disturbset)
close(10)
!初期化
pz_VelX = 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
XMinMPI = XMin !領域の最小値
XMaxMPI = nprocs * (Xmax - Xmin) + XMinMPI !領域の最大値
s_Xmpi = s_X + myrank * (Xmax - Xmin) !当該 CPU で担当する水平領域
RegXMinMPI = RegXMin
RegXMaxMPI = nprocs * (RegXMax - RegXMin) + RegXMinMPI
DimXMinMPI = RegXMinMPI - MarginX
DimXMaxMPI = RegXMaxMPI + MarginX
allocate(RandomNum(DimXMinMPI:DimXmaxMPI))
allocate(RandomNum2(DimXMinMPI:DimXmaxMPI))
Xr = XMaxMPI * XrRate
Xc = XMaxMPI * XcRate
! Xr = minval( s_X, 1, s_X > (XMax - XMin) * XrRate )
! Xc = minval( s_X, 1, s_X > (XMax - XMin) * XcRate )
Zr = minval( s_Z, 1, s_Z > (ZMax - ZMin) * ZrRate )
Zc = minval( s_Z, 1, s_Z > (ZMax - ZMin) * ZcRate )
! write(*,*) "DisturbEnvMPI: ", minval(s_Xmpi), maxval(s_Xmpi)
! write(*,*) "DisturbEnvMPI: ", nprocs, XmaxMPI
! write(*,*) "DisturbEnvMPI: ", Xr, Xc
!-------------------------------------------------------------
! 温位の擾乱
!-------------------------------------------------------------
select case(Type)
case ("Thermal-KW1978")
! 指定された領域内に温位擾乱を与える (Klemp and Wilhelmson, 1978)
do k = DimZMin, DimZMax
do i = DimXMin, DimXMax
beta(i,k) = ( ( ( s_Xmpi(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_Xmpi(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 ("Thermal-MixRt-Gauss")
! 指定された座標を中心としたガウス分布の温位擾乱と混合比擾乱を与える.
! ただし, 混合比擾乱の最大値は, 基本場の混合比の 0.01 倍とする.
do k = DimZMin, DimZMax
do i = DimXMin, DimXMax
xz_PotTemp(i,k) = DelMax * dexp( - ( (s_X(i) - Xc) / Xr )**3.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 = DimXMinMPI, DimXmaxMPI
read(20,*) random
RandomNum(i) = random
end do
close(20)
do i = DimXMinMPI, DimXmaxMPI
!擾乱が全体としてはゼロとなるように調整
RandomNum2(i) = RandomNum(i) - sum( RandomNum(RegXMinMPI+1:RegXMaxMPI) ) / real(RegXMaxMPI - RegXMinMPI, 8)
end do
do i = DimXMin, DimXMax
k = i + (RegXMax - RegXMin) * myrank
xz_PotTemp(i, maxloc(s_Z, s_Z <= Zpos) - MarginZ - 1) = DelMax * RandomNum2(k) / xz_ExnerBasicZ(i, maxloc(s_Z, s_Z <= Zpos) - MarginZ - 1)
write(*,*) "myrank", myrank, i, RandomNum2(k)
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
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 DisturbEnvMPI