Linux R2.6.5-7.282-sn2 FORTRAN90/SX Rev.360 Tue Oct 11 12:33:51 2011 FILE NAME: damping.f90 PROGRAM NAME: damping DIAGNOSTIC LIST LINE LEVEL( NO.): DIAGNOSTIC MESSAGE 114 vec ( 4): Vectorized array expression. 115 vec ( 4): Vectorized array expression. 116 vec ( 4): Vectorized array expression. 117 vec ( 4): Vectorized array expression. 137 vec ( 3): Unvectorized loop. 140 vec ( 4): Vectorized array expression. 145 vec ( 4): Vectorized array expression. 150 vec ( 4): Vectorized array expression. 156 vec ( 4): Vectorized array expression. 164 opt ( 11): Fused array assignments. :line 164 - 165 164 vec ( 4): Vectorized array expression. 164 vec ( 4): Vectorized array expression. 178 vec ( 3): Unvectorized loop. 181 vec ( 4): Vectorized array expression. 181 vec ( 4): Vectorized array expression. 186 vec ( 4): Vectorized array expression. 186 vec ( 4): Vectorized array expression. 191 vec ( 4): Vectorized array expression. 191 vec ( 4): Vectorized array expression. 197 vec ( 4): Vectorized array expression. 197 vec ( 4): Vectorized array expression. 205 opt ( 11): Fused array assignments. :line 205 - 206 205 vec ( 4): Vectorized array expression. 205 vec ( 4): Vectorized array expression. 215 vec ( 3): Unvectorized loop. 218 vec ( 4): Vectorized array expression. 225 vec ( 4): Vectorized array expression. 234 opt ( 11): Fused array assignments. :line 234 - 235 234 vec ( 4): Vectorized array expression. 234 vec ( 4): Vectorized array expression. 330 opt ( 11): Fused array assignments. :line 330 - 348 330 vec ( 4): Vectorized array expression. Linux R2.6.5-7.282-sn2 FORTRAN90/SX Rev.360 Tue Oct 11 12:33:51 2011 FILE NAME: damping.f90 PROGRAM NAME: damping TRANSFORMATION LIST LINE FORTRAN STATEMENT 1 != Module Damping 2 ! 3 ! Authors:: SUGIYAMA Ko-ichiro, ODAKA Masatsugu 4 ! Version:: $Id: damping.f90,v 1.7 2011-10-05 03:35:23 sugiyama Exp $ 5 ! Tag Name:: $Name: arare5-20111010 $ 6 ! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved. 7 ! License:: See COPYRIGHT[link:../../COPYRIGHT] 8 ! 9 !== Overview 10 ! 11 !減衰率とその計算を行うためのパッケージ型モジュール 12 ! * 音波減衰項の係数 13 ! * スポンジ層の設定(境界付近で波の反射を抑え吸収するための層) 14 ! 15 !== Error Handling 16 ! 17 !== Bugs 18 ! 19 !== Note 20 ! 21 ! * この関数は, 基本場が零な変数(速度, エクスナー関数の擾乱)に 22 ! 適用することを想定している. 23 ! * 各格子点に対する関数を定義する必要がある 24 ! 25 !== Future Plans 26 ! 27 ! 28 29 module Damping 30 ! 31 !減衰率とその計算を行うためのパッケージ型モジュール 32 ! * 音波減衰項の係数 33 ! * スポンジ層の設定(境界付近で波の反射を抑え吸収するための層) 34 ! 35 36 !モジュール読み込み 37 use dc_types, only : DP 38 39 !暗黙の型宣言禁止 40 implicit none 41 42 !private 属性を指定 43 private 44 45 !関数には public 属性を指定 46 public Damping_Init 47 public SpongeLayer_forcing 48 49 !変数定義 50 real(DP), save :: EFTime = 100.0d0 !スポンジ層の e-folding time 51 real(DP), save :: DampDepthH = 0.0d0 !スポンジ層の厚さ(水平方向) 52 real(DP), save :: DampDepthV = 0.0d0 !スポンジ層の厚さ(鉛直方向) 53 real(DP), allocatable, save :: xyz_Gamma(:,:,:) !xyz 格子減衰係数(水平方向) 54 real(DP), allocatable, save :: pyz_Gamma(:,:,:) !pyz 格子減衰係数(鉛直方向) 55 real(DP), allocatable, save :: xqz_Gamma(:,:,:) !xqz 格子減衰係数(鉛直方向) 56 real(DP), allocatable, save :: xyr_Gamma(:,:,:) !xyr 格子減衰係数(鉛直方向) 57 58 contains 59 60 !!!------------------------------------------------------------------------!!! 61 subroutine Damping_Init 62 ! 63 ! 音波減衰項とスポンジ層の減衰係数の初期化 64 ! 65 use dc_iounit, only: FileOpen 66 use dc_message, only: MessageNotify 67 use gtool_historyauto, only: HistoryAutoAddVariable 68 use mpi_wrapper, only: myrank 69 use gridset, only: imin, &! x 方向の配列の下限 70 & imax, &! x 方向の配列の上限 71 & jmin, &! y 方向の配列の下限 72 & jmax, &! y 方向の配列の上限 73 & kmin, &! z 方向の配列の下限 74 & kmax, &! z 方向の配列の上限 75 & nx, &! x 方向の物理領域の上限 76 & ny, &! y 方向の物理領域の上限 77 & nz ! z 方向の物理領域の上限 78 use axesset, only: x_X, &!X 座標軸(スカラー格子点) 79 & y_Y, &!Y 座標軸(スカラー格子点) 80 & z_Z, &!Z 座標軸(スカラー格子点) 81 & p_X, &!X 座標軸(フラックス格子点) 82 & q_Y, &!Y 座標軸(フラックス格子点) 83 & r_Z, &!Z 座標軸(フラックス格子点) 84 & x_dx, y_dy, z_dz, &! 格子間隔 85 & XMax, &!X 座標の最大値 86 & YMax, &!Y 座標の最大値 87 & ZMax !Z 座標の最大値 88 use namelist_util, only: namelist_filename 89 90 !暗黙の型宣言禁止 91 implicit none 92 93 !変数定義 94 real(DP) :: Time ! 95 real(DP) :: DepthH !スポンジ層の厚さ(水平方向) 96 real(DP) :: DepthV !スポンジ層の厚さ(鉛直方向) 97 real(DP), parameter :: Pi =3.1415926535897932385d0 !円周率 98 integer :: unit ! 装置番号 99 integer :: i, j, k 100 101 !NAMELIST から取得 102 NAMELIST /damping_nml/ Time, DepthH, DepthV 103 104 call FileOpen(unit, file=namelist_filename, mode='r') 105 read(unit, NML=damping_nml) 106 close(unit) 107 108 !初期化 109 allocate( & 110 & xyz_Gamma(imin:imax,jmin:jmax,kmin:kmax), & 111 & pyz_Gamma(imin:imax,jmin:jmax,kmin:kmax), & 112 & xqz_Gamma(imin:imax,jmin:jmax,kmin:kmax), & 113 & xyr_Gamma(imin:imax,jmin:jmax,kmin:kmax) ) 114 xyz_Gamma = 0.0d0 . !CDIR NODEP . !CDIR NOASSUME . do t719 = 1, (t4 - t5 + 1)*(t2 - t3 + 1)*(t0 - t1 + 1) . xyz_gamma(t1+t719-1,t3,t5) = 0.0000000000000000e+000 . end do 115 pyz_Gamma = 0.0d0 . !CDIR NODEP . !CDIR NOASSUME . do t728 = 1, (t15 - t16 + 1)*(t13 - t14 + 1)*(t11 - t12 + 1) . pyz_gamma(t12+t728-1,t14,t16) = 0.0000000000000000e+000 . end do 116 xqz_Gamma = 0.0d0 . !CDIR NODEP . !CDIR NOASSUME . do t737=1,(t26-t27+1)*(xqz_gamma.DSC.U2-t25+1)*(t22-t23+1) . xqz_gamma(t23+t737-1,t25,t27) = 0.0000000000000000e+000 . end do 117 xyr_Gamma = 0.0d0 . !CDIR NODEP . !CDIR NOASSUME . do t746 = 1, (t37 - t38 + 1)*(t35 - t36 + 1)*(t33 - t34 + 1) . xyr_gamma(t34+t746-1,t36,t38) = 0.0000000000000000e+000 . end do 118 119 !値の入力 120 EFTime = Time 121 DampDepthH = DepthH 122 DampDepthV = DepthV 123 124 !----------------------------------------------------------------- 125 ! スポンジ層の減衰率 126 ! 127 !水平方向の東側・西側境界 128 if ( DampDepthH < x_dx(1) ) then 129 if (myrank == 0) & 130 & call MessageNotify( "W", "Damping_init", "DampDepthH is too thin. DelX is %f", d=(/x_dx(1)/)) 131 132 else if ( DampDepthH < x_dx(nx) ) then 133 if (myrank == 0) & 134 & call MessageNotify( "W", "Damping_init", "DampDepthH is too thin. DelX is %f", d=(/x_dx(nx)/)) 135 136 else 137 do i = imin, imax 138 !スカラー格子点の西側境界 139 if ( x_X(i) < DampDepthH) then 140 xyz_Gamma(i,:,:) = ((1.0d0 - x_X(i) / DampDepthH) ** 3.0d0) / EFTime . D1 = 1.D0/dampdepthh . D2 = 1.D0/eftime . !CDIR NODEP . !CDIR NOASSUME . do t863 = 1, (t4 - t5 + 1)*(t2 - t3 + 1) . xyz_gamma(i,t3+t863-1,t5) = (1.00000000000000e+000 - x_x(i)*D1) . 1 **3.00000000000000e+000*D2 . end do 141 end if 142 143 !フラックス格子点の西側境界 144 if ( p_X(i) < DampDepthH) then 145 pyz_Gamma(i,:,:) = ((1.0d0 - p_X(i) / DampDepthH) ** 3.0d0) / EFTime . D3 = 1.D0/dampdepthh . D4 = 1.D0/eftime . !CDIR NODEP . !CDIR NOASSUME . do t857 = 1, (t15 - t16 + 1)*(t13 - t14 + 1) . pyz_gamma(i,t14+t857-1,t16) = (1.00000000000000e+000 - p_x(i)* . 1 D3)**3.00000000000000e+000*D4 . end do 146 end if 147 148 !スカラー格子点の東側境界 149 if ( x_X(i) > ( XMax - DampDepthH ) ) then 150 xyz_Gamma(i,:,:) = & . D5 = 1.D0/dampdepthh . D6 = 1.D0/eftime . !CDIR NODEP . !CDIR NOASSUME . do t851 = 1, (t4 - t5 + 1)*(t2 - t3 + 1) . xyz_gamma(i,t3+t851-1,t5) = (1.00000000000000e+000 - (xmax - . 1 x_x(i))*D5)**3.00000000000000e+000*D6 . end do 151 & ((1.0d0 - (XMax - x_X(i)) / DampDepthH) ** 3.0d0) / EFTime 152 end if 153 154 !フラックス格子点の東側境界 155 if ( p_X(i) > ( XMax - DampDepthH ) ) then 156 pyz_Gamma(i,:,:) = & . D7 = 1.D0/dampdepthh . D8 = 1.D0/eftime . !CDIR NODEP . !CDIR NOASSUME . do t845 = 1, (t15 - t16 + 1)*(t13 - t14 + 1) . pyz_gamma(i,t14+t845-1,t16) = (1.00000000000000e+000 - (xmax - . 1 p_x(i))*D7)**3.00000000000000e+000*D8 . end do 157 & ((1.0d0 - (XMax - p_X(i)) / DampDepthH) ** 3.0d0) / EFTime 158 end if 159 end do 160 end if 161 162 ! x 方向には同じ 163 ! 164 xyr_Gamma = xyz_Gamma . if (t35 - t36 + 1 .gt. 0) then . J1 = and(t35 - t36 + 1,1) . !CDIR NODEP . do t757 = 1, J1 . !CDIR NODEP . do t759 = 1, t33 - t34 + 1 . xyr_gamma(t34+t759-1,t757-1+t36,t755+t38) = xyz_gamma(t1+ . 1 t759-1,t757-1+t3,t755+t5) . xqz_gamma(t23+t759-1,t757-1+t25,t755+t27) = xyz_gamma(t1+ . 1 t759-1,t757-1+t3,t755+t5) . end do . end do . !CDIR NODEP . do t757 = J1 + 1, t35 - t36 + 1, 2 . !CDIR NODEP . do t759 = 1, t33 - t34 + 1 . xyr_gamma(t34+t759-1,t757-1+t36,t755+t38) = xyz_gamma(t1+ . 1 t759-1,t757-1+t3,t755+t5) . xyr_gamma(t34+t759-1,t757+t36,t755+t38) = xyz_gamma(t1+ . 1 t759-1,t757+t3,t755+t5) . xqz_gamma(t23+t759-1,t757-1+t25,t755+t27) = xyz_gamma(t1+ . 1 t759-1,t757-1+t3,t755+t5) . xqz_gamma(t23+t759-1,t757+t25,t755+t27) = xyz_gamma(t1+ . 1 t759-1,t757+t3,t755+t5) . end do . end do . endif 165 xqz_Gamma = xyz_Gamma 166 167 168 !水平方向の南側・北側境界 169 if ( DampDepthH < y_dy(1) ) then 170 if (myrank == 0) & 171 & call MessageNotify( "W", "Damping_init", "DampDepthH is too thin. DelY is %f", d=(/x_dx(1)/)) 172 173 else if ( DampDepthH < y_dy(ny) ) then 174 if (myrank == 0) & 175 & call MessageNotify( "W", "Damping_init", "DampDepthH is too thin. DelY is %f", d=(/y_dy(ny)/)) 176 177 else 178 do j = jmin, jmax 179 !スカラー格子点の西側境界 180 if ( y_Y(j) < DampDepthH) then 181 xyz_Gamma(:,j,:) = ((1.0d0 - y_Y(j) / DampDepthH) ** 3.0d0) / EFTime . if (t4 - t5 + 1 .gt. 0) then . J2 = and(t4 - t5 + 1,3) . !CDIR NODEP . do t839 = 1, J2 . D13 = 1.D0/dampdepthh . D14 = 1.D0/eftime . !CDIR NODEP . do t841 = 1, t0 - t1 + 1 . xyz_gamma(t1+t841-1,j,t839-1+t5) = (1.00000000000000e+000 . 1 - y_y(j)*D13)**3.00000000000000e+000*D14 . end do . end do . !CDIR NODEP . do t839 = J2 + 1, t4 - t5 + 1, 4 . D15 = 1.D0/dampdepthh . D16 = 1.D0/eftime . D17 = 1.D0/dampdepthh . D18 = 1.D0/eftime . D19 = 1.D0/dampdepthh . D20 = 1.D0/eftime . D21 = 1.D0/dampdepthh . D22 = 1.D0/eftime . D9 = y_y(j) . !CDIR NODEP . do t841 = 1, t0 - t1 + 1 . xyz_gamma(t1+t841-1,j,t839-1+t5) = (1.00000000000000e+000 . 1 - D9*D15)**3.00000000000000e+000*D16 . xyz_gamma(t1+t841-1,j,t839+t5) = (1.00000000000000e+000 . 1 - D9*D17)**3.00000000000000e+000*D18 . xyz_gamma(t1+t841-1,j,t839+1+t5) = (1.00000000000000e+000 . 1 - D9*D19)**3.00000000000000e+000*D20 . xyz_gamma(t1+t841-1,j,t839+2+t5) = (1.00000000000000e+000 . 1 - D9*D21)**3.00000000000000e+000*D22 . end do . end do . endif 182 end if 183 184 !フラックス格子点の西側境界 185 if ( q_Y(j) < DampDepthH) then 186 xqz_Gamma(:,j,:) = ((1.0d0 - q_Y(j) / DampDepthH) ** 3.0d0) / EFTime . if (t26 - t27 + 1 .gt. 0) then . J3 = and(t26 - t27 + 1,3) . !CDIR NODEP . do t833 = 1, J3 . D23 = 1.D0/dampdepthh . D24 = 1.D0/eftime . !CDIR NODEP . do t835 = 1, t22 - t23 + 1 . xqz_gamma(t23+t835-1,j,t833-1+t27) = ( . 1 1.00000000000000e+000 - q_y(j)*D23)** . 2 3.00000000000000e+000*D24 . end do . end do . !CDIR NODEP . do t833 = J3 + 1, t26 - t27 + 1, 4 . D25 = 1.D0/dampdepthh . D26 = 1.D0/eftime . D27 = 1.D0/dampdepthh . D28 = 1.D0/eftime . D29 = 1.D0/dampdepthh . D30 = 1.D0/eftime . D31 = 1.D0/dampdepthh . D32 = 1.D0/eftime . D10 = q_y(j) . !CDIR NODEP . do t835 = 1, t22 - t23 + 1 . xqz_gamma(t23+t835-1,j,t833-1+t27) = ( . 1 1.00000000000000e+000 - D10*D25)** . 2 3.00000000000000e+000*D26 . xqz_gamma(t23+t835-1,j,t833+t27) = (1.00000000000000e+000 . 1 - D10*D27)**3.00000000000000e+000*D28 . xqz_gamma(t23+t835-1,j,t833+1+t27) = ( . 1 1.00000000000000e+000 - D10*D29)** . 2 3.00000000000000e+000*D30 . xqz_gamma(t23+t835-1,j,t833+2+t27) = ( . 1 1.00000000000000e+000 - D10*D31)** . 2 3.00000000000000e+000*D32 . end do . end do . endif 187 end if 188 189 !スカラー格子点の東側境界 190 if ( y_Y(j) > ( YMax - DampDepthH ) ) then 191 xyz_Gamma(:,j,:) = & . if (t4 - t5 + 1 .gt. 0) then . J4 = and(t4 - t5 + 1,3) . !CDIR NODEP . do t827 = 1, J4 . D33 = 1.D0/dampdepthh . D34 = 1.D0/eftime . !CDIR NODEP . do t829 = 1, t0 - t1 + 1 . xyz_gamma(t1+t829-1,j,t827-1+t5) = (1.00000000000000e+000 . 1 - (ymax - y_y(j))*D33)**3.00000000000000e+000*D34 . end do . end do . !CDIR NODEP . do t827 = J4 + 1, t4 - t5 + 1, 4 . D35 = 1.D0/dampdepthh . D36 = 1.D0/eftime . D37 = 1.D0/dampdepthh . D38 = 1.D0/eftime . D39 = 1.D0/dampdepthh . D40 = 1.D0/eftime . D41 = 1.D0/dampdepthh . D42 = 1.D0/eftime . D11 = y_y(j) . !CDIR NODEP . do t829 = 1, t0 - t1 + 1 . xyz_gamma(t1+t829-1,j,t827-1+t5) = (1.00000000000000e+000 . 1 - (ymax - D11)*D35)**3.00000000000000e+000*D36 . xyz_gamma(t1+t829-1,j,t827+t5) = (1.00000000000000e+000 . 1 - (ymax - D11)*D37)**3.00000000000000e+000*D38 . xyz_gamma(t1+t829-1,j,t827+1+t5) = (1.00000000000000e+000 . 1 - (ymax - D11)*D39)**3.00000000000000e+000*D40 . xyz_gamma(t1+t829-1,j,t827+2+t5) = (1.00000000000000e+000 . 1 - (ymax - D11)*D41)**3.00000000000000e+000*D42 . end do . end do . endif 192 & ((1.0d0 - (YMax - y_Y(j)) / DampDepthH) ** 3.0d0) / EFTime 193 end if 194 195 !フラックス格子点の東側境界 196 if ( q_Y(j) > ( YMax - DampDepthH ) ) then 197 xqz_Gamma(:,j,:) = & . if (t26 - t27 + 1 .gt. 0) then . J5 = and(t26 - t27 + 1,3) . !CDIR NODEP . do t821 = 1, J5 . D43 = 1.D0/dampdepthh . D44 = 1.D0/eftime . !CDIR NODEP . do t823 = 1, t22 - t23 + 1 . xqz_gamma(t23+t823-1,j,t821-1+t27) = ( . 1 1.00000000000000e+000 - (ymax - q_y(j))*D43)** . 2 3.00000000000000e+000*D44 . end do . end do . !CDIR NODEP . do t821 = J5 + 1, t26 - t27 + 1, 4 . D45 = 1.D0/dampdepthh . D46 = 1.D0/eftime . D47 = 1.D0/dampdepthh . D48 = 1.D0/eftime . D49 = 1.D0/dampdepthh . D50 = 1.D0/eftime . D51 = 1.D0/dampdepthh . D52 = 1.D0/eftime . D12 = q_y(j) . !CDIR NODEP . do t823 = 1, t22 - t23 + 1 . xqz_gamma(t23+t823-1,j,t821-1+t27) = ( . 1 1.00000000000000e+000 - (ymax - D12)*D45)** . 2 3.00000000000000e+000*D46 . xqz_gamma(t23+t823-1,j,t821+t27) = (1.00000000000000e+000 . 1 - (ymax - D12)*D47)**3.00000000000000e+000*D48 . xqz_gamma(t23+t823-1,j,t821+1+t27) = ( . 1 1.00000000000000e+000 - (ymax - D12)*D49)** . 2 3.00000000000000e+000*D50 . xqz_gamma(t23+t823-1,j,t821+2+t27) = ( . 1 1.00000000000000e+000 - (ymax - D12)*D51)** . 2 3.00000000000000e+000*D52 . end do . end do . endif 198 & ((1.0d0 - (YMax - q_Y(j)) / DampDepthH) ** 3.0d0) / EFTime 199 end if 200 end do 201 end if 202 203 ! y 方向には同じ 204 ! 205 pyz_Gamma = xyz_Gamma . if (t13 - t14 + 1 .gt. 0) then . J6 = and(t13 - t14 + 1,1) . !CDIR NODEP . do t775 = 1, J6 . !CDIR NODEP . do t777 = 1, t11 - t12 + 1 . pyz_gamma(t12+t777-1,t775-1+t14,t773+t16) = xyz_gamma(t1+ . 1 t777-1,t775-1+t3,t773+t5) . xyr_gamma(t34+t777-1,t775-1+t36,t773+t38) = xyz_gamma(t1+ . 1 t777-1,t775-1+t3,t773+t5) . end do . end do . !CDIR NODEP . do t775 = J6 + 1, t13 - t14 + 1, 2 . !CDIR NODEP . do t777 = 1, t11 - t12 + 1 . pyz_gamma(t12+t777-1,t775-1+t14,t773+t16) = xyz_gamma(t1+ . 1 t777-1,t775-1+t3,t773+t5) . pyz_gamma(t12+t777-1,t775+t14,t773+t16) = xyz_gamma(t1+ . 1 t777-1,t775+t3,t773+t5) . xyr_gamma(t34+t777-1,t775-1+t36,t773+t38) = xyz_gamma(t1+ . 1 t777-1,t775-1+t3,t773+t5) . xyr_gamma(t34+t777-1,t775+t36,t773+t38) = xyz_gamma(t1+ . 1 t777-1,t775+t3,t773+t5) . end do . end do . endif 206 xyr_Gamma = xyz_Gamma 207 208 209 !鉛直方向の上部境界 210 if ( DampDepthV < z_dz(nz) ) then 211 if (myrank == 0) & 212 & call MessageNotify( "W", "Damping_init", "DampDepthV is too thin. DelZ is %f", d=(/z_dz(nz)/) ) 213 214 else 215 do k = kmin, kmax 216 !スカラー格子点 217 if ( z_Z(k) >= ( ZMax - DampDepthV ) ) then 218 xyz_Gamma(:,:,k) = & . D53 = 3.14159265358979e+000/dampdepthv . D54 = 1.D0/eftime . !CDIR NODEP . !CDIR NOASSUME . do t815 = 1, (t2 - t3 + 1)*(t0 - t1 + 1) . xyz_gamma(t1+t815-1,t3,k) = (1.00000000000000e+000 - dcos((z_z( . 1 k)-zmax+dampdepthv)*D53))*D54 . end do 219 & (1.0d0 - dcos(Pi * (z_Z(k) - ZMax + DampDepthV) / DampDepthV)) & 220 & / EFTime 221 end if 222 223 !フラックス格子点 224 if ( r_Z(k) >= ( ZMax - DampDepthV ) ) then 225 xyr_Gamma(:,:,k) = & . D55 = 3.14159265358979e+000/dampdepthv . D56 = 1.D0/eftime . !CDIR NODEP . !CDIR NOASSUME . do t809 = 1, (t35 - t36 + 1)*(t33 - t34 + 1) . xyr_gamma(t34+t809-1,t36,k) = (1.00000000000000e+000 - dcos(( . 1 r_z(k)-zmax+dampdepthv)*D55))*D56 . end do 226 & (1.0d0 - dcos(Pi * (r_Z(k) - ZMax + DampDepthV)/ DampDepthV)) & 227 & / EFTime 228 end if 229 end do 230 end if 231 232 ! z 方向には同じ 233 ! 234 pyz_Gamma = xyz_Gamma . if (t13 - t14 + 1 .gt. 0) then . J7 = and(t13 - t14 + 1,1) . !CDIR NODEP . do t793 = 1, J7 . !CDIR NODEP . do t795 = 1, t11 - t12 + 1 . pyz_gamma(t12+t795-1,t793-1+t14,t791+t16) = xyz_gamma(t1+ . 1 t795-1,t793-1+t3,t791+t5) . xqz_gamma(t23+t795-1,t793-1+t25,t791+t27) = xyz_gamma(t1+ . 1 t795-1,t793-1+t3,t791+t5) . end do . end do . !CDIR NODEP . do t793 = J7 + 1, t13 - t14 + 1, 2 . !CDIR NODEP . do t795 = 1, t11 - t12 + 1 . pyz_gamma(t12+t795-1,t793-1+t14,t791+t16) = xyz_gamma(t1+ . 1 t795-1,t793-1+t3,t791+t5) . pyz_gamma(t12+t795-1,t793+t14,t791+t16) = xyz_gamma(t1+ . 1 t795-1,t793+t3,t791+t5) . xqz_gamma(t23+t795-1,t793-1+t25,t791+t27) = xyz_gamma(t1+ . 1 t795-1,t793-1+t3,t791+t5) . xqz_gamma(t23+t795-1,t793+t25,t791+t27) = xyz_gamma(t1+ . 1 t795-1,t793+t3,t791+t5) . end do . end do . endif 235 xqz_Gamma = xyz_Gamma 236 237 238 !----------------------------------------------------------------- 239 ! 値の確認 240 ! 241 if (myrank == 0) then 242 call MessageNotify( "M", "Damping_init", "EFTime = %f", d=(/EFTime/) ) 243 call MessageNotify( "M", "Damping_init", "DampDepthH = %f", d=(/DampDepthH/) ) 244 call MessageNotify( "M", "Damping_init", "DampDepthV = %f", d=(/DampDepthV/) ) 245 end if 246 247 248 !----------------------------------------------------------------- 249 ! 出力 250 ! 251 call HistoryAutoAddVariable( & 252 & varname='PTempDamp', & 253 & dims=(/'x','y','z','t'/), & 254 & longname='Damping term of potential temperature', & 255 & units='K.s-1', & 256 & xtype='float') 257 258 call HistoryAutoAddVariable( & 259 & varname='ExnerDamp', & 260 & dims=(/'x','y','z','t'/), & 261 & longname='Damping term of exner function', & 262 & units='s-1', & 263 & xtype='float') 264 265 call HistoryAutoAddVariable( & 266 & varname='VelXDamp', & 267 & dims=(/'x','y','z','t'/), & 268 & longname='Damping term of VelX', & 269 & units='m.s-1', & 270 & xtype='float') 271 272 call HistoryAutoAddVariable( & 273 & varname='VelYDamp', & 274 & dims=(/'x','y','z','t'/), & 275 & longname='Damping term of VelY', & 276 & units='m.s-1', & 277 & xtype='float') 278 279 call HistoryAutoAddVariable( & 280 & varname='VelZDamp', & 281 & dims=(/'x','y','z','t'/), & 282 & longname='Damping term of VelZ', & 283 & units='m.s-1', & 284 & xtype='float') 285 286 end subroutine Damping_Init 287 288 289 290 subroutine SpongeLayer_forcing( & 291 & pyz_VelXBl, xqz_VelYBl, xyr_VelZBl, xyz_PTempBl, xyz_ExnerBl, & !(in) 292 & pyz_VelXAl, xqz_VelYAl, xyr_VelZAl, xyz_PTempAl, xyz_ExnerAl ) !(inout) 293 294 use gtool_historyauto, only: HistoryAutoPut 295 use timeset, only: TimeN, DelTimeLong 296 use gridset, only: imin, &! x 方向の配列の下限 297 & imax, &! x 方向の配列の上限 298 & jmin, &! y 方向の配列の下限 299 & jmax, &! y 方向の配列の上限 300 & kmin, &! z 方向の配列の下限 301 & kmax, &! z 方向の配列の上限 302 & nx, &! x 方向の物理領域の上限 303 & ny, &! y 方向の物理領域の上限 304 & nz ! z 方向の物理領域の上限 305 306 !暗黙の型宣言禁止 307 implicit none 308 309 real(8), intent(in) :: pyz_VelXBl(imin:imax, jmin:jmax, kmin:kmax) 310 real(8), intent(in) :: xqz_VelYBl(imin:imax, jmin:jmax, kmin:kmax) 311 real(8), intent(in) :: xyr_VelZBl(imin:imax, jmin:jmax, kmin:kmax) 312 real(8), intent(in) :: xyz_PTempBl(imin:imax, jmin:jmax, kmin:kmax) 313 real(8), intent(in) :: xyz_ExnerBl(imin:imax, jmin:jmax, kmin:kmax) 314 real(8), intent(inout) :: pyz_VelXAl(imin:imax, jmin:jmax, kmin:kmax) 315 real(8), intent(inout) :: xqz_VelYAl(imin:imax, jmin:jmax, kmin:kmax) 316 real(8), intent(inout) :: xyr_VelZAl(imin:imax, jmin:jmax, kmin:kmax) 317 real(8), intent(inout) :: xyz_PTempAl(imin:imax, jmin:jmax, kmin:kmax) 318 real(8), intent(inout) :: xyz_ExnerAl(imin:imax, jmin:jmax, kmin:kmax) 319 real(8) :: pyz_VelXAl0(imin:imax, jmin:jmax, kmin:kmax) 320 real(8) :: xqz_VelYAl0(imin:imax, jmin:jmax, kmin:kmax) 321 real(8) :: xyr_VelZAl0(imin:imax, jmin:jmax, kmin:kmax) 322 real(8) :: xyz_PTempAl0(imin:imax, jmin:jmax, kmin:kmax) 323 real(8) :: xyz_ExnerAl0(imin:imax, jmin:jmax, kmin:kmax) 324 real(8) :: pyz_DampVelX(imin:imax, jmin:jmax, kmin:kmax) 325 real(8) :: xqz_DampVelY(imin:imax, jmin:jmax, kmin:kmax) 326 real(8) :: xyr_DampVelZ(imin:imax, jmin:jmax, kmin:kmax) 327 real(8) :: xyz_DampPTemp(imin:imax, jmin:jmax, kmin:kmax) 328 real(8) :: xyz_DampExner(imin:imax, jmin:jmax, kmin:kmax) 329 330 pyz_VelXAl0 = pyz_VelXAl . !CDIR NODEP . do t539 = 1, pyz_velxal0.DSC.U1 + 1 - pyz_velxal0.DSC.L1 . pyz_velxal0(pyz_velxal0.DSC.L1+t539-1,t537+pyz_velxal0.DSC.L2, . 1 t535+pyz_velxal0.DSC.L3) = pyz_velxal(t184+t539-1,t537+t186, . 2 t535+t188) . pyz_dampvelx(pyz_dampvelx.DSC.L1+t539-1,t537+ . 1 pyz_dampvelx.DSC.L2,t535+pyz_dampvelx.DSC.L3) = -xyz_gamma( . 2 xyz_gamma.DSC.L1+t539-1,t537+xyz_gamma.DSC.L2,t535+ . 3 xyz_gamma.DSC.L3)*pyz_velxbl(t134+t539-1,t537+t136,t535+t138 . 4 ) . pyz_velxal(t184+t539-1,t537+t186,t535+t188) = pyz_velxal0( . 1 pyz_velxal0.DSC.L1+t539-1,t537+pyz_velxal0.DSC.L2,t535+ . 2 pyz_velxal0.DSC.L3) + 2.00000000000000e+000*deltimelong* . 3 pyz_dampvelx(pyz_dampvelx.DSC.L1+t539-1,t537+ . 4 pyz_dampvelx.DSC.L2,t535+pyz_dampvelx.DSC.L3) . xqz_velyal0(xqz_velyal0.DSC.L1+t539-1,t537+xqz_velyal0.DSC.L2, . 1 t535+xqz_velyal0.DSC.L3) = xqz_velyal(t194+t539-1,t537+t196, . 2 t535+t198) . xqz_dampvely(xqz_dampvely.DSC.L1+t539-1,t537+ . 1 xqz_dampvely.DSC.L2,t535+xqz_dampvely.DSC.L3) = -xqz_gamma( . 2 xqz_gamma.DSC.L1+t539-1,t537+xqz_gamma.DSC.L2,t535+ . 3 xqz_gamma.DSC.L3)*xqz_velybl(t144+t539-1,t537+t146,t535+t148 . 4 ) . xqz_velyal(t194+t539-1,t537+t196,t535+t198) = xqz_velyal0( . 1 xqz_velyal0.DSC.L1+t539-1,t537+xqz_velyal0.DSC.L2,t535+ . 2 xqz_velyal0.DSC.L3) + 2.00000000000000e+000*deltimelong* . 3 xqz_dampvely(xqz_dampvely.DSC.L1+t539-1,t537+ . 4 xqz_dampvely.DSC.L2,t535+xqz_dampvely.DSC.L3) . xyr_velzal0(xyr_velzal0.DSC.L1+t539-1,t537+xyr_velzal0.DSC.L2, . 1 t535+xyr_velzal0.DSC.L3) = xyr_velzal(t204+t539-1,t537+t206, . 2 t535+t208) . xyr_dampvelz(xyr_dampvelz.DSC.L1+t539-1,t537+ . 1 xyr_dampvelz.DSC.L2,t535+xyr_dampvelz.DSC.L3) = -xyr_gamma( . 2 xyr_gamma.DSC.L1+t539-1,t537+xyr_gamma.DSC.L2,t535+ . 3 xyr_gamma.DSC.L3)*xyr_velzbl(t154+t539-1,t537+t156,t535+t158 . 4 ) . xyr_velzal(t204+t539-1,t537+t206,t535+t208) = xyr_velzal0( . 1 xyr_velzal0.DSC.L1+t539-1,t537+xyr_velzal0.DSC.L2,t535+ . 2 xyr_velzal0.DSC.L3) + 2.00000000000000e+000*deltimelong* . 3 xyr_dampvelz(xyr_dampvelz.DSC.L1+t539-1,t537+ . 4 xyr_dampvelz.DSC.L2,t535+xyr_dampvelz.DSC.L3) . xyz_ptempal0(xyz_ptempal0.DSC.L1+t539-1,t537+ . 1 xyz_ptempal0.DSC.L2,t535+xyz_ptempal0.DSC.L3) = xyz_ptempal( . 2 t214+t539-1,t537+t216,t535+t218) . xyz_dampptemp(xyz_dampptemp.DSC.L1+t539-1,t537+ . 1 xyz_dampptemp.DSC.L2,t535+xyz_dampptemp.DSC.L3) = -xyz_gamma . 2 (xyz_gamma.DSC.L1+t539-1,t537+xyz_gamma.DSC.L2,t535+ . 3 xyz_gamma.DSC.L3)*xyz_ptempbl(t164+t539-1,t537+t166,t535+ . 4 t168) . xyz_ptempal(t214+t539-1,t537+t216,t535+t218) = xyz_ptempal0( . 1 xyz_ptempal0.DSC.L1+t539-1,t537+xyz_ptempal0.DSC.L2,t535+ . 2 xyz_ptempal0.DSC.L3) + 2.00000000000000e+000*deltimelong* . 3 xyz_dampptemp(xyz_dampptemp.DSC.L1+t539-1,t537+ . 4 xyz_dampptemp.DSC.L2,t535+xyz_dampptemp.DSC.L3) . xyz_exneral0(xyz_exneral0.DSC.L1+t539-1,t537+ . 1 xyz_exneral0.DSC.L2,t535+xyz_exneral0.DSC.L3) = xyz_exneral( . 2 t224+t539-1,t537+t226,t535+t228) . xyz_dampexner(xyz_dampexner.DSC.L1+t539-1,t537+ . 1 xyz_dampexner.DSC.L2,t535+xyz_dampexner.DSC.L3) = -xyz_gamma . 2 (xyz_gamma.DSC.L1+t539-1,t537+xyz_gamma.DSC.L2,t535+ . 3 xyz_gamma.DSC.L3)*xyz_exnerbl(t174+t539-1,t537+t176,t535+ . 4 t178) . xyz_exneral(t224+t539-1,t537+t226,t535+t228) = xyz_exneral0( . 1 xyz_exneral0.DSC.L1+t539-1,t537+xyz_exneral0.DSC.L2,t535+ . 2 xyz_exneral0.DSC.L3) + 2.00000000000000e+000*deltimelong* . 3 xyz_dampexner(xyz_dampexner.DSC.L1+t539-1,t537+ . 4 xyz_dampexner.DSC.L2,t535+xyz_dampexner.DSC.L3) . end do 331 pyz_DampVelX = - xyz_Gamma * pyz_VelXBl 332 pyz_VelXAl = pyz_VelXAl0 + (2.0d0 * DelTimeLong) * pyz_DampVelX 333 334 xqz_VelYAl0 = xqz_VelYAl 335 xqz_DampVelY = - xqz_Gamma * xqz_VelYBl 336 xqz_VelYAl = xqz_VelYAl0 + (2.0d0 * DelTimeLong) * xqz_DampVelY 337 338 xyr_VelZAl0 = xyr_VelZAl 339 xyr_DampVelZ = - xyr_Gamma * xyr_VelZBl 340 xyr_VelZAl = xyr_VelZAl0 + (2.0d0 * DelTimeLong) * xyr_DampVelZ 341 342 xyz_PTempAl0 = xyz_PTempAl 343 xyz_DampPTemp = - xyz_Gamma * xyz_PTempBl 344 xyz_PTempAl = xyz_PTempAl0 + (2.0d0 * DelTimeLong) * xyz_DampPTemp 345 346 xyz_ExnerAl0 = xyz_ExnerAl 347 xyz_DampExner = - xyz_Gamma * xyz_ExnerBl 348 xyz_ExnerAl = xyz_ExnerAl0 + (2.0d0 * DelTimeLong) * xyz_DampExner 349 350 call HistoryAutoPut(TimeN, 'VelXDamp', pyz_DampVelX(1:nx,1:ny,1:nz)) 351 call HistoryAutoPut(TimeN, 'VelYDamp', xqz_DampVelY(1:nx,1:ny,1:nz)) 352 call HistoryAutoPut(TimeN, 'VelZDamp', xyr_DampVelZ(1:nx,1:ny,1:nz)) 353 call HistoryAutoPut(TimeN, 'PTempDamp', xyz_DampPTemp(1:nx,1:ny,1:nz)) 354 call HistoryAutoPut(TimeN, 'ExnerDamp', xyz_DampExner(1:nx,1:ny,1:nz)) 355 356 end subroutine SpongeLayer_forcing 357 358 end module Damping Linux R2.6.5-7.282-sn2 FORTRAN90/SX Rev.360 Tue Oct 11 12:33:51 2011 FILE NAME: damping.f90 PROGRAM NAME: damping FORMAT LIST LINE LOOP FORTRAN STATEMENT 1: != Module Damping 2: ! 3: ! Authors:: SUGIYAMA Ko-ichiro, ODAKA Masatsugu 4: ! Version:: $Id: damping.f90,v 1.7 2011-10-05 03:35:23 sugiyama Exp $ 5: ! Tag Name:: $Name: arare5-20111010 $ 6: ! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved. 7: ! License:: See COPYRIGHT[link:../../COPYRIGHT] 8: ! 9: !== Overview 10: ! 11: !減衰率とその計算を行うためのパッケージ型モジュール 12: ! * 音波減衰項の係数 13: ! * スポンジ層の設定(境界付近で波の反射を抑え吸収するための層) 14: ! 15: !== Error Handling 16: ! 17: !== Bugs 18: ! 19: !== Note 20: ! 21: ! * この関数は, 基本場が零な変数(速度, エクスナー関数の擾乱)に 22: ! 適用することを想定している. 23: ! * 各格子点に対する関数を定義する必要がある 24: ! 25: !== Future Plans 26: ! 27: ! 28: 29: module Damping 30: ! 31: !減衰率とその計算を行うためのパッケージ型モジュール 32: ! * 音波減衰項の係数 33: ! * スポンジ層の設定(境界付近で波の反射を抑え吸収するための層) 34: ! 35: 36: !モジュール読み込み 37: use dc_types, only : DP 38: 39: !暗黙の型宣言禁止 40: implicit none 41: 42: !private 属性を指定 43: private 44: 45: !関数には public 属性を指定 46: public Damping_Init 47: public SpongeLayer_forcing 48: 49: !変数定義 50: real(DP), save :: EFTime = 100.0d0 !スポンジ層の e-folding time 51: real(DP), save :: DampDepthH = 0.0d0 !スポンジ層の厚さ(水平方向) 52: real(DP), save :: DampDepthV = 0.0d0 !スポンジ層の厚さ(鉛直方向) 53: real(DP), allocatable, save :: xyz_Gamma(:,:,:) !xyz 格子減衰係数(水平方向) 54: real(DP), allocatable, save :: pyz_Gamma(:,:,:) !pyz 格子減衰係数(鉛直方向) 55: real(DP), allocatable, save :: xqz_Gamma(:,:,:) !xqz 格子減衰係数(鉛直方向) 56: real(DP), allocatable, save :: xyr_Gamma(:,:,:) !xyr 格子減衰係数(鉛直方向) 57: 58: contains 59: 60: !!!------------------------------------------------------------------------!!! 61: subroutine Damping_Init 62: ! 63: ! 音波減衰項とスポンジ層の減衰係数の初期化 64: ! 65: use dc_iounit, only: FileOpen 66: use dc_message, only: MessageNotify 67: use gtool_historyauto, only: HistoryAutoAddVariable 68: use mpi_wrapper, only: myrank 69: use gridset, only: imin, &! x 方向の配列の下限 70: & imax, &! x 方向の配列の上限 71: & jmin, &! y 方向の配列の下限 72: & jmax, &! y 方向の配列の上限 73: & kmin, &! z 方向の配列の下限 74: & kmax, &! z 方向の配列の上限 75: & nx, &! x 方向の物理領域の上限 76: & ny, &! y 方向の物理領域の上限 77: & nz ! z 方向の物理領域の上限 78: use axesset, only: x_X, &!X 座標軸(スカラー格子点) 79: & y_Y, &!Y 座標軸(スカラー格子点) 80: & z_Z, &!Z 座標軸(スカラー格子点) 81: & p_X, &!X 座標軸(フラックス格子点) 82: & q_Y, &!Y 座標軸(フラックス格子点) 83: & r_Z, &!Z 座標軸(フラックス格子点) 84: & x_dx, y_dy, z_dz, &! 格子間隔 85: & XMax, &!X 座標の最大値 86: & YMax, &!Y 座標の最大値 87: & ZMax !Z 座標の最大値 88: use namelist_util, only: namelist_filename 89: 90: !暗黙の型宣言禁止 91: implicit none 92: 93: !変数定義 94: real(DP) :: Time ! 95: real(DP) :: DepthH !スポンジ層の厚さ(水平方向) 96: real(DP) :: DepthV !スポンジ層の厚さ(鉛直方向) 97: real(DP), parameter :: Pi =3.1415926535897932385d0 !円周率 98: integer :: unit ! 装置番号 99: integer :: i, j, k 100: 101: !NAMELIST から取得 102: NAMELIST /damping_nml/ Time, DepthH, DepthV 103: 104: call FileOpen(unit, file=namelist_filename, mode='r') 105: read(unit, NML=damping_nml) 106: close(unit) 107: 108: !初期化 109: allocate( & 110: & xyz_Gamma(imin:imax,jmin:jmax,kmin:kmax), & 111: & pyz_Gamma(imin:imax,jmin:jmax,kmin:kmax), & 112: & xqz_Gamma(imin:imax,jmin:jmax,kmin:kmax), & 113: & xyr_Gamma(imin:imax,jmin:jmax,kmin:kmax) ) 114: WWW==== xyz_Gamma = 0.0d0 115: W**==== pyz_Gamma = 0.0d0 116: ***==== xqz_Gamma = 0.0d0 117: ***==== xyr_Gamma = 0.0d0 118: 119: !値の入力 120: EFTime = Time 121: DampDepthH = DepthH 122: DampDepthV = DepthV 123: 124: !----------------------------------------------------------------- 125: ! スポンジ層の減衰率 126: ! 127: !水平方向の東側・西側境界 128: if ( DampDepthH < x_dx(1) ) then 129: if (myrank == 0) & 130: & call MessageNotify( "W", "Damping_init", "DampDepthH is too thin. DelX is %f", d=(/x_dx(1)/)) 131: 132: else if ( DampDepthH < x_dx(nx) ) then 133: if (myrank == 0) & 134: & call MessageNotify( "W", "Damping_init", "DampDepthH is too thin. DelX is %f", d=(/x_dx(nx)/)) 135: 136: else 137: +------> do i = imin, imax 138: | !スカラー格子点の西側境界 139: | if ( x_X(i) < DampDepthH) then 140: |WW==== xyz_Gamma(i,:,:) = ((1.0d0 - x_X(i) / DampDepthH) ** 3.0d0) / EFTime 141: | end if 142: | 143: | !フラックス格子点の西側境界 144: | if ( p_X(i) < DampDepthH) then 145: |WW==== pyz_Gamma(i,:,:) = ((1.0d0 - p_X(i) / DampDepthH) ** 3.0d0) / EFTime 146: | end if 147: | 148: | !スカラー格子点の東側境界 149: | if ( x_X(i) > ( XMax - DampDepthH ) ) then 150: |**==== xyz_Gamma(i,:,:) = & 151: | & ((1.0d0 - (XMax - x_X(i)) / DampDepthH) ** 3.0d0) / EFTime 152: | end if 153: | 154: | !フラックス格子点の東側境界 155: | if ( p_X(i) > ( XMax - DampDepthH ) ) then 156: |**==== pyz_Gamma(i,:,:) = & 157: | & ((1.0d0 - (XMax - p_X(i)) / DampDepthH) ** 3.0d0) / EFTime 158: | end if 159: +------ end do 160: end if 161: 162: ! x 方向には同じ 163: ! 164: **V----> xyr_Gamma = xyz_Gamma 165: **V---- xqz_Gamma = xyz_Gamma 166: 167: 168: !水平方向の南側・北側境界 169: if ( DampDepthH < y_dy(1) ) then 170: if (myrank == 0) & 171: & call MessageNotify( "W", "Damping_init", "DampDepthH is too thin. DelY is %f", d=(/x_dx(1)/)) 172: 173: else if ( DampDepthH < y_dy(ny) ) then 174: if (myrank == 0) & 175: & call MessageNotify( "W", "Damping_init", "DampDepthH is too thin. DelY is %f", d=(/y_dy(ny)/)) 176: 177: else 178: +------> do j = jmin, jmax 179: | !スカラー格子点の西側境界 180: | if ( y_Y(j) < DampDepthH) then 181: |+V==== xyz_Gamma(:,j,:) = ((1.0d0 - y_Y(j) / DampDepthH) ** 3.0d0) / EFTime 182: | end if 183: | 184: | !フラックス格子点の西側境界 185: | if ( q_Y(j) < DampDepthH) then 186: |+V==== xqz_Gamma(:,j,:) = ((1.0d0 - q_Y(j) / DampDepthH) ** 3.0d0) / EFTime 187: | end if 188: | 189: | !スカラー格子点の東側境界 190: | if ( y_Y(j) > ( YMax - DampDepthH ) ) then 191: |+V==== xyz_Gamma(:,j,:) = & 192: | & ((1.0d0 - (YMax - y_Y(j)) / DampDepthH) ** 3.0d0) / EFTime 193: | end if 194: | 195: | !フラックス格子点の東側境界 196: | if ( q_Y(j) > ( YMax - DampDepthH ) ) then 197: |+V==== xqz_Gamma(:,j,:) = & 198: | & ((1.0d0 - (YMax - q_Y(j)) / DampDepthH) ** 3.0d0) / EFTime 199: | end if 200: +------ end do 201: end if 202: 203: ! y 方向には同じ 204: ! 205: **V----> pyz_Gamma = xyz_Gamma 206: **V---- xyr_Gamma = xyz_Gamma 207: 208: 209: !鉛直方向の上部境界 210: if ( DampDepthV < z_dz(nz) ) then 211: if (myrank == 0) & 212: & call MessageNotify( "W", "Damping_init", "DampDepthV is too thin. DelZ is %f", d=(/z_dz(nz)/) ) 213: 214: else 215: +------> do k = kmin, kmax 216: | !スカラー格子点 217: | if ( z_Z(k) >= ( ZMax - DampDepthV ) ) then 218: |WW==== xyz_Gamma(:,:,k) = & 219: | & (1.0d0 - dcos(Pi * (z_Z(k) - ZMax + DampDepthV) / DampDepthV)) & 220: | & / EFTime 221: | end if 222: | 223: | !フラックス格子点 224: | if ( r_Z(k) >= ( ZMax - DampDepthV ) ) then 225: |**==== xyr_Gamma(:,:,k) = & 226: | & (1.0d0 - dcos(Pi * (r_Z(k) - ZMax + DampDepthV)/ DampDepthV)) & 227: | & / EFTime 228: | end if 229: +------ end do 230: end if 231: 232: ! z 方向には同じ 233: ! 234: **V----> pyz_Gamma = xyz_Gamma 235: **V---- xqz_Gamma = xyz_Gamma 236: 237: 238: !----------------------------------------------------------------- 239: ! 値の確認 240: ! 241: if (myrank == 0) then 242: call MessageNotify( "M", "Damping_init", "EFTime = %f", d=(/EFTime/) ) 243: call MessageNotify( "M", "Damping_init", "DampDepthH = %f", d=(/DampDepthH/) ) 244: call MessageNotify( "M", "Damping_init", "DampDepthV = %f", d=(/DampDepthV/) ) 245: end if 246: 247: 248: !----------------------------------------------------------------- 249: ! 出力 250: ! 251: call HistoryAutoAddVariable( & 252: & varname='PTempDamp', & 253: & dims=(/'x','y','z','t'/), & 254: & longname='Damping term of potential temperature', & 255: & units='K.s-1', & 256: & xtype='float') 257: 258: call HistoryAutoAddVariable( & 259: & varname='ExnerDamp', & 260: & dims=(/'x','y','z','t'/), & 261: & longname='Damping term of exner function', & 262: & units='s-1', & 263: & xtype='float') 264: 265: call HistoryAutoAddVariable( & 266: & varname='VelXDamp', & 267: & dims=(/'x','y','z','t'/), & 268: & longname='Damping term of VelX', & 269: & units='m.s-1', & 270: & xtype='float') 271: 272: call HistoryAutoAddVariable( & 273: & varname='VelYDamp', & 274: & dims=(/'x','y','z','t'/), & 275: & longname='Damping term of VelY', & 276: & units='m.s-1', & 277: & xtype='float') 278: 279: call HistoryAutoAddVariable( & 280: & varname='VelZDamp', & 281: & dims=(/'x','y','z','t'/), & 282: & longname='Damping term of VelZ', & 283: & units='m.s-1', & 284: & xtype='float') 285: 286: end subroutine Damping_Init 287: 288: 289: 290: subroutine SpongeLayer_forcing( & 291: & pyz_VelXBl, xqz_VelYBl, xyr_VelZBl, xyz_PTempBl, xyz_ExnerBl, & !(in) 292: & pyz_VelXAl, xqz_VelYAl, xyr_VelZAl, xyz_PTempAl, xyz_ExnerAl ) !(inout) 293: 294: use gtool_historyauto, only: HistoryAutoPut 295: use timeset, only: TimeN, DelTimeLong 296: use gridset, only: imin, &! x 方向の配列の下限 297: & imax, &! x 方向の配列の上限 298: & jmin, &! y 方向の配列の下限 299: & jmax, &! y 方向の配列の上限 300: & kmin, &! z 方向の配列の下限 301: & kmax, &! z 方向の配列の上限 302: & nx, &! x 方向の物理領域の上限 303: & ny, &! y 方向の物理領域の上限 304: & nz ! z 方向の物理領域の上限 305: 306: !暗黙の型宣言禁止 307: implicit none 308: 309: real(8), intent(in) :: pyz_VelXBl(imin:imax, jmin:jmax, kmin:kmax) 310: real(8), intent(in) :: xqz_VelYBl(imin:imax, jmin:jmax, kmin:kmax) 311: real(8), intent(in) :: xyr_VelZBl(imin:imax, jmin:jmax, kmin:kmax) 312: real(8), intent(in) :: xyz_PTempBl(imin:imax, jmin:jmax, kmin:kmax) 313: real(8), intent(in) :: xyz_ExnerBl(imin:imax, jmin:jmax, kmin:kmax) 314: real(8), intent(inout) :: pyz_VelXAl(imin:imax, jmin:jmax, kmin:kmax) 315: real(8), intent(inout) :: xqz_VelYAl(imin:imax, jmin:jmax, kmin:kmax) 316: real(8), intent(inout) :: xyr_VelZAl(imin:imax, jmin:jmax, kmin:kmax) 317: real(8), intent(inout) :: xyz_PTempAl(imin:imax, jmin:jmax, kmin:kmax) 318: real(8), intent(inout) :: xyz_ExnerAl(imin:imax, jmin:jmax, kmin:kmax) 319: real(8) :: pyz_VelXAl0(imin:imax, jmin:jmax, kmin:kmax) 320: real(8) :: xqz_VelYAl0(imin:imax, jmin:jmax, kmin:kmax) 321: real(8) :: xyr_VelZAl0(imin:imax, jmin:jmax, kmin:kmax) 322: real(8) :: xyz_PTempAl0(imin:imax, jmin:jmax, kmin:kmax) 323: real(8) :: xyz_ExnerAl0(imin:imax, jmin:jmax, kmin:kmax) 324: real(8) :: pyz_DampVelX(imin:imax, jmin:jmax, kmin:kmax) 325: real(8) :: xqz_DampVelY(imin:imax, jmin:jmax, kmin:kmax) 326: real(8) :: xyr_DampVelZ(imin:imax, jmin:jmax, kmin:kmax) 327: real(8) :: xyz_DampPTemp(imin:imax, jmin:jmax, kmin:kmax) 328: real(8) :: xyz_DampExner(imin:imax, jmin:jmax, kmin:kmax) 329: 330: **V----> pyz_VelXAl0 = pyz_VelXAl 331: ||| pyz_DampVelX = - xyz_Gamma * pyz_VelXBl 332: ||| pyz_VelXAl = pyz_VelXAl0 + (2.0d0 * DelTimeLong) * pyz_DampVelX 333: ||| 334: ||| xqz_VelYAl0 = xqz_VelYAl 335: ||| xqz_DampVelY = - xqz_Gamma * xqz_VelYBl 336: ||| xqz_VelYAl = xqz_VelYAl0 + (2.0d0 * DelTimeLong) * xqz_DampVelY 337: ||| 338: ||| xyr_VelZAl0 = xyr_VelZAl 339: ||| xyr_DampVelZ = - xyr_Gamma * xyr_VelZBl 340: ||| xyr_VelZAl = xyr_VelZAl0 + (2.0d0 * DelTimeLong) * xyr_DampVelZ 341: ||| 342: ||| xyz_PTempAl0 = xyz_PTempAl 343: ||| xyz_DampPTemp = - xyz_Gamma * xyz_PTempBl 344: ||| xyz_PTempAl = xyz_PTempAl0 + (2.0d0 * DelTimeLong) * xyz_DampPTemp 345: ||| 346: ||| xyz_ExnerAl0 = xyz_ExnerAl 347: ||| xyz_DampExner = - xyz_Gamma * xyz_ExnerBl 348: **V---- xyz_ExnerAl = xyz_ExnerAl0 + (2.0d0 * DelTimeLong) * xyz_DampExner 349: 350: call HistoryAutoPut(TimeN, 'VelXDamp', pyz_DampVelX(1:nx,1:ny,1:nz)) 351: call HistoryAutoPut(TimeN, 'VelYDamp', xqz_DampVelY(1:nx,1:ny,1:nz)) 352: call HistoryAutoPut(TimeN, 'VelZDamp', xyr_DampVelZ(1:nx,1:ny,1:nz)) 353: call HistoryAutoPut(TimeN, 'PTempDamp', xyz_DampPTemp(1:nx,1:ny,1:nz)) 354: call HistoryAutoPut(TimeN, 'ExnerDamp', xyz_DampExner(1:nx,1:ny,1:nz)) 355: 356: end subroutine SpongeLayer_forcing 357: 358: end module Damping