subroutine DisturbEnvMPI( myrank, nprocs, 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, RegZMin, RegZMax, MarginX, MarginZ, DelZ, SpcNum, XMin, XMax, ZMin, ZMax, s_X, s_Z ! Z 座標軸(スカラー格子点)
use fileset, only:RandomFile ! 乱数ファイル
use basicset, only: SpcWetMolFr, MolWtWet, MolWtDry, xz_TempBasicZ, xz_PressBasicZ, xz_ExnerBasicZ, xza_MixRtBasicZ ! 基本場の混合比
use Boundary, only: BoundaryXCyc_xza , BoundaryZSym_xza
use ECCM, only: ECCM_MolFr
!暗黙の型宣言禁止
implicit none
!変数定義
integer, intent(in) :: myrank, nprocs
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 /disturbenv/ Type, Humidity
NAMELIST /disturbenv_thermal/ DelMax, XrRate, XcRate, ZrRate, ZcRate
NAMELIST /disturbenv_random/ DelMax, Zpos
NAMELIST /disturbenv_windshear/ Us, ZposMin, ZposMax
NAMELIST /disturbenv_dryregion/ XposMin, XposMax, ZposMin, ZposMax
!初期化
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
! namelist の読み込み
open (10, FILE=cfgfile)
read(10, NML=disturbset)
close(10)
! 領域を決める
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 - XMinMPI) * XrRate
Xc = (XMaxMPI - XMinMPI) * 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 )
!-------------------------------------------------------------
! 温位の擾乱
!-------------------------------------------------------------
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 )**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 = 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)
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
! DryHeight2 以上の水平領域の半分の初期の湿度をゼロにするために
! 基本場と逆符号の水蒸気擾乱を与える
! (水平方向の格子点数/2)をif文内に直接書いているので格子点数が変わったときは
! 注意が必要
do k = RegZMin+1,DimZMax
do i = DimXMin,DimXMax
if (s_Z(k).gt.DryHeight2.and.i.ge.DimXMin.and.i.le.512) then
xza_MixRt(i,k,1) = - xza_MixRtBasicZ(i,k,1)
end if
end do
end do
!境界条件
call BoundaryXCyc_xza( xza_MixRt )
call BoundaryZSym_xza( xza_MixRt )
! xza_MixRt = xza_BoundaryXCyc_xza( xza_MixRt )
! xza_MixRt = xza_BoundaryZSym_xza( xza_MixRt )
!-------------------------------------------------------------!
! シアーの設定 (Takemi,2007) !
!-------------------------------------------------------------!
!
!= 概要
!* case "Takemi2007" での計算時に鉛直シアーのある風を与える時に使用する
!* 風の与え方には, 以下のようなバリエーションがある
! (1) シアーを与える高度を変える
! (2) シアーのある風の最大風速 (U_s) を変える
!
! (1) については, (a) 0 - 2.5 km, (b) 2.5 - 5.0 km, (c) 5.0 - 7.5 km の
! 三パターンがある
! (2) については, Takemi (2007) では熱帯場と中緯度場の温度場毎に
! 異なる値を設定している
!
! その強度(Us)は, 以下の通り
! <熱帯場> (1) 5 m/s, (2) 10 m/s, (3) 15 m/s
! <中緯度場> (1) 10 m/s, (2) 15 m/s, (3) 20 m/s
!
!* シアーの形の模式図 (Takemi, 2007) |
! /| 7.5 km
! / |
! / |
! / ←|
! ｜ /| 5.0 km
! ｜ / |
! ｜/ |
! / ←|
! ｜ /| 2.5 km
! ｜ / |
! ｜/ |
! / ←|
!---------------------------------+------------- 0.0 km
! Us (m/s)
!
! シアーを与える高度を決める
! 傾き始めの高度の格子点 (B) と 傾きの終わりの格子点 (U) で決定
do k = RegZMin+1, DimZMax
if (s_Z(k).ge.Hb.and.s_Z(k).le.(Hb + DelZ)) then
B = k
end if
if (s_Z(k).ge.Hu.and.s_Z(k).le.(Hu + DelZ)) then
U = k
end if
end do
! 上記で決めた格子点BからUの間に鉛直シアーのある風を与える
do k = RegZMin+1, U
do i = DimXMin, DimXMax
if (k.lt.B) then
pz_VelX(i,k) = - Us
else
pz_VelX(i,k) = - Us + (Us / (s_Z(U) - s_Z(B))) * (s_Z(k) - s_Z(B))
end if
end do
end do
end subroutine DisturbEnvMPI