: 4. 参考文献 : 2 次元非静力学モデルの離散化 : 2. 空間方向の離散化

3. 時間方向の離散化

3.1 運動方程式と圧力方程式

空間離散化された運動方程式(2.28), (2.29)と圧力方程式 (2.30)を時間方向に離散化する. 音波に関連する項は短いタイムステップ $\Delta \tau$ で離散化し, その他の項は長いタイムステップ $\Delta t$ で離散化する. 音波に関連する項の離散化には HE-VI 法を採用し, の式は前進差分, $w, \pi$ の式は後退差分 (クランク・ニコルソン法)で離散化する. その他の項の離散化にはリープフロッグ法を用いる. 離散化した式の計算はまずの式から行う. 得られた $\tau +\Delta \tau$ のを用いて $\pi$ を計算し, $u, \pi$ を用いてを計算する.

運動方程式の各項のうち, 音波に関係しない項をとしてまとめると, 運動方程式と圧力方程式は以下のように書ける.

		$\displaystyle \DP{u_{i(u),k}}{t} = - \left[\bar{c_{p}} \bar{\theta}_{v} \DP{(\pi - \alpha Div )}{x}\right]_{i(u),k} + [F_{u}]_{i(u),k}^{t},$	(3.1)
		$\displaystyle \DP{w_{i,k(w)}}{t} = - \left[\bar{c_{p}} \bar{\theta}_{v} \DP{(\pi - \alpha Div )}{z}\right]_{i,k(w)} + [F_{w}]_{i,k(w)}^{t},$	(3.2)
		$\displaystyle \DP{\pi_{i,k}}{t} + \left[\frac{\bar{c}^{2}}{\bar{c_{p}} \bar{\rh... ... \left[\frac{\bar{c}^{2}}{\bar{c_{p}} \bar{\theta}_{v}} \DP{u}{x}\right]_{i,k}.$	(3.3)

ただし

の式には音波減衰項 $\alpha Div$ を加えてある (Skamarock and Klemp, 1992). 音波に関連しない項 $F_{u}, F_{w}$ は,

		$\displaystyle [F_{u}]_{i(u),k}^{t} = - \left[{\rm Adv}.{u}\right]_{i(u),k}^{t} ... ...ht]_{i(u),k}^{t - \Delta t} + \left[{\rm Diff}.u\right]_{i(u),k}^{t - \Delta t}$	(3.4)
		$\displaystyle \left[F_{w}\right]_{i,k(w)}^{t} = - \left[{\rm Adv}.{w}\right]_{i... ..._{i,k(w)}^{t - \Delta t}. + \left[{\rm Diff}.w \right]_{i,k(w)}^{t - \Delta t}.$	(3.5)

であり, 時刻

で評価することにする. 但し, 中心差分でリープフロッグ法を用いるため, 数値粘性項 Diff を追加してある.

3.1.1 音波に関連する項の時間方向の離散化

3.1.1.1 水平方向の運動方程式の離散化

(3.1)を時間方向に離散化すると以下のようになる.

$\displaystyle u^{\tau + \Delta \tau}_{i(u),k} = u^{\tau}_{i(u), k} - \left[ \ba... ...\alpha Div)^{\tau}}{x} \right\} \right]_{i(u),k} + F_{u,i(u),k}^{t} \Delta \tau$

(3.6)

3.1.1.2 鉛直方向の運動方程式と圧力方程式の離散化

HE-VI 法を用いるので, と $\pi$ の式を連立して解く. の式において音波減衰項は前進差分, 圧力項は後退差分で離散化する. $\pi$ の式において水平微分項は(3.6)で求めた $u^{\tau +\Delta \tau}$ を用いて離散化し, 鉛直微分項は後退差分で離散化する.

$\displaystyle w^{\tau + \Delta \tau} = w^{\tau} - \bar{c_{p}} \bar{\theta}_{v} ... ...{\pi^{\tau}}{z} - \DP{(\alpha Div)^{\tau}}{z} \right\} + F_{w}^{t} \Delta \tau.$

(3.7)

$\displaystyle \pi^{\tau + \Delta \tau} + \beta \frac{\bar{c}^{2}\Delta \tau}{\b... ...^{2} \Delta \tau}{\bar{c_{p}} \bar{\theta}_{v}} \DP{u^{\tau + \Delta \tau}}{x}.$

ここでは簡単のため格子点位置を表す添字は省略した. (3.8) 式に (3.7) を代入して $w^{\tau + \Delta \tau}$ を消去する.

$\displaystyle \pi^{\tau + \Delta \tau}$	$\textstyle -$	$\displaystyle \beta^{2} \frac{\bar{c}^{2}{\Delta \tau}^{2}}{\bar{c_{p}} \bar{\r... ...\theta}_{v}^{2}\right) \left( \DP{\pi^{\tau + \Delta \tau}}{z} \right) \right\}$
	$\textstyle =$	$\displaystyle \pi^{\tau} -(1 - \beta) \frac{\bar{c}^{2}\Delta \tau}{\bar{c_{p}}... ...}^{2} \Delta \tau}{\bar{c_{p}} \bar{\theta}_{v}} \DP{u^{\tau + \Delta \tau}}{x}$
		$\displaystyle - \beta \frac{\bar{c}^{2}\Delta \tau}{\bar{c_{p}} \bar{\rho} \bar... ... \DP{(\alpha Div)^{\tau}}{z} \right\} + F_{w}^{t} \Delta \tau \right\} \right].$

(3.9) 式右辺を空間方向に離散化し, 格子点位置を表す添字を付けて表すと以下のようになる (計算の詳細は第A章参照).

		$\displaystyle \left\{ - \beta^{2} \left( \frac{\bar{c}^{2}{\Delta \tau}^{2}}{\b... ...} \bar{\theta}_{v}^{2} \right)_{k(w)} \right\} \pi^{\tau + \Delta \tau}_{i,k+1}$
		$\displaystyle \hspace{10mm}+ \left[ 1 + \beta^{2} \left( \frac{\bar{c}^{2}{\Del... ...theta}_{v}^{2} \right)_{k-1(w)} \right\} \right] \pi^{\tau + \Delta \tau}_{i,k}$
		$\displaystyle \hspace{10mm}+ \left\{ - \beta^{2} \left( \frac{\bar{c}^{2}{\Delt... ...\bar{\theta}_{v}^{2} \right)_{k-1(w)} \right\} \pi^{\tau + \Delta \tau}_{i,k-1}$
		$\displaystyle = \pi^{\tau}_{i,k} - (1 - \beta) \left( \frac{\bar{c}^{2}\Delta \... ...ar{\theta}_{v}} \right)_{k} \left( \DP{u^{\tau + \Delta \tau}}{x} \right)_{i,k}$
		$\displaystyle \hspace{2mm} - \beta \left( \frac{\bar{c}^{2}\Delta \tau}{\bar{c_... ...ho} \bar{\theta}_{v} \right)_{i,k(w)} \left\{ w^{\tau}_{i,k(w)} \right. \right.$
		$\displaystyle \hspace{10mm} \left. \left. - \left( \bar{c_{p}} \bar{\theta}_{v}... ...i,k(w)} + \left( F_{w}^{t} \right)_{i,k(w)} \Delta \tau \right\} \right]_{i,k}.$

但し平均場の量は鉛直方向にしか依存しないので

方向の添字のみ付けてある.

3.1.1.3 境界条件

上下境界を固定壁とする場合, 境界条件は上部下部境界で,

		$\displaystyle w(i,0(w)) = 0,$	(3.8)
		$\displaystyle w(i,km(w)) = 0$	(3.9)

である.

下部境界:

下部境界()について考える. この時 (3.7) 式に添字を付けて書き下すと,

$\displaystyle \beta \left( \DP{\pi^{\tau + \Delta \tau}}{z} \right)_{i,0(w)}$	$\textstyle =$	$\displaystyle \left( \DP{(\alpha Div)^{\tau}}{z} \right)_{i,0(w)} - (1 - \beta)... ..._{i,0(w)} + \left(\Dinv{\bar{c_{p}} \bar{\theta}_{v}} F_{w}^{t}\right)_{i,0(w)}$
	$\textstyle \equiv$	$\displaystyle E_{i,0(w)}$	(3.10)

となる. したがって (3.10) 式は以下のようになる.

	$\displaystyle \left\{ - \beta^{2} \left( \frac{\bar{c}^{2}{\Delta \tau}^{2}}{\b... ...ho} \bar{\theta}_{v}^{2} \right)_{1(w)} \right\} \pi^{\tau + \Delta \tau}_{i,1}$
$\textstyle =$	$\displaystyle \pi^{\tau}_{i,1} -(1 - \beta) \left( \frac{\bar{c}^{2}\Delta \tau... ...ar{\theta}_{v}} \right)_{1} \left( \DP{u^{\tau + \Delta \tau}}{x} \right)_{i,1}$
	$\displaystyle - \beta \left( \frac{\bar{c}^{2}\Delta \tau}{\bar{c_{p}} \bar{\rh... ...(\alpha Div)^{\tau}}{z} \right\} + F_{w}^{t} \Delta \tau \right\} \right]_{i,1}$
	$\displaystyle - \beta \left( \frac{\bar{c}^{2}{\Delta \tau}^{2}}{\bar{c_{p}} \b... ...\left( \bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2} \right)_{i,0(w)} E_{i,0(w)}.$	(3.11)

上部境界:

上部境界()について考える. この時 (3.7) 式を添字を付けて書き下すと,

$\displaystyle \beta \left( \DP{\pi^{\tau + \Delta \tau}}{z} \right)_{i,km(w)}$	$\textstyle =$	$\displaystyle \left( \DP{(\alpha Div)^{\tau}}{z} \right)_{i,km(w)} - (1 - \beta... ...i,km(w)} + \left(\Dinv{\bar{c_{p}} \bar{\theta}_{v}} F_{w}^{t}\right)_{i,km(w)}$
	$\textstyle \equiv$	$\displaystyle E_{i,km(w)}$	(3.12)

となる. したがって (3.10) 式は以下のようになる.

	$\displaystyle \left\{ 1 + \beta^{2} \left( \frac{\bar{c}^{2}{\Delta \tau}^{2}}{... ...\bar{\theta}_{v}^{2} \right)_{km-1(w)} \right\} \pi^{\tau + \Delta \tau}_{i,km}$
	$\displaystyle + \left\{ - \beta^{2} \left( \frac{\bar{c}^{2}{\Delta \tau}^{2}}{... ...ar{\theta}_{v}^{2} \right)_{km-1(w)} \right\} \pi^{\tau + \Delta \tau}_{i,km-1}$
$\textstyle =$	$\displaystyle \pi^{\tau}_{i,km} -(1 - \beta) \left( \frac{\bar{c}^{2}\Delta \ta... ...{\theta}_{v}} \right)_{km} \left( \DP{u^{\tau + \Delta \tau}}{x} \right)_{i,km}$
	$\displaystyle - \beta \left( \frac{\bar{c}^{2}\Delta \tau}{\bar{c_{p}} \bar{\rh... ...\alpha Div)^{\tau}}{z} \right\} + F_{w}^{t} \Delta \tau \right\} \right]_{i,km}$
	$\displaystyle + \frac{\beta}{\Delta z} \left( \frac{\bar{c}^{2}{\Delta \tau}^{2... ...\left( \bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2} \right)_{km(w)} E_{i,km(w)}.$	(3.13)

3.1.1.4 圧力方程式の時間積分方法

(3.10), (3.14), (3.16) 式を連立すると, 以下のような行列式の形式で書くことができる.

$\displaystyle \left(\begin{array}{cccc} A_{1} & B_{2} & & 0 \\ C_{1} & \ddots ... ...m} & \cdots & \cdots & \pi_{im, km} \\ \end{array}\right)^{\tau + \Delta \tau}$
$\displaystyle = \left(\begin{array}{cccc} D_{1,1} & D_{2,1} & \cdots & D_{im,1}... ...dots \\ D_{1,km} & \cdots & \cdots & D_{im,km} \\ \end{array}\right)^{\tau} .$			(3.14)

この連立方程式を解くことで $\pi_{i, k}$ を求める. この連立方程式の係数は以下のように書ける.

$\displaystyle A_{k}$	$\textstyle =$	$\displaystyle 1 + \beta^{2} \left( \frac{\bar{c}^{2}{\Delta \tau}^{2}}{\bar{c_{... ... + \left( \bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2} \right)_{k-1(w)} \right\}$
		$\displaystyle (k = 2, 3, \cdots km-1),$
$\displaystyle A_{1}$	$\textstyle =$	$\displaystyle 1 + \beta^{2} \left( \frac{\bar{c}^{2}{\Delta \tau}^{2}}{\bar{c_{... ...Delta z^{2}} \left( \bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2} \right)_{1(w)},$
$\displaystyle A_{km}$	$\textstyle =$	$\displaystyle 1 + \beta^{2} \left( \frac{\bar{c}^{2}{\Delta \tau}^{2}}{\bar{c_{... ...ta z^{2}} \left( \bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2} \right)_{km-1(w)},$
$\displaystyle B_{k}$	$\textstyle =$	$\displaystyle - \beta^{2} \left( \frac{\bar{c}^{2}{\Delta \tau}^{2}}{\bar{c_{p}... ...lta z^{2}} \left( \bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2} \right)_{k-1(w)},$
		$\displaystyle (k = 2, 3, \cdots km),$
$\displaystyle C_{k}$	$\textstyle =$	$\displaystyle - \beta^{2} \left( \frac{\bar{c}^{2}{\Delta \tau}^{2}}{\bar{c_{p}... ...Delta z^{2}} \left( \bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2} \right)_{k(w)},$
		$\displaystyle (k = 1, 2, \cdots km-1),$
$\displaystyle D_{i,k}$	$\textstyle =$	$\displaystyle \pi^{\tau}_{i,k} -(1 - \beta) \left( \frac{\bar{c}^{2}\Delta \tau... ..._{v}} \right)_{k} \left( \DP{u^{\tau + \Delta \tau}}{x} \right)_{i,k} + F_{i,k}$
		$\displaystyle (k = 2, 3, \cdots km-1),$
$\displaystyle D_{i,1}$	$\textstyle =$	$\displaystyle \pi^{\tau}_{i,1} -(1 - \beta) \left( \frac{\bar{c}^{2}\Delta \tau... ..._{v}} \right)_{1} \left( \DP{u^{\tau + \Delta \tau}}{x} \right)_{i,1} + F_{i,1}$
		$\displaystyle - \beta \left( \frac{\bar{c}^{2}{\Delta \tau}^{2}}{\bar{c_{p}} \b... ...\left( \bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2} \right)_{i,0(w)} E_{i,0(w)},$
$\displaystyle D_{i,km}$	$\textstyle =$	$\displaystyle \pi^{\tau}_{i,km} -(1 - \beta) \left( \frac{\bar{c}^{2}\Delta \ta... ...}} \right)_{km} \left( \DP{u^{\tau + \Delta \tau}}{x} \right)_{i,km} + F_{i,km}$
		$\displaystyle + \beta \left( \frac{\bar{c}^{2}{\Delta \tau}^{2}}{\bar{c_{p}} \b... ...\left( \bar{c_{p}} \bar{\rho} \bar{\theta}_{v}^{2} \right)_{km(w)} E_{i,km(w)}.$

ただし,

$\displaystyle E_{i,k(w)} \equiv \left( \DP{(\alpha Div)^{\tau}}{z} \right)_{i,k... ..._{i,k(w)} + \left(\Dinv{\bar{c_{p}} \bar{\theta}_{v}} F_{w}^{t}\right)_{i,k(w)}$

$\displaystyle F_{i,k}$	$\textstyle \equiv$	$\displaystyle \hspace{2mm} - \beta \left( \frac{\bar{c}^{2}\Delta \tau}{\bar{c_... ...ho} \bar{\theta}_{v} \right)_{i,k(w)} \left\{ w^{\tau}_{i,k(w)} \right. \right.$
		$\displaystyle \hspace{10mm} \left. \left. - \left( \bar{c_{p}} \bar{\theta}_{v}... ...i,k(w)} + \left( F_{w}^{t} \right)_{i,k(w)} \Delta \tau \right\} \right]_{i,k}.$

である.

3.1.2 音波に関連しない項の時間方向の離散化

運動方程式の音波に関連しない項 (3.1), (3.2) 式を離散化する.

		$\displaystyle F_{u,i(u),k}^{t} = - \left[ {\rm Adv}.u \right]_{i(u),k}^{t} + \l... ...ight]_{i(u),k}^{t-\Delta t} + \left[{\rm Diff}.u \right]_{i(u),k}^{t-\Delta t},$	(3.15)
		$\displaystyle F_{w,i,k(w)} = - \left[ {\rm Adv}.w \right]_{i(u),k}^{t} + \left[... ..._{i,k(w)}^{t - \Delta t} + \left[ {\rm Diff}.w \right]_{i,k(w)}^{t - \Delta t}.$	(3.16)

ここで, Adv は移流項, D は粘性拡散項, Buoy は浮力項, Diff は数値粘性項である. それぞれの項を書き下すと,

$\displaystyle \left[ {\rm Adv}.{u} \right]_{i(u),k}^{t}$	$\textstyle =$	$\displaystyle u_{i(u),k}^{t} \left[\DP{u}{x}\right]_{i(u),k}^{t} + w_{i(u),k}^{t} \left[\DP{u}{z}\right]_{i(u),k}^{t}$	(3.17)
$\displaystyle \left[ {\rm Adv}.{w} \right]_{i,k(w)}^{t}$	$\textstyle =$	$\displaystyle u_{i,k(w)}\left[\DP{w}{x}\right]_{i,k(w)}^{t} + w_{i,k(w)}\left[\DP{w}{z}\right]_{i,k(w)}^{t}$	(3.18)

であり, 浮力項は,

$\displaystyle [{\rm Buoy}]^{t}_{i,k(w)}$	$\textstyle =$	$\displaystyle g \frac{\theta_{i,k(w)}^{t}}{\overline{\theta}_{i,k(w)}}$
		$\displaystyle + g \frac{\sum [q_{v}]_{i,k(w)}^{t}/M_{v}}{1/M_{d} + \sum [\bar{q_{v}}]_{i,k(w)}/M_{v}}$
		$\displaystyle - g \frac{\sum [q_{v}]_{i,k(w)}^{t} + \sum [q_{c}]_{i,k(w)}^{t} + \sum [q_{r}]_{i,k(w)}^{t}} {1 + \sum [\bar{q_{v}}]_{i,k(w)}}$	(3.19)

であり, 粘性拡散項は,

$\displaystyle \left[ {\rm Turb}.{u} \right]_{i(u),k}^{t - \Delta t}$	$\textstyle =$	$\displaystyle 2 \left[ \DP{}{x}\left\{ \left( K_{m} \right)_{i,k} \left( \DP{u}{x} \right)_{i,k} \right\} \right]_{i(u),k}^{t - \Delta t}$
		$\displaystyle +\left[ \DP{}{z}\left\{ \left( K_{m} \right)_{i(u),k(w)} \left( \... ...} \left( \DP{u}{z} \right)_{i(u),k(w)} \right\} \right]_{i(u),k}^{t - \Delta t}$
		$\displaystyle - \frac{2}{3 C_{m}^{2} l^{2}} \left( \DP{ K_{m}^{2} }{x} \right)_{i(u),k}^{t - \Delta t}$	(3.20)
$\displaystyle \left[ {\rm Turb}.{w} \right]_{i,k(w)}^{t - \Delta t}$	$\textstyle =$	$\displaystyle 2 \left[ \DP{}{z}\left\{ \left( K_{m} \right)_{i,k} \left( \DP{w}{z} \right)_{i,k} \right\} \right]_{i,k(w)}^{t - \Delta t}$
		$\displaystyle +\left[ \DP{}{x}\left\{ \left( K_{m} \right)_{i(u),k(w)} \left( \... ...} \left( \DP{u}{z} \right)_{i(u),k(w)} \right\} \right]_{i,k(w)}^{t - \Delta t}$
		$\displaystyle - \frac{2}{3 C_{m}^{2} l^{2}} \left( \DP{ K_{m}^{2} }{z} \right)_{i,k(w)}^{t - \Delta t}$	(3.21)

である. 数値粘性項は,

$\displaystyle \left[ {\rm Diff}.u \right]_{i(u),k}^{t - \Delta t}$	$\textstyle =$	$\displaystyle \nu_{h} \left\{ \DP{}{x} \left(\DP{u}{x}\right)_{i,k} \right\}_{i... ...\{ \DP{}{z} \left(\DP{u}{z}\right)_{i(u),k(w)} \right\}_{i(u),k}^{t - \Delta t}$	(3.22)
$\displaystyle \left[ {\rm Diff}.w \right]_{i,k(w)}^{t - \Delta t}$	$\textstyle =$	$\displaystyle \nu_{h} \left\{ \DP{}{x} \left(\DP{w}{x}\right)_{i(u),k(w)} \righ... ... \left\{ \DP{}{z} \left(\DP{w}{z}\right)_{i,k} \right\}_{i,k(w)}^{t - \Delta t}$	(3.23)

である. $K_{m}$ は乱流エネルギーの時間発展方程式から計算し(詳細は後述), $\nu_{h}, \nu_{v}$ は以下のように定める.

$\displaystyle \nu_{h} = \frac{\alpha_{h} \Delta x^{2}}{\Delta t}$			(3.24)
$\displaystyle \nu_{v} = \frac{\alpha_{v} \Delta z^{2}}{\Delta t}$			(3.25)

ここで $\Delta x, \Delta z$ は水平・鉛直方向の格子間隔を意味し, $\alpha_{h}, \alpha_{v}$ はそれぞれ,

$\displaystyle \alpha_{h} \le \Dinv{8}, \hspace{3em} \alpha_{v} \le \Dinv{8}$

(3.26)

とする.

3.2 熱力学の式と混合比の保存式の離散化

熱の式と混合比の保存式の右辺をまとめてで表し, 時間方向にリープフロッグ法を用いて離散化する.

$\displaystyle \theta_{i,k}^{t + \Delta t}$	$\textstyle =$	$\displaystyle \theta_{i,k}^{t - \Delta t} + 2 \Delta t [F_{\theta}]_{i,k}^{t}$	(3.27)
$\displaystyle \left[ q_{v} \right]_{i,k}^{t+\Delta t}$	$\textstyle =$	$\displaystyle \left[ q_{v} \right]_{i,k}^{t-\Delta t} + 2 \Delta t [F_{q_{v}}]_{i,k}^{t}$	(3.28)
$\displaystyle \left[ q_{c} \right]_{i,k}^{t+\Delta t}$	$\textstyle =$	$\displaystyle \left[ q_{c} \right]_{i,k}^{t-\Delta t} + 2 \Delta t [F_{q_{c}}]_{i,k}^{t}$	(3.29)
$\displaystyle \left[ q_{r} \right]_{i,k}^{t+\Delta t}$	$\textstyle =$	$\displaystyle \left[ q_{r} \right]_{i,k}^{t-\Delta t} + 2 \Delta t [F_{q_{r}}]_{i,k}^{t}$	(3.30)

ここで,

$\displaystyle [F_{\theta}]_{i,k}$	$\textstyle =$	$\displaystyle - \left[{\rm Adv}.{\theta}\right]_{i,k}^{t} - \left[{\rm Adv}.{\b... ...]_{i,k}^{t - \Delta t} + \left[{\rm Diff}.{\theta} \right]_{i,k}^{t - \Delta t}$
		$\displaystyle + [Q_{cnd}]_{i,k}^{t} + [Q_{rad}]_{i,k}^{t-\Delta t} + [Q_{dis}]_{i,k}^{t-\Delta t}$	(3.31)
$\displaystyle \left[F_{q_{v}}\right]_{i,k}^{t}$	$\textstyle =$	$\displaystyle - \left[ {\rm Adv}.q_{v} \right]_{i,k}^{t} - \left[ {\rm Adv}.\ba... ...,k}^{t - \Delta t} + \left[ {\rm Turb}.\bar{q_{v}} \right]_{i,k}^{t - \Delta t}$
		$\displaystyle + \left[ {\rm Diff}.q_{v} \right]_{i,k}^{t - \Delta t} + \left[ EV_{rv} \right]_{i,k}^{t}$	(3.32)
$\displaystyle \left[F_{q_{c}}\right]_{i,k}^{t}$	$\textstyle =$	$\displaystyle - \left[ {\rm Adv}.q_{c} \right]_{i,k}^{t} + \left[ {\rm Turb}.q_{c} \right]^{t-\Delta t} + \left[ {\rm Diff}.q_{c} \right]^{t-\Delta t}$
		$\displaystyle - \left[ CN_{cr} + CL_{cr} \right]_{i,k}^{t}$	(3.33)
$\displaystyle \left[F_{q_{r}}\right]_{i,k}^{t}$	$\textstyle =$	$\displaystyle - \left[ {\rm Adv}.q_{r} \right]_{i,k}^{t} + \left[ {\rm Turb}.q_{r} \right]_{i,k}^{t-\Delta t} + \left[ {\rm Diff}.q_{c} \right]^{t-\Delta t}$
		$\displaystyle + \left[CN_{cr} + CL_{cr} - EV_{rv} \right]_{i,k}^{t} + \left[ PR_{r} \right]_{i,k}^{t}$	(3.34)

である. 移流を中心差分で安定して解くために, 数値粘性項 Diff を追加してある. また, $CN_{vc}, EV_{cv}$ 項は湿潤飽和調節法より決めるため, それらの項を含めない.

$\theta$ , $q_{v}$ , $q_{c}$ , $q_{r}$ をまとめて $\phi$ で表し, それぞれの項を書き下す. 移流項は,

$\displaystyle \left[{\rm Adv}.{\phi}\right]_{i,k}^{t}$

$\textstyle =$

$\displaystyle \left[ u_{i(u),k} \left[ \DP{\phi}{x} \right]_{i(u),k} \right]_{i,k}^{t} + \left[ w_{i,k(w)} \left[ \DP{\phi}{z} \right]_{i,k(w)} \right]_{i,k}^{t}$

(3.35)

であり, 基本場の移流項は,

$\displaystyle \left[{\rm Adv}.{\bar{\phi}}\right]_{i,k}^{t} = \left[ w_{i,k(w)} \left[ \DP{\overline{\phi}}{z} \right]_{i,k(w)} \right]_{i,k}^{t}$

(3.36)

である. 粘性拡散項は CReSS と同様に 1.5 次のクロージャーを用いることで,

$\displaystyle \left[{\rm Turb}.{\phi} \right]_{i,k}^{t - \Delta t}$	$\textstyle =$	$\displaystyle \left[ \DP{}{x} \left\{ \left( K_{h} \right)_{i(u),k} \left( \DP{\phi}{x} \right)_{i(u),k} \right\} \right]_{i,k}^{t - \Delta t}$
		$\displaystyle + \left[ \DP{}{z}\left\{ \left( K_{h} \right)_{i,k(w)} \left( \DP{\phi }{z} \right)_{i,k(w)} \right\} \right]_{i,k}^{t - \Delta t}$	(3.37)

となり, 基本場の粘性拡散項は,

$\displaystyle \left[{\rm Turb}.{\bar{\phi}} \right]_{i,k}^{t - \Delta t}$

$\textstyle =$

$\displaystyle \left[ \DP{}{z}\left\{ \left( K_{h} \right)_{i,k(w)} \left( \DP{\overline{\phi}}{z} \right)_{i,k(w)} \right\} \right]_{i,k}^{t - \Delta t}$

(3.38)

となる. 数値粘性項は,

$\displaystyle \left[ {\rm Diff}_{\phi} \right]_{i,k}^{t - \Delta t}$

$\textstyle =$

$\displaystyle \nu_{h} \left\{ \DP{}{x} \left(\DP{\phi}{x}\right)_{i(u),k} \righ... ...eft\{ \DP{}{z} \left(\DP{\phi}{z}\right)_{i,k(w)} \right\}_{i,k}^{t - \Delta t}$

(3.39)

である. $K_{h}$ は乱流エネルギーの時間発展方程式から計算する(詳細は後述). $\nu_{h}, \nu_{v}$ は (3.29) 式を利用する.

凝縮加熱項 $Q_{cnd}$ は

$\begin{displaymath}[Q_{cnd}]_{i,k}^{t} = - \left[ \frac{L}{{c_{p} \bar{\pi}}_{... ...{i,k}^{t}) (\bar{\rho}_{i,k} [q_{r}]_{i,k})^{0.65} \right\} \end{displaymath}$

(3.40)

である.

散逸加熱項 $Q_{dis}$ は

$\begin{displaymath}[Q_{dis}]_{i,k}^{t-\Delta t} = \frac{1}{{c_{p}}_{d} \bar{\p... ...{\pi}} \frac{(K_{m,i,k}^{t-\Delta t})^{3}}{{C_{m}}^{2} l^{4}} \end{displaymath}$

(3.41)

と与える. ここで $l=(\Delta x\Delta z)^{1/2}$ である.

放射強制 $[Q_{rad}]_{i,k}$ は計算設定ごとに与える.

雲水から雨水への変換を表す $CN_{cr}$ , $CL_{cr}$ は以下のようになる.

		$\displaystyle [CN_{cr}]_{i,k}^{t} = (q_{c,i,k}^{t} - q_{c0})/\tau _{ac}$	(3.42)
		$\displaystyle \left[ CL_{cr} \right]_{i,k}^{t} = 2.2 [q_{c}]_{i,k}^{t} \left( [\bar{\rho}]_{i,k} \left[ q_{r} \right]_{i,k}^{t} \right)^{0.875}$	(3.43)

雨水の蒸発を表す $EV_{rv}$ は以下のようになる.

$\displaystyle \left[ EV_{rv} \right]_{i,k}^{t} = 4.85 \times 10^{-2} ([q_{vsw}]_{i,k}^{t} - [q_{v}]_{i,k}^{t}) ([\bar{\rho}]_{i,k} [q_{r}]_{i,k}^{t})^{0.65}$

(3.44)

降水による雨水フラックスを表す $PR_{r}$ は以下のように書ける.

		$\displaystyle \left[ PR_{r} \right]_{i,k}^{t} = \Dinv{[\bar{\rho}]_{i,k}} \DP{}{z}([\bar{\rho}]_{i,k} [U_{r}]_{i,k}^{t} [q_{r}]_{i,k}^{t}).$	(3.45)
		$\displaystyle [U_{r}]_{i,k}^{t} = 12.2 ([q_{r}]_{i,k}^{t})^{0.125}$	(3.46)

3.2.1 湿潤飽和調節法

Klemp and Wilhelmson (1983), CReSS ユーザーマニュアル(坪木と榊原, 2001) では, 水蒸気と雲水の間の変換を表す $-CN_{vc} + EV_{cv}$ は, Soong and Ogura (1973) において開発された湿潤飽和調節法を用いる. この方法はの断熱線と, $\mu_{vapar} = \mu_{condensed phase}$ の平衡条件( $\mu$ は化学ポテンシャル)の交わる温度・圧力・組成を反復的に求める数値解法である. 以下ではそのやり方を解説する.

3.2.1.1 飽和蒸気圧を用いる場合

湿潤飽和調節法を用いる場合, まず始めに (3.30) - (3.37) 式から求まる量にを添付し, $[\theta]^{*}$ , $[q_{v}]^{*}$ , $[q_{c}]^{*}$ , $[q_{r}]^{*}$ とする. 水に対する過飽和混合比

$\displaystyle \Delta q_{c} = MAX\{0, [q_{v}]^{*} - q_{vsw}([\theta]^{*})\}$

(3.47)

が $\Delta q_{c} > 0$ , もしくは雲粒混合比が $q_{c}^{*} > 0$ ならば, 次式を用いて暫定的に $\theta$ , $q_{v}$ , $q_{c}$ を求める.

$\displaystyle \left[ \theta \right]^{t + \Delta t}$	$\textstyle =$	$\displaystyle \theta^{} + \frac { \gamma ( [q_{v}]^{} - q_{vsw}([\theta]^{})) } { 1 + \gamma \DP{q_{vsw}([\theta]^{})}{\theta} }$	(3.48)
$\displaystyle \left[ q_{v} \right]^{t + \Delta t}$	$\textstyle =$	$\displaystyle [q_{v}]^{} + \frac{[\theta]^{} - [\theta]^{t + \Delta t}}{\gamma},$	(3.49)
$\displaystyle \left[ q_{c} \right]^{t + \Delta t}$	$\textstyle =$	$\displaystyle [q_{v}]^{} + [q_{c}]^{} - [q_{v}]^{t + \Delta t} .$	(3.50)

ただし, $\gamma = L_{v}/(c_{p} \Pi)$ である. もしも $[q_c]^{t + \Delta t} > 0$ ならば, 暫定的に得られた値を

付きのものに置き換え, (3.51) - (3.53) 式の値が収束するまで繰り返し適用する. 普通, 高々数回繰り返せば収束し, 調整後の値が得られるそうである.

もしも $q_{c}^{t + \Delta t} < 0$ の場合には,

		$\displaystyle \left[ \theta \right]^{t + \Delta t} = [\theta]^{} - \gamma [q_{c}]^{},$	(3.51)
		$\displaystyle \left[ q_{v} \right]^{t + \Delta t}, = [q_{v}]^{} + [q_{c}]^{}$	(3.52)
		$\displaystyle \left[ q_{c} \right]^{t + \Delta t} = 0$	(3.53)

とし, 繰り返しを中止する.

3.2.1.2 圧平衡定数を用いる場合

硫化アンモニウムの声齔瘤�竚癈鷭∂焜聨纃瘟赧漓�籬�㏍聽轣蛹就恭汲桶岩霞⊂桿轣蛹Γ蔚飴頏阡繝�籟鹿齔瘤梢応

$\displaystyle {\rm NH_{3}} + {\rm H_{2}S} \rightarrow {\rm NH_{4}SH}$

(3.54)

のような, 2 種類の気体 1 モルづつから凝縮物質 1 モルが生成されるような生成反応の場合の, 湿潤飽和調節法を考える.

硫化アンモニウムの声齔瘤�竚癈鷭∂焜聨纃瘟赧漓�籬�㏍聽轣蛹就宴牡泣害鬼⊂桿轣蛹Γ蔚飴頏阡繝�籟鹿齔瘤梢応の圧平衡定数は,

$\displaystyle K_{p} \equiv \ln(p_{\rm NH_{3}} \cdot p_{\rm H_{2}S}) = 61.781 - \frac{10834}{T} - \ln{10^{2}}$

(3.55)

である. 圧平衡定数を用いることで, 任意の温度に対するアンモニアと硫化水素のモル比の積を求めることができる.

任意の温度における NH $_{4}$ SH の声齔瘤�竚癈鷭∂焜聨纃瘟赧漓�籬�㏍聽轣蛹就聿齋砌箟皹盂⊂桿轣蛹Γ蔚飴頏阡繝�籟鹿齔瘤召� とすると, 圧平衡定数の式は以下のように書ける.

		$\displaystyle (p_{\rm NH_{3}} - X) ( p_{\rm H_{2}S} - X ) = e^{k_{p}}$
		$\displaystyle X^{2} - (p_{\rm NH_{3}} + p_{\rm H_{2}S}) X + p_{\rm NH_{3}} \cdot p_{\rm H_{2}S} - e^{k_{p}} = 0$	(3.56)

解の公式を使うと, 声齔瘤�竚癈鷭∂焜聨纃瘟赧漓�籬�㏍聽轣蛹就号畄粭皹箟箜⊂桿轣蛹Γ蔚飴頏阡繝�籟鹿齔瘤� X は以下となる.

		$\displaystyle X = \Dinv{2} \left\{ (p_{\rm NH_{3}} + p_{\rm H_{2}S}) \pm \sqrt{... ...m H_{2}S})^{2} - 4 (p_{\rm NH_{3}} \cdot p_{\rm H_{2}S} - e^{K_{p}}) } \right\}$
		$\displaystyle X = \Dinv{2} \left\{ (p_{\rm NH_{3}} + p_{\rm H_{2}S}) \pm \sqrt{ (p_{\rm NH_{3}} - p_{\rm H_{2}S})^{2} + 4 e^{K_{p}} } \right\}$	(3.57)

根号の符号は $\exp{(K_{p})} \approx 0$ の場合にとりうる

の値を仮定することで決める. $\exp{(K_{p})} \approx 0$ の場合, 明らかに

$\displaystyle X = {\rm min}(P_{\rm NH_3}, P_{\rm H_{2}S} )$

(3.58)

である. ここで木星大気を想定し, $P_{\rm NH_3} > P_{\rm H_{2}S}$ であることを仮定すると $X = P_{\rm H_{2}S}$ である. そこで (3.60)の根号の符号は $\exp{(K_{p})} \approx 0$ のとき $X = P_{\rm H_{2}S}$ となるよう, 負を選択する.

$\displaystyle X = \Dinv{2} \left\{ (p_{\rm NH_{3}} + p_{\rm H_{2}S}) - \sqrt{ (p_{\rm NH_{3}} - p_{\rm H_{2}S})^{2} + 4 e^{K_{p}} }. \right\}$

(3.59)

の満たすべき条件は,

$\displaystyle 0 \leq X \leq {\rm min}(P_{\rm NH_{3}}, P_{\rm H_{2}S})$

(3.60)

である. 上記の条件を満たさない場合には

とする.

が (3.63) 式の条件を満たすならば, 次式を用いて暫定的に $\theta$ , $q_{v}$ , $q_{c}$ を求める.

		$\displaystyle \left[ q_{\rm NH_3} \right]^{t + \Delta t} = [q_{\rm NH_3}]^{*} + \Delta q_{\rm NH_3},$	(3.61)
		$\displaystyle \left[ q_{\rm H_2S} \right]^{t + \Delta t} = [q_{\rm H_2S}]^{*} + \Delta q_{\rm H_2S},$	(3.62)
		$\displaystyle \left[ q_{\rm NH_4SH} \right]^{t + \Delta t} = [q_{\rm NH_3}]^{}... ...rm H_2S}]^{} - [q_{\rm NH_3}]^{t + \Delta t} - [q_{\rm H_2S}]^{t + \Delta t} ,$	(3.63)
		$\displaystyle \left[ \theta \right]^{t + \Delta t} = \theta^{} + \gamma \left(... ...}]^{} - [q_{\rm NH_3}]^{t + \Delta t} - [q_{\rm H_2S}]^{t + \Delta t} \right).$	(3.64)

ただし, $\gamma = L_{\rm NH_4SH}/({c_{p}}_{d} \Pi)$ であり, $\Delta q_{\rm NH_3}$ と $\Delta q_{\rm H_2S}$ はそれぞれ, 声齔瘤�竚癈鷭∂焜聨纃瘟赧漓�籬�㏍聽轣蛹就鍵盡祕痳竇窿⊂桿轣蛹Γ蔚飴頏阡繝�籟鹿齔瘤�

に対応する NH

と H

S の混合比である. 温位が収束するまで反復改良を行う.

3.3 乱流運動エネルギーの式

Klemp and Wilhelmson (1978) および CReSS (坪木と榊原篤志, 2001) と同様に, 1.5 次のクロージャーを用いる. 乱流エネルギーの時間発展方程式をリープフロッグ法を用いて時間方向に離散化すると, 以下のようになる.

$\displaystyle [K_{m}]_{i,k}^{t + \Delta t} = [K_{m}]_{i,k}^{t - \Delta t} + 2 \Delta t [F_{K_m}]_{i,k}^{t}$

(3.65)

ここで,

$\displaystyle [F_{K_m}]_{i,k}^{t}$	$\textstyle =$	$\displaystyle - [{\rm Adv}.K_m]_{i,k}^{t} + [{\rm Buoy}.K_m]_{i,k}^{t - \Delta t} + [{\rm Shear}.K_m]_{i,k}^{t - \Delta t}$
		$\displaystyle + [{\rm Turb}.K_m]_{i,k}^{t - \Delta t} + [{\rm Disp}.K_m]_{i,k}^{t - \Delta t}$	(3.66)

である. CReSS にならい, 移流項を

で, 移流項以外を $t - \Delta t$ で評価した.

$F_{K_m}$ に含まれる各項は以下のように書き下すことができる.

$\displaystyle [{\rm Adv}.K_m]_{i,k}^{t}$	$\textstyle =$	$\displaystyle \left\{ u_{i(u),k} \left( \DP{K_{m}}{x} \right)_{i(u), k} \right\... ... + \left\{ w_{i,k(w)} \left( \DP{K_{m}}{z} \right)_{i, k(w)} \right\}_{i,k}^{t}$	(3.67)
$\displaystyle \left[{\rm Buoy}.K_m\right]_{i,k}^{t - \Delta t}$	$\textstyle =$	$\displaystyle - \left\{ \frac{3 g C_{m}^{2} l^{2}}{ 2 \overline{\theta}} \left(\DP{\theta_{el}}{z} \right)_{i,k(w)} \right\}_{i,k}^{t-\Delta t}$	(3.68)
$\displaystyle \left[{\rm Shear}.K_m\right]_{i,k}^{t - \Delta t}$	$\textstyle =$	$\displaystyle \left( C_{m}^{2} l^{2} \right)_{i,k} \left[ \left( \DP{u}{x} \rig... ... \right)_{i,k} \left[ \left( \DP{w}{z} \right)^{2} \right]_{i,k}^{t - \Delta t}$
		$\displaystyle + \left( \frac{ C_{m}^{2} l^{2} }{2} \right)_{i,k} \left[ \left\{... ... \left( \DP{w}{x} \right)_{i(u),k(w)} \right\}_{i,k}^{t - \Delta t} \right]^{2}$
		$\displaystyle - \left( \frac{K_{m}}{3} \right)_{i,k}^{t - \Delta t} \left\{ \le... ...ight)_{i,k}^{t-\Delta t} + \left( \DP{w}{z} \right)_{i,k}^{t-\Delta t} \right\}$	(3.69)
$\displaystyle \left[{\rm Turb}.K_m\right]_{i,k}^{t - \Delta t}$	$\textstyle =$	$\displaystyle \Dinv{2} \left[ \left\{ \DP{}{x} \left( \DP{K_{m}^{2}}{x} \right)... ...} \left( \DP{K_{m}^{2}}{z} \right)_{i,k(w)} \right\}_{i,k}^{t-\Delta t} \right]$
		$\displaystyle + \left[ \left\{ \left( \DP{K_{m}}{x}\right)^{2} \right\}_{i(u),k... ...{ \left(\DP{K_{m}}{z}\right)^{2} \right\}_{i,k(w)} \right]_{i,k}^{t - \Delta t}$	(3.70)
$\displaystyle \left[{\rm Disp}.K_m\right]_{i,k}^{t - \Delta t}$	$\textstyle =$	$\displaystyle - \Dinv{2 l^{2}} \left( K_{m}^{2} \right)_{i,k}^{t - \Delta t}$	(3.71)

ここで $C_{\varepsilon} = C_{m} = 0.2$ , 混合距離 $l = \left(\Delta x \Delta z \right)^{1/2}$ とする. また $\theta_{el}$ は以下で与えられる.

$\displaystyle \theta_{el}$	$\textstyle =$	$\displaystyle \overline{ \theta_{v}} + \theta_{v}^{'} \;\;\; (for \;\; q_{c} = 0)$	(3.72)
$\displaystyle \theta_{el}$	$\textstyle =$	$\displaystyle \overline{\theta_{v}} + \theta_{v}^{'} + \frac{ \sum L q_{v}}{{c_p}_{d} \bar{\pi}} \;\;\; (for \;\; q_{c} > 0)$	(3.73)

ただし,

$\displaystyle \overline{\theta_{v}} + \theta_{v}^{'}$

$\textstyle =$

$\displaystyle \bar{\theta_{v}} \left\{ 1 + \frac{\theta}{\bar{\theta}} + \frac{... ...} - \frac{\sum q_{v} + \sum q_{c} + \sum q_{r}} {1 + \sum \bar{q_{v}}} \right\}$

(3.74)

である.

3.4 時間フィルター

リープフロッグ法を用いたことによって生じる計算モードの増幅を抑制するため, Asselin (1972) の時間フィルターを長い時間刻みで 1 ステップ計算する毎に(実際には短い時間刻みの計算を $N_{\tau}\equiv 2\Delta t/\Delta \tau$ ステップ計算する毎に)適用する.

たとえば(3.6)を用いて $u^{t + \Delta t}_{i(u),k}$ を計算する場合, 以下のように時間フィルターを適用する.

$\displaystyle u^{*}_{i(u),k}$	$\textstyle =$	$\displaystyle u^{\tau + (N_{\tau}-1)\Delta \tau}_{i(u), k} - \left[ \bar{c_{p}}... ...\DP{(\alpha Div)^{\tau + (N_{\tau}-1)\Delta \tau}}{x} \right\} \right]_{i(u),k}$
		$\displaystyle + F_{u,i(u),k}^{t} \Delta \tau,$
$\displaystyle u^{t+\Delta t}_{i(u),k}$	$\textstyle =$	$\displaystyle (1-2 \gamma)u^{t}_{i(u),k} + \gamma (u^{*}_{i(u),k} + u^{t -\Delta t}_{i(u),k})$	(3.75)

ここで $\gamma$ はフィルターの係数であり, その値は 0.05 を用いる. (3.7), (3.8)の計算に対しても同様に時間フィルターを適用する.

3.5 スポンジ層

境界面付近での波の反射を抑えるために, 基礎方程式の付加的な項を付け加える.

$\displaystyle \DP{\phi}{t} = -{\rm Adv}.\phi + \cdots + \gamma_{h}(x) (\phi - \phi_{e}) + \gamma_{v}(z) (\phi - \phi_{e})$

(3.76)

ただし, $\phi$ は任意の予報変数であり, $\phi_{e}$ は客観解析値等の既知の値である. この項は1 つ前のタイムステップ $t - \Delta t$ で計算され, 小さいタイムステップで扱われる予報変数に対しても, 移流項や数値粘性項と同様に $2 \Delta t$ の大きなタイムステップ間の値として評価される。具体的には,

		$\displaystyle [\pi]^{t + \Delta t} = 2 \Delta t \left\{ [{\rm Adv}.\pi]^{t} + \... ...gamma_{h}(x) + \gamma_{v}(z) \right\} (\pi - \bar{\pi})^{t - \Delta t} \right\}$	(3.77)
		$\displaystyle [u]^{t + \Delta t} = 2 \Delta t \left\{ [{\rm Adv}.u]^{t} + \cdots + \left\{ \gamma_{h}(x) + \gamma_{v}(z) \right\} [u]^{t - \Delta t} \right\}$	(3.78)
		$\displaystyle [w]^{t + \Delta t} = 2 \Delta t \left\{ [{\rm Adv}.w]^{t} + \cdots + \left\{ \gamma_{h}(x) + \gamma_{v}(z) \right\} [w]^{t - \Delta t} \right\}$	(3.79)
		$\displaystyle [\theta]^{t + \Delta t} = 2 \Delta t \left\{ [{\rm Adv}.\theta]^{... ...\left\{ \gamma_{h}(x) + \gamma_{v}(z) \right\} [\theta]^{t - \Delta t} \right\}$	(3.80)

とする. 但し $\bar{\pi}$ はエクスナー関数の基本場である.

$\gamma_{h}, \gamma_{v}$ はそれぞれ水平方向には各境界面に向かって, 鉛直方向には上境界面に向かって小さくなる減衰係数である. これらの減衰係数は, 水平方向には吸収層の厚みを $d_{h}$ とし, の範囲を $0 \leq x \leq x_{max}$ とすれば,

		$\displaystyle \gamma_{h} = \alpha_{h} \left( 1 - \frac{x}{d_{h}}\right)^{3} \hspace{5em} (x < d_{h}),$
		$\displaystyle \gamma_{h} = 0 \hspace{10em} ( d_{h} \leq x \leq x_{max} - d_{h}),$
		$\displaystyle \gamma_{h} = \alpha_{h} \left( 1 - \frac{(x_{max} - x)}{d_{h}}\right)^{3} \;\; (x > x_{max} - d_{h}),$	(3.81)

であり, 鉛直方向には吸収層の厚さを $d_{v}$ とし,

の範囲を $0 \leq z \leq z_{max}$ とすれば,

		$\displaystyle \gamma_{v} = 0 \hspace{12em} ( z \leq z_{max} - d_{v}),$
		$\displaystyle \gamma_{v} = \alpha_{v} \left ( 1 - \cos{\frac{\pi (z - z_{max} - d_{v})}{d_{v}}} \right)^{3} \;\; (z > z_{max} - d_{v}),$	(3.82)

である. ここで, $\alpha_h, \alpha_v$ はそれぞれ水平・鉛直方向の減衰定数である. $\alpha_h, \alpha_v$ は時間の逆数の次元を持ち, それらの逆数 $1/\alpha_{h}, 1/\alpha_{v}$ は e-folding time と呼ばれる. e-folding time は通常 100 - 300 s に設定する. また吸収層の厚み $d_{h}, d_{v}$ はそれぞれ, 水平方向には数格子分, 鉛直方向には上面から1/3 程度設定すれば良い.

: 4. 参考文献 : 2 次元非静力学モデルの離散化 : 2. 空間方向の離散化

Odaka Masatsugu 平成19年9月11日