!= Module Damping_3D ! ! Authors:: SUGIYAMA Ko-ichiro, ODAKA Masatsugu ! Version:: $Id: damping_3d.f90,v 1.1 2007/08/03 06:56:39 odakker Exp $ ! Tag Name:: $Name: arare4-20071012 $ ! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved. ! License:: See COPYRIGHT[link:../../COPYRIGHT] ! !== Overview ! !減衰率とその計算を行うためのパッケージ型モジュール ! * 音波減衰項の係数 ! * スポンジ層の設定(境界付近で波の反射を抑え吸収するための層) ! !== Error Handling ! !== Bugs ! !== Note ! ! * この関数は, 基本場が零な変数(速度, エクスナー関数の擾乱)に ! 適用することを想定している. ! * 各格子点に対する関数を定義する必要がある ! !== Future Plans ! ! module Damping_3d ! !減衰率とその計算を行うためのパッケージ型モジュール ! * 音波減衰項の係数 ! * スポンジ層の設定(境界付近で波の反射を抑え吸収するための層) ! !モジュール読み込み 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 方向の配列の上限 & RegXMin, &! x 方向の物理領域の下限 & RegXMax, &! x 方向の物理領域の上限 & RegYMin, &! y 方向の物理領域の下限 & RegYMax, &! y 方向の物理領域の上限 & RegZMin, &! z 方向の物理領域の下限 & RegZMax, &! z 方向の物理領域の上限 & x_X, &!X 座標軸(スカラー格子点) & y_Y, &!Y 座標軸(スカラー格子点) & z_Z, &!Z 座標軸(スカラー格子点) & p_X, &!X 座標軸(フラックス格子点) & q_Y, &!Y 座標軸(フラックス格子点) & r_Z, &!Z 座標軸(フラックス格子点) & x_dx, y_dy, z_dz, &! 格子間隔 & XMax, &!X 座標の最大値 & YMax, &!Y 座標の最大値 & ZMax !Z 座標の最大値 use timeset, only: DelTimeShort !短い時間ステップ !暗黙の型宣言禁止 implicit none !private 属性を指定 private !関数には public 属性を指定 public Damping_Init public DampSound_Init public DampSponge_Init public xyz_DampSponge_xyz public xyr_DampSponge_xyr public pyz_DampSponge_pyz public xqz_DampSponge_xqz public xyz_DampSponge public xyr_DampSponge public pyz_DampSponge public xqz_DampSponge public DampSound !変数定義 real(DP) :: DampSound = 0.0d0 !音波減衰項の減衰係数 real(DP) :: EFTime = 100.0d0 !スポンジ層の e-folding time real(DP) :: DampDepthH = 0.0d0 !スポンジ層の厚さ(水平方向) real(DP) :: DampDepthV = 0.0d0 !スポンジ層の厚さ(鉛直方向) real(DP), allocatable :: xyz_DampRateH(:,:,:) !xyz 格子減衰係数(水平方向) real(DP), allocatable :: xyz_DampRateV(:,:,:) !xyz 格子減衰係数(鉛直方向) real(DP), allocatable :: pyz_DampRateH(:,:,:) !pyz 格子減衰係数(水平方向) real(DP), allocatable :: pyz_DampRateV(:,:,:) !pyz 格子減衰係数(鉛直方向) real(DP), allocatable :: xqz_DampRateH(:,:,:) !xqz 格子減衰係数(水平方向) real(DP), allocatable :: xqz_DampRateV(:,:,:) !xqz 格子減衰係数(鉛直方向) real(DP), allocatable :: xyr_DampRateH(:,:,:) !xyr 格子減衰係数(水平方向) real(DP), allocatable :: xyr_DampRateV(:,:,:) !xyr 格子減衰係数(鉛直方向) integer :: unit ! 装置番号 !値を save する save DampSound, EFTime, DampDepthH, DampDepthV save xyz_DampRateH, xyz_DampRateV save pyz_DampRateH, pyz_DampRateV save xqz_DampRateH, xqz_DampRateV save xyr_DampRateH, xyr_DampRateV contains !!!------------------------------------------------------------------------!!! subroutine Damping_Init( cfgfile ) ! ! 音波減衰項とスポンジ層の減衰係数の初期化 ! !暗黙の型宣言禁止 implicit none !変数定義 character(*), intent(in) :: cfgfile real(DP) :: Alpha !音波減衰項の係数 real(DP) :: Time ! real(DP) :: DepthH !スポンジ層の厚さ(水平方向) real(DP) :: DepthV !スポンジ層の厚さ(鉛直方向) !NAMELIST から取得 NAMELIST /damping/ Alpha, Time, DepthH, DepthV call FileOpen(unit, file=cfgfile, mode='r') read(unit, NML=damping) close(unit) !初期化 call DampSound_Init( Alpha ) call DampSponge_Init( Time, DepthH, DepthV ) end subroutine Damping_Init !!!------------------------------------------------------------------------!!! subroutine DampSound_Init( damp ) ! ! 音波減衰項の減衰係数の初期化 ! !暗黙の型宣言禁止 implicit none !変数定義 real(DP), intent(in) :: damp real(DP) :: DelXMin, DelZMin DelXMin = minval( x_dx ) DelZMin = minval( z_dz ) !------------------------------------------------------------------- ! 音波減衰項の減衰率 Min(DelX, DelZ) ** 2.0 に比例 !------------------------------------------------------------------- DampSound = Damp * ( Min(DelXMin, DelZMin) ** 2.0d0 ) / DelTimeShort !----------------------------------------------------------------- ! 値の確認 !----------------------------------------------------------------- call MessageNotify( "M", & & "DampSound_init", "DampSound = %f", d=(/DampSound/) ) end subroutine DampSound_Init !!!------------------------------------------------------------------------!!! subroutine DampSponge_Init( Time, DepthH, DepthV ) ! ! スポンジ層の減衰係数を決める ! !暗黙の型宣言禁止 implicit none !入力変数 real(DP), intent(in) :: Time !スポンジ層の e-folding time real(DP), intent(in) :: DepthH !スポンジ層の厚さ(水平方向) real(DP), intent(in) :: DepthV !スポンジ層の厚さ(鉛直方向) real(DP), parameter :: Pi =3.1415926535897932385d0 !円周率 integer :: i, j, k !初期化 allocate( & & xyz_DampRateH(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax), & & xyz_DampRateV(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax), & & pyz_DampRateH(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax), & & pyz_DampRateV(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax), & & xqz_DampRateH(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax), & & xqz_DampRateV(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax), & & xyr_DampRateH(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax), & & xyr_DampRateV(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) ) xyz_DampRateH = 0.0d0 xyz_DampRateV = 0.0d0 pyz_DampRateH = 0.0d0 pyz_DampRateV = 0.0d0 xqz_DampRateH = 0.0d0 xqz_DampRateV = 0.0d0 xyr_DampRateH = 0.0d0 xyr_DampRateV = 0.0d0 !値の入力 EFTime = Time DampDepthH = DepthH DampDepthV = DepthV !----------------------------------------------------------------- ! スポンジ層の減衰率 !----------------------------------------------------------------- !水平方向の東側・西側境界 if ( DampDepthH < x_dx(RegXMin) ) then call MessageNotify( "W", & & "DampSponge_Init", "DampDepthH is too thin. DelX is %f", d=(/x_dx(RegXMin)/)) else if ( DampDepthH < x_dx(RegXMax) ) then call MessageNotify( "W", & & "DampSponge_Init", "DampDepthH is too thin. DelX is %f", d=(/x_dx(RegXMax)/)) else do i = DimXMin, DimXMax !スカラー格子点の西側境界 if ( x_X(i) < DampDepthH) then xyz_DampRateH(i,:,:) = ((1.0d0 - x_X(i) / DampDepthH) ** 3.0d0) / EFTime end if !フラックス格子点の西側境界 if ( p_X(i) < DampDepthH) then pyz_DampRateH(i,:,:) = ((1.0d0 - p_X(i) / DampDepthH) ** 3.0d0) / EFTime end if !スカラー格子点の東側境界 if ( x_X(i) > ( XMax - DampDepthH ) ) then xyz_DampRateH(i,:,:) = & & ((1.0d0 - (XMax - x_X(i)) / DampDepthH) ** 3.0d0) / EFTime end if !フラックス格子点の東側境界 if ( p_X(i) > ( XMax - DampDepthH ) ) then pyz_DampRateH(i,:,:) = & & ((1.0d0 - (XMax - p_X(i)) / DampDepthH) ** 3.0d0) / EFTime end if end do end if !sf と ss は X 方向に関しては同じ xyr_DampRateH = xyz_DampRateH !水平方向の南側・北側境界 if ( DampDepthH < y_dy(RegYMin) ) then call MessageNotify( "W", & & "DampSponge_Init", "DampDepthH is too thin. DelY is %f", d=(/x_dx(RegYMin)/)) else if ( DampDepthH < y_dy(RegYMax) ) then call MessageNotify( "W", & & "DampSponge_Init", "DampDepthH is too thin. DelY is %f", d=(/y_dy(RegYMax)/)) else do j = DimYMin, DimYMax !スカラー格子点の西側境界 if ( y_Y(j) < DampDepthH) then xyz_DampRateH(:,j,:) = ((1.0d0 - y_Y(j) / DampDepthH) ** 3.0d0) / EFTime end if !フラックス格子点の西側境界 if ( q_Y(j) < DampDepthH) then pyz_DampRateH(:,j,:) = ((1.0d0 - q_Y(j) / DampDepthH) ** 3.0d0) / EFTime end if !スカラー格子点の東側境界 if ( y_Y(j) > ( YMax - DampDepthH ) ) then xyz_DampRateH(:,j,:) = & & ((1.0d0 - (YMax - y_Y(j)) / DampDepthH) ** 3.0d0) / EFTime end if !フラックス格子点の東側境界 if ( q_Y(j) > ( YMax - DampDepthH ) ) then pyz_DampRateH(:,j,:) = & & ((1.0d0 - (YMax - q_Y(j)) / DampDepthH) ** 3.0d0) / EFTime end if end do end if !sf と ss は X 方向に関しては同じ xqz_DampRateH = xyz_DampRateH !鉛直方向の上部境界 if ( DampDepthV < z_dz(RegZMax) ) then call MessageNotify( "W", & & "DampSponge_Init", "DampDepthV is too thin. DelZ is %f", d=(/z_dz(RegZMax)/) ) else do k = DimZMin, DimZMax !スカラー格子点 if ( z_Z(k) >= ( ZMax - DampDepthV ) ) then xyz_DampRateV(:,:,k) = & & (1.0d0 - dcos(Pi * (z_Z(k) - ZMax + DampDepthV) / DampDepthV)) & & / EFTime end if !フラックス格子点 if ( r_Z(k) >= ( ZMax - DampDepthV ) ) then xyr_DampRateV(:,:,k) = & & (1.0d0 - dcos(Pi * (r_Z(k) - ZMax + DampDepthV)/ DampDepthV)) & & / EFTime end if end do end if !fs と ss は Z 方向に関しては同じ pyz_DampRateV = xyz_DampRateV !----------------------------------------------------------------- ! 値の確認 !----------------------------------------------------------------- call MessageNotify( "M", "DampSponge_Init", "EFTime = %f", d=(/EFTime/) ) call MessageNotify( "M", "DampSponge_Init", "DampDepthH = %f", d=(/DampDepthH/) ) call MessageNotify( "M", "DampSponge_Init", "DampDepthV = %f", d=(/DampDepthV/) ) end subroutine DampSponge_Init !!!------------------------------------------------------------------------!!! function xyz_DampSponge_xyz(xyz_Var) ! ! xyz 格子点に対するスポンジ層 ! !暗黙の型宣言禁止 implicit none !変数定義 real(DP), intent(in) :: xyz_Var(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) real(DP) :: xyz_DampSponge_xyz(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !スポンジ層によるダンピングを計算 xyz_DampSponge_xyz = ( xyz_DampRateH + xyz_DampRateV ) * xyz_Var end function xyz_DampSponge_xyz !!!------------------------------------------------------------------------!!! function xyz_DampSponge(xyz_VarA, xyz_VarB, DelTime) ! ! xyz 格子点に対するスポンジ層 ! !暗黙の型宣言禁止 implicit none !変数定義 real(DP), intent(in) :: xyz_VarA(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) real(DP), intent(in) :: xyz_VarB(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) real(DP), intent(in) :: DelTime real(DP) :: xyz_DampSponge(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !スポンジ層によるダンピングを計算 xyz_DampSponge = xyz_VarA - ( xyz_DampRateH + xyz_DampRateV ) * xyz_VarB * DelTime end function xyz_DampSponge !!!------------------------------------------------------------------------!!! function xyr_DampSponge_xyr(xyr_Var) ! ! xyr 格子点に対するスポンジ層 ! !暗黙の型宣言禁止 implicit none !変数定義 real(DP), intent(in) :: xyr_Var(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) real(DP) :: xyr_DampSponge_xyr(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !スポンジ層によるダンピングを計算 xyr_DampSponge_xyr = ( xyr_DampRateH + xyr_DampRateV )* xyr_Var end function xyr_DampSponge_xyr !!!------------------------------------------------------------------------!!! function xyr_DampSponge(xyr_VarA, xyr_VarB, DelTime) ! ! xyr 格子点に対するスポンジ層 ! !暗黙の型宣言禁止 implicit none !変数定義 real(DP), intent(in) :: xyr_VarA(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) real(DP), intent(in) :: xyr_VarB(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) real(DP), intent(in) :: DelTime real(DP) :: xyr_DampSponge(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !スポンジ層によるダンピングを計算 xyr_DampSponge = xyr_VarA - ( xyr_DampRateH + xyr_DampRateV )* xyr_VarB * DelTime end function xyr_DampSponge !!!------------------------------------------------------------------------!!! function pyz_DampSponge_pyz(pyz_Var) ! ! pyz 格子点に対するスポンジ層 ! !暗黙の型宣言禁止 implicit none !変数定義 real(DP), intent(in) :: pyz_Var(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) real(DP) :: pyz_DampSponge_pyz(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !スポンジ層によるダンピングを計算 pyz_DampSponge_pyz = ( pyz_DampRateH + pyz_DampRateV ) * pyz_Var end function pyz_DampSponge_pyz !!!------------------------------------------------------------------------!!! function pyz_DampSponge(pyz_VarA, pyz_VarB, DelTime) ! ! pyz 格子点に対するスポンジ層 ! !暗黙の型宣言禁止 implicit none !変数定義 real(DP), intent(in) :: pyz_VarA(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) real(DP), intent(in) :: pyz_VarB(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) real(DP), intent(in) :: DelTime real(DP) :: pyz_DampSponge(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !スポンジ層によるダンピングを計算 pyz_DampSponge = pyz_VarA - ( pyz_DampRateH + pyz_DampRateV ) * pyz_VarB * DelTime end function pyz_DampSponge !!!------------------------------------------------------------------------!!! function xqz_DampSponge_xqz(xqz_Var) ! ! xqz 格子点に対するスポンジ層 ! !暗黙の型宣言禁止 implicit none !変数定義 real(DP), intent(in) :: xqz_Var(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) real(DP) :: xqz_DampSponge_xqz(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !スポンジ層によるダンピングを計算 xqz_DampSponge_xqz = ( xqz_DampRateH + xqz_DampRateV ) * xqz_Var end function xqz_DampSponge_xqz !!!------------------------------------------------------------------------!!! function xqz_DampSponge(xqz_VarA, xqz_VarB, DelTime) ! ! xqz 格子点に対するスポンジ層 ! !暗黙の型宣言禁止 implicit none !変数定義 real(DP), intent(in) :: xqz_VarA(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) real(DP), intent(in) :: xqz_VarB(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) real(DP), intent(in) :: DelTime real(DP) :: xqz_DampSponge(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !スポンジ層によるダンピングを計算 xqz_DampSponge = xqz_VarA - ( xqz_DampRateH + xqz_DampRateV ) * xqz_VarB * DelTime end function xqz_DampSponge end module Damping_3d