!= Module Damping ! ! Authors:: SUGIYAMA Ko-ichiro, ODAKA Masatsugu ! Version:: $Id: damping.f90,v 1.7 2011-10-05 03:35:23 sugiyama Exp $ ! Tag Name:: $Name: arare5-20111010 $ ! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved. ! License:: See COPYRIGHT[link:../../COPYRIGHT] ! !== Overview ! !減衰率とその計算を行うためのパッケージ型モジュール ! * 音波減衰項の係数 ! * スポンジ層の設定(境界付近で波の反射を抑え吸収するための層) ! !== Error Handling ! !== Bugs ! !== Note ! ! * この関数は, 基本場が零な変数(速度, エクスナー関数の擾乱)に ! 適用することを想定している. ! * 各格子点に対する関数を定義する必要がある ! !== Future Plans ! ! module Damping ! !減衰率とその計算を行うためのパッケージ型モジュール ! * 音波減衰項の係数 ! * スポンジ層の設定(境界付近で波の反射を抑え吸収するための層) ! !モジュール読み込み use dc_types, only : DP !暗黙の型宣言禁止 implicit none !private 属性を指定 private !関数には public 属性を指定 public Damping_Init public SpongeLayer_forcing !変数定義 real(DP), save :: EFTime = 100.0d0 !スポンジ層の e-folding time real(DP), save :: DampDepthH = 0.0d0 !スポンジ層の厚さ(水平方向) real(DP), save :: DampDepthV = 0.0d0 !スポンジ層の厚さ(鉛直方向) real(DP), allocatable, save :: xyz_Gamma(:,:,:) !xyz 格子減衰係数(水平方向) real(DP), allocatable, save :: pyz_Gamma(:,:,:) !pyz 格子減衰係数(鉛直方向) real(DP), allocatable, save :: xqz_Gamma(:,:,:) !xqz 格子減衰係数(鉛直方向) real(DP), allocatable, save :: xyr_Gamma(:,:,:) !xyr 格子減衰係数(鉛直方向) contains !!!------------------------------------------------------------------------!!! 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, &! x 方向の配列の下限 & imax, &! x 方向の配列の上限 & jmin, &! y 方向の配列の下限 & jmax, &! y 方向の配列の上限 & kmin, &! z 方向の配列の下限 & kmax, &! z 方向の配列の上限 & nx, &! x 方向の物理領域の上限 & ny, &! y 方向の物理領域の上限 & nz ! z 方向の物理領域の上限 use axesset, only: 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 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,:,:) = ((1.0d0 - x_X(i) / DampDepthH) ** 3.0d0) / EFTime end if !フラックス格子点の西側境界 if ( p_X(i) < DampDepthH) then pyz_Gamma(i,:,:) = ((1.0d0 - p_X(i) / DampDepthH) ** 3.0d0) / EFTime end if !スカラー格子点の東側境界 if ( x_X(i) > ( XMax - DampDepthH ) ) then 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,:,:) = & & ((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,:) = ((1.0d0 - y_Y(j) / DampDepthH) ** 3.0d0) / EFTime end if !フラックス格子点の西側境界 if ( q_Y(j) < DampDepthH) then xqz_Gamma(:,j,:) = ((1.0d0 - q_Y(j) / DampDepthH) ** 3.0d0) / EFTime end if !スカラー格子点の東側境界 if ( y_Y(j) > ( YMax - DampDepthH ) ) then 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,:) = & & ((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) = & & (1.0d0 - dcos(Pi * (z_Z(k) - ZMax + DampDepthV) / DampDepthV)) & & / EFTime end if !フラックス格子点 if ( r_Z(k) >= ( ZMax - DampDepthV ) ) then 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 SpongeLayer_forcing( & & pyz_VelXBl, xqz_VelYBl, xyr_VelZBl, xyz_PTempBl, xyz_ExnerBl, & !(in) & pyz_VelXAl, xqz_VelYAl, xyr_VelZAl, xyz_PTempAl, xyz_ExnerAl ) !(inout) use gtool_historyauto, only: HistoryAutoPut use timeset, only: TimeN, DelTimeLong use gridset, only: imin, &! x 方向の配列の下限 & imax, &! x 方向の配列の上限 & jmin, &! y 方向の配列の下限 & jmax, &! y 方向の配列の上限 & kmin, &! z 方向の配列の下限 & kmax, &! z 方向の配列の上限 & nx, &! x 方向の物理領域の上限 & ny, &! y 方向の物理領域の上限 & 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_VelXAl(imin:imax, jmin:jmax, kmin:kmax) real(8), intent(inout) :: xqz_VelYAl(imin:imax, jmin:jmax, kmin:kmax) real(8), intent(inout) :: xyr_VelZAl(imin:imax, jmin:jmax, kmin:kmax) real(8), intent(inout) :: xyz_PTempAl(imin:imax, jmin:jmax, kmin:kmax) real(8), intent(inout) :: xyz_ExnerAl(imin:imax, jmin:jmax, kmin:kmax) real(8) :: pyz_VelXAl0(imin:imax, jmin:jmax, kmin:kmax) real(8) :: xqz_VelYAl0(imin:imax, jmin:jmax, kmin:kmax) real(8) :: xyr_VelZAl0(imin:imax, jmin:jmax, kmin:kmax) real(8) :: xyz_PTempAl0(imin:imax, jmin:jmax, kmin:kmax) real(8) :: xyz_ExnerAl0(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_VelXAl0 = pyz_VelXAl pyz_DampVelX = - xyz_Gamma * pyz_VelXBl pyz_VelXAl = pyz_VelXAl0 + (2.0d0 * DelTimeLong) * pyz_DampVelX xqz_VelYAl0 = xqz_VelYAl xqz_DampVelY = - xqz_Gamma * xqz_VelYBl xqz_VelYAl = xqz_VelYAl0 + (2.0d0 * DelTimeLong) * xqz_DampVelY xyr_VelZAl0 = xyr_VelZAl xyr_DampVelZ = - xyr_Gamma * xyr_VelZBl xyr_VelZAl = xyr_VelZAl0 + (2.0d0 * DelTimeLong) * xyr_DampVelZ xyz_PTempAl0 = xyz_PTempAl xyz_DampPTemp = - xyz_Gamma * xyz_PTempBl xyz_PTempAl = xyz_PTempAl0 + (2.0d0 * DelTimeLong) * xyz_DampPTemp xyz_ExnerAl0 = xyz_ExnerAl xyz_DampExner = - xyz_Gamma * xyz_ExnerBl xyz_ExnerAl = xyz_ExnerAl0 + (2.0d0 * DelTimeLong) * 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 end module Damping