3. 時間方向の離散化

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

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

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

$$\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)

$$\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)

$$\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)

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

$$[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)

$$\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)

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

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

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

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

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

$$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)

$$\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}$ を消去する.

$$\pi^{\tau + \Delta \tau} - \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\}$$

$$= \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}$$

$$- \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章 参照).

$$\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}$$

$$\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}$$

$$\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}$$

$$= \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}$$

$$\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.$$

$$\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}.$$

但し平均場の量は鉛直方向にしか依存しないので $z$ 方向の添字のみ 付けてある.

3.1.1.3 境界条件

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

$$w(i,0(w)) = 0,$$ (3.8)

$$w(i,km(w)) = 0$$ (3.9)

である.

下部境界:

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

$$\beta \left( \DP{\pi^{\tau + \Delta \tau}}{z} \right)_{i,0(w)} = \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)}$$

$$\equiv E_{i,0(w)}$$ (3.10)

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

$$\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}$$

$$= \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}$$

$$- \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}$$

$$- \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)

上部境界:

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

$$\beta \left( \DP{\pi^{\tau + \Delta \tau}}{z} \right)_{i,km(w)} = \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)}$$

$$\equiv E_{i,km(w)}$$ (3.12)

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

$$\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}$$

$$+ \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}$$

$$= \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}$$

$$- \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}$$

$$+ \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) 式を連立すると, 以下のような行列式の形式で書く ことができる.

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

$$= \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}$ を求める. この連立方程式の係数は以下の ように書ける.

$$A_{k} = 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\}$$

$$(k = 2, 3, \cdots km-1),$$

$$A_{1} = 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)},$$

$$A_{km} = 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)},$$

$$B_{k} = - \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)},$$

$$(k = 2, 3, \cdots km),$$

$$C_{k} = - \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)},$$

$$(k = 1, 2, \cdots km-1),$$

$$D_{i,k} = \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}$$

$$(k = 2, 3, \cdots km-1),$$

$$D_{i,1} = \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}$$

$$- \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)},$$

$$D_{i,km} = \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}$$

$$+ \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)}.$$

ただし,

$$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)}$$

$$F_{i,k} \equiv \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.$$

$$\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) 式を 離散化する.

$$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)

$$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 は数値粘性項である. それぞれの項を書き下すと,

$$\left[ {\rm Adv}.{u} \right]_{i(u),k}^{t} = 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)

$$\left[ {\rm Adv}.{w} \right]_{i,k(w)}^{t} = 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)

であり, 浮力項は,

$$[{\rm Buoy}]^{t}_{i,k(w)} = g \frac{\theta_{i,k(w)}^{t}}{\overline{\theta}_{i,k(w)}}$$

$$+ g \frac{\sum [q_{v}]_{i,k(w)}^{t}/M_{v}}{1/M_{d} + \sum [\bar{q_{v}}]_{i,k(w)}/M_{v}}$$

$$- 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)

であり, 粘性拡散項は,

$$\left[ {\rm Turb}.{u} \right]_{i(u),k}^{t - \Delta t} = 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}$$

$$+\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}$$

$$- \frac{2}{3 C_{m}^{2} l^{2}} \left( \DP{ K_{m}^{2} }{x} \right)_{i(u),k}^{t - \Delta t}$$ (3.20)

$$\left[ {\rm Turb}.{w} \right]_{i,k(w)}^{t - \Delta t} = 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}$$

$$+\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}$$

$$- \frac{2}{3 C_{m}^{2} l^{2}} \left( \DP{ K_{m}^{2} }{z} \right)_{i,k(w)}^{t - \Delta t}$$ (3.21)

である. 数値粘性項は,

$$\left[ {\rm Diff}.u \right]_{i(u),k}^{t - \Delta t} = \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)

$$\left[ {\rm Diff}.w \right]_{i,k(w)}^{t - \Delta t} = \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}$ は以下のように定める.

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

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

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

$$\alpha_{h} \le \Dinv{8}, \hspace{3em} \alpha_{v} \le \Dinv{8}$$ (3.26)

とする.

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

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

$$\theta_{i,k}^{t + \Delta t} = \theta_{i,k}^{t - \Delta t} + 2 \Delta t [F_{\theta}]_{i,k}^{t}$$ (3.27)

$$\left[ q_{v} \right]_{i,k}^{t+\Delta t} = \left[ q_{v} \right]_{i,k}^{t-\Delta t} + 2 \Delta t [F_{q_{v}}]_{i,k}^{t}$$ (3.28)

$$\left[ q_{c} \right]_{i,k}^{t+\Delta t} = \left[ q_{c} \right]_{i,k}^{t-\Delta t} + 2 \Delta t [F_{q_{c}}]_{i,k}^{t}$$ (3.29)

$$\left[ q_{r} \right]_{i,k}^{t+\Delta t} = \left[ q_{r} \right]_{i,k}^{t-\Delta t} + 2 \Delta t [F_{q_{r}}]_{i,k}^{t}$$ (3.30)

ここで,

$$[F_{\theta}]_{i,k} = - \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}$$

$$+ [Q_{cnd}]_{i,k}^{t} + [Q_{rad}]_{i,k}^{t-\Delta t} + [Q_{dis}]_{i,k}^{t-\Delta t}$$ (3.31)

$$\left[F_{q_{v}}\right]_{i,k}^{t} = - \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}$$

$$+ \left[ {\rm Diff}.q_{v} \right]_{i,k}^{t - \Delta t} + \left[ EV_{rv} \right]_{i,k}^{t}$$ (3.32)

$$\left[F_{q_{c}}\right]_{i,k}^{t} = - \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}$$

$$- \left[ CN_{cr} + CL_{cr} \right]_{i,k}^{t}$$ (3.33)

$$\left[F_{q_{r}}\right]_{i,k}^{t} = - \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}$$

$$+ \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$ で表し, それぞれの項を書き下す. 移流項は,

$$\left[{\rm Adv}.{\phi}\right]_{i,k}^{t} = \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)

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

$$\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 次のクロージャーを用いることで,

$$\left[{\rm Turb}.{\phi} \right]_{i,k}^{t - \Delta t} = \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}$$

$$+ \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)

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

$$\left[{\rm Turb}.{\bar{\phi}} \right]_{i,k}^{t - \Delta t} = \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)

となる. 数値粘性項は,

$$\left[ {\rm Diff}_{\phi} \right]_{i,k}^{t - \Delta t} = \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}$ は

$$[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\}$$ (3.40)

である.

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

$$[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}}$$ (3.41)

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

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

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

$$CN_{cr} = \frac{10^{6} [\bar{\rho}]_{i,k} \left( [q_{c}]_{i,k}^{t... ...60 (2 [q_{c}]_{i,k}^{t} + 2.66 \times 10^{-8} N_{0}/ [\bar{\rho}]_{i,k} D_{0})}$$ (3.42)

$$\hspace{2em} N_{0} = 5.0 \times 10^{7}, \hspace{1em} D_{0} = 0.366$$ (3.43)

$$\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.44)

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

$$\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.45)

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

$$\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.46)

$$[U_{r}]_{i,k}^{t} = 12.2 ([q_{r}]_{i,k}^{t})^{0.125}$$ (3.47)

3.2.1 湿潤飽和調節法

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

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

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

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

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

$$\left[ \theta \right]^{t + \Delta t} = \theta^{*} + \frac { \gamma ( [q_{v}]^{*} - q_{vsw}([\theta]^{*})) } { 1 + \gamma \DP{q_{vsw}([\theta]^{*})}{\theta} }$$ (3.49)

$$\left[ q_{v} \right]^{t + \Delta t} = [q_{v}]^{*} + \frac{[\theta]^{*} - [\theta]^{t + \Delta t}}{\gamma},$$ (3.50)

$$\left[ q_{c} \right]^{t + \Delta t} = [q_{v}]^{*} + [q_{c}]^{*} - [q_{v}]^{t + \Delta t} .$$ (3.51)

ただし, $\gamma = L_{v}/(c_{p} \Pi)$ である. もしも $[q_c]^{t + \Delta t} > 0$ ならば, 暫定的に得られた値を $*$ 付き のものに置き換え, (3.52) - (3.54) 式 の値が収束するまで繰り返し適用する. 普通, 高々数回繰り返せば収束し, 調整後の値が得られるそうである.

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

$$\left[ \theta \right]^{t + \Delta t} = [\theta]^{*} - \gamma [q_{c}]^{*},$$ (3.52)

$$\left[ q_{v} \right]^{t + \Delta t}, = [q_{v}]^{*} + [q_{c}]^{*}$$ (3.53)

$$\left[ q_{c} \right]^{t + \Delta t} = 0$$ (3.54)

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

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

硫化アンモニウムの生成反応

$${\rm NH_{3}} + {\rm H_{2}S} \rightarrow {\rm NH_{4}SH}$$ (3.55)

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

硫化アンモニウムの生成反応の圧平衡定数は,

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

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

任意の温度 $T$ における NH$_{4}$SH の生成量を $X$ とすると, 圧平衡定数の式は以下のように書ける.

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

$$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.57)

解の公式を使うと, 生成量 X は以下となる.

$$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\}$$

$$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.58)

ただし $\exp{(K_{p})} \approx 0$ の場合には, 明らかに

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

なために, $X$ は符号は負しか取り得ない.

$$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.60)

$X$ の満たすべき条件は,

$${\rm min}(P_{\rm NH_{3}}, P_{\rm H_{2}S}) > X, \hspace{1em} X < -... ...M_{\rm NH_4SH} \cdot {\rm min}(P_{\rm NH_{3}}, P_{\rm H_{2}S})}{M_{d} P_{all} }$$ (3.61)

である. 上記の条件を満たさない場合には $X = 0$ とする.

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

$$\left[ q_{\rm NH_3} \right]^{t + \Delta t} = [q_{\rm NH_3}]^{*} + \Delta q_{\rm NH_3},$$ (3.62)

$$\left[ q_{\rm H_2S} \right]^{t + \Delta t} = [q_{\rm H_2S}]^{*} + \Delta q_{\rm H_2S},$$ (3.63)

$$\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.64)

$$\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.65)

ただし, $\gamma = L_{\rm NH_4SH}/({c_{p}}_{d} \Pi)$ であり, $\Delta q_{\rm NH_3}$ と $\Delta q_{\rm H_2S}$ はそれぞれ, 生成量 $X$ に対応する NH$_3$ と H$_2$S の混合比である. 温位が収束するまで反復改良を行う.

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

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

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

ここで,

$$[F_{K_m}]_{i,k}^{t} = - [{\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}$$

$$+ [{\rm Turb}.K_m]_{i,k}^{t - \Delta t} + [{\rm Disp}.K_m]_{i,k}^{t - \Delta t}$$ (3.67)

である. CReSS にならい, 移流項を $t$ で, 移流項以外を $t - \Delta t$ で評価した.

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

$$[{\rm Adv}.K_m]_{i,k}^{t} = \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.68)

$$\left[{\rm Buoy}.K_m\right]_{i,k}^{t - \Delta t} = - \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.69)

$$\left[{\rm Shear}.K_m\right]_{i,k}^{t - \Delta t} = \left( C_{m}^{2} l^{2} \right)_{i,k} \left[ \left\{ \left( \DP{u}... ...t\{ \left( \DP{w}{z} \right)^{2} \right\}_{i,k(w)} \right]_{i,k}^{t - \Delta t}$$

$$+ \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}$$

$$- \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.70)

$$\left[{\rm Turb}.K_m\right]_{i,k}^{t - \Delta t} = \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]$$

$$+ \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.71)

$$\left[{\rm Disp}.K_m\right]_{i,k}^{t - \Delta t} = - \Dinv{2 l^{2}} \left( K_{m}^{2} \right)_{i,k}^{t - \Delta t}$$ (3.72)

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

$$\theta_{el} = \overline{ \theta_{v}} + \theta_{v}^{'} \;\;\; (for \;\; q_{c} = 0)$$ (3.73)

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

ただし,

$$\overline{\theta_{v}} + \theta_{v}^{'} = \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.75)

である.

3.4 時間フィルター

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

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

$$u^{*}_{i(u),k} = 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}$$

$$+ F_{u,i(u),k}^{t} \Delta \tau,$$

$$u^{t+\Delta t}_{i(u),k} = (1-2 \gamma)u^{t}_{i(u),k} + \gamma (u^{*}_{i(u),k} + u^{t -\Delta t}_{i(u),k})$$ (3.76)

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

3.5 スポンジ層

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

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

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

$$[\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.78)

$$[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.79)

$$[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.80)

$$[\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.81)

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

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

$$\gamma_{h} = \alpha_{h} \left( 1 - \frac{x}{d_{h}}\right)^{3} \hspace{5em} (x < d_{h}),$$

$$\gamma_{h} = 0 \hspace{10em} ( d_{h} \leq x \leq x_{max} - d_{h}),$$

$$\gamma_{h} = \alpha_{h} \left( 1 - \frac{(x - x_{max})}{d_{h}}\right)^{3} \;\; (x > x_{max} - d_{h}),$$ (3.82)

であり, 鉛直方向には吸収層の厚さを $d_{v}$ とし, $z$ の範囲を $0 \leq z \leq z_{max}$ とすれば,

$$\gamma_{v} = 0 \hspace{12em} ( z \leq z_{max} - d_{v}),$$

$$\gamma_{v} = \alpha_{v} \left ( 1 - \cos{\frac{\pi (z - z_{max})}{d_{v}}} \right)^{3} \;\; (z > z_{max} - d_{v}),$$ (3.83)

である. ここで, $\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 程度設定すれば良い.