Class Damping
In: ../src/utils/damping.f90

減衰率とその計算を行うためのパッケージ型モジュール

 * 音波減衰項の係数
 * スポンジ層の設定(境界付近で波の反射を抑え吸収するための層)

Methods

Included Modules

dc_types dc_iounit dc_message gtool_historyauto mpi_wrapper gridset axesset namelist_util timeset

Public Instance methods

Subroutine :

音波減衰項とスポンジ層の減衰係数の初期化

This procedure input/output NAMELIST#damping_nml .

[Source]

  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)

[Source]

  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