!--------------------------------------------------------------------- ! Copyright (C) GFD Dennou Club, 2004, 2005. All rights reserved. !--------------------------------------------------------------------- != Subroutine DisturbEnv_3d ! ! * Developer: SUGIYAMA Ko-ichiro, ODAKA Masatsugu ! * Version: $Id: disturbenv_3d.f90,v 1.6 2008-06-24 20:16:40 odakker2 Exp $ ! * Tag Name: $Name: arare4-20120911 $ ! * Change History: ! !== Overview ! !擾乱のデフォルト値を与えるためのルーチン. ! !== Error Handling ! !== Known Bugs ! !== Note ! !== Future Plans ! ! * 設定可能な擾乱のタイプを増やす ! * エラー処理に gf4f90io を利用 ! subroutine DisturbEnv_3d( & & cfgfile, xyz_PotTemp, xyz_Exner, & & pyz_VelX, xqz_VelY, xyr_VelZ, & & xyza_MixRt, xyz_Km, xyz_Kh & & ) ! !擾乱のデフォルト値を与えるためのルーチン. ! !モジュール読み込み use dc_types, only : DP use dc_iounit, only: FileOpen use dc_message, only: MessageNotify use gridset_3d, only:DimXMin, &! 配列の X 方向の下限 & DimXMax, &! 配列の X 方向の上限 & DimYMin, &! 配列の Y 方向の下限 & DimYMax, &! 配列の Y 方向の上限 & DimZMin, &! 配列の Z 方向の下限 & DimZMax, &! 配列の Z 方向の上限 & Zmargin, &! & SpcNum, &! 凝縮性化学種の数 & XMin, XMax, &! & YMin, YMax, &! & ZMin, ZMax, &! & x_X, &! X 座標軸(スカラー格子点) & y_Y, &! X 座標軸(スカラー格子点) & z_Z, &! Z 座標軸(スカラー格子点) & xyz_X, xyz_Y, xyz_Z, & & x_dx, y_dy, z_dz use fileset_3d, only:RandomFile ! 乱数ファイル use basicset_3d,only: & & SpcWetMolFr, &!凝縮成分の初期モル比 & MolWtWet, &!凝縮成分の分子量 & MolWtDry, &!乾燥成分の分子量 & xyz_TempBasicZ, &! 基本場の温度 & xyz_PressBasicZ, &! 基本場の圧力 & xyz_ExnerBasicZ ! 無次元圧力 ! & xyza_MixRtBasicZ ! 基本場の混合比 use xyz_module, only: BoundaryXCyc_xyz, & & BoundaryYCyc_xyz, & & BoundaryZSym_xyz ! use ECCM_3d, only: ECCM_MolFr !暗黙の型宣言禁止 implicit none !変数定義 character(*), intent(in) :: cfgfile real(DP), intent(out) :: pyz_VelX(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !水平風速の擾乱成分 real(DP), intent(out) :: xqz_VelY(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !水平風速の擾乱成分 real(DP), intent(out) :: xyr_VelZ(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !鉛直風速の擾乱成分 real(DP), intent(out) :: xyz_Exner(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !エクスナー関数の擾乱成分 real(DP), intent(out) :: xyz_PotTemp(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !温位の擾乱成分 real(DP), intent(out) :: xyza_MixRt(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax,SpcNum) !凝縮成分の混合比(擾乱成分) real(DP), intent(out) :: xyz_Km(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !運動量に対する拡散係数 real(DP), intent(out) :: xyz_Kh(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !熱に対する拡散係数 real(DP), parameter :: Pi = 3.1415926535897932385d0 !円周率 real(DP) :: Humidity !相対湿度 real(DP) :: XcRate !擾乱の中心位置(水平方向)の領域に対する割合 real(DP) :: XrRate !擾乱の半径(水平方向)の領域に対する割合 real(DP) :: YcRate !擾乱の中心位置(水平方向)の領域に対する割合 real(DP) :: YrRate !擾乱の半径(水平方向)の領域に対する割合 real(DP) :: ZcRate !擾乱の中心位置(鉛直方向)の領域に対する割合 real(DP) :: ZrRate !擾乱の半径(鉛直方向)の領域に対する割合 real(DP) :: Xc !擾乱の中心位置(水平方向) real(DP) :: Xr !擾乱の半径(水平方向) real(DP) :: Yc !擾乱の中心位置(水平方向) real(DP) :: Yr !擾乱の半径(水平方向) real(DP) :: Zc !擾乱の中心位置(鉛直方向) real(DP) :: Zr !擾乱の半径(鉛直方向) real(DP) :: Xpos !擾乱の X 座標 [m] (Therma-Random 用) real(DP) :: Ypos !擾乱の Y 座標 [m] (Therma-Random 用) real(DP) :: Zpos !擾乱の Z 座標 [m] (Therma-Random 用) real(DP) :: beta(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !擾乱の最大値に対する割合 real(DP) :: DelMax !温位擾乱の最大値 real(DP) :: Random !ファイルから取得した乱数 real(DP) :: RandomNum(DimXMin:DimXMax,DimYMin:DimYMax) real(DP) :: RandomNum2(DimXMin:DimXMax,DimYMin:DimYMax) real(DP) :: RandomMean character(20) :: Type !温位擾乱のタイプ ! real(DP) :: xyza_MolFr(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum) !モル比 integer :: i, j, k, s, n ! DO ループ変数 integer :: unit ! 装置番号 !------------------------------------------------------------- ! 初期化 !------------------------------------------------------------- !NAMELIST ファイルの読み込み NAMELIST /disturbset/ & & Type, DelMax, XrRate, XcRate, YrRate, YcRate, ZrRate, ZcRate, & & Humidity, Xpos, Ypos, Zpos call FileOpen(unit, file=cfgfile, mode='r') read(unit, NML=disturbset) close(unit) !初期化 pyz_VelX = 0.0d0 xqz_VelY = 0.0d0 xyr_VelZ = 0.0d0 xyz_Exner = 0.0d0 xyz_PotTemp = 0.0d0 xyza_MixRt = 0.0d0 xyz_Km = 0.0d0 xyz_Kh = 0.0d0 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 ("Thermal-KW1978") ! 指定された領域内に温位擾乱を与える (Klemp and Wilhelmson, 1978) beta = & & ( & & ( ( xyz_X - Xc ) / Xr ) ** 2.0d0 & & + ( ( xyz_Y - Yc ) / Yr ) ** 2.0d0 & & + ( ( xyz_Z - Zc ) / Zr ) ** 2.0d0 & & ) ** 5.0d-1 where ( beta < 1.0d0 ) xyz_PotTemp = DelMax * ( dcos( Pi * 5.0d-1 * beta ) ** 2.0d0 ) & & / xyz_ExnerBasicZ end where case ("Thermal-Gauss") ! 指定された座標を中心としたガウス分布の温位擾乱を与える. xyz_PotTemp = & & DelMax & & * dexp( - ( (xyz_X - Xc) / Xr )**2.0d0 * 5.0d-1 & & - ( (xyz_Y - Yc) / Yr )**2.0d0 * 5.0d-1 & & - ( (xyz_Z - Zc) / Yr )**2.0d0 * 5.0d-1 ) & & / xyz_ExnerBasicZ where ( sign(1.0d0, DelMax) * xyz_PotTemp < DelMax * 1.0d-4 ) xyz_PotTemp = 0.0d0 end where case ("Thermal-GaussXZ") ! 指定された座標を中心としたガウス分布の温位擾乱を与える. xyz_PotTemp = & & DelMax & & * dexp( - ( (xyz_X - Xc) / Xr )**2.0d0 * 5.0d-1 & & - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) & & / xyz_ExnerBasicZ where ( sign(1.0d0, DelMax) * xyz_PotTemp < DelMax * 1.0d-4 ) xyz_PotTemp = 0.0d0 end where case ("Thermal-GaussYZ") ! 指定された座標を中心としたガウス分布の温位擾乱を与える. xyz_PotTemp = & & DelMax & & * dexp( - ( (xyz_Y - Yc) / Yr )**2.0d0 * 5.0d-1 & & - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) & & / xyz_ExnerBasicZ where ( sign(1.0d0, DelMax) * xyz_PotTemp < DelMax * 1.0d-4 ) xyz_PotTemp = 0.0d0 end where case ("Adv-XZ-X") ! X 方向移流のテスト ! 指定された座標を中心としたガウス分布の温位擾乱を与える. xyz_PotTemp = & & DelMax & & * dexp( - ( (xyz_X - Xc) / Xr )**2.0d0 * 5.0d-1 & & - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) where ( sign(1.0d0, DelMax) * xyz_PotTemp < DelMax * 1.0d-4 ) xyz_PotTemp = 0.0d0 end where where ( xyz_PotTemp > DelMax ) xyz_PotTemp = DelMax end where pyz_VelX = 20.0d0 case ("Adv-XZ-Z") ! Z 方向移流のテスト ! 指定された座標を中心としたガウス分布の温位擾乱を与える. xyz_PotTemp = & & DelMax & & * dexp( - ( (xyz_X - Xc) / Xr )**2.0d0 * 5.0d-1 & & - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) where ( sign(1.0d0, DelMax) * xyz_PotTemp < DelMax * 1.0d-4 ) xyz_PotTemp = 0.0d0 end where where ( xyz_PotTemp > DelMax ) xyz_PotTemp = DelMax end where xyr_VelZ = 20.0d0 case ("Adv-YZ-Y") ! Y 方向移流のテスト ! 指定された座標を中心としたガウス分布の温位擾乱を与える. xyz_PotTemp = & & DelMax & & * dexp( - ( (xyz_Y - Yc) / Yr )**2.0d0 * 5.0d-1 & & - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) where ( sign(1.0d0, DelMax) * xyz_PotTemp < DelMax * 1.0d-4 ) xyz_PotTemp = 0.0d0 end where where ( xyz_PotTemp > DelMax ) xyz_PotTemp = DelMax end where xqz_VelY = 20.0d0 case ("Adv-YZ-Z") ! Z 方向移流のテスト ! 指定された座標を中心としたガウス分布の温位擾乱を与える. xyz_PotTemp = & & DelMax & & * dexp( - ( (xyz_Y - Yc) / Yr )**2.0d0 * 5.0d-1 & & - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) where ( sign(1.0d0, DelMax) * xyz_PotTemp < DelMax * 1.0d-4 ) xyz_PotTemp = 0.0d0 end where where ( xyz_PotTemp > DelMax ) xyz_PotTemp = DelMax end where xyr_VelZ = 20.0d0 case ("Adv-XZ-XZ") ! XZ 方向移流のテスト ! 指定された座標を中心としたガウス分布の温位擾乱を与える. xyz_PotTemp = & & DelMax & & * dexp( - ( (xyz_X - Xc) / Xr )**2.0d0 * 5.0d-1 & & - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) where ( sign(1.0d0, DelMax) * xyz_PotTemp < DelMax * 1.0d-4 ) xyz_PotTemp = 0.0d0 end where where ( xyz_PotTemp > DelMax ) xyz_PotTemp = DelMax end where pyz_VelX = 20.0d0 xyr_VelZ = 20.0d0 case ("Adv-YZ-YZ") ! YZ 方向移流のテスト ! 指定された座標を中心としたガウス分布の温位擾乱を与える. xyz_PotTemp = & & DelMax & & * dexp( - ( (xyz_Y - Yc) / Yr )**2.0d0 * 5.0d-1 & & - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) where ( sign(1.0d0, DelMax) * xyz_PotTemp < DelMax * 1.0d-4 ) xyz_PotTemp = 0.0d0 end where where ( xyz_PotTemp > DelMax ) xyz_PotTemp = DelMax end where xqz_VelY = 20.0d0 xyr_VelZ = 20.0d0 case ("Exner-Gauss") ! 指定された場所に, ガウシアンな圧力の擾乱を与える. xyz_Exner = & & DelMax & & * dexp( - ( (xyz_X - Xc) / Xr )**2.0d0 * 5.0d-1 ) & & * dexp( - ( (xyz_Y - Yc) / Yr )**2.0d0 * 5.0d-1 ) & & * dexp( - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) where ( xyz_Exner < DelMax * 1.0d-4) xyz_Exner = 0.0d0 end where case ("Exner-GaussX") ! 指定された場所に, ガウシアンな圧力の擾乱を与える. ! (X 方向一様) xyz_Exner = & & DelMax & & * dexp( - ( (xyz_Y - Yc) / Yr )**2.0d0 * 5.0d-1 ) & & * dexp( - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) where ( xyz_Exner < DelMax * 1.0d-4) xyz_Exner = 0.0d0 end where case ("Exner-GaussY") ! 指定された場所に, ガウシアンな圧力の擾乱を与える. ! (Y 方向一様) xyz_Exner = & & DelMax & & * dexp( - ( (xyz_X - Xc) / Xr )**2.0d0 * 5.0d-1 ) & & * dexp( - ( (xyz_Z - Zc) / Zr )**2.0d0 * 5.0d-1 ) where ( xyz_Exner < DelMax * 1.0d-4) xyz_Exner = 0.0d0 end where case ("Exner-GaussZ") ! 指定された場所に, ガウシアンな圧力の擾乱を与える. ! (Z 方向一様) xyz_Exner = & & DelMax & & * dexp( - ( (xyz_X - Xc) / Xr )**2.0d0 * 5.0d-1 ) & & * dexp( - ( (xyz_Y - Yc) / Yr )**2.0d0 * 5.0d-1 ) where ( xyz_Exner < DelMax * 1.0d-4) xyz_Exner = 0.0d0 end where case ("Thermal-Random") ! 下層にランダムな擾乱を与える call FileOpen(unit, file=RandomFile, mode='r') do j = DimYMin, DimYMax do i = DimXMin, DimXMax read(unit,*) random RandomNum(i,j) = random end do end do close(unit) ! 乱数の平均値を計算 RandomMean = 0.0d0 do j = DimYMin, DimYMax do i = DimXMin, DimXMax RandomMean = RandomMean + & & (RandomNum(i,j)*x_dx(i)*y_dy(j))/((XMax-Xmin)*(YMax-YMin)) end do end do do j = DimYMin, DimYMax do i = DimXMin, DimXMax !擾乱が全体としてはゼロとなるように調整 RandomNum2(i,j) = RandomNum(i,j) - RandomMean xyz_PotTemp(i,j,maxloc(z_Z, z_Z <= Zpos) - Zmargin ) = & & DelMax * RandomNum2(i,j) / xyz_ExnerBasicZ(i,j,maxloc(z_Z, z_Z <= Zpos) - Zmargin - 1) end do end do case ("HS2001") ! Hueso and Sanchez-Lavega を模した設定 i = ( DimXMax - DimXMin - 10) / 2 j = ( DimYMax - DimYMin - 10) / 2 k = minloc( z_Z, 1, z_Z > 2.5d4 ) - Zmargin n = int( 5.0d3 / z_dz(1) ) xyz_PotTemp(i-n:i,j-n:j,k-n:k) = DelMax case ("SK1989") ! Skamarock and Klemp (1989) の Cold-bubble 実験 xyz_PotTemp = 0.0d0 beta = & & sqrt( & & ( ( xyz_X - Xc ) / Xr ) ** 2.0d0 & & + ( ( xyz_Y - Yc ) / Yr ) ** 2.0d0 & & + ( ( xyz_Z - Zc ) / Zr ) ** 2.0d0 & & ) where ( beta < 1.0d0 ) xyz_PotTemp = 0.5d0*DelMax*(cos(pi*beta) + 1.0d0) end where end select !------------------------------------------------------------- ! 蒸気圧の設定. !------------------------------------------------------------- ! if (Humidity /= 0.0d0) then ! do i = DimXMin, DimXMax ! do j = DimYMin, DimYMax ! call ECCM_MolFr( SpcWetMolFr(1:SpcNum), Humidity, & ! & xyz_TempBasicZ(i,j,:), & ! & xyz_PressBasicZ(i,j,:), xyza_MolFr(i,j,:,:) ) ! end do ! end do !気相のモル比を混合比に変換 ! do s = 1, SpcNum ! xyza_MixRt(:,:,:,s) = & ! & xyza_MolFr(:,:,:,s) & ! & * MolWtWet(s) / MolWtDry - xyza_MixRtBasicZ(:,:,:,s) ! end do ! end if !値が小さくなりすぎないように最低値を与える ! where (xza_MixRt <= 1.0d-20 ) ! xza_MixRt = 1.0d-20 ! end where !境界条件 ! do s =1, SpcNum ! call BoundaryXCyc_xyz( xyza_MixRt(:,:,:,s) ) ! call BoundaryYCyc_xyz( xyza_MixRt(:,:,:,s) ) ! call BoundaryZSym_xyz( xyza_MixRt(:,:,:,s) ) ! end do end subroutine DisturbEnv_3d