!--------------------------------------------------------------------- ! Copyright (C) GFD Dennou Club, 2004, 2005. All rights reserved. !--------------------------------------------------------------------- != Subroutine DisturbEnv ! ! * Developer: SUGIYAMA Ko-ichiro, ODAKA Masatsugu ! * Version: $Id: disturbenv.f90,v 1.19 2009-03-10 16:17:24 sugiyama Exp $ ! * Tag Name: $Name: arare4-20100306 $ ! * Change History: ! !== Overview ! !擾乱のデフォルト値を与えるためのルーチン. ! !== Error Handling ! !== Known Bugs ! !== Note ! !== Future Plans ! ! * 設定可能な擾乱のタイプを増やす ! * エラー処理に gf4f90io を利用 ! 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, &! 配列の X 方向の下限 & DimXMax, &! 配列の X 方向の上限 & RegXMin, &! 配列の X 方向の下限 & RegXMax, &! 配列の X 方向の上限 & RegZMin, &! 配列の Z 方向の下限 & DimZMin, &! 配列の Z 方向の下限 & DimZMax, &! 配列の Z 方向の上限 & MarginZ, &! 計算領域のマージン & DelZ, &! Z 方向の刻み幅 & SpcNum, &! 凝縮性化学種の数 & XMin, XMax, &! & ZMin, ZMax, &! & s_X, &! 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) :: 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) !モル比 real(8) :: Us ! 水平風の最大値 real(8) :: pz_XVel(DimXMin:DimXMax, DimZMin:DimZMax) integer :: i, k, s, n, B, U ! 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 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) & 0.0d0 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 では鉛直シアーを持つ西向き風を与えている !! シアーは以下の 3 種類の高度分布, !! 中緯度場・熱帯場にそれぞれ 3 種類ずつの強さ (U_s) を用意する !! !! <高度分布> (1) 0 - 2.5 km (低層), (2) 2.5 - 5.0 km (中層), (3) 5.0 - 7.5 km (高層). !! <シアー強度> !! <熱帯> (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 !! U_s (m/s) !! !! * 風速は深さ 2.5 km の層内だけに与える !! !!!---問題---!!! !! deepconv ではどっちが西でどっちが東なんだろう??? !! ひとまず, pz_VelX > 0 が西風 (西から東へ吹く) 逆を東風と考える !! そのため, 描画したときの絵の左が西, 右が東と思うことにする !!! <熱帯のシアー強度> !!! Us = 5.0d0 ! Us = 10.0d0 ! Us = 15.0d0 !!! <中緯度のシアー強度> !!! ! Us = 10.0d0 ! Us = 15.0d0 ! Us = 20.0d0 !!! <シアーを置く高度と深さの設定> !!! do k = RegZMin+1, DimZMax if (s_Z(k).ge.2500.and.s_Z(k).le.(2500 + DelZ)) then B = k end if if (s_Z(k).ge.5000.and.s_Z(k).le.(5000 + DelZ)) then U = k end if end do !!! <シアーの計算> !!! do k = RegZMin+1, DimZMax do i = DimXMin, DimXMax if (k.ge.B.and.k.le.U) then 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