| Class | Damping |
| In: |
../src/utils/damping.f90
|
減衰率とその計算を行うためのパッケージ型モジュール
* 音波減衰項の係数 * スポンジ層の設定(境界付近で波の反射を抑え吸収するための層)
| Subroutine : |
音波減衰項とスポンジ層の減衰係数の初期化
This procedure input/output NAMELIST#damping_nml .
subroutine Damping_Init
!
! 音波減衰項とスポンジ層の減衰係数の初期化
!
use dc_iounit, only: FileOpen
use dc_message, only: MessageNotify
use gtool_historyauto, only: HistoryAutoAddVariable
use mpi_wrapper, only: myrank
use gridset, only: imin, imax, jmin, jmax, kmin, kmax, nx, ny, nz ! z 方向の物理領域の上限
use axesset, only: x_X, y_Y, z_Z, p_X, q_Y, r_Z, x_dx, y_dy, z_dz, XMax, YMax, ZMax !Z 座標の最大値
use namelist_util, only: namelist_filename
!暗黙の型宣言禁止
implicit none
!変数定義
real(DP) :: Time !
real(DP) :: DepthH !スポンジ層の厚さ(水平方向)
real(DP) :: DepthV !スポンジ層の厚さ(鉛直方向)
real(DP), parameter :: Pi =3.1415926535897932385d0 !円周率
integer :: unit ! 装置番号
integer :: i, j, k
!NAMELIST から取得
NAMELIST /damping_nml/ Time, DepthH, DepthV
call FileOpen(unit, file=namelist_filename, mode='r')
read(unit, NML=damping_nml)
close(unit)
!初期化
allocate( xyz_Gamma(imin:imax,jmin:jmax,kmin:kmax), pyz_Gamma(imin:imax,jmin:jmax,kmin:kmax), xqz_Gamma(imin:imax,jmin:jmax,kmin:kmax), xyr_Gamma(imin:imax,jmin:jmax,kmin:kmax) )
xyz_Gamma = 0.0d0
pyz_Gamma = 0.0d0
xqz_Gamma = 0.0d0
xyr_Gamma = 0.0d0
!値の入力
EFTime = Time
DampDepthH = DepthH
DampDepthV = DepthV
!-----------------------------------------------------------------
! スポンジ層の減衰率
!
!水平方向の東側・西側境界
!
if ( DampDepthH < x_dx(1) ) then
if (myrank == 0) call MessageNotify( "W", "Damping_init", "DampDepthH is too thin. DelX is %f", d=(/x_dx(1)/))
else if ( DampDepthH < x_dx(nx) ) then
if (myrank == 0) call MessageNotify( "W", "Damping_init", "DampDepthH is too thin. DelX is %f", d=(/x_dx(nx)/))
else
do i = imin, imax
!スカラー格子点の西側境界
if ( x_X(i) < DampDepthH) then
xyz_Gamma(i,:,:) = xyz_Gamma(i,:,:) + ((1.0d0 - x_X(i) / DampDepthH) ** 3.0d0) / EFTime
end if
!フラックス格子点の西側境界
if ( p_X(i) < DampDepthH) then
pyz_Gamma(i,:,:) = pyz_Gamma(i,:,:) + ((1.0d0 - p_X(i) / DampDepthH) ** 3.0d0) / EFTime
end if
!スカラー格子点の東側境界
if ( x_X(i) > ( XMax - DampDepthH ) ) then
xyz_Gamma(i,:,:) = xyz_Gamma(i,:,:) + ((1.0d0 - (XMax - x_X(i)) / DampDepthH) ** 3.0d0) / EFTime
end if
!フラックス格子点の東側境界
if ( p_X(i) > ( XMax - DampDepthH ) ) then
pyz_Gamma(i,:,:) = pyz_Gamma(i,:,:) + ((1.0d0 - (XMax - p_X(i)) / DampDepthH) ** 3.0d0) / EFTime
end if
end do
end if
! x 方向には同じ
!
xyr_Gamma = xyz_Gamma
xqz_Gamma = xyz_Gamma
!水平方向の南側・北側境界
!
if ( DampDepthH < y_dy(1) ) then
if (myrank == 0) call MessageNotify( "W", "Damping_init", "DampDepthH is too thin. DelY is %f", d=(/x_dx(1)/))
else if ( DampDepthH < y_dy(ny) ) then
if (myrank == 0) call MessageNotify( "W", "Damping_init", "DampDepthH is too thin. DelY is %f", d=(/y_dy(ny)/))
else
do j = jmin, jmax
!スカラー格子点の西側境界
if ( y_Y(j) < DampDepthH) then
xyz_Gamma(:,j,:) = xyz_Gamma(:,j,:) + ((1.0d0 - y_Y(j) / DampDepthH) ** 3.0d0) / EFTime
end if
!フラックス格子点の西側境界
if ( q_Y(j) < DampDepthH) then
xqz_Gamma(:,j,:) = xqz_Gamma(:,j,:) + ((1.0d0 - q_Y(j) / DampDepthH) ** 3.0d0) / EFTime
end if
!スカラー格子点の東側境界
if ( y_Y(j) > ( YMax - DampDepthH ) ) then
xyz_Gamma(:,j,:) = xyz_Gamma(:,j,:) + ((1.0d0 - (YMax - y_Y(j)) / DampDepthH) ** 3.0d0) / EFTime
end if
!フラックス格子点の東側境界
if ( q_Y(j) > ( YMax - DampDepthH ) ) then
xqz_Gamma(:,j,:) = xqz_Gamma(:,j,:) + ((1.0d0 - (YMax - q_Y(j)) / DampDepthH) ** 3.0d0) / EFTime
end if
end do
end if
! y 方向には同じ
!
pyz_Gamma = xyz_Gamma
xyr_Gamma = xyz_Gamma
!鉛直方向の上部境界
!
if ( DampDepthV < z_dz(nz) ) then
if (myrank == 0) call MessageNotify( "W", "Damping_init", "DampDepthV is too thin. DelZ is %f", d=(/z_dz(nz)/) )
else
do k = kmin, kmax
!スカラー格子点
if ( z_Z(k) >= ( ZMax - DampDepthV ) ) then
xyz_Gamma(:,:,k) = xyz_Gamma(:,:,k) + (1.0d0 - dcos(Pi * (z_Z(k) - ZMax + DampDepthV) / DampDepthV)) / EFTime
end if
!フラックス格子点
if ( r_Z(k) >= ( ZMax - DampDepthV ) ) then
xyr_Gamma(:,:,k) = xyr_Gamma(:,:,k) + (1.0d0 - dcos(Pi * (r_Z(k) - ZMax + DampDepthV)/ DampDepthV)) / EFTime
end if
end do
end if
! z 方向には同じ
!
pyz_Gamma = xyz_Gamma
xqz_Gamma = xyz_Gamma
!-----------------------------------------------------------------
! 値の確認
!
if (myrank == 0) then
call MessageNotify( "M", "Damping_init", "EFTime = %f", d=(/EFTime/) )
call MessageNotify( "M", "Damping_init", "DampDepthH = %f", d=(/DampDepthH/) )
call MessageNotify( "M", "Damping_init", "DampDepthV = %f", d=(/DampDepthV/) )
end if
!-----------------------------------------------------------------
! 出力
!
call HistoryAutoAddVariable( varname='PTempDamp', dims=(/'x','y','z','t'/), longname='Damping term of potential temperature', units='K.s-1', xtype='float')
call HistoryAutoAddVariable( varname='ExnerDamp', dims=(/'x','y','z','t'/), longname='Damping term of exner function', units='s-1', xtype='float')
call HistoryAutoAddVariable( varname='VelXDamp', dims=(/'x','y','z','t'/), longname='Damping term of VelX', units='m.s-1', xtype='float')
call HistoryAutoAddVariable( varname='VelYDamp', dims=(/'x','y','z','t'/), longname='Damping term of VelY', units='m.s-1', xtype='float')
call HistoryAutoAddVariable( varname='VelZDamp', dims=(/'x','y','z','t'/), longname='Damping term of VelZ', units='m.s-1', xtype='float')
end subroutine Damping_Init
| Subroutine : | |
| pyz_VelXBl(imin:imax, jmin:jmax, kmin:kmax) : | real(8), intent(in) |
| xqz_VelYBl(imin:imax, jmin:jmax, kmin:kmax) : | real(8), intent(in) |
| xyr_VelZBl(imin:imax, jmin:jmax, kmin:kmax) : | real(8), intent(in) |
| xyz_PTempBl(imin:imax, jmin:jmax, kmin:kmax) : | real(8), intent(in) |
| xyz_ExnerBl(imin:imax, jmin:jmax, kmin:kmax) : | real(8), intent(in) |
| pyz_DVelXDt(imin:imax, jmin:jmax, kmin:kmax) : | real(8), intent(inout) |
| xqz_DVelYDt(imin:imax, jmin:jmax, kmin:kmax) : | real(8), intent(inout) |
| xyr_DVelZDt(imin:imax, jmin:jmax, kmin:kmax) : | real(8), intent(inout) |
| xyz_DPTempDt(imin:imax, jmin:jmax, kmin:kmax) : | real(8), intent(inout) |
| xyz_DExnerDt(imin:imax, jmin:jmax, kmin:kmax) : | real(8), intent(inout) |
(inout)
subroutine SpongeLayer_forcing( pyz_VelXBl, xqz_VelYBl, xyr_VelZBl, xyz_PTempBl, xyz_ExnerBl, pyz_DVelXDt, xqz_DVelYDt, xyr_DVelZDt, xyz_DPTempDt, xyz_DExnerDt ) !(inout)
use gtool_historyauto, only: HistoryAutoPut
use timeset, only: TimeN ! 現在の時刻
use gridset, only: imin, imax, jmin, jmax, kmin, kmax, nx, ny, nz ! z 方向の物理領域の上限
!暗黙の型宣言禁止
implicit none
real(8), intent(in) :: pyz_VelXBl(imin:imax, jmin:jmax, kmin:kmax)
real(8), intent(in) :: xqz_VelYBl(imin:imax, jmin:jmax, kmin:kmax)
real(8), intent(in) :: xyr_VelZBl(imin:imax, jmin:jmax, kmin:kmax)
real(8), intent(in) :: xyz_PTempBl(imin:imax, jmin:jmax, kmin:kmax)
real(8), intent(in) :: xyz_ExnerBl(imin:imax, jmin:jmax, kmin:kmax)
real(8), intent(inout) :: pyz_DVelXDt(imin:imax, jmin:jmax, kmin:kmax)
real(8), intent(inout) :: xqz_DVelYDt(imin:imax, jmin:jmax, kmin:kmax)
real(8), intent(inout) :: xyr_DVelZDt(imin:imax, jmin:jmax, kmin:kmax)
real(8), intent(inout) :: xyz_DPTempDt(imin:imax, jmin:jmax, kmin:kmax)
real(8), intent(inout) :: xyz_DExnerDt(imin:imax, jmin:jmax, kmin:kmax)
real(8) :: pyz_DampVelX(imin:imax, jmin:jmax, kmin:kmax)
real(8) :: xqz_DampVelY(imin:imax, jmin:jmax, kmin:kmax)
real(8) :: xyr_DampVelZ(imin:imax, jmin:jmax, kmin:kmax)
real(8) :: xyz_DampPTemp(imin:imax, jmin:jmax, kmin:kmax)
real(8) :: xyz_DampExner(imin:imax, jmin:jmax, kmin:kmax)
!--------------------------------------------------------
! 擾乱場の値をゼロに戻すような強制を与える
!
pyz_DampVelX = - pyz_Gamma * pyz_VelXBl
xqz_DampVelY = - xqz_Gamma * xqz_VelYBl
xyr_DampVelZ = - xyr_Gamma * xyr_VelZBl
xyz_DampPTemp = - xyz_Gamma * xyz_PTempBl
xyz_DampExner = - xyz_Gamma * xyz_ExnerBl
pyz_DVelXDt = pyz_DVelXDt + pyz_DampVelX
xqz_DVelYDt = xqz_DVelYDt + xqz_DampVelY
xyr_DVelZDt = xyr_DVelZDt + xyr_DampVelZ
xyz_DPTempDt = xyz_DPTempDt + xyz_DampPTemp
xyz_DExnerDt = xyz_DExnerDt + xyz_DampExner
call HistoryAutoPut(TimeN, 'VelXDamp', pyz_DampVelX(1:nx,1:ny,1:nz))
call HistoryAutoPut(TimeN, 'VelYDamp', xqz_DampVelY(1:nx,1:ny,1:nz))
call HistoryAutoPut(TimeN, 'VelZDamp', xyr_DampVelZ(1:nx,1:ny,1:nz))
call HistoryAutoPut(TimeN, 'PTempDamp', xyz_DampPTemp(1:nx,1:ny,1:nz))
call HistoryAutoPut(TimeN, 'ExnerDamp', xyz_DampExner(1:nx,1:ny,1:nz))
end subroutine SpongeLayer_forcing