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 write(*,*) "++++++++ damping.f90 +++++++++", Time, DepthH, DepthV !値の入力 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='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) |
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_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) real(8) :: xyz_SpngExner(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 xyz_SpngExner = - xyz_Gamma * xyz_ExnerBl pyz_DVelXDt = pyz_DVelXDt + pyz_SpngVelX xqz_DVelYDt = xqz_DVelYDt + xqz_SpngVelY xyr_DVelZDt = xyr_DVelZDt + xyr_SpngVelZ xyz_DPTempDt = xyz_DPTempDt + xyz_SpngPTemp xyz_DExnerDt = xyz_DExnerDt + xyz_SpngExner 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)) call HistoryAutoPut(TimeN, 'ExnerSpng', xyz_SpngExner(1:nx,1:ny,1:nz)) end subroutine SpongeLayer_forcing