| 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, ZMin, ZMax !Z 座標の最大値
use namelist_util, only: namelist_filename
!暗黙の型宣言禁止
implicit none
!変数定義
real(DP), parameter :: Pi =3.1415926535897932385d0 !円周率
integer :: unit ! 装置番号
integer :: i, j, k
!NAMELIST から取得
NAMELIST /damping_nml/ EFTime, DepthH, DepthV, DepthVb
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
!-----------------------------------------------------------------
! 厚さのチェック
!
if ( DepthH > 0.0d0 ) then
if ( DepthH < x_dx(1) ) then
if (myrank == 0) call MessageNotify( "E", "Damping_init", "DepthH is too thin. DelX is %f", d=(/x_dx(1)/))
else if ( DepthH < x_dx(nx) ) then
if (myrank == 0) call MessageNotify( "E", "Damping_init", "DepthH is too thin. DelX is %f", d=(/x_dx(nx)/))
end if
if ( DepthH < y_dy(1) ) then
if (myrank == 0) call MessageNotify( "E", "Damping_init", "DepthH is too thin. DelY is %f", d=(/x_dx(1)/))
else if ( DepthH < y_dy(ny) ) then
if (myrank == 0) call MessageNotify( "E", "Damping_init", "DepthH is too thin. DelY is %f", d=(/y_dy(ny)/))
end if
end if
if ( DepthV > 0.0d0 ) then
if ( DepthV < z_dz(1) ) then
if (myrank == 0) call MessageNotify( "E", "Damping_init", "DepthV is too thin. DelZ is %f", d=(/z_dz(nz)/) )
elseif ( DepthV < z_dz(nz) ) then
if (myrank == 0) call MessageNotify( "E", "Damping_init", "DepthV is too thin. DelZ is %f", d=(/z_dz(nz)/) )
end if
end if
if ( DepthVb > 0.0d0 ) then
if ( DepthVb < z_dz(1) ) then
if (myrank == 0) call MessageNotify( "E", "Damping_init", "DepthVb is too thin. DelZ is %f", d=(/z_dz(nz)/) )
elseif ( DepthVb < z_dz(nz) ) then
if (myrank == 0) call MessageNotify( "E", "Damping_init", "DepthVb is too thin. DelZ is %f", d=(/z_dz(nz)/) )
end if
end if
!-----------------------------------------------------------------
! スポンジ層の減衰率
!
!水平方向の東側・西側境界
!
if ( DepthH > 0.0d0 ) then
do i = imin, imax
!スカラー格子点の西側境界
if ( x_X(i) < DepthH) then
xyz_Gamma(i,:,:) = xyz_Gamma(i,:,:) + ((1.0d0 - x_X(i) / DepthH) ** 3.0d0) / EFTime
end if
!フラックス格子点の西側境界
if ( p_X(i) < DepthH) then
pyz_Gamma(i,:,:) = pyz_Gamma(i,:,:) + ((1.0d0 - p_X(i) / DepthH) ** 3.0d0) / EFTime
end if
!スカラー格子点の東側境界
if ( x_X(i) > ( XMax - DepthH ) ) then
xyz_Gamma(i,:,:) = xyz_Gamma(i,:,:) + ((1.0d0 - (XMax - x_X(i)) / DepthH) ** 3.0d0) / EFTime
end if
!フラックス格子点の東側境界
if ( p_X(i) > ( XMax - DepthH ) ) then
pyz_Gamma(i,:,:) = pyz_Gamma(i,:,:) + ((1.0d0 - (XMax - p_X(i)) / DepthH) ** 3.0d0) / EFTime
end if
end do
! x 方向には同じ
!
xyr_Gamma = xyz_Gamma
xqz_Gamma = xyz_Gamma
end if
!水平方向の南側・北側境界
!
if ( DepthH > 0.0d0 ) then
do j = jmin, jmax
!スカラー格子点の南側境界
if ( y_Y(j) < DepthH) then
xyz_Gamma(:,j,:) = xyz_Gamma(:,j,:) + ((1.0d0 - y_Y(j) / DepthH) ** 3.0d0) / EFTime
end if
!フラックス格子点の南側境界
if ( q_Y(j) < DepthH) then
xqz_Gamma(:,j,:) = xqz_Gamma(:,j,:) + ((1.0d0 - q_Y(j) / DepthH) ** 3.0d0) / EFTime
end if
!スカラー格子点の北側境界
if ( y_Y(j) > ( YMax - DepthH ) ) then
xyz_Gamma(:,j,:) = xyz_Gamma(:,j,:) + ((1.0d0 - (YMax - y_Y(j)) / DepthH) ** 3.0d0) / EFTime
end if
!フラックス格子点の北側境界
if ( q_Y(j) > ( YMax - DepthH ) ) then
xqz_Gamma(:,j,:) = xqz_Gamma(:,j,:) + ((1.0d0 - (YMax - q_Y(j)) / DepthH) ** 3.0d0) / EFTime
end if
end do
! y 方向には同じ
!
pyz_Gamma = xyz_Gamma
xyr_Gamma = xyz_Gamma
end if
!鉛直方向の上部境界
!
if ( DepthV > 0.0d0 ) then
do k = kmin, kmax
!スカラー格子点
if ( z_Z(k) >= ( ZMax - DepthV ) ) then
xyz_Gamma(:,:,k) = xyz_Gamma(:,:,k) + (1.0d0 - dcos(Pi * (z_Z(k) - ZMax + DepthV) / DepthV)) / EFTime
end if
!フラックス格子点
if ( r_Z(k) >= ( ZMax - DepthV ) ) then
xyr_Gamma(:,:,k) = xyr_Gamma(:,:,k) + (1.0d0 - dcos(Pi * (r_Z(k) - ZMax + DepthV)/ DepthV)) / EFTime
end if
end do
! z 方向には同じ
!
pyz_Gamma = xyz_Gamma
xqz_Gamma = xyz_Gamma
end if
!鉛直方向の下部境界
!
if ( DepthVb > 0.0d0 ) then
do k = kmin, kmax
!スカラー格子点
if ( z_Z(k) <= (ZMin + DepthV) ) then
xyz_Gamma(:,:,k) = xyz_Gamma(:,:,k) + (1.0d0 + dcos(Pi * ( z_Z(k) - ZMin ) / DepthV) ) / EFTime
end if
!フラックス格子点
if ( r_Z(k) <= (ZMin + DepthV) ) then
xyr_Gamma(:,:,k) = xyr_Gamma(:,:,k) + (1.0d0 + dcos(Pi * ( r_Z(k) - ZMin )/ DepthV ) ) / EFTime
end if
end do
! z 方向には同じ
!
pyz_Gamma = xyz_Gamma
xqz_Gamma = xyz_Gamma
end if
! 確認
!
! i = 1; j = 1
! write(*,*) 'z_Z, xyz_Gamma'
! do k = kmin, kmax
! write(*,*) z_Z(k), '=>', xyz_Gamma(i,j,k)
! end do
! write(*,*) 'r_Z, xyr_Gamma'
! do k = kmin, kmax
! write(*,*) r_Z(k), '=>', xyr_Gamma(i,j,k)
! end do
!-----------------------------------------------------------------
! 値の確認
!
if (myrank == 0) then
call MessageNotify( "M", "Damping_init", "EFTime = %f", d=(/EFTime/) )
call MessageNotify( "M", "Damping_init", "DepthH = %f", d=(/DepthH/) )
call MessageNotify( "M", "Damping_init", "DepthV = %f", d=(/DepthV/) )
call MessageNotify( "M", "Damping_init", "DepthVb= %f", d=(/DepthVb/) )
end if
!-----------------------------------------------------------------
! 出力
!
call HistoryAutoAddVariable( varname='PTempSpng', dims=(/'x','y','z','t'/), longname='Damping term of potential temperature', units='K.s-1', xtype='float')
! call HistoryAutoAddVariable( &
! & varname='ExnerSpng', &
! & dims=(/'x','y','z','t'/), &
! & longname='Damping term of exner function', &
! & units='s-1', &
! & xtype='float')
call HistoryAutoAddVariable( varname='VelXSpng', dims=(/'x','y','z','t'/), longname='Damping term of VelX', units='m.s-1', xtype='float')
call HistoryAutoAddVariable( varname='VelYSpng', dims=(/'x','y','z','t'/), longname='Damping term of VelY', units='m.s-1', xtype='float')
call HistoryAutoAddVariable( varname='VelZSpng', 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) |
| 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) |
(inout)
subroutine SpongeLayer_forcing( pyz_VelXBl, xqz_VelYBl, xyr_VelZBl, xyz_PTempBl, pyz_DVelXDt, xqz_DVelYDt, xyr_DVelZDt, xyz_DPTempDt ) !(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(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) :: pyz_SpngVelX(imin:imax, jmin:jmax, kmin:kmax)
real(8) :: xqz_SpngVelY(imin:imax, jmin:jmax, kmin:kmax)
real(8) :: xyr_SpngVelZ(imin:imax, jmin:jmax, kmin:kmax)
real(8) :: xyz_SpngPTemp(imin:imax, jmin:jmax, kmin:kmax)
!--------------------------------------------------------
! 擾乱場の値をゼロに戻すような強制を与える
!
pyz_SpngVelX = - pyz_Gamma * pyz_VelXBl
xqz_SpngVelY = - xqz_Gamma * xqz_VelYBl
xyr_SpngVelZ = - xyr_Gamma * xyr_VelZBl
xyz_SpngPTemp = - xyz_Gamma * xyz_PTempBl
pyz_DVelXDt = pyz_DVelXDt + pyz_SpngVelX
xqz_DVelYDt = xqz_DVelYDt + xqz_SpngVelY
xyr_DVelZDt = xyr_DVelZDt + xyr_SpngVelZ
xyz_DPTempDt = xyz_DPTempDt + xyz_SpngPTemp
call HistoryAutoPut(TimeN, 'VelXSpng', pyz_SpngVelX(1:nx,1:ny,1:nz))
call HistoryAutoPut(TimeN, 'VelYSpng', xqz_SpngVelY(1:nx,1:ny,1:nz))
call HistoryAutoPut(TimeN, 'VelZSpng', xyr_SpngVelZ(1:nx,1:ny,1:nz))
call HistoryAutoPut(TimeN, 'PTempSpng', xyz_SpngPTemp(1:nx,1:ny,1:nz))
end subroutine SpongeLayer_forcing