概要
スペクトル法を用いた数値解析法に興味を持っています。特に、昔から気象予報において使われております。今回の記事では、気象分野で使われるピリティブ方程式について紹介したいと思います。
乾燥気体における運動方程式、連続の式、熱力学方程式に基づいて紹介します。水蒸気の方程式は勉強したいです。
支配方程式
乾燥気体における運動方程式、連続の式、熱力学方程式は、
\begin{align}
&\frac{\mathrm{d} u_i }{\mathrm{d}t } + 2\varepsilon_{ijk}\Omega_ju_k+ \frac{1}{\rho}\frac{\partial P}{\partial x_i} + \frac{\partial \phi}{\partial x_i} =F_i \\
&\frac{\mathrm{d} \rho }{\mathrm{d}t } + \rho \frac{\partial u_i}{\partial x_i} =0 \\
&C_p \frac{T}{\theta}\frac{\mathrm{d} \theta}{\mathrm{d}t } = Q
\end{align}
と書ける。$u_i$は風速であり、$\rho$は密度、$\theta$は熱ポテンシャル（温位）である。また下付き添字は$x,y,z$の座標をとる。 時間微分には移流項が含まれている。
\begin{align}
\frac{\mathrm{d} }{\mathrm{d}t } = \frac{\partial }{\partial t } + u_i \frac{\partial }{\partial x_i }
\end{align}
運動方程式の右辺２項目は、地球の回転による力（コリオリ力）で$\Omega$は地球の角速度である。３項目は圧力$P$による傾度力、４項目はGeopotentialであり重力加速度を$g$とすれば
\begin{align}
\frac{\partial \phi}{\partial x_i} = g k_{i}
\end{align}
と書ける。また$k$は$k =(0,0,1)$のベクトルである。熱力学方程式にでてくる温位は、圧力$P$の空気を乾燥断熱的に基準圧力$\bar{P}$まで下げた（上げた）ときの空気の温度を表し
\begin{align}
\theta = T\left(\frac{\bar{P}}{P} \right)^{R/C_p}
\end{align}
と書ける。温位は、乾燥断熱変化をしているならば保存される量である。
以上の方程式から$u_i$,$\rho$,$\theta$は求まるが、圧力$P$と温度$T$が未知数である。そのため、理想気体の状態方程式を使う。
理想気体においては、
\begin{align}
P = \rho R T
\end{align}
が成立するため、定圧比熱$C_p$と定積比熱$C_v$の関係
\begin{align}
C_p &= C_v + R \\
1-\frac{R}{C_p}&=\frac{C_p-R}{C_p}= \frac{C_v}{C_p}\\
\frac{C_p}{C_v} &=1 +\frac{R}{C_v}
\end{align}
を用いれば、温位の式は
\begin{align}
&\theta = \frac{P}{\rho R }\left(\frac{\bar{P}}{P} \right)^{R/C_p}
= \frac{\bar{P}}{\rho R }\left(\frac{\bar{P}}{P} \right)^{R/C_p-1}
=\frac{\bar{P}}{\rho R }\left(\frac{\bar{P}}{P} \right)^{-C_v/C_p} \\
&\therefore P = \bar{P} \left(\frac{\rho R \theta }{\bar{P}}\right)^{C_v/C_p}
\end{align}
そして、
\begin{align}
T &= \frac{P}{R \rho} = \frac{\bar{P}}{R \rho} \left(\frac{\rho R \theta }{\bar{P}}\right)^{C_v/C_p}
= \theta \left(\frac{\rho R \theta }{\bar{P}}\right)^{C_v/C_p-1} \\
&= \theta \left(\frac{\rho R \theta }{\bar{P}}\right)^{R/C_v}
\end{align}
と書けるので、未知数である圧力$P$と温度$T$が温位と密度から求まる。
シグマ座標系
プリミティブ方程式は静力学平衡を仮定するため、気圧$p$は高さ$z$とともに指数関数的に減少する。そのため、鉛直座標は$z$ではなく気圧座標系$p$が用いられる。ただし、ヒマラヤ山脈やロッキー山脈などは対流圏中層の高さがあり、気圧座標系を用いるのは難しい。
そのため、$P^*$を表面圧力としてシグマ座標系
\begin{align}
\sigma=\frac{P}{P^*}
\end{align}
と設定する。$\sigma=1$のとき地表面に、$\sigma=0$は大気上端$p=0$に対応する。この時間微分の境界条件は$\sigma=1$と$\sigma=0$のとき
\begin{align}
\dot{\sigma}=0
\end{align}
と設定する。
プリミティブ方程式
プリミティブ方程式は、静力学平衡を仮定して
\begin{align}
\frac{1}{\rho}\frac{\partial P}{\partial \sigma} = -\frac{\partial \phi}{\partial \sigma} = -g
\end{align}
支配方程式から緯度に依存しないコリオリ力を無視すると以下の方程式が得られる。
\begin{align}
&\frac{\partial u_i }{\partial t} + u_j\frac{\partial u_i}{\partial x_j}\Bigg|_{\sigma}+ \dot{\sigma}\frac{\partial u_i}{\partial \sigma}
+ f\varepsilon_{ijl}k_{j}u_l+ \frac{1}{\rho}\frac{\partial P}{\partial x_i}\Bigg|_{\sigma} + \frac{\partial \phi }{\partial x_i}\Bigg|_{\sigma}=F_i \\
&\frac{\partial \rho}{\partial t} + u_j\frac{\partial \rho}{\partial x_j}\Bigg|_{\sigma}+ \dot{\sigma}\frac{\partial \rho}{\partial \sigma}+\rho \frac{\partial u_i}{\partial x_i}\Bigg|_{\sigma}+\rho \frac{\partial \dot{\sigma}}{\partial \sigma}
=0 \\
&C_p \frac{T}{\theta}\left(\frac{\partial \theta}{\partial t} + u_j\frac{\partial \theta}{\partial x_j}\Bigg|_{\sigma}+ \dot{\sigma}\frac{\partial \theta}{\partial \sigma}\right)= Q
\end{align}
下付き添字は$x,y,z$の座標をとるが、風速の$z$成分$u_z$は0とする。また、風速の$z$成分の方程式は、静力学平衡の式になっている。$\big|_{\sigma}$ は$xy$平面における微分を表している。また、$f$はコリオリパラメータである。
\begin{align}
f = 2\Omega\sin\varphi
\end{align}
$\varphi$は緯度である。以下では便利上
\begin{align}
\frac{\mathrm{d}}{\mathrm{d} t} = \frac{\partial }{\partial t} + u_j\frac{\partial }{\partial x_j}\Bigg|_{\sigma}+ \dot{\sigma}\frac{\partial }{\partial \sigma}
\end{align}
と記述する。
シグマ座標を使い書き直すと静力学平衡の式は、理想気体の状態方程式を使い
\begin{align}
\frac{1}{\rho}\frac{\partial P}{\partial \sigma} + \frac{\partial \phi}{\partial \sigma} &= 0 \\
\therefore \frac{RT}{\sigma P^*}\frac{\partial (\sigma P^*)}{\partial \sigma} + \frac{\partial \phi}{\partial \sigma} &= 0 \\
\therefore \frac{RT}{\sigma} + \frac{\partial \phi}{\partial \sigma} &= 0 \\
\end{align}
圧力勾配は、
\begin{align}
\frac{1}{\rho}\frac{\partial P}{\partial x_i}\Bigg|_{\sigma} =
\frac{RT}{\sigma P^*}\frac{\partial (\sigma P^*)}{\partial x_i}\Bigg|_{\sigma} = \frac{RT}{P^*}\frac{\partial P^*}{\partial x_i}\Bigg|_{\sigma}
\end{align}
温位の微分は、
\begin{align}
\frac{\mathrm{d} \theta}{\mathrm{d} t} =& \frac{\mathrm{d} }{\mathrm{d} t} \left\{ T\left(\frac{\bar{P}}{P} \right)^{R/C_p} \right\}\\
=& \frac{\theta}{T} \frac{\mathrm{d} T}{\mathrm{d} t} +
T\bar{P}^{R/C_p} \left(-\frac{R}{C_p}\right) P^{-R/C_p-1}\frac{\mathrm{d} P}{\mathrm{d} t} \\
=& \frac{\theta}{T} \frac{\mathrm{d} T}{\mathrm{d} t}-\frac{R\theta}{P C_p} \frac{\mathrm{d} P}{\mathrm{d} t}
\end{align}
とかけ、新しい変数とパラメータを
\begin{align}
\omega &= \frac{\mathrm{d} P}{\mathrm{d} t} \ , \ \kappa = \frac{R}{C_p}
\end{align}
とすれば、温位の方程式は温度に関する方程式
\begin{align}
\frac{\mathrm{d} T}{\mathrm{d} t}-\kappa T\frac{\omega}{P} = \frac{Q}{C_p}
\end{align}
となる。さらに、密度は
\begin{align}
\rho = - P^*\left(\frac{\partial \phi}{\partial \sigma}\right)^{-1}
\end{align}
と書ける。$\phi$は時間依存しないので連続の式は
\begin{align}
&\frac{1}{\rho}\frac{\mathrm{d} \rho}{\mathrm{d} t}+ \frac{\partial u_i}{\partial x_i}\Bigg|_{\sigma}+ \frac{\partial \dot{\sigma}}{\partial \sigma} = 0 \\
\therefore & \frac{1}{P^*}\frac{\mathrm{d} P^*}{\mathrm{d} t} + \frac{\partial u_i}{\partial x_i}\Bigg|_{\sigma}+ \frac{\partial \dot{\sigma}}{\partial \sigma}= 0 \\
\therefore & \frac{\partial P^*}{\partial t}+u_j\frac{\partial P^*}{\partial x_j}\Bigg|_{\sigma}+ \dot{\sigma}\frac{\partial P^*}{\partial \sigma}+P^*\frac{\partial u_i}{\partial x_i}\Bigg|_{\sigma}+P^*\frac{\partial \dot{\sigma}}{\partial \sigma} = 0
\\
\therefore & \frac{\partial P^*}{\partial t}+\frac{\partial u_jP^*}{\partial x_j}\Bigg|_{\sigma}+P^*\frac{\partial \dot{\sigma}}{\partial \sigma} = 0
\end{align}
となる。上記より３つの方程式は、
\begin{align}
\frac{\partial u_i }{\partial t} &= -u_j\frac{\partial u_i}{\partial x_j}\Bigg|_{\sigma}- \dot{\sigma}\frac{\partial u_i}{\partial \sigma}
- f\varepsilon_{ijl}k_{j}u_l- \frac{RT}{P^*}\frac{\partial P^*}{\partial x_i}\Bigg|_{\sigma} - \frac{\partial \phi }{\partial x_i}\Bigg|_{\sigma}+F_i \\
\frac{\partial P^*}{\partial t} &= -\frac{\partial u_jP^*}{\partial x_j}\Bigg|_{\sigma}-P^*\frac{\partial \dot{\sigma}}{\partial \sigma}
\\
\frac{\partial T}{\partial t} &=-u_j\frac{\partial T}{\partial x_j}\Bigg|_{\sigma}- \dot{\sigma}\frac{\partial T}{\partial \sigma}+\kappa T\frac{\omega}{P} + \frac{Q}{C_p}
\end{align}
となり、水平風速と表面圧力、温度の方程式となる。
鉛直方向への積分
連続の式と静力学平衡の式を鉛直方向に積分を行う。まず、連続の式を$P^*$で割り、鉛直方向に積分すると
\begin{align}
& \int_0^{\sigma} \mathrm{d}\hat{\sigma} \left\{\frac{1}{P^*}\frac{\partial P^*}{\partial t}+\frac{1}{P^*}\frac{\partial u_jP^*}{\partial x_j}\Bigg|_{\sigma}+\frac{\partial \dot{\sigma}}{\partial \hat{\sigma}} \right\}= 0 \\
\therefore & \frac{\sigma}{P^*} \frac{\partial P^*}{\partial t}+\int_0^{\sigma} \mathrm{d}\hat{\sigma}\frac{1}{P^*}\frac{\partial u_jP^*}{\partial x_j}\Bigg|_{\sigma}+\left(\dot{\sigma}(\sigma) -\dot{\sigma}(0) \right) = 0 \\
\end{align}
となる。境界条件から$\dot{\sigma}(0)=0$なので最後の項は消える。$\sigma=1$のとき、境界条件から$\dot{\sigma}(1)=0$で、表面圧力傾度方程式
\begin{align}
\frac{1}{P^*}\frac{\partial P^*}{\partial t}+\int_0^{1} \mathrm{d}\hat{\sigma}\frac{1}{P^*}\frac{\partial u_jP^*}{\partial x_j}\Bigg|_{\sigma} = 0 \\
\end{align}
が成立する。シグマ座標における鉛直速度は、
\begin{align}
\dot{\sigma}(\sigma) &= -\frac{\sigma}{P^*} \frac{\partial P^*}{\partial t}-\int_0^{\sigma} \mathrm{d}\hat{\sigma}\frac{1}{P^*}\frac{\partial u_jP^*}{\partial x_j}\Bigg|_{\sigma} \\
&= \sigma\int_0^{1} \mathrm{d}\hat{\sigma}\frac{1}{P^*}\frac{\partial u_jP^*}{\partial x_j}\Bigg|_{\sigma} -\int_0^{\sigma} \mathrm{d}\hat{\sigma}\frac{1}{P^*}\frac{\partial u_jP^*}{\partial x_j}\Bigg|_{\sigma}
\end{align}
と書ける。 また、一般的な鉛直速度は、上記の式を使い
\begin{align}
\omega &= \frac{\mathrm{d} P}{\mathrm{d} t} = \frac{\mathrm{d} \sigma P^*}{\mathrm{d} t} = \sigma\frac{\mathrm{d} P^*}{\mathrm{d} t}+ P^*\frac{\mathrm{d} \sigma}{\mathrm{d} t} \\
&= \sigma\frac{\partial P^*}{\partial t}+\sigma u_j \frac{\partial P^*}{\partial x_j}+P^*\dot{\sigma}\\
&= -P^*\int_0^{\sigma} \mathrm{d}\hat{\sigma}\frac{1}{P^*}\frac{\partial u_jP^*}{\partial x_j}\Bigg|_{\sigma} +\sigma u_j \frac{\partial P^*}{\partial x_j} \\
\therefore \frac{\omega}{P} &= \frac{\omega}{\sigma P^*} = -\frac{1}{\sigma} \int_0^{\sigma} \mathrm{d}\hat{\sigma}\frac{1}{P^*}\frac{\partial u_jP^*}{\partial x_j}\Bigg|_{\sigma} +\frac{1}{P^*} u_j \frac{\partial P^*}{\partial x_j}
\end{align}
と書ける。
次に、静力学平衡の式を鉛直方向に積分すると
\begin{align}
& \int_{\sigma}^1 \mathrm{d}\hat{\sigma} \left\{ \frac{RT}{\hat{\sigma}} + \frac{\partial \phi}{\partial \hat{\sigma}}\right\}= 0 \\
\therefore & \int^{\sigma}_1 \mathrm{d}\hat{\sigma} \frac{RT}{\hat{\sigma}} +\phi(1)- \phi(\sigma) = 0
\end{align}
が得られる。
最後に、表面圧力を基準圧力$\bar{P}$を使い
\begin{align}
\tilde{\pi} = \ln \frac{P^*}{\bar{P}}
\end{align}
と表記すれば
\begin{align}
\frac{\partial \tilde{\pi}}{\partial x_j} = \frac{1}{P^*}\frac{\partial P^*}{\partial x_j}
\end{align}
と書ける。上記をまとめると、以下の方程式が得られる。
予報方程式
時間微分を含んだ方程式を予報方程式と呼び、プリミティブ方程式系では以下となる。しかし、添字$i$に関しては$z$を含まない。
\begin{align}
\frac{\partial u_i }{\partial t} &= -u_j\frac{\partial u_i}{\partial x_j}\Bigg|_{\sigma}- \dot{\sigma}\frac{\partial u_i}{\partial \sigma}
- f\varepsilon_{ijl}k_{j}u_l- RT\frac{\partial \tilde{\pi}}{\partial x_i} \Bigg|_{\sigma} - \frac{\partial \phi }{\partial x_i}\Bigg|_{\sigma}+F_i \\
\frac{\partial T}{\partial t} &=-u_j\frac{\partial T}{\partial x_j}\Bigg|_{\sigma}- \dot{\sigma}\frac{\partial T}{\partial \sigma}+\kappa T\frac{\omega}{P} + \frac{Q}{C_p}\\
\frac{\partial \tilde{\pi}}{\partial t} &= -\int_0^{1} \mathrm{d}\hat{\sigma}A
\end{align}
診断方程式
時間微分を含まない方程式を予報方程式と呼び、プリミティブ方程式系では以下となる。
\begin{align}
A &= \frac{1}{P^*}\frac{\partial u_jP^*}{\partial x_j}\Bigg|_{\sigma}\\
\phi(\sigma) &= \phi(1)+\int_{\sigma}^1 \mathrm{d}\hat{\sigma} \frac{RT}{\hat{\sigma}} \\
\frac{\omega}{P} & = -\frac{1}{\sigma} \int_0^{\sigma} \mathrm{d}\hat{\sigma}A+ u_j \frac{\partial \tilde{\pi}}{\partial x_j} \\
\dot{\sigma}(\sigma) &= \sigma\int_0^{1} \mathrm{d}\hat{\sigma}A -\int_0^{\sigma} \mathrm{d}\hat{\sigma}A
\end{align}
渦度方程式と発散方程式
上記の予報方程式を解けば、乾燥大気のシミュレーションが可能であるが、座標はデカルト座標ではなく（地球を球体としているので）球座標を用いるため、物理量は座標依存しない量が望ましい。従って、座標系に依存する水平風速$u_i$を求めるのではなく、スカラー量である渦度や発散を求める。
相対渦度$\zeta$は、水平風速$u_i$によって
\begin{align}
\zeta = \varepsilon_{ijl} k_i\frac{\partial u_l}{\partial x_j}\Bigg|_{\sigma}
\end{align}
発散$\delta$は、
\begin{align}
\delta = \frac{\partial u_i}{\partial x_i}\Bigg|_{\sigma}
\end{align}
と書ける。
ここで、水平風速の方程式を渦度と発散の方程式に書き換えるため、ベクトル計算の準備を行う。まず、3次元ベクトル${\bf a}$ に対して、$i$成分は
\begin{align}
\{(\nabla \times {\bf a})\times {\bf a} \}_i &= \varepsilon_{ijl}\{(\nabla \times {\bf a})\}_j a_l = \varepsilon_{ijl}\varepsilon_{jmn}\frac{\partial a_n}{\partial x_m} a_l \\
& = (\delta_{in}\delta_{lm}-\delta_{im}\delta_{ln})\frac{\partial a_n}{\partial x_m} a_l \\
& = \frac{\partial a_i}{\partial x_l} a_l - \frac{\partial a_l}{\partial x_i} a_l \\
& = a_l\frac{\partial a_i}{\partial x_l} -\frac{1}{2}\frac{\partial (a_l a_l)}{\partial x_i}
\end{align}
である。ベクトル${\bf a}$の$z$成分がゼロの場合、この式と
\begin{align}
\{{\bf k}\cdot (\nabla \times {\bf a})({\bf k} \times {\bf a})\}_i = \varepsilon_{lmn} k_l \frac{\partial a_n}{\partial x_m} \varepsilon_{ijp} k_ja_p
\end{align}
は等価になる。2成分しかないので直接代入して確認すると、$x$成分は
\begin{align}
\{(\nabla \times {\bf a})\times {\bf a} \}_x
& = a_x\frac{\partial a_x}{\partial x}+a_y\frac{\partial a_x}{\partial y} -\frac{\partial a_x }{\partial x} a_x-\frac{\partial a_y }{\partial x} a_y\\
& =a_y\frac{\partial a_x}{\partial y} -\frac{\partial a_y }{\partial x} a_y \\
\{{\bf k}\cdot (\nabla \times {\bf a})({\bf k} \times {\bf a})\}_x &= \left( \varepsilon_{zxy} \frac{\partial a_y}{\partial x} + \varepsilon_{zyx} \frac{\partial a_x}{\partial y}\right)\varepsilon_{xzy} a_y \\
& = -\frac{\partial a_y}{\partial x}a_y + \frac{\partial a_x}{\partial y} a_y
\end{align}
$y$成分は
\begin{align}
\{(\nabla \times {\bf a})\times {\bf a} \}_y
& = a_x\frac{\partial a_y}{\partial x}+a_y\frac{\partial a_y}{\partial y} -\frac{\partial a_x }{\partial y} a_x-\frac{\partial a_y }{\partial y} a_y\\
& =a_x\frac{\partial a_y}{\partial x}-\frac{\partial a_x }{\partial y} a_x \\
\{{\bf k}\cdot (\nabla \times {\bf a})({\bf k} \times {\bf a})\}_y &= \left( \varepsilon_{zxy} \frac{\partial a_y}{\partial x} + \varepsilon_{zyx} \frac{\partial a_x}{\partial y}\right)\varepsilon_{yzx} a_x \\
& = \frac{\partial a_y}{\partial x}a_x - \frac{\partial a_x}{\partial y} a_x
\end{align}
となり同じになる。
上記の関係式を使うと
\begin{align}
u_j\frac{\partial u_i}{\partial x_j} \Bigg|_{\sigma} &= \frac{1}{2}\frac{\partial (u_j u_j)}{\partial x_i}\Bigg|_{\sigma} +\varepsilon_{lmn} k_l \frac{\partial u_n}{\partial x_m}\Bigg|_{\sigma} \varepsilon_{ijp} k_ju_p \\
&=\frac{1}{2}\frac{\partial (u_j u_j)}{\partial x_i}\Bigg|_{\sigma} +\zeta \varepsilon_{ijp} k_ju_p
\end{align}
と書ける。従って風速の方程式は、
\begin{align}
\frac{\partial u_i }{\partial t} &= -\frac{1}{2}\frac{\partial (u_j u_j)}{\partial x_i}\Bigg|_{\sigma} - \dot{\sigma}\frac{\partial u_i}{\partial \sigma}
- (\zeta +f)\varepsilon_{ijl}k_{j}u_l- RT\frac{\partial \tilde{\pi}}{\partial x_i} \Bigg|_{\sigma} - \frac{\partial \phi }{\partial x_i}\Bigg|_{\sigma} \\
\end{align}
と書き直せる。相対渦度$\zeta$とコリオリパラメータ（プラネタリー渦度）$f$の和
\begin{align}
\eta = \zeta +f
\end{align}
を絶対渦度と言う。
以下の演算のスカラー$\psi$に対して、微分の順序は入れ替え、添字を交換すると
\begin{align}
\varepsilon_{lmi} k_l \frac{\partial }{\partial x_m}\frac{\partial }{\partial x_i} \psi\Bigg|_{\sigma} &= \varepsilon_{lmi} k_l \frac{\partial }{\partial x_i}\frac{\partial }{\partial x_m} \psi\Bigg|_{\sigma} = -\varepsilon_{lim} k_l \frac{\partial }{\partial x_i}\frac{\partial }{\partial x_m} \psi\Bigg|_{\sigma}\\
&= -\varepsilon_{lmi} k_l \frac{\partial }{\partial x_m}\frac{\partial }{\partial x_i} \psi\Bigg|_{\sigma} = 0
\end{align}
となる。$k$は$k =(0,0,1)$であり、風速は水平成分しかないので$k_i u_i=0$となることから、以下の計算は
\begin{align}
\varepsilon_{lmi} k_l \frac{\partial }{\partial x_m} \eta\varepsilon_{ijn} k_{j}u_n \Bigg|_{\sigma}&= (\delta_{lj}\delta_{mn} -\delta_{ln}\delta_{mj}) k_l k_{j} \frac{\partial }{\partial x_m} \eta u_n\Bigg|_{\sigma} \\
&= k_j k_{j} \frac{\partial\eta u_m }{\partial x_m}\Bigg|_{\sigma} - k_{m} \frac{\partial }{\partial x_m} \eta k_n u_n\Bigg|_{\sigma}\\
&= \frac{\partial\eta u_m }{\partial x_m}\Bigg|_{\sigma}
\end{align}
となる。
従って、風速の方程式に
\begin{align}
\varepsilon_{lmi} k_l \frac{\partial }{\partial x_m}
\end{align}
を作用させれば、渦度方程式
\begin{align}
\frac{\partial \zeta }{\partial t} &= -\frac{\partial\eta u_m }{\partial x_m}\Bigg|_{\sigma} - \varepsilon_{lmi} k_l \frac{\partial }{\partial x_m}\left\{ \dot{\sigma}\frac{\partial u_i}{\partial \sigma}
+ RT\frac{\partial \tilde{\pi}}{\partial x_i} \Bigg|_{\sigma} \right\} \\
\end{align}
が得られる。
また、添字を交換すると
\begin{align}
\frac{\partial }{\partial x_i} \eta\varepsilon_{ijn} k_{j}u_n \Bigg|_{\sigma}&= k_{j}\varepsilon_{ijn}\frac{\partial }{\partial x_i} \eta u_n \Bigg|_{\sigma} = -k_{i}\varepsilon_{ijn}\frac{\partial }{\partial x_j} \eta u_n \Bigg|_{\sigma}
\end{align}
となるので、風速の方程式に
\begin{align}
\frac{\partial }{\partial x_i}
\end{align}
を作用させれば、発散方程式
\begin{align}
\frac{\partial \delta}{\partial t} &= k_{i}\varepsilon_{ijn}\frac{\partial \eta u_n }{\partial x_j} \Bigg|_{\sigma}
-\frac{\partial^2 }{\partial x_i\partial x_i}\left(\frac{1}{2}u_j u_j +\phi \right)\Bigg|_{\sigma} - \frac{\partial }{\partial x_i}\left\{ \dot{\sigma}\frac{\partial u_i}{\partial \sigma}
+ RT\frac{\partial \tilde{\pi}}{\partial x_i} \Bigg|_{\sigma} \right\}
\end{align}
が得られる。
Reference Profile Temperature
温度$T$を鉛直方向のみ依存する基準温度${\bar T}$ とその摂動部分$T'$に分ける。
\begin{align}
T(x,y,\sigma,t) = \bar{T}(\sigma) +T'(x,y,\sigma,t)
\end{align}
渦度や発散、温度の方程式を書き直すと
\begin{align}
R\bar{T}\varepsilon_{lmi} k_l \frac{\partial }{\partial x_m}\frac{\partial \tilde{\pi}}{\partial x_i} \Bigg|_{\sigma} = 0
\end{align}
であることと、温度移流の項を
\begin{align}
u_j\frac{\partial T'}{\partial x_j}\Bigg|_{\sigma} = \frac{\partial u_jT'}{\partial x_j}\Bigg|_{\sigma} - T'\frac{\partial u_j}{\partial x_j}\Bigg|_{\sigma}
\end{align}
と書けば、最終的に解きたい方程式は、
\begin{align}
\frac{\partial \eta }{\partial t} &= -\frac{\partial\eta u_m }{\partial x_m}\Bigg|_{\sigma} - \varepsilon_{lmi} k_l \frac{\partial }{\partial x_m}\left\{ \dot{\sigma}\frac{\partial u_i}{\partial \sigma}
+ RT'\frac{\partial \tilde{\pi}}{\partial x_i} \Bigg|_{\sigma} \right\} \\
\frac{\partial \delta}{\partial t} &= k_{i}\varepsilon_{ijn}\frac{\partial \eta u_n }{\partial x_j} \Bigg|_{\sigma}
-\frac{\partial^2 }{\partial x_i\partial x_i}\left(\frac{1}{2}u_j u_j+R\bar{T} \tilde{\pi} +\phi \right)\Bigg|_{\sigma} \\
& \qquad \qquad \qquad - \frac{\partial }{\partial x_i}\left\{ \dot{\sigma}\frac{\partial u_i}{\partial \sigma}
+ RT'\frac{\partial \tilde{\pi}}{\partial x_i} \Bigg|_{\sigma} \right\} \\
\frac{\partial T'}{\partial t} &=-\frac{\partial u_jT'}{\partial x_j}\Bigg|_{\sigma} + T'\frac{\partial u_j}{\partial x_j}\Bigg|_{\sigma}- \dot{\sigma}\frac{\partial T}{\partial \sigma}+\kappa T\frac{\omega}{P} + \frac{Q}{C_p}\\
\end{align}
となる。またコリオリパラメータ$f$は定数であることから$\zeta$の時間微分を書き直した。
上記の方程式において水平風速を求める必要があるが、流線関数$\psi$と速度ポテンシャル$\chi$ は、
\begin{align}
\zeta = \frac{\partial^2 \psi}{\partial x_i\partial x_i} \ ,\ \delta = \frac{\partial^2 \chi}{\partial x_i\partial x_i}
\end{align}
を満たすように定義され、回転風と発散風を
\begin{align}
u_{\psi i} = \varepsilon_{ijl} k_l \frac{\partial \psi}{\partial x_j}
\ ,\
u_{\chi i} = \frac{\partial \chi}{\partial x_j}
\end{align}
とずれば、水平風速が
\begin{align}
u_{i} = u_{\psi i}+u_{\chi i}
\end{align}
と求まる。
まとめ
プリミティブ方程式で解くべき方程式は以下である。
予報方程式
\begin{align}
\frac{\partial \eta }{\partial t} &= -\frac{\partial\eta u_m }{\partial x_m}\Bigg|_{\sigma} - \varepsilon_{lmi} k_l \frac{\partial }{\partial x_m}\left\{ \dot{\sigma}\frac{\partial u_i}{\partial \sigma}
+ RT'\frac{\partial \tilde{\pi}}{\partial x_i} \Bigg|_{\sigma} \right\} \\
\frac{\partial \delta}{\partial t} &= k_{i}\varepsilon_{ijn}\frac{\partial \eta u_n }{\partial x_j} \Bigg|_{\sigma}
-\frac{\partial^2 }{\partial x_i\partial x_i}\left(\frac{1}{2}u_j u_j+R\bar{T} \tilde{\pi} +\phi \right)\Bigg|_{\sigma} \\
& \qquad \qquad \qquad - \frac{\partial }{\partial x_i}\left\{ \dot{\sigma}\frac{\partial u_i}{\partial \sigma}
+ RT'\frac{\partial \tilde{\pi}}{\partial x_i} \Bigg|_{\sigma} \right\} \\
\frac{\partial \tilde{\pi}}{\partial t} &= -\int_0^{1} \mathrm{d}\hat{\sigma}A \\
\frac{\partial T'}{\partial t} &=-\frac{\partial u_jT'}{\partial x_j}\Bigg|_{\sigma} + T'\frac{\partial u_j}{\partial x_j}\Bigg|_{\sigma}- \dot{\sigma}\frac{\partial T}{\partial \sigma}+\kappa T\frac{\omega}{P} + \frac{Q}{C_p}\\
\end{align}
診断方程式
\begin{align}
A &= \frac{1}{P^*}\frac{\partial u_jP^*}{\partial x_j}\Bigg|_{\sigma}\\
\phi(\sigma) &= \phi(1)+\int_{\sigma}^1 \mathrm{d}\hat{\sigma} \frac{RT}{\hat{\sigma}} \\
\frac{\omega}{P} & = -\frac{1}{\sigma} \int_0^{\sigma} \mathrm{d}\hat{\sigma}A+ u_j \frac{\partial \tilde{\pi}}{\partial x_j} \\
\dot{\sigma}(\sigma) &= \sigma\int_0^{1} \mathrm{d}\hat{\sigma}A -\int_0^{\sigma} \mathrm{d}\hat{\sigma}A
\end{align}
上記の方程式を解くためには、球座標において微分を計算する必要がある。次回は、球座標における方程式について説明したい。