!--------------------------------------------------------------------- ! Copyright (C) GFD Dennou Club, 2004, 2005. All rights reserved. !--------------------------------------------------------------------- != Module DisturbEnv ! ! * Developer: SUGIYAMA Ko-ichiro, ODAKA Masatsugu ! * Version: $Id: distset.f90,v 1.1 2011-03-30 13:53:06 sugiyama Exp $ ! * Tag Name: $Name: arare4-20120911 $ ! * Change History: ! !== Overview ! !擾乱のデフォルト値を与えるためのルーチン. ! !== Error Handling ! !== Known Bugs ! !== Note ! !== Future Plans ! ! * 設定可能な擾乱のタイプを増やす ! * エラー処理に gf4f90io を利用 ! * randomfile はこのサブルーチンで定義 ! module DistSet ! !擾乱のデフォルト値を与えるためのルーチン. ! !モジュール読み込み use dc_message, only: MessageNotify use gridset, only:DimXMin, &! 配列の X 方向の下限 & DimXMax, &! 配列の X 方向の上限 & RegXMin, &! 配列の X 方向の下限 & RegXMax, &! 配列の X 方向の上限 & DimZMin, &! 配列の Z 方向の下限 & DimZMax, &! 配列の Z 方向の上限 & RegZMin, &! 配列の Z 方向の下限 & RegZMax, &! 配列の Z 方向の上限 & MarginX, &! 計算領域のマージン & MarginZ, &! 計算領域のマージン & DelZ, &! Z 方向の刻み幅 & SpcNum, &! 凝縮性化学種の数 & XMin, XMax, &! & ZMin, ZMax, &! & x_X, &! X 座標軸(スカラー格子点) & z_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, & & BoundaryXCyc_xz , &! & BoundaryZSym_xz , & & BoundaryXCyc_pz , &! & BoundaryZSym_pz use ECCM, only: ECCM_Wet !暗黙の型宣言禁止 implicit none !属性 private !変数定義 real(8), parameter, save :: Pi = 3.1415926535897932385d0 !円周率 !public public distset_pottemp public distset_exner public distset_mixrt public distset_velx contains subroutine distset_pottemp(cfgfile, xz_PotTemp) ! 暗黙の型宣言禁止 implicit none !変数定義 character(*), intent(in) :: cfgfile real(8), intent(out) :: xz_PotTemp(DimXMin:DimXMax, DimZMin:DimZMax) character(20) :: Type ="" real(8) :: DelMax = 0.0d0 !擾乱の最大値 real(8) :: XcRate = 0.0d0 !擾乱の中心位置(水平方向)の領域に対する割合 real(8) :: XrRate = 0.0d0 !擾乱の半径(水平方向)の領域に対する割合 real(8) :: YcRate = 0.0d0 !擾乱の中心位置(水平方向)の領域に対する割合 real(8) :: YrRate = 0.0d0 !擾乱の半径(水平方向)の領域に対する割合 real(8) :: ZcRate= 0.0d0 !擾乱の中心位置(鉛直方向)の領域に対する割合 real(8) :: ZrRate = 0.0d0 !擾乱の半径(鉛直方向)の領域に対する割合 real(8) :: Zpos = 0.0d0 !擾乱の Z 座標 [m] (Therma-Random 用) real(8) :: Xc !擾乱の中心位置(X方向) real(8) :: Yc !擾乱の中心位置(Y方向) real(8) :: Zc !擾乱の中心位置(鉛直方向) real(8) :: Xr !擾乱の半径(X方向) real(8) :: Yr !擾乱の半径(Y方向) real(8) :: Zr !擾乱の半径(鉛直方向) real(8) :: Random !ファイルから取得した乱数 real(8) :: beta(DimXMin:DimXMax, DimZMin:DimZMax) !擾乱の最大値に対する割合 real(8) :: RandomNum(DimXMin:DimXMax) real(8) :: RandomNum2(DimXMin:DimXMax) integer :: i, j, k, n ! NAMELIST ファイルの定義 NAMELIST /disturbenv_pottemp/ Type, & & DelMax, XrRate, XcRate, YrRate, YcRate, ZrRate, ZcRate, Zpos ! namelist の読み込み open (10, FILE=cfgfile) read(10, NML=disturbenv_pottemp) close(10) ! 位置の決定 ! XMax には MPI を含めた全領域が入っている Xr = minval( x_X, 1, x_X > (XMax - XMin) * XrRate ) Xc = minval( x_X, 1, x_X > (XMax - XMin) * XcRate ) ! Yr = minval( y_Y, 1, y_Y > (YMax - YMin) * YrRate ) ! Yc = minval( y_Y, 1, y_Y > (YMax - YMin) * YcRate ) Zr = minval( z_Z, 1, z_Z > (ZMax - ZMin) * ZrRate ) Zc = minval( z_Z, 1, z_Z > (ZMax - ZMin) * ZcRate ) ! 擾乱のタイプの設定 select case(Type) ! 指定された高度にランダムな擾乱を与える case ("Random") ! 乱数発生 do i = RegXMin, RegXMax call random_number(random) RandomNum(i) = random end do ! 擾乱が全体としてはゼロとなるように調整. 平均からの差にする. do i = RegXMin, RegXmax RandomNum2(i) = RandomNum(i) - sum( RandomNum(RegXMin+1:RegXMax) ) / real(RegXMax - RegXMin, 8) end do ! 指定された高度に擾乱を置く. DelMax は温度なので, 温位に変換する. do i = RegXMin, RegXMax xz_PotTemp(i, maxloc(z_Z, z_Z <= Zpos) - MarginZ - 1) = & & DelMax * RandomNum2(i) / xz_ExnerBasicZ(i, maxloc(z_Z, z_Z <= Zpos) - MarginZ - 1) end do ! 指定された座標を中心としたガウス分布の温位擾乱を与える. case ("Gauss") do k = DimZMin, DimZMax do i = DimXMin, DimXMax xz_PotTemp(i,k) = & & DelMax * dexp( - ( (x_X(i) - Xc) / Xr )**2.0d0 * 5.0d-1 & & - ( (z_Z(k) - Zc) / Zr )**2.0d0 * 5.0d-1 ) & & / xz_ExnerBasicZ(i,k) end do end do ! 指定された領域内に温位擾乱を与える (Klemp and Wilhelmson, 1978) case ("KW1978") do k = DimZMin, DimZMax do i = DimXMin, DimXMax beta(i,k) = & & ( & & ( ( x_X(i) - Xc ) / Xr ) ** 2.0d0 & & + ( ( z_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 ! Hueso and Sanchez-Lavega を模した設定 case ("HS2001") i = ( DimXMax - DimXMin - 10) / 2 k = minloc( z_Z, 1, z_Z > 2.5d4 ) - MarginZ n = int( 5.0d3 / DelZ ) xz_PotTemp(i-n:i,k-n:k) = DelMax ! Skamarock and Klemp (1989) の Cold-bubble 実験 case ("SK1989") do k = DimZMin, DimZMax do i = DimXMin, DimXMax beta(i,k) = & & sqrt( & & ( ( x_X(i) - Xc ) / Xr ) ** 2.0d0 & & + ( ( z_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 ! 境界条件 call BoundaryXCyc_xz( xz_PotTemp ) call BoundaryZSym_xz( xz_PotTemp ) end subroutine distset_pottemp subroutine distset_exner(cfgfile, xz_Exner) ! 暗黙の型宣言禁止 implicit none !変数定義 character(*), intent(in) :: cfgfile real(8), intent(out) :: xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax) character(20) :: Type ="" real(8) :: DelMax = 0.0d0 !擾乱の最大値 real(8) :: XcRate = 0.0d0 !擾乱の中心位置(水平方向)の領域に対する割合 real(8) :: XrRate = 0.0d0 !擾乱の半径(水平方向)の領域に対する割合 real(8) :: YcRate = 0.0d0 !擾乱の中心位置(水平方向)の領域に対する割合 real(8) :: YrRate = 0.0d0 !擾乱の半径(水平方向)の領域に対する割合 real(8) :: ZcRate= 0.0d0 !擾乱の中心位置(鉛直方向)の領域に対する割合 real(8) :: ZrRate = 0.0d0 !擾乱の半径(鉛直方向)の領域に対する割合 real(8) :: Xc !擾乱の中心位置(X方向) real(8) :: Yc !擾乱の中心位置(Y方向) real(8) :: Zc !擾乱の中心位置(鉛直方向) real(8) :: Xr !擾乱の半径(X方向) real(8) :: Yr !擾乱の半径(Y方向) real(8) :: Zr !擾乱の半径(鉛直方向) integer :: i, j, k ! NAMELIST ファイルの定義 NAMELIST /disturbenv_exner/ Type, & & DelMax, XrRate, XcRate, YrRate, YcRate, ZrRate, ZcRate ! namelist の読み込み open (10, FILE=cfgfile) read(10, NML=disturbenv_exner) close(10) ! 位置の決定 ! XMax には MPI を含めた全領域が入っている Xr = minval( x_X, 1, x_X > (XMax - XMin) * XrRate ) Xc = minval( x_X, 1, x_X > (XMax - XMin) * XcRate ) ! Yr = minval( y_Y, 1, y_Y > (YMax - YMin) * YrRate ) ! Yc = minval( y_Y, 1, y_Y > (YMax - YMin) * YcRate ) Zr = minval( z_Z, 1, z_Z > (ZMax - ZMin) * ZrRate ) Zc = minval( z_Z, 1, z_Z > (ZMax - ZMin) * ZcRate ) ! エクスナー関数の設定. select case (Type) ! 指定された場所に, ガウシアンなエクスナー関数の擾乱を与える. case ("Gauss") do k = DimZMin, DimZMax do i = DimXMin, DimXMax xz_Exner(i,k) = & & DelMax * dexp( - ( (x_X(i) - Xc) / Xr )**2.0d0 * 5.0d-1 & & - ( (z_Z(k) - Zc) / Zr )**2.0d0 * 5.0d-1 ) end do end do end select ! 境界条件 call BoundaryXCyc_xz( xz_Exner ) call BoundaryZSym_xz( xz_Exner ) end subroutine distset_exner subroutine distset_mixrt(cfgfile, xza_MixRt) ! 暗黙の型宣言禁止 implicit none !変数定義 character(*), intent(in) :: cfgfile real(8), intent(out) :: xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, 1:SpcNum) character(20) :: Type ="" real(8) :: DelMax = 0.0d0 !擾乱の最大値 real(8) :: XcRate = 0.0d0 !擾乱の中心位置(水平方向)の領域に対する割合 real(8) :: XrRate = 0.0d0 !擾乱の半径(水平方向)の領域に対する割合 real(8) :: YcRate = 0.0d0 !擾乱の中心位置(水平方向)の領域に対する割合 real(8) :: YrRate = 0.0d0 !擾乱の半径(水平方向)の領域に対する割合 real(8) :: ZcRate= 0.0d0 !擾乱の中心位置(鉛直方向)の領域に対する割合 real(8) :: ZrRate = 0.0d0 !擾乱の半径(鉛直方向)の領域に対する割合 real(8) :: XposMin = 0.0d0 ! real(8) :: XposMax = 0.0d0 ! real(8) :: ZposMin = 0.0d0 ! real(8) :: ZposMax = 0.0d0 ! real(8) :: Humidity = 0.0d0 !相対湿度 real(8) :: Xc !擾乱の中心位置(X方向) real(8) :: Yc !擾乱の中心位置(Y方向) real(8) :: Zc !擾乱の中心位置(鉛直方向) real(8) :: Xr !擾乱の半径(X方向) real(8) :: Yr !擾乱の半径(Y方向) real(8) :: Zr !擾乱の半径(鉛直方向) integer :: i, j, k, s real(8) :: MolFrIni(SpcNum) real(8) :: za_MolFr(DimZMin:DimZMax, SpcNum) ! NAMELIST ファイルの定義 NAMELIST /disturbenv_mixrt/ Type, & & DelMax, XrRate, XcRate, YrRate, YcRate, ZrRate, ZcRate, & & XposMin, XposMax, ZposMin, ZposMax, Humidity ! namelist の読み込み open (10, FILE=cfgfile) read(10, NML=disturbenv_mixrt) close(10) ! 位置の決定 ! XMax には MPI を含めた全領域が入っている Xr = minval( x_X, 1, x_X > (XMax - XMin) * XrRate ) Xc = minval( x_X, 1, x_X > (XMax - XMin) * XcRate ) ! Yr = minval( y_Y, 1, y_Y > (YMax - YMin) * YrRate ) ! Yc = minval( y_Y, 1, y_Y > (YMax - YMin) * YcRate ) Zr = minval( z_Z, 1, z_Z > (ZMax - ZMin) * ZrRate ) Zc = minval( z_Z, 1, z_Z > (ZMax - ZMin) * ZcRate ) ! 擾乱の設定 select case (Type) ! 指定された座標を中心としたガウス分布の混合比擾乱を与える. case ("Gauss") do k = DimZMin, DimZMax do i = DimXMin, DimXMax do s = 1, SpcNum xza_MixRt(i,:,s) = & & DelMax & & * dexp( - ( (x_X(i) - Xc) / Xr )**2.0d0 * 5.0d-1 & & - ( (z_Z(k) - Zc) / Zr )**2.0d0 * 5.0d-1 ) end do end do end do ! XposMin:XposMax,ZposMin:ZposMax で囲まれた領域の初期の湿度をゼロにするために ! 基本場と逆符号の水蒸気擾乱を与える case ("Dry") do s = 1, SpcNum do k = DimZMin,DimZMax do i = DimXMin,DimXMax if (z_Z(k) >= ZposMin .AND. z_Z(k) < ZposMax & & .AND. x_X(i) >= XposMin .AND. x_X(i) < XposMax) then xza_MixRt(i,k,s) = - xza_MixRtBasicZ(i,k,s) end if end do end do end do ! 与えられた相対湿度に設定する. case ("Humidity") MolFrIni = SpcWetMolFr(1:SpcNum) ! 水平一様なので, i=0 だけ計算. call ECCM_Wet( MolFrIni, Humidity, xz_TempBasicZ(1,:), xz_PressBasicZ(1,:), za_MolFr ) !気相のモル比を混合比に変換 do s = 1, SpcNum do k = RegZMin, RegZMax do i = RegXMin, RegXMax xza_MixRt(i,k,s) = za_MolFr(k,s) * MolWtWet(s) / MolWtDry - xza_MixRtBasicZ(i,k,s) end do end do end do end select ! 境界条件 call BoundaryXCyc_xza( xza_MixRt ) call BoundaryZSym_xza( xza_MixRt ) end subroutine distset_mixrt subroutine distset_VelX(cfgfile, pz_VelX) ! 暗黙の型宣言禁止 implicit none !変数定義 character(*), intent(in) :: cfgfile real(8), intent(out) :: pz_VelX(DimXMin:DimXMax, DimZMin:DimZMax) character(20) :: Type ="" real(8) :: Umax = 0.0d0 !擾乱の中心位置(水平方向)の領域に対する割合 real(8) :: ZposMin = 0.0d0 !擾乱の Z 座標 [m] (下端) real(8) :: ZposMax = 0.0d0 !擾乱の Z 座標 [m] (上端) integer :: i, j, k, b, u ! namelist の定義 NAMELIST /disturbenv_VelX/ Type, Umax, ZposMin, ZposMax ! 値を読み出す open (10, FILE=cfgfile) read(10, NML=disturbenv_velx) close(10) select case(Type) case ("Takemi2007") ! != 概要 !* 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 (z_Z(k) <= ZposMin .AND. ZposMin < z_Z(k+1)) then B = k end if if (z_Z(k) <= ZposMax .AND. ZposMax < z_Z(k+1)) then U = k end if end do ! 上記で決めた格子点BからUの間に鉛直シアーのある風を与える do k = RegZMin+1, U do i = DimXMin, DimXMax if (k < B) then pz_VelX(i,k) = Umax else pz_VelX(i,k) = Umax - (Umax / (z_Z(U) - z_Z(B)) * (z_Z(k) - z_Z(B))) end if end do end do end select !境界条件 call BoundaryXCyc_pz( pz_VelX ) call BoundaryZSym_pz( pz_VelX ) end subroutine Distset_VelX end module DistSet