3. 力学過程

3.1 はじめに

この章では力学過程の支配方程式を記し, その支配方程式の離散化を 行う.

ここで述べる力学過程とは, 流体の支配方程式における外力項を除いた部分を指す. 外力項である放射や鉛直乱流拡散や雲などに関する過程については 別紙を参照のこと.

離散化については, 空間に関する離散化である鉛直離散化と, 水平離散化の方法ならびに時間に関する離散化を行う.

3.2 支配方程式

ここでは力学過程の支配方程式系を示す. この方程式系の詳細に関しては, Haltiner and Williams (1980) もしくは 別紙『 支配方程式系の導出に関する参考資料』 の『力学過程の支配方程式系の導出』を参照せよ.

3.2.1 連続の式

$$\DP{\pi}{t} + \Dvect{v}_H \cdot \nabla_{\sigma} \pi = - D - \DP{\dot{\sigma}}{\sigma}.$$ (3.1)

3.2.2 静水圧の式

$$\DP{\Phi}{\sigma} = - \frac{RT_v}{\sigma}.$$ (3.2)

3.2.3 運動方程式

$$\DP{\zeta}{t} \ = \ \Dinv{a} \left( \Dinv{1 - \mu^2} \DP{V_A}{\lambda} - \DP{U_A}{\mu} \right) + {\cal D}(\zeta),$$ (3.3)

$$\DP{D}{t} \ = \ \Dinv{a} \left( \Dinv{1 - \mu^2} \DP{U_A}{\lambda} + \DP{V_A}{\mu} \right) - \nabla^{2}_{\sigma} ( \Phi + R \overline{T} \pi + ) + {\cal D}(D).$$ (3.4)

3.2.4 熱力学の式

$$\begin{align*}\begin{split}\DP{T}{t} \ &= \ - \Dinv{a} \left( \Dinv{1 - \mu^2} \... ...ac{Q}{C_p} + {\cal D}(T) + {\cal D}^{\prime}(\Dvect{v}). \end{split}\end{align*}$$ (3.5)

3.2.5 水蒸気の式

$$\begin{align*}\begin{split}\DP{q}{t} \ &= \ - \Dinv{a} \left( \Dinv{1 - \mu^2} \... ...uad - \dot{\sigma} \DP{q}{\sigma} + S_{q} + {\cal D}(q). \end{split}\end{align*}$$ (3.6)

ここで, 独立変数は以下の通りである.

$$\varphi : \quad 緯度 [\mathrm{deg.}],$$ (3.7)

$$\lambda : \quad 経度 [\mathrm{deg.}],$$ (3.8)

$$\sigma \equiv p/p_s,$$ (3.9)

$$t : \quad 時間 [\mathrm{s}].$$ (3.10)

ここで, $p$は気圧, $p_s$ は地表面気圧である. また $\mu \equiv \sin \varphi$ である.

モデルで時間発展を計算することとなる予報変数は以下の通りである.

$$\pi\ (\varphi, \lambda) \equiv \ln p_s,$$ (3.11)

$$T\ (\varphi, \lambda, \sigma) : \quad 気温 [\mathrm{K}],$$ (3.12)

$$q\ (\varphi, \lambda, \sigma) : \quad 比湿 [\mathrm{kg}\ \mathrm{kg}^{-1}],$$ (3.13)

$$\zeta\ (\varphi, \lambda, \sigma) \equiv \Dinv{a} \left( \Dinv{1 - \mu^2} \DP{V}{\lambda} - \DP{U}{\mu} \right) : \quad 渦度 [\mathrm{s}^{-1}],$$ (3.14)

$$D\ (\varphi, \lambda, \sigma) \equiv \Dinv{a} \left( \Dinv{1 - \mu^2} \DP{U}{\lambda} + \DP{V}{\mu} \right) : \quad 発散 [\mathrm{s}^{-1}].$$ (3.15)

ここで,

$$U(\varphi, \lambda, \sigma) \equiv u(\varphi, \lambda, \sigma) \cos \varphi,$$ (3.16)

$$V(\varphi, \lambda, \sigma) \equiv v(\varphi, \lambda, \sigma) \cos \varphi,$$ (3.17)

$$u : ,$$ (3.18)

$$v :$$ (3.19)

である. 流線関数$\psi$と速度ポテンシャル$\chi$を導入すると, $U$, $V$, $\zeta$, $D$はそれぞれ以下のように表わされる.

$$U = \Dinv{a} \left( \DP{\chi}{\lambda} - (1-\mu^2) \DP{\psi}{\mu} \right),$$ (3.20)

$$V = \Dinv{a} \left( \DP{\psi}{\lambda} + (1-\mu^2) \DP{\chi}{\mu} \right),$$ (3.21)

$$\zeta = \Dlapla \psi,$$ (3.22)

$$D = \Dlapla \chi.$$ (3.23)

各時間ステップで診断的に求められる変数は以下の通りである.

$$\Phi \equiv gz : \quad ジオポテンシャル高度 [\mathrm{m}^{2}\ \mathrm{s}^{-2}],$$ (3.24)

$$\dot{\sigma} \equiv \DD{\sigma}{t} \ \equiv \ \DP{\sigma}{t} + \frac{u}{a \cos... ... \DP{\sigma}{\lambda} + \frac{v}{a} \DP{\sigma}{\varphi} + \DP{\sigma}{\sigma},$$ (3.25)

$$\overline{T}\ (\sigma) : \quad 基準温度 [\mathrm{K}],$$ (3.26)

$$T^{\prime}\ (\varphi, \lambda, \sigma) \equiv T - \overline{T},$$ (3.27)

$$T_v\ (\varphi, \lambda, \sigma) \equiv T \left\{ 1 + \left(\epsilon_v^{-1} - 1\right) q \right\},$$ (3.28)

$$T_v^{\prime}\ (\varphi, \lambda, \sigma) \equiv T_v - \overline{T},$$ (3.29)

$$U_A\ (\varphi, \lambda, \sigma) \equiv ( \zeta + f ) V - \dot{\sigma} \DP{U}{\sigma} - \frac{R T_v^{\prime}}{a} \DP{\pi}{\lambda} + {\cal F}_{\lambda} \cos \varphi,$$ (3.30)

$$V_A\ (\varphi, \lambda, \sigma) \equiv - ( \zeta + f ) U - \dot{\sigma} \DP{V}{\sigma} - \frac{R T_v^{\prime}}{a} (1-\mu^2) \DP{\pi}{\mu} + {\cal F}_{\varphi} \cos \varphi,$$ (3.31)

$$\Dvect{v}_H \cdot \nabla_{\sigma} \pi \equiv \frac{U}{a (1 - \mu^2)} \DP{\pi}{\lambda} + \frac{V}{a} \DP{\pi}{\mu}$$ (3.32)

$$\nabla^{2}_{\sigma} \equiv \frac{1}{a^{2} (1-\mu^2)} \DP[2]{}{\lambda} + \frac{1}{a^{2}} \DP{}{\mu} \left[ (1-\mu^2) \DP{}{\mu} \right],$$ (3.33)

$$\ (\varphi, \lambda, \sigma) \equiv \frac{U^{2}+V^{2}}{2 (1-\mu^2) }$$ (3.34)

$${\cal D}(\zeta) : \qquad 渦度の水平拡散とスポンジ層における散逸,$$ (3.35)

$${\cal D}(D) : \qquad 発散の水平拡散とスポンジ層における散逸,$$ (3.36)

$${\cal D}(T) : \qquad 熱の水平拡散,$$ (3.37)

$${\cal D}(q) : \qquad 水蒸気の水平拡散,$$ (3.38)

$${\cal F}_\lambda \ (\varphi, \lambda, \sigma) : \qquad 小規模運動過程 (経度方向),$$ (3.39)

$${\cal F}_\varphi \ (\varphi, \lambda, \sigma) : \qquad 小規模運動過程 (緯度方向),$$ (3.40)

$$Q \ (\varphi, \lambda, \sigma) : \qquad 放射, 凝結, 小規模運動過程等による加熱・温度変化,$$ (3.41)

$$S_q \ (\varphi, \lambda, \sigma) : \qquad 凝結, 小規模運動過程等による水蒸気ソース,$$ (3.42)

$${\cal D}' \ (\Dvect{v}) : \qquad 摩擦熱.$$ (3.43)

各水平拡散(3.35)〜(3.38) に関しては3.2.7節で説明される. 定数は以下の通りである.

$$a : \quad 惑星半径 [\mathrm{m}],$$ (3.44)

$$R : \quad 乾燥大気の気体定数 [\mathrm{J\ kg}^{-1}\ \mathrm{K}^{-1}],$$ (3.45)

$$C_p : \quad 乾燥大気の大気定圧比熱 [\mathrm{J\ kg}^{-1}\ \mathrm{K}^{-1}],$$ (3.46)

$$f : \quad コリオリパラメータ [\mathrm{s}^{-1}],$$ (3.47)

$$\kappa \equiv R/C_p,$$ (3.48)

$$\epsilon_v : \quad 水蒸気分子量比.$$ (3.49)

3.2.6 境界条件

鉛直流に関する境界条件は

$$\dot{\sigma} = 0 \ \ \ at \ \ \sigma = 0 , \ 1 .$$ (3.50)

である. よって(3.1) から, 地表気圧の時間変化式と $\sigma$系での鉛直速度 $\dot{\sigma}$を求める診断式

$$\DP{\pi}{t} = - \int_{0}^{1} \Dvect{v}_{H} \cdot \nabla_{\sigma} \pi d \sigma - \int_{0}^{1} D d \sigma ,$$ (3.51)

$$\dot{\sigma} = - \sigma \DP{\pi}{t} - \int_{0}^{\sigma} D d \sigma - \int_{0}^{\sigma} \Dvect{v}_{H} \cdot \nabla_{\sigma} \pi d \sigma ,$$ (3.52)

が導かれる.

3.2.7 水平拡散とスポンジ層

水平拡散とスポンジ層における渦度と発散の散逸は次のように表現する.

$${\cal D}(\zeta) = {\cal D_{HD}}(\zeta) + {\cal D_{SL}}(\zeta)$$ (3.53)

$${\cal D}(D) = {\cal D_{HD}}(D) + {\cal D_{SL}}(D)$$ (3.54)

$${\cal D}(T) = {\cal D_{HD}}(T) + {\cal D_{SL}}(T)$$ (3.55)

$${\cal D}(q) = {\cal D_{HD}}(q)$$ (3.56)

ここで, ${\cal D_{HD}}$, ${\cal D_{SL}}$ はそれぞれ水平拡散とスポンジ層における 散逸を表す.

水平拡散項は, 次のように $\nabla^{N_D}$ の形で計算する.

$${\cal D_{HD}}(\zeta) = - K_{HD} \left[ (-1)^{N_D/2} \nabla^{N_D} - \left( \frac{2}{a^2} \right)^{N_D/2} \right] \zeta ,$$ (3.57)

$${\cal D_{HD}}(D) = - K_{HD} \left[ (-1)^{N_D/2} \nabla^{N_D} - \left( \frac{2}{a^2} \right)^{N_D/2} \right] D ,$$ (3.58)

$${\cal D_{HD}}(T) = - (-1)^{N_D/2} K_{HD} \nabla^{N_D} T ,$$ (3.59)

$${\cal D_{HD}}(q) = - (-1)^{N_D/2} K_{HD} \nabla^{N_D} q .$$ (3.60)

小さなスケールに選択的な水平拡散を表すため, 慣例として $N_D$ には 4$\sim$16 を用いることが多い.

スポンジ層における運動量の散逸項は, 東西平均成分を減衰させる場合とさせない場合の 2 通りの 計算法を導入する. 東西平均成分も減衰させる場合には,

$${\cal D_{SL}}(\zeta) = - \gamma_M \zeta,$$ (3.61)

$${\cal D_{SL}}(D) = - \gamma_M D,$$ (3.62)

となる. ここで, $\gamma_M$ はスポンジ層における運動量の減衰係数である. 東西平均成分を減衰させない場合には,

$${\cal D_{SL}}(\zeta) = - \gamma_M ( \zeta - \bar{\zeta} ),$$ (3.63)

$${\cal D_{SL}}(D) = - \gamma_M ( D - \bar{D} ),$$ (3.64)

となる. ここで, $\bar{}$ は, 東西平均を表す.

スポンジ層内の温度擾乱の減衰には以下の項を導入する.

$${\cal D_{SL}}(T) = - \gamma_H ( T - \bar{T} ),$$ (3.65)

ここで, $\gamma_H$ はスポンジ層における温度擾乱の減衰係数である.

減衰係数 $\gamma_M$, $\gamma_H$ の $\sigma$ 依存性に一般形はないが, dcpam では 下のような $\sigma$ 依存性を考慮する.

$$\gamma_M = \left\{ \begin{array}{ll} \gamma_{M,0} \left( \frac{\s... ...e \sigma_{lim}$)} \\ 0 . & \text{($\sigma > \sigma_{lim}$)} \end{array} \right.$$ (3.66)

$$\gamma_H = \left\{ \begin{array}{ll} \gamma_{H,0} \left( \frac{\s... ...e \sigma_{lim}$)} \\ 0 . & \text{($\sigma > \sigma_{lim}$)} \end{array} \right.$$ (3.67)

ここで, $\gamma_{M,0}$, $\gamma_{H,0}$, $N_{SL}$, $\sigma_{lim}$ はそれぞれ, $\sigma = \sigma_0$ における減衰係数, $\sigma$ 依存性の指数, スポンジ層の 下限の $\sigma$ である. dcpam では, $\sigma_0$ はモデル最上層の $\sigma$ としている.

3.3 鉛直離散化

ここでは支配方程式を鉛直方向に離散化する. Arakawa and Suarez(1983) に従って, (3.1)〜(3.6) を鉛直方向に差分によって離散化する. 各方程式の離散化表現は次のようになる.

3.3.1 連続の式, 鉛直速度

$$\DP{\pi}{t} = - \sum_{k=1}^{K} ( D_k + \Dvect{v}_k \cdot \nabla \pi ) \Delta \sigma_k,$$ (3.68)

$$\dot{\sigma}_{k-1/2} = - \sigma_{k-1/2} \DP{\pi}{t} - \sum_{l=k}^{K} ( D_l + \Dvect{v}_l \cdot \nabla \pi ) \Delta \sigma_l \qquad (k = 2, \cdots, K),$$ (3.69)

$$\dot{\sigma}_{1/2} = \dot{\sigma}_{K+1/2} = 0.$$ (3.70)

ここで,

$$\Dvect{v}_k \cdot \nabla \pi = \frac{U_k}{a (1-\mu^2)} \DP{\pi}{\lambda} + \frac{V_k}{a (1-\mu^2)} (1-\mu^2) \DP{\pi}{\mu}.$$ (3.71)

3.3.2 静水圧の式

$$\begin{align*}\begin{split}\Phi_{1} & = \Phi_{s} + C_{p} ( \sigma_{1}^{-\kappa} - 1 ) T_{v,1} \\ & = \Phi_{s} + C_{p} \alpha_{1} T_{v,1}. \end{split}\end{align*}$$ (3.72)

$$\begin{align*}\begin{split}\Phi_k - \Phi_{k-1} & = C_{p} \left[ \left( \frac{ \s... ... = C_{p} \alpha_k T_{v,k} + C_{p} \beta_{k-1} T_{v,k-1}. \end{split}\end{align*}$$ (3.73)

ここで,

$$\alpha_k = \left( \frac{ \sigma_{k-1/2} } { \sigma_k } \right)^{\kappa} -1 ,$$ (3.74)

$$\beta_k = 1- \left( \frac{ \sigma_{k+1/2} } { \sigma_k } \right)^{\kappa} ,$$ (3.75)

$$\Phi_{s} = gz_{s}$$ (3.76)

であり, $z_{s}$は地表面高度である.

3.3.3 運動方程式

$$\DP{\zeta_k}{t} = \Dinv{a} \left( \Dinv{1 - \mu^2} \DP{{V_A}_{,k}}{\lambda} - \DP{{U_A}_{,k}}{\mu} \right) + {\cal D}(\zeta_k),$$ (3.77)

$$\DP{D_k}{t} = \Dinv{a} \left( \Dinv{1 - \mu^2} \DP{{U_A}_{,k}}{\lambda} + \DP... ...t) - \nabla^{2}_{\sigma} ( \Phi_k + C_{p} \hat{\kappa}_k \overline{T}_k \pi + ( )_k ) + {\cal D}(D_k).$$ (3.78)

ここで,

$$\begin{align*}\begin{split}{U_A}_{,1} & = ( \zeta_1 + f ) V_1 - \frac{1}{2 \Delt... ... \DP{\pi}{\lambda} + {\cal F}_{\lambda, K} \cos \varphi, \end{split}\end{align*}$$ (3.79)

$$\begin{align*}\begin{split}{V_A}_{,1} & = - ( \zeta_1 + f ) U_1 - \frac{1}{2 \De... ...u^2) \DP{\pi}{\mu} + {\cal F}_{\varphi, K} \cos \varphi, \end{split}\end{align*}$$ (3.80)

$$\begin{align*}\begin{split}\hat{\kappa}_k & = \frac{ \sigma_{k-1/2}( \sigma^{\ka... ...\alpha_k + \sigma_{k+1/2} \beta_k } { \Delta \sigma_k }, \end{split}\end{align*}$$ (3.81)

$$% T_{v,k}' = T_{v,k} - \overline{T}_k,$$ (3.82)

$$% ( )_k = \frac{U^{2}_k + V^{2}_k}{2 (1-\mu^2)}.$$ (3.83)

3.3.4 熱力学の式

$$\begin{align*}\begin{split}\DP{T_k}{t} & = - \Dinv{a \cos \varphi} \left( \Dinv{... ...frac{Q_k}{C_{p}} + {\cal D}(T_k) + {\cal D}'(\Dvect{v}). \end{split}\end{align*}$$ (3.84)

ここで,

$$\begin{align*}\begin{split}H_k & \equiv T_k' D_k - \frac{1}{\Delta \sigma_k} [ \... ...la \pi ) \Delta \sigma_K \frac{T_{v,K}}{\Delta \sigma_K} \end{split}\end{align*}$$ (3.85)

であり,

$$\begin{align*}\begin{split}\hat{T}_{k-1/2} & = \frac{ \left[ \left( \displaystyl... ... K), \\ \hat{T}_{1/2} & = 0, \\ \hat{T}_{K + 1/2} & = 0, \end{split}\end{align*}$$ (3.86)

$$a_k = \alpha_k \left[ 1- \left( \frac{ \sigma_k }{ \sigma_{k-1} } \right)^{\kappa} \right]^{-1},$$ (3.87)

$$b_k = \beta_k \left[ \left( \frac{ \sigma_k }{ \sigma_{k+1} } \right)^{\kappa} - 1 \right]^{-1} .$$ (3.88)

3.3.5 水蒸気の式

$$\DP{q_k}{t} = - \Dinv{a} \left( \Dinv{1 - \mu^2} \DP{U_k q_k}{\lambda} + \DP{V_k q_k}{\mu} \right) + R_k + S_{q,k} + {\cal D}(q_k).$$ (3.89)

ここで,

$$\begin{align*}\begin{split}R_1 &= q_1 D_1 - \frac{1}{2 \Delta \sigma_1} \dot{\si... ...Delta \sigma_K} \dot{\sigma}_{K-1/2} ( q_{K-1} - q_K ) . \end{split}\end{align*}$$ (3.90)

3.4 水平離散化

ここでは支配方程式を水平離散化する. 水平方向の離散化はスペクトル変換法を用いる (Bourke, 1988). 非線形項は格子点上で計算する. 各方程式のスペクトル表現は以下のようになる. スペクトル表現に関する記号の意味については 2.5節を参照されたい. その詳細については第A章 を参照せよ. なお, 簡単化のため, 部分的に鉛直方向添字$k$を省略する.

3.4.1 連続の式

$$\DP{\tilde{\pi}_n^m}{t} = - \sum_{k=1}^{K} (\tilde{D}_n^m)_k \Delta \sigma_k + \frac{1}{I} \sum_{i=1}^{I} \sum_{j=1}^{J} Z_{ij} Y_n^{m *} ( \lambda_i, \mu_j ) w_j.$$ (3.91)

ここで,

$$Z \equiv - \sum_{k=1}^{K} \Dvect{v}_k \cdot \nabla \pi \Delta \sigma_k.$$ (3.92)

3.4.2 運動方程式

$$\begin{align*}\begin{split}\DP{\tilde{\zeta}_n^m}{t} & = \frac{1}{I} \sum_{i=1}^... ... \\ & \quad + \tilde{\cal D}_{M,n}^m \tilde{\zeta}_n^m , \end{split}\end{align*}$$ (3.93)

$$\begin{align*}\begin{split}\DP{\tilde{D}_n^m}{t} & = \frac{1}{I} \sum_{i=1}^{I} ... ...{T}_k \pi_n^m ) + \tilde{\cal D}_{M,n}^m \tilde{D}_n^m . \end{split}\end{align*}$$ (3.94)

ここで,

$$\tilde{\cal D}_{M,n}^m = - K_{HD} \left[ \left( \frac{-n(n+1)}{a^{2}} \right)^{N_D/2} - \left( \frac{2}{a^2} \right)^{N_D/2} \right] - \tilde{\gamma}_{M,k,n}^m ,$$ (3.95)

$$\tilde{\gamma}_{M,k,n}^m = \left\{ \begin{array}{ll} \tilde{\gamma}_{M,0,n}^m \left( \frac... ...text{($k \ge k_{SLlim}$)} \\ 0 . & \text{($k < k_{SLlim}$)} \end{array} \right.$$ (3.96)

ここで, $k_{SLlim}$ はスポンジ層を適応する下限の $k$ である. また, スポンジ層において東西平均成分も減衰させる場合には, $\tilde{\gamma}_{M,0,n}^m = \gamma_{M,0}$ であり, 東西平均成分を 減衰させない場合には,

$$\tilde{\gamma}_{M,0,n}^m = \left\{ \begin{array}{ll} \gamma_{M,0}, & \text{($m \ne 0$)} \\ 0 , & \text{($m = 0$)} \end{array} \right.$$ (3.97)

である.

3.4.3 熱力学の式

$$\begin{align*}\begin{split}\DP{\tilde{T}_n^m}{t} & = - \frac{1}{I} \sum_{i=1}^{I... ...D}'_{ij}(\Dvect{v}) Y_n^{m *} ( \lambda_i, \mu_j ) w_j . \end{split}\end{align*}$$ (3.98)

ここで,

$$\tilde{\cal D}_{H,n}^m = - K_{HD} \left( \frac{-n(n+1)}{a^{2}} \right)^{N_D/2} - \tilde{\gamma}_{H,k,n}^m .$$ (3.99)

$$\tilde{\gamma}_{H,k,n}^m = \left\{ \begin{array}{ll} \tilde{\gamma}_{H,0,n}^m \left( \frac... ...text{($k \ge k_{SLlim}$)} \\ 0 , & \text{($k < k_{SLlim}$)} \end{array} \right.$$ (3.100)

$$\tilde{\gamma}_{H,0,n}^m = \left\{ \begin{array}{ll} \gamma_{H,0}, & \text{($m \ne 0$)} \\ 0 , & \text{($m = 0$)} \end{array} \right.$$ (3.101)

である.

3.4.4 水蒸気の式

$$\begin{align*}\begin{split}\DP{\tilde{q}_n^m}{t} & = - \frac{1}{I} \sum_{i=1}^{I... ... w_j \\ & \quad + \tilde{\cal D}_{q,n}^m \tilde{q}_n^m . \end{split}\end{align*}$$ (3.102)

ここで,

$$\tilde{\cal D}_{q,n}^m = - K_{HD} \left( \frac{-n(n+1)}{a^{2}} \right)^{N_D/2}$$ (3.103)

である.

3.5 時間積分

ここでは時間積分スキームについて記す.

時間差分には, 複数の方法を組み合わせて用いる. 用いる方法の 概要を以下に示す.

• 力学過程
  • 水平拡散およびスポンジ層における減衰項には, 後方差分を用いる.
  • その他の項には, leap frog 法と Crank-Nicolson 法を組み合わせた semi-implicit 法 (Bourke, 1988) を用いる.
• 物理過程
  • 予報型の物理過程には, 前方差分を用いる.
  • 調節型の物理過程は, semi-implicit 法での力学過程積分後に計算された値を 用いて計算する.
• 時間フィルタ
  • 力学過程, 物理過程のすべての計算後に, 力学過程で用いている leap frog 法を起源とする計算モード抑制のための時間フィルター (Asselin, 1972) を適応する.

この方法は, 予報変数を ${\cal A}$ と表すと, 以下の 3 式で表現される.

$$\frac{ \hat{\cal A}^{t+\Delta t} - \bar{\cal A}^{t-\Delta t} }{ 2... ...+\Delta t} \right) \right\} + \dot{\cal A}_{dyn,NG} \left( {\cal A}^{t} \right)$$

$$+ \dot{\cal A}_{dyn,dis }\left( \hat{\cal A}^{t+\Delta t} \right) + \dot{\cal A}_{phy,pred}\left( \bar{\cal A}^{t-\Delta t} \right) ,$$ (3.104)

$${\cal A}^{t+\Delta t} = \hat{\cal A}^{t+\Delta t} + 2 \Delta t \d... ...) + 2 \Delta t \dot{\cal A}_{phy,adj}\left( \hat{\cal A}^{t+\Delta t} \right) ,$$ (3.105)

$$\bar{\cal A}^{t} = {\cal A}^{t} + \epsilon_f \left( \bar{\cal A}^{t-\Delta t} - 2 {\cal A}^{t} + {\cal A}^{t+\Delta t} \right).$$ (3.106)

ここで, $\dot{\cal A}_{dyn,G}$, $\dot{\cal A}_{dyn,NG}$ はそれぞれ, 力学過程において semi-implicit 法で分離された重力波項 (線型項) と非重力波項 (非線型項), $\dot{\cal A}_{dyn,dis}$ は水平拡散とスポンジ層における減衰項, $\dot{\cal A}_{phy,pred}$ は予報型の物理過程項である. $\dot{\cal A}_{fric}$, $\dot{\cal A}_{phy,adj}$ は, それぞれ 摩擦熱による加熱項および調節型の物理過程項である. $\epsilon_f$ は時間フィルタの係数であり, dcpam での標準値は 0.05 としている.

3.5.1 力学過程の方程式系の時間差分式

まず, semi-implicit 法を用いるために, 方程式系を $T=\overline{T}_k$ である 静止場に基づいて線形重力波項とそれ以外の項に分離する. 鉛直方向のベクトル表現 $\Dvect{A}=\{ A_{k} \}$, および行列表現 $\underline{A}=\{ A_{kl} \}$ を用いると, 連続の式, 発散方程式, 熱力学の式は,

$$\DP{\tilde{\pi}^{m}_{n}}{t} = \left( \DP{\tilde{\pi}^{m}_{n}}{t} \right)^{\rm NG} - \Dvect{C} \cdot \tilde{\Dvect{D}}^{m}_{n} ,$$ (3.107)

$$\DP{\tilde{\Dvect{D}}^{m}_{n}}{t} = \left( \DP{\tilde{\Dvect{D}}^... ...^{m}_{n} ) + \underline{ \tilde{\cal D}_M }_{n}^{m} \tilde{\Dvect{D}}^{m}_{n} ,$$ (3.108)

$$\DP{\tilde{\Dvect{T}}^{m}_{n}}{t} = \left( \DP{\tilde{\Dvect{T}}^{m}_{n}}{t} \right)^{\rm NG} - \un... ...}}^{m}_{n} + \underline{ \tilde{\cal D}_{H} }_{n}^{m} \tilde{\Dvect{T}}^{m}_{n}$$ (3.109)

となる ^3.1. $\widetilde{( \hspace{0.5cm} )}^{m}_{n}$や $\widetilde{[ \hspace{0.5cm} ]}^{m}_{n}$ といった表記については 2.5節の (2.10), (2.15), (2.17) を参照のこと. ここで, 添字 NG の付いた項は, 非重力波項であり, 以下のように表される.

$$\left( \DP{\tilde{\pi}^{m}_{n}}{t} \right)^{\rm NG} = \tilde{Z}^{m}_{n},$$ (3.110)

$$\begin{align*}\begin{split}\left( \DP{\tilde{D}^{m}_{k,n}}{t} \right)^{\rm NG} &... ...um_{l=1}^{K} W_{kl} ( T_{v,l}-T_{l} ) \right] }^{m}_{n}, \end{split}\end{align*}$$ (3.111)

$$\begin{align*}\begin{split}\left( \DP{\tilde{T_{k}}^{m}_{,n}}{t} \right)^{\rm NG... ...& \qquad + \widetilde{ \left[ H_{ijk} \right]^{m}_{n} }. \end{split}\end{align*}$$ (3.112)

各項は以下の通りである. 簡単化のため経度, 緯度方向添字 $i,j$ の 表記を省略する.

$$Z = - \sum_{k=1}^{K} \Dvect{v}_{k} \cdot \nabla \pi \Delta \sigma_{k},$$ (3.113)

$$\begin{align*}\begin{split}H_k & = T_{k}^{\prime} D_{k} \\ & \quad - \frac{1}{\D... ...lta \sigma_{K} + T'_{v,K} D_K \Delta \sigma_{K} \right], \end{split}\end{align*}$$ (3.114)

$$\begin{align*}\begin{split}\dot{\sigma}^{\rm NG}_{k-1/2} &= - \sigma_{k-1/2} \le... ...k}^{K} \Dvect{v}_{l} \cdot \nabla \pi \Delta \sigma_{l}, \end{split}\end{align*}$$ (3.115)

$$\hat{T}_{k-1/2}' = \left\{ \begin{array}{ll} 0 , & \text{($k = 1$... ..., & \text{($k = 2, \cdots, K$)} \\ 0 , & \text{($k = K+1$)} \end{array} \right.$$ (3.116)

$$\hat{\overline{T}}_{k-1/2} = \left\{ \begin{array}{ll} 0 , & \tex... ..., & \text{($k = 2, \cdots, K$)} \\ 0 . & \text{($k = K+1$)} \end{array} \right.$$ (3.117)

また, 重力波項のベクトルおよび行列は以下のとおりである.

$$C_{k} = \Delta \sigma_{k} ,$$ (3.118)

$$W_{kl} = C_{p} \alpha_{l} \delta_{k \geq l} + C_{p} \beta_{l} \delta_{k-1 \geq l} ,$$ (3.119)

$$G_{k} = \hat{\kappa}_{k} C_{p} \overline{T}_{k} ,$$ (3.120)

$$\underline{h} = \underline{Q}\underline{S} - \underline{R} ,$$ (3.121)

$$Q_{kl} = \frac{1}{\Delta \sigma_{k}} ( \hat{\overline{T}}_{k-1/2} - \ove... ... \sigma_{k}} ( \overline{T}_{k} - \hat{\overline{T}}_{k+1/2} ) \delta_{k+1=l} ,$$ (3.122)

$$S_{kl} = \sigma_{k-1/2} \Delta \sigma_{l} - \Delta \sigma_{l} \delta_{k \leq l } ,$$ (3.123)

$$R_{kl} = - \left( \frac{ \alpha_{k} }{ \Delta \sigma_{k} } \Delta \sigma... ...a \sigma_{k} } \Delta \sigma_{l} \delta_{k+1 \leq l} \right) \overline{T}_{k} ,$$ (3.124)

$$(\tilde{ {\cal D}_{M,kl} } )_n^m = - K_{HD} \left[ \left( \frac{-n(n+1)}{a^{2}} \right)^{N_D/2} - \left( \frac{2}{a^2} \right)^{N_D/2} \right] \delta_{k=l}$$

$$\hspace{1cm} - \gamma_{M,0,n}^m \left( \frac{\sigma_k}{\sigma_K} \right)^{N_{SL}} \delta_{k=l} \delta_{k \ge k_{SLlim}} .$$ (3.125)

$$(\tilde{ {\cal D}_{H,kl} } )_n^m = - K_{HD} \left( \frac{-n(n+1)}{a^{2}} \right)^{N_D/2} \delta_{k=l}$$

$$\hspace{1cm} - \gamma_{H,0,n}^m \left( \frac{\sigma_k}{\sigma_K} \right)^{N_{SL}} \delta_{k=l} \delta_{k \ge k_{SLlim}} .$$ (3.126)

$\delta_{k \leq l}$ は, $k \leq l$ が成り立つとき 1, そうでないとき 0 となる関数である.

なお, 渦度方程式には線型重力波項がないため, ここでは示さない. ^3.2

これらの方程式に,

• 水平拡散とスポンジ層における減衰項には後退差分
• その他の項には, leap frog 法と中心差分を組み合わせた semi-implicit 法

を適応すると,

$$\delta_{t} \tilde{\pi}^{m}_{n} = \left( \DP{\tilde{\pi}^{m}_{n}}{t} \right)^{\rm NG} - \Dvect{C} \cdot \overline{ \tilde{\Dvect{D}}^{m}_{n} }^{t} ,$$ (3.127)

$$\delta_{t} \tilde{\Dvect{D}}^{m}_{n} = \left( \DP{\tilde{\Dvect{D}}^{m}_{n}}{t} \right)^{\rm NG} - \le... ...\underline{ \tilde{\cal D}_{M} }_{n}^{m} \tilde{\Dvect{D}}^{m,t+\Delta t}_{n} ,$$ (3.128)

$$\delta_{t} \tilde{\Dvect{T}}^{m}_{n} = \left( \DP{\tilde{\Dvect{T... ... \underline{ \tilde{\cal D}_{H} }_{n}^{m} \tilde{\Dvect{T}}^{m,t+\Delta t}_{n}.$$ (3.129)

となる. ただし,

$$\delta_{t} {\cal A} \equiv \frac{1}{2 \Delta t} \left( {\cal A}^{t+\Delta t} - {\cal A}^{t-\Delta t} \right) ,$$ (3.130)

$$\overline{\cal A}^{t} \equiv \frac{1}{2} \left( {\cal A}^{t+\Delta t} + {\cal A}^{t-\Delta t} \right) = {\cal A}^{t-\Delta t} + \delta_{t} {\cal A} \Delta t .$$ (3.131)

である.

(3.127), (3.128), (3.129) より, $\overline{\tilde{\Dvect{D}}^{m}_{n}}^{t}$ について整理すると,

$$\begin{align*}\begin{split}& \Biggl[ ( \underline{I}-2\Delta t \underline{ \tild... ...tilde{\pi}^{m}_{n}}{t} \right)^{\rm NG} \right\} \right] \end{split}\end{align*}$$ (3.132)

となる. ここで $\underline{I}$は単位行列, $\Dvect{C}^{T}$は$\Dvect{C}$の 転置ベクトルである. (3.132) を $\overline{\tilde{\Dvect{D}}^{m}_{n}}^{t}$ について解き,

$$\tilde{\Dvect{D}}^{m,t+\Delta t}_{n} = 2\overline{\tilde{\Dvect{D}}^{m}_{n}}^{t} - \tilde{\Dvect{D}}^{m,t-\Delta t}_{n}$$ (3.133)

および, (3.127), (3.129) により $\hat{\cal A}^{t+\Delta t}$ が求められる.

3.6 参考文献

Arakawa, A., Suarez, M. J., 1983: Vertical differencing of the primitive equations in sigma coordinates. Mon. Wea. Rev., 111, 34-35.

Asselin, R. A., 1972: Frequency filter for time integrations. Mon. Wea. Rev., 100, 487-490.

Bourke, W.P., 1988: Spectral methods in global climate and weather prediction models. Physically-Based Modelling and Simulation of Climates and Climatic Change. Part I., M.E. Schlesinger (ed.), Kluwer Academic Publishers, Dordrecht, 169-220.

Haltiner, G.J., Williams, R.T., 1980: Numerical Prediction and Dynamic Meteorology (2nd ed.). John Wiley & Sons, 477pp.

石岡 圭一, 2004: スペクトル法による数値計算入門 . 東京大学出版会, 232pp.

... となる^3.1

念のため注記しておくと, $\tilde{\Dvect{\Phi}}^{m}_{s,n} = \left( \tilde{\Phi}^{m}_{s,n}, \tilde{\Phi}^{m}_{s,n}, \cdots, \tilde{\Phi}^{m}_{s,n}\right)$ である.

... ここでは示さない.^3.2

ここは本当は方程式を書くべきだろう. 後で書く. (YOT, 2009/10/11)