subroutine DisturbEnv( 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, 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
!変数定義
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) :: DryHeight2 !乾燥を与える高度の下限
real(8) :: Us !シアーの強度
real(8) :: Hb !シアーを与える高度の下限
real(8) :: Hu !シアーを与える高度の上限
real(8) :: beta(DimXMin:DimXMax, DimZMin:DimZMax)
!擾乱の最大値に対する割合
real(8) :: DelMax !温位擾乱の最大値
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, B, U ! DO ループ変数
!-------------------------------------------------------------
! 初期化
!-------------------------------------------------------------
!NAMELIST ファイルの読み込み
NAMELIST /disturbset/ Type, DelMax, XrRate, XcRate, ZrRate, ZcRate, Humidity, Xpos, Zpos, DryHeight2, Us, Hu, Hb
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
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_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 ("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
end do
close(20)
do i = DimXMin, DimXMax
!擾乱が全体としてはゼロとなるように調整
RandomNum2(i) = RandomNum(i) - sum( RandomNum(RegXMin+1:RegXMax) ) / real(RegXMax - RegXMin, 8)
! 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)
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 )
!-------------------------------------------------------------!
! シアーの設定 (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 DisturbEnv