プリミティブ方程式における水平方向の微分

# 概要
　スペクトル法を用いた数値解析法に興味を持っています。特に、昔から気象予報において使われております。
 
 https://qiita.com/fujitagodai4/items/177a096ee270bff4722f
 
 今回の記事では、球座標における微分について紹介したいと思います。


 
# 直交曲線座標
　デカルト座標を$x^i$と表現すると変換された新しい座標系$q^j$は
```math
\begin{align}
q^j = q^j(x^i)
\end{align}
```
と書ける。逆の関係も成立する。
```math
\begin{align}
x^i = x^i(q^j)
\end{align}
```

## 共変・反共変基底ベクトル

 一般的なベクトルは、共変・反共変基底ベクトル
```math
\begin{align}
{\bf g}_j = \frac{\partial {\bf x} }{\partial q^i} \ \ , \ \ {\bf g}^i =  {\bf \nabla} q^i
\end{align}
```
を使い
```math
\begin{align}
{\bf A} = A^j {\bf g}_j =  A_m {\bf g}^m 
\end{align}
```
と書ける。
　共変・反共変基底ベクトルの積は
```math
\begin{align}
 g_{ij} \equiv {\bf g}_i \cdot {\bf g}_j  \ \ , \ \  g^{ij} \equiv {\bf g}^i \cdot {\bf g}^j
\end{align}
```
と定義され、
```math
\begin{align}
{\bf g}^i \cdot {\bf g}_j = {\bf \nabla} q^i \frac{\partial {\bf x} }{\partial q^j} = \delta^i_j {\bf \nabla} {\bf x} = \delta^i_j
\end{align}
```
ヤコビアンは
```math
\begin{align}
J^2 = \Pi_{i=1}^{3} g_{ii} 
\end{align}
```
で与えられる。ベクトルに共変・反共変基底ベクトルを作用させると、
```math
\begin{align}
{\bf A}\cdot{\bf g}^k  &= A^j {\bf g}_j \cdot{\bf g}^k  =A^k \\
& = A_j {\bf g}^j \cdot{\bf g}^k =  A_j g^{jk}  \\
{\bf A}\cdot{\bf g}_k  &= A_j {\bf g}^j \cdot{\bf g}_k  =A_k \\
& = A^j {\bf g}_j \cdot{\bf g}_k =  A^j g_{jk}  \\
\end{align}
```
となる。共変・反共変基底ベクトルを使うと、恒等テンソルは
```math
\begin{align}
{\bf I} = {\bf g}_j {\bf g}^j = {\bf g}^j {\bf g}_j
\end{align}
```
とあらわせ、外積の計算は
```math
\begin{align}
{\bf A}\times {\bf B} = J\varepsilon_{ijk}A^i B^j {\bf g}^k 
\end{align}
```
となる。
　
## Christoffel 記号
　
　Christoffel 記号は、
```math
\begin{align}
 \Gamma^{p}_{ki} \equiv {\bf g}^p \cdot \frac{\partial {\bf g}_k }{\partial q^i}
\end{align}
```
と定義する。両辺から反共変ベクトルを作用させると
```math
\begin{align}
 \Gamma^{p}_{ki} {\bf g}_p =  {\bf g}_p {\bf g}^p  \frac{\partial {\bf g}_k }{\partial q^i} = {\bf I} \frac{\partial {\bf g}_k }{\partial q^i} = \frac{\partial {\bf g}_k }{\partial q^i}
\end{align}
```
の関係式が得られる。また、$k\leftrightarrow j$ に関して
```math
\begin{align}
 \Gamma^{p}_{ki} &= {\bf g}^p \cdot \frac{\partial {\bf g}_k }{\partial q^i} = {\bf g}^p \cdot \frac{\partial {\bf x} }{\partial q^i \partial q^k} \\
 &= {\bf g}^p \cdot \frac{\partial {\bf g}_i }{\partial q^k} =   \Gamma^{p}_{ik} 
\end{align}
```
となるので対称であり、さらに反共変基底ベクトルの定義から
```math
\begin{align}
 \Gamma^{p}_{ki} &= {\bf g}^p \cdot \frac{\partial {\bf g}_k }{\partial q^i} ={\bf \nabla} q^i\cdot \frac{\partial {\bf x} }{\partial q^i \partial q^k}  = \frac{\partial q^i}{\partial x_p}\frac{\partial x_p}{\partial q^i \partial q^k} 
\end{align}
```
と書ける。


## 微分演算子
　
 　微分は、
```math
\begin{align}
{\bf \nabla } \equiv  {\bf g}^k \frac{\partial }{\partial q^k}
\end{align}
```
と書ける。従って、発散・回転・勾配は、
```math
\begin{align}
{\bf \nabla }\cdot {\bf A} &= {\bf g}^k \cdot \frac{\partial {\bf A} }{\partial q^k} \\
{\bf \nabla }\times {\bf A} &= {\bf g}^k \times  \frac{\partial {\bf A} }{\partial q^k} \\
{\bf \nabla } {\bf A} &= {\bf g}^k \frac{\partial {\bf A} }{\partial q^k} \\
\end{align}
```
と書ける。ラプラシアンは、
```math
\begin{align}
{\bf \nabla }^2 a &=  {\bf \nabla }\cdot{\bf \nabla } a = 
{\bf g}^i \cdot \frac{\partial }{\partial q^i} \left( {\bf g}^k \cdot \frac{\partial a}{\partial q^k}\right) \\
&= g^{ik} \frac{\partial^2 a}{\partial q^i \partial q^k}
 + {\bf g}^i \cdot \frac{\partial {\bf g}^k}{\partial q^i}\frac{\partial a}{\partial q^k}\\
\end{align}
```
となるが、
```math
\begin{align}
{\bf g}^k = g^{kj} {\bf g}_j 
\end{align}
```
使うと
```math
\begin{align}
{\bf g}^i \cdot \frac{\partial {\bf g}^k}{\partial q^i} &= 
{\bf g}^i \cdot \frac{\partial (g^{kj} {\bf g}_j )}{\partial q^i} = {\bf g}^i \cdot {\bf g}_j \frac{\partial g^{kj}}{\partial q^i} + g^{kj} {\bf g}^i \cdot \frac{\partial  {\bf g}_j }{\partial q^i} \\
&= \delta^i_j \frac{\partial g^{kj}}{\partial q^i} + g^{kj} {\bf g}^i \cdot \frac{\partial  {\bf g}_j }{\partial q^i} \\
&= \frac{\partial g^{ki}}{\partial q^i} + g^{kj} \Gamma^i_{ji} \\
\end{align}
```
なるので、
```math
\begin{align}
{\bf \nabla }^2 a &= g^{ik} \frac{\partial^2 a}{\partial q^i \partial q^k}
 + \frac{\partial g^{ki}}{\partial q^i}\frac{\partial a}{\partial q^k} + g^{kj} \Gamma^i_{ji} \frac{\partial a}{\partial q^k}\\
\end{align}
```
と書ける。


# 球極座標

 　通常の球極座標は、北極と南極で特異点になることが知られている。従って、修正された球極座標$(\lambda,\mu,r)$は、デカルト座標を用いた
```math
\begin{align}
x &= r\sqrt{1-\mu^2 } \cos\lambda \\
y &= r\sqrt{1-\mu^2 }  \sin\lambda \\
z &= r \mu
\end{align}
``` 
を用いる。通常の極座標は、
```math
\begin{align}
\mu = \sin \varphi
\end{align}
``` 
である。$\lambda$は経度で$\varphi$は緯度に相当する。$r$は地球の半径または高度で$r=a$と固定する。

　この場合、共変基底ベクトルは、
```math
\begin{align}
{\bf g}_{\lambda} &= \frac{\partial {\bf x}}{\partial \lambda} = \left(-r\sqrt{1-\mu^2 } \sin\lambda \ , \ r\sqrt{1-\mu^2 } \cos\lambda  \ , \ 0 \right)  \\ 
{\bf g}_{\mu} &= \frac{\partial {\bf x}}{\partial \mu} = \left(\frac{-r \mu}{\sqrt{1-\mu^2} } \cos\lambda \ , \ \frac{-r \mu}{\sqrt{1-\mu^2} }\sin\lambda  \ , \ r \right) \\ 
{\bf g}_{r} &= \frac{\partial {\bf x}}{\partial r} = \left(\sqrt{1-\mu^2 } \cos\lambda \ , \ \sqrt{1-\mu^2 }  \sin\lambda  \ , \ \mu \right) \\ 
\end{align}
```
となるので
```math
\begin{align}
g_{\lambda\lambda} &= {\bf g}_{\lambda} \cdot {\bf g}_{\lambda} = r^2(1-\mu^2) \\
g_{\mu\mu} &= {\bf g}_{\mu} \cdot {\bf g}_{\mu} = \frac{r^2\mu^2}{(1-\mu^2) } + r^2 = \frac{r^2}{(1-\mu^2) } \\
g_{rr} &= {\bf g}_{r} \cdot {\bf g}_{r} =1
\end{align}
```
となる。非対角成分はゼロであるので、逆行列は
```math
\begin{align}
g^{\lambda\lambda} =  \frac{1}{r^2(1-\mu^2)} \ \ , \ \
g^{\mu\mu}  = \frac{(1-\mu^2) }{r^2} \ \ , \ \
g^{rr} =1
\end{align}
```
となる。また、共変基底ベクトルの外積は、
```math
\begin{align}
{\bf g}_{\lambda} \times {\bf g}_{\mu} &=  r^2\sqrt{1-\mu^2}\left( 
\cos\lambda  \ , \  \sin\lambda \ , \ \frac{\mu}{\sqrt{1-\mu^2}}  \right) = r^2 {\bf g}_{r} \\
{\bf g}_{\mu} \times {\bf g}_{r} &=  \frac{r}{\sqrt{1-\mu^2}}\left( 
-\sin\lambda  \ , \  \cos\lambda \ , \ 0 \right) = \frac{1}{1-\mu^2} {\bf g}_{\lambda} \\
{\bf g}_{r} \times {\bf g}_{\lambda} &=  r\sqrt{1-\mu^2}\left( 
-\mu\cos\lambda  \ , \  -\mu\sin\lambda \ , \ \sqrt{1-\mu^2}  \right) = (1-\mu^2) {\bf g}_{\mu} \\
\end{align}
```
と書ける。Christoffel 記号の非ゼロ成分は、
```math
\begin{align}
\Gamma^{r}_{\mu\mu} &=  {\bf g}^{r} \cdot \frac{\partial {\bf g}_{\mu}}{\partial \mu}
=  g^{rr}{\bf g}_{r} \cdot \frac{\partial {\bf g}_{\mu}}{\partial \mu}\\
&= -r + \frac{r\mu^2}{1-\mu^2} = -\frac{r}{1-\mu^2} \\
\Gamma^{r}_{\lambda\lambda} &= -r(1-\mu^2) \\
\Gamma^{\mu}_{r\mu} &=\Gamma^{\mu}_{\mu r}= \frac{1}{r} \\
\Gamma^{\mu}_{\mu\mu} &= \frac{\mu}{1-\mu^2} \\
\Gamma^{\mu}_{\lambda\lambda} &= \mu(1-\mu^2) \\
\Gamma^{\lambda}_{r\lambda} &=\Gamma^{\lambda}_{\lambda r}= \frac{1}{r} \\
\Gamma^{\lambda}_{\mu\lambda} &= \Gamma^{\lambda}_{\lambda\mu}  =-\frac{\mu}{1-\mu^2} \\
\end{align}
```
となる。

 ## 球極座標における水平勾配
 
　修正された球極座標において水平勾配は、$\phi$をスカラーとすると
```math
\begin{align}
\nabla \phi &= {\bf g}^{\lambda} \frac{\partial \phi}{ \partial \lambda}+ {\bf g}^{\mu} \frac{\partial \phi}{ \partial \mu} \\
& =  g^{\lambda\lambda} {\bf g}_{\lambda} \frac{\partial \phi}{ \partial \lambda}+ g^{\mu\mu}{\bf g}_{\mu} \frac{\partial \phi}{ \partial \mu} \\
& =\frac{1}{a^2(1-\mu^2)} {\bf g}_{\lambda} \frac{\partial \phi}{ \partial \lambda}+ \frac{(1-\mu^2) }{a^2} {\bf g}_{\mu} \frac{\partial \phi}{ \partial \mu} 
\end{align}
```
と書ける。また$a$は地球の半径で$r=a$と置いた。

 ## 球極座標における回転
 
　回転は、
```math
\begin{align}
{\bf k}\cdot {\bf \nabla }\times {\bf A} &= {\bf g}^r\cdot \left\{ \left({\bf g}^{\lambda} \frac{\partial }{ \partial \lambda}+ {\bf g}^{\mu} \frac{\partial }{ \partial \mu}\right)\times\left( {\bf g}_{\lambda} A^{\lambda}+ {\bf g}_{\mu} A^{\mu}\right)\right\}
\end{align}
```
と書けるが、
```math
\begin{align}
{\bf g}^r\cdot\left\{ {\bf g}^{\lambda} \frac{\partial }{ \partial \lambda}\times ({\bf g}_{\lambda} A^{\lambda})\right\} &= {\bf g}^r\cdot\left\{ {\bf g}^{\lambda} \times \frac{\partial {\bf g}_{\lambda} }{ \partial \lambda} A^{\lambda} + {\bf g}^{\lambda} \times {\bf g}_{\lambda}  \frac{\partial A^{\lambda}}{ \partial \lambda}\right\} \\
&= g^{\lambda\lambda}{\bf g}^r \cdot\left\{ {\bf g}_{\lambda} \times  \frac{\partial {\bf g}^{\lambda} }{ \partial \lambda} A^{\lambda} + {\bf g}_{\lambda} \times {\bf g}_{\lambda}  \frac{\partial A^{\lambda}}{ \partial \lambda}\right\} \\
&= g^{\lambda\lambda}{\bf g}^r \cdot \left\{ \Gamma^{p}_{\lambda\lambda} {\bf g}_{\lambda} \times {\bf g}_p  \right\} A^{\lambda}\\
&= g^{\lambda\lambda}{\bf g}^r \cdot \left\{ \Gamma^{\mu}_{\lambda\lambda} {\bf g}_{\lambda} \times {\bf g}_{\mu} 
+\Gamma^{r}_{\lambda\lambda} {\bf g}_{\lambda} \times {\bf g}_r 
\right\} A^{\lambda}\\
& = r^2 g^{\lambda\lambda}\Gamma^{\mu}_{\lambda\lambda}A^{\lambda}
\end{align}
```
```math
\begin{align}
{\bf g}^r\cdot\left\{ {\bf g}^{\lambda} \frac{\partial }{ \partial \lambda}\times ({\bf g}_{\mu} A^{\mu})\right\} &= {\bf g}^r\cdot\left\{ {\bf g}^{\lambda} \times \frac{\partial {\bf g}_{\mu} }{\partial \lambda} A^{\mu} + {\bf g}^{\lambda} \times {\bf g}_{\mu}  \frac{\partial A^{\mu}}{ \partial \lambda}\right\} \\
&= {\bf g}^r \cdot\left\{ {\bf g}^{\lambda} \times \frac{\partial {\bf g}_{\mu} }{\partial \lambda} A^{\mu} + {\bf g}^{\lambda} \times {\bf g}_{\mu}  \frac{\partial A^{\mu}}{ \partial \lambda}\right\} \\
&= {\bf g}^r \cdot \left\{ \Gamma^{p}_{\mu\lambda} {\bf g}^{\lambda} \times {\bf g}_p A^{\mu} + {\bf g}^{\lambda} \times {\bf g}_{\mu}  \frac{\partial A^{\mu}}{ \partial \lambda} \right\} \\
&= g^{\lambda\lambda}{\bf g}^r \cdot \left\{ \Gamma^{p}_{\mu\lambda} {\bf g}_{\lambda} \times {\bf g}_p A^{\mu}+ {\bf g}_{\lambda} \times {\bf g}_{\mu}  \frac{\partial A^{\mu}}{ \partial \lambda} \right\} \\
& =r^2 g^{\lambda\lambda} \left\{ \Gamma^{\mu}_{\mu\lambda}A^{\mu} +  \frac{\partial A^{\mu}}{ \partial \lambda} \right\} 
\end{align}
```
同様に、
```math
\begin{align}
{\bf g}^r\cdot\left\{ {\bf g}^{\mu} \frac{\partial }{ \partial \mu}\times ({\bf g}_{\lambda} A^{\lambda})\right\} &= {\bf g}^r\cdot\left\{ {\bf g}^{\mu} \times \frac{\partial {\bf g}_{\lambda} }{ \partial \mu} A^{\lambda} + {\bf g}^{\mu} \times {\bf g}_{\lambda}  \frac{\partial A^{\lambda}}{ \partial \mu}\right\} \\
 &= g^{\mu\mu} {\bf g}^r\cdot\left\{ \Gamma^{p}_{\mu\lambda} {\bf g}_{\mu} \times {\bf g}_{p} A^{\lambda} + {\bf g}_{\mu} \times {\bf g}_{\lambda}  \frac{\partial A^{\lambda}}{ \partial \mu}\right\} \\
 &= g^{\mu\mu} {\bf g}^r\cdot\left\{ \Gamma^{p}_{\mu\lambda} {\bf g}_{\mu} \times {\bf g}_{p} A^{\lambda} + {\bf g}_{\mu} \times {\bf g}_{\lambda} \frac{\partial A^{\lambda}}{ \partial \mu}\right\} \\
 &= -r^2 g^{\mu\mu} \left\{ \Gamma^{\lambda}_{\mu\lambda} A^{\lambda} + \frac{\partial A^{\lambda}}{ \partial \mu}  \right\}
\end{align}
```
```math
\begin{align}
{\bf g}^r\cdot\left\{ {\bf g}^{\mu} \frac{\partial }{ \partial \mu}\times ({\bf g}_{\mu} A^{\mu})\right\} &= {\bf g}^r\cdot\left\{ {\bf g}^{\mu} \times \frac{\partial {\bf g}_{\mu} }{ \partial \mu} A^{\mu} + {\bf g}^{\mu} \times {\bf g}_{\mu}  \frac{\partial A^{\mu}}{ \partial \mu}\right\} \\
 &= g^{\mu\mu} {\bf g}^r\cdot\left\{ \Gamma^{p}_{\mu\mu} {\bf g}_{\mu} \times {\bf g}_{p} A^{\mu} \right\} \\
 &= -r^2 g^{\mu\mu}\Gamma^{\lambda}_{\mu\mu} A^{\mu} 
\end{align}
```
となるので、まとめると
```math
\begin{align}
{\bf k}\cdot {\bf \nabla }\times {\bf A} &= r^2 g^{\lambda\lambda}\Gamma^{\mu}_{\lambda\lambda}A^{\lambda} +r^2 g^{\lambda\lambda} \left\{ \Gamma^{\mu}_{\mu\lambda}A^{\mu} +  \frac{\partial A^{\mu}}{ \partial \lambda} \right\} \\
&-r^2 g^{\mu\mu} \left\{ \Gamma^{\lambda}_{\mu\lambda} A^{\lambda} + \frac{\partial A^{\lambda}}{ \partial \mu}  \right\} -r^2 g^{\mu\mu}\Gamma^{\lambda}_{\mu\mu} A^{\mu}\\
&= \mu A^{\lambda} +  \frac{1}{(1-\mu^2)}  \frac{\partial A^{\mu}}{ \partial \lambda} +\mu A^{\lambda} -(1-\mu^2)\frac{\partial A^{\lambda}}{ \partial \mu}  
\\
&=\frac{1}{1-\mu^2}\frac{\partial A^{\mu}}{\partial \lambda} + \frac{\partial }{\partial \mu} \left\{ (\mu^2-1)A^{\lambda}\right\} 
\end{align}
```
と書ける。

 ## 球極座標における発散
 
　発散は、
```math
\begin{align}
{\bf \nabla }\cdot {\bf A} &= \left({\bf g}^{\lambda} \frac{\partial }{ \partial \lambda}+ {\bf g}^{\mu} \frac{\partial }{ \partial \mu}\right)\cdot \left( {\bf g}_{\lambda} A^{\lambda}+ {\bf g}_{\mu} A^{\mu}\right) \\
&= \frac{\partial A^{\lambda}}{ \partial \lambda} + A^{\lambda}{\bf g}^{\lambda} \frac{\partial {\bf g}_{\lambda} }{ \partial \lambda} +A^{\mu}{\bf g}^{\lambda} \frac{\partial {\bf g}_{\mu} }{ \partial \lambda} \\
&\qquad	\qquad	+ A^{\lambda}{\bf g}^{\mu} \frac{\partial {\bf g}_{\lambda} }{ \partial \mu} + A^{\mu}{\bf g}^{\mu} \frac{\partial {\bf g}_{\mu} }{ \partial \mu} + \frac{\partial A^{\mu} }{ \partial \mu} \\
&= \frac{\partial A^{\lambda}}{ \partial \lambda} + A^{\lambda}\Gamma^{\lambda}_{\lambda\lambda} +A^{\mu}\Gamma^{\lambda}_{\mu\lambda} \\
&\qquad	 \qquad	 + A^{\lambda}\Gamma^{\mu}_{\lambda\mu}  + A^{\mu}\Gamma^{\mu}_{\mu\mu} + \frac{\partial A^{\mu} }{ \partial \mu} \\
& = \frac{\partial A^{\lambda}}{ \partial \lambda}  + \frac{\partial A^{\mu} }{ \partial \mu} 
\end{align}
```
と書ける。

 ## 球極座標におけるラプラシアン

　ラプラシアンは、
```math
\begin{align}
{\bf \nabla }^2 \phi &= g^{ik} \frac{\partial^2 \phi}{\partial q^i \partial q^k}
 + \frac{\partial g^{ki}}{\partial q^i}\frac{\partial \phi}{\partial q^k} + g^{kj} \Gamma^i_{ji} \frac{\partial \phi}{\partial q^k}\\
\end{align}
```
と書けるが、1項目は

 ```math
\begin{align}
g^{ik} \frac{\partial^2 \phi}{\partial q^i \partial q^k}
&=g^{rr} \frac{\partial^2 \phi}{\partial r^2} +g^{\lambda\lambda} \frac{\partial^2 \phi}{\partial \lambda^2} + g^{\mu\mu} \frac{\partial^2 \phi}{\partial \mu^2} \\
&=\frac{\partial^2 \phi}{\partial r^2} +\frac{1}{r^2(1-\mu^2)}\frac{\partial^2 \phi}{\partial \lambda^2} + \frac{(1-\mu^2) }{r^2}  \frac{\partial^2 \phi}{\partial \mu^2}
\end{align}
```
2項目は
```math
\begin{align}
\frac{\partial g^{ki}}{\partial q^i}\frac{\partial \phi}{\partial q^k} 
= \frac{\partial g^{rr}}{\partial r}\frac{\partial \phi}{\partial r} 
+ \frac{\partial g^{\lambda\lambda}}{\partial \lambda}\frac{\partial \phi}{\partial \lambda} + \frac{\partial g^{\mu\mu}}{\partial \mu}\frac{\partial \phi}{\partial \mu} = -\frac{2\mu}{r^2} \frac{\partial \phi}{\partial \mu} 
\end{align}
```

3項目は
```math
\begin{align}
g^{kj} \Gamma^i_{ji} \frac{\partial \phi}{\partial q^k} &= 
g^{rr} \Gamma^r_{rr} \frac{\partial \phi}{\partial r} +
g^{rr} \Gamma^{\lambda}_{\lambda r} \frac{\partial \phi}{\partial r} +
g^{rr} \Gamma^{\mu}_{\mu r} \frac{\partial \phi}{\partial r} \\
&+ g^{\lambda\lambda} \Gamma^r_{\lambda r} \frac{\partial \phi}{\partial \lambda} +
g^{\lambda\lambda} \Gamma^{\lambda}_{\lambda \lambda} \frac{\partial \phi}{\partial \lambda} +
g^{\lambda\lambda} \Gamma^{\mu}_{\mu \lambda} \frac{\partial \phi}{\partial \lambda} \\
&+ g^{\mu\mu} \Gamma^r_{\mu r} \frac{\partial \phi}{\partial \mu} +
g^{\mu\mu} \Gamma^{\lambda}_{\mu \lambda} \frac{\partial \phi}{\partial \mu} +
g^{\mu\mu} \Gamma^{\mu}_{\mu \mu} \frac{\partial \phi}{\partial \mu}\\
&=\frac{2}{r}\frac{\partial \phi}{\partial r}
\end{align}
```
となるので、最終的に
```math
\begin{align}
{\bf \nabla }^2 \phi &= \frac{1}{a^2} \left\{\frac{1}{1-\mu^2} \frac{\partial^2 }{\partial \lambda^2} + \frac{\partial}{\partial \mu} 
\left\{\left(1-\mu^2 \right)\frac{\partial}{\partial \mu}\right\}\right\}\phi
\end{align}
```
と書ける。

# Helmholtz 分解

　Helmholtz 分解を行うと水平風速は、流線関数$\psi$と速度ポテンシャル$\chi$を用いると
 ```math
\begin{align}
{\bf v} &\equiv  {\dot \lambda} {\bf g}_{\lambda} + {\dot \mu} {\bf g}_{\mu}\\
&={\bf k} \times \nabla \psi + \nabla \chi
\end{align}
```
と書ける。球極座標においては
 ```math
\begin{align}
{\bf v} 
&={\bf k} \times \nabla \psi + \nabla \chi\\
&={\bf g}_{r} \times \left({\bf g}^{\lambda} \frac{\partial \psi}{ \partial \lambda}+ {\bf g}^{\mu} \frac{\partial \psi}{ \partial \mu}\right) + {\bf g}^{\lambda} \frac{\partial \chi}{ \partial \lambda}+ {\bf g}^{\mu} \frac{\partial \chi}{ \partial \mu}\\
&=  g^{\lambda\lambda} {\bf g}_{r} \times  {\bf g}_{\lambda} \frac{\partial \psi}{ \partial \lambda} + g^{\mu\mu}{\bf g}_{r} \times {\bf g}_{\mu} \frac{\partial \psi}{ \partial \mu}
+  g^{\lambda\lambda}{\bf g}_{\lambda} \frac{\partial \chi}{ \partial \lambda}+ g^{\mu\mu}{\bf g}_{\mu} \frac{\partial \chi}{ \partial \mu}\\
&=  \frac{1}{a^2} \frac{\partial \psi}{ \partial \lambda}{\bf g}_{\mu}  - \frac{1}{a^2}\frac{\partial \psi}{ \partial \mu}{\bf g}_{\lambda}
+  \frac{1}{a^2(1-\mu^2)}\frac{\partial \chi}{ \partial \lambda}{\bf g}_{\lambda} + \frac{(1-\mu^2)}{a^2} \frac{\partial \chi}{ \partial \mu}{\bf g}_{\mu} \\
&= \frac{1}{a^2}\left\{-\frac{\partial \psi}{ \partial \mu}+\frac{1}{1-\mu^2}\frac{\partial \chi}{ \partial \lambda} \right\}{\bf g}_{\lambda}
+ \frac{1}{a^2}\left\{\frac{\partial \psi}{ \partial \lambda}+ (1-\mu^2) \frac{\partial \chi}{ \partial \mu}\right\}{\bf g}_{\mu}
\end{align}
```
なので
 ```math
\begin{align}
{\dot \lambda}  &=  \frac{1}{a^2}\left\{-\frac{\partial \psi}{ \partial \mu}+\frac{1}{1-\mu^2}\frac{\partial \chi}{ \partial \lambda} \right\} \\
{\dot \mu} &=  \frac{1}{a^2}\left\{\frac{\partial \psi}{ \partial \lambda}+ (1-\mu^2)\frac{\partial \chi}{ \partial \mu}\right\}
\end{align}
```
の関係式が得られる。


 ## 移流項

　物理的な速度を
 ```math
\begin{align}
U &\equiv  {\dot \lambda}a (1-\mu^2) = a^{-1}\left\{(\mu^2-1) \frac{\partial \psi}{ \partial \mu}+\frac{\partial \chi}{ \partial \lambda} \right\} \\
V &\equiv {\dot \mu}a = a^{-1} \left\{\frac{\partial \psi}{ \partial \lambda}+ (1-\mu^2)\frac{\partial \chi}{ \partial \mu}\right\}
\end{align}
```
と定義する。すると移流項は、任意の物理量$\phi$に対して
 ```math
\begin{align}
{\bf v}\cdot\nabla \phi &= \left( {\dot \lambda}{\bf g}_{\lambda} + {\dot \mu}{\bf g}_{\mu} \right) \cdot  \left({\bf g}^{\lambda} \frac{\partial \phi}{ \partial \lambda}+ {\bf g}^{\mu} \frac{\partial \phi}{ \partial \mu}\right) \\
&= {\dot \lambda} \frac{\partial \phi}{ \partial \lambda}+ {\dot \mu}\frac{\partial \phi}{ \partial \mu} \\
&=\frac{U}{a (1-\mu^2)} \frac{\partial \phi}{ \partial \lambda}+ \frac{V}{a}\frac{\partial \phi}{ \partial \mu}
\end{align}
```
となる。

## 渦度と発散

　渦度は、極座標で表すと
 ```math
\begin{align}
\zeta = {\bf k}\cdot {\bf \nabla }\times {\bf v} &=\frac{1}{1-\mu^2}\frac{\partial {\dot \mu}}{\partial \lambda} + \frac{\partial }{\partial \mu} \left\{ (\mu^2-1){\dot \lambda}\right\} \\
&=\frac{1}{a^2(1-\mu^2)}\frac{\partial }{\partial \lambda} \left\{\frac{\partial \psi}{ \partial \lambda}+ (1-\mu^2)\frac{\partial \chi}{ \partial \mu}\right\} \\
&\qquad	\qquad	+ \frac{1}{a^2}\frac{\partial }{\partial \mu} \left\{ (1-\mu^2)\frac{\partial \psi}{ \partial \mu} - \frac{\partial \chi}{ \partial \lambda}\right\} \\
 &= \frac{1}{a^2} \left\{\frac{1}{1-\mu^2} \frac{\partial^2 }{\partial \lambda^2} + \frac{\partial}{\partial \mu} 
\left\{\left(1-\mu^2 \right)\frac{\partial}{\partial \mu}\right\}\right\}\psi \\
&={\bf \nabla }^2 \psi
\end{align}
```
となりラプラシアンで表現できる。発散も同様に
```math
\begin{align}
\delta ={\bf \nabla }\cdot {\bf v} 
 &= \frac{\partial {\dot \lambda}}{ \partial \lambda}  + \frac{\partial {\dot \mu} }{ \partial \mu} \\
&=\frac{1}{a^2(1-\mu^2)} \frac{\partial }{ \partial \lambda} \left\{ (\mu^2+1)\frac{\partial \psi}{ \partial \mu} + \frac{\partial \chi}{ \partial \lambda}\right\} \\
&\qquad	\qquad	+\frac{1}{a^2}\frac{\partial }{ \partial \mu}\left\{\frac{\partial \psi}{ \partial \lambda}+ (1-\mu^2)\frac{\partial \chi}{ \partial \mu}\right\}\\
 &= \frac{1}{a^2} \left\{\frac{1}{1-\mu^2} \frac{\partial^2 }{\partial \lambda^2} + \frac{\partial}{\partial \mu} 
\left\{\left(1-\mu^2 \right)\frac{\partial}{\partial \mu}\right\}\right\}\chi \\
&={\bf \nabla }^2 \chi
\end{align}
```
ラプラシアンで表現できる。

## α-演算子

　$\alpha$-演算子を
```math
\begin{align}
 \alpha(A,B) &\equiv\frac{1}{1-\mu^2}\frac{\partial A}{\partial \lambda} + \frac{\partial B}{\partial \mu}\\
\end{align}
```
と定義する。例えば、渦度と発散は$U,V$を用いて
 ```math
\begin{align}
\zeta = {\bf k}\cdot {\bf \nabla }\times {\bf v} &=\frac{1}{1-\mu^2}\frac{\partial {\dot \mu}}{\partial \lambda} + \frac{\partial }{\partial \mu} \left\{ (\mu^2-1){\dot \lambda}\right\} \\
&=\frac{1}{a(1-\mu^2)}\frac{\partial V}{\partial \lambda} + \frac{1}{a} \frac{\partial U}{\partial \mu} \\
&= \frac{1}{a}\alpha(V,-U) \\
\delta ={\bf \nabla }\cdot {\bf v} 
 &= \frac{\partial {\dot \lambda}}{ \partial \lambda}  + \frac{\partial {\dot \mu} }{ \partial \mu} \\
 &= \frac{1}{a(1-\mu^2)}\frac{\partial U}{\partial \lambda} + \frac{1}{a} \frac{\partial V}{\partial \mu} \\
 &= \frac{1}{a}\alpha(U,V) \\
\end{align}
```
と書ける。
　$\alpha$-演算子は、プリミティブ方程式を書き換えるときに用いる。

# プリミティブ方程式

　プリミティブ方程式を上記の極座標系の微分に置き換える。プリミティブ方程式は以下であった。
 
https://qiita.com/fujitagodai4/items/177a096ee270bff4722f#%E3%81%BE%E3%81%A8%E3%82%81
 
 前回の記事ではテンソル表記であったが、便利上ベクトル表記をしている。
### 予報方程式
 ```math
\begin{align}
\frac{\partial \eta }{\partial t} &= - \nabla\cdot(\eta {\bf u} )  - {\bf k}\cdot\nabla \times\left\{ \dot{\sigma}\frac{\partial {\bf u} }{\partial \sigma}
+ RT' \nabla\tilde{\pi}  \right\} \\

\frac{\partial \delta}{\partial t} &= {\bf k}\cdot\nabla \times(\eta  {\bf u})
-\nabla^2 \left(\frac{1}{2}{\bf u}^2+R\bar{T} \tilde{\pi} +\phi \right) \\
& \qquad　\qquad \qquad - \nabla \cdot\left\{ \dot{\sigma}\frac{\partial {\bf u}}{\partial \sigma}
+ RT' \nabla \tilde{\pi} \right\} \\

\frac{\partial \tilde{\pi}}{\partial t} &= -\int_0^{1} \mathrm{d}\hat{\sigma}A \\
\frac{\partial T'}{\partial t} &=- \nabla \cdot({\bf u}T') + T'\delta - \dot{\sigma}\frac{\partial T}{\partial \sigma}+\kappa T\frac{\omega}{P} + \frac{Q}{C_p}\\
\end{align}
```

### 診断方程式

```math
\begin{align}
A &= \frac{1}{P^*} \nabla \cdot ({\bf u} P^*)= \delta + {\bf u}\cdot\nabla \tilde{\pi} \\
 \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+{\bf u}\cdot\nabla \tilde{\pi} \\
\dot{\sigma}(\sigma) &= \sigma\int_0^{1} \mathrm{d}\hat{\sigma}A -\int_0^{\sigma} \mathrm{d}\hat{\sigma}A
\end{align}
```


　渦度方程式は、上記の発散・回転・勾配の式を用いて
 ```math
\begin{align}
\frac{\partial \eta }{\partial t} &= - \nabla\cdot(\eta {\bf u} )  - {\bf k}\cdot\nabla \times\left\{ \dot{\sigma}\frac{\partial {\bf u} }{\partial \sigma}
+ RT' \nabla\tilde{\pi}  \right\} \\
&= -\frac{\partial }{ \partial \lambda}(\eta\dot{\lambda})  - \frac{\partial }{ \partial \mu}(\eta\dot{\mu}) \\
&\qquad- \frac{1}{1-\mu^2}\frac{\partial }{\partial \lambda}\left\{ \dot{\sigma}\frac{\partial \dot{\mu} }{\partial \sigma}
+ RT' \frac{(1-\mu^2) }{a^2}\frac{\partial \tilde{\pi}}{\partial \mu}  \right\} \\
&\qquad- \frac{\partial }{\partial \mu} \left[ (\mu^2-1)\left\{ \dot{\sigma}\frac{\partial \dot{\lambda} }{\partial \sigma}
+ \frac{RT' }{a^2(1-\mu^2)}\frac{\partial \tilde{\pi}}{\partial \lambda}  \right\} \right] \\
 
&= -\frac{1}{a(1-\mu^2)}\frac{\partial }{ \partial \lambda}\left(\eta U \right)  - \frac{1}{a}\frac{\partial }{ \partial \mu}\left(\eta V \right) \\
&\qquad- \frac{1}{a(1-\mu^2)}\frac{\partial }{\partial \lambda}\left\{ \dot{\sigma}\frac{\partial V}{\partial \sigma}
+ RT' \frac{(1-\mu^2) }{a}\frac{\partial \tilde{\pi}}{\partial \mu}  \right\} \\
&\qquad + \frac{1}{a}\frac{\partial }{\partial \mu} \left\{\dot{\sigma}\frac{\partial U }{\partial \sigma}
+ \frac{RT'}{a}\frac{\partial \tilde{\pi}}{\partial \lambda}  \right\}  \\
&=  -\frac{1}{a(1-\mu^2)}\frac{\partial }{\partial \lambda}\left\{ \eta U + \dot{\sigma}\frac{\partial V}{\partial \sigma}
+ RT' \frac{(1-\mu^2) }{a}\frac{\partial \tilde{\pi}}{\partial \mu} \right\} \\
&\qquad - \frac{1}{a}\frac{\partial }{\partial \mu} \left\{ \eta V- \dot{\sigma}\frac{\partial U }{\partial \sigma}
- \frac{RT'}{a}\frac{\partial \tilde{\pi}}{\partial \lambda} \right\}  \\
\end{align}
```

発散方程式も同様に
 ```math
\begin{align}
\frac{\partial \delta}{\partial t} &= {\bf k}\cdot\nabla \times(\eta  {\bf u})
-\nabla^2 \left(\frac{1}{2}{\bf u}^2+R\bar{T} \tilde{\pi} +\phi \right) - \nabla \cdot\left\{ \dot{\sigma}\frac{\partial {\bf u}}{\partial \sigma}
+ RT' \nabla \tilde{\pi} \right\} \\
&=\frac{1}{1-\mu^2}\frac{\partial }{\partial \lambda}\left(\eta\dot{\mu} \right) + \frac{\partial }{\partial \mu}  \left\{(\mu^2-1)\eta\dot{\lambda} \right\} \\
&- \nabla^2 \left\{\frac{1}{2}\left( {\dot \lambda} g_{\lambda\lambda} + {\dot \mu} g_{\mu\mu} \right)+R\bar{T} \tilde{\pi} +\phi \right\}\\
& - \frac{\partial }{ \partial \lambda}\left\{ \dot{\sigma}\frac{\partial \dot{\lambda} }{\partial \sigma}
+ \frac{RT' }{a^2(1-\mu^2)}\frac{\partial \tilde{\pi}}{\partial \lambda} \right\} -\frac{\partial }{ \partial \mu}\left\{ \dot{\sigma}\frac{\partial \dot{\mu} }{\partial \sigma}
+ RT' \frac{(1-\mu^2) }{a^2}\frac{\partial \tilde{\pi}}{\partial \mu}\right\}  \\
&=\frac{1}{a(1-\mu^2)}\frac{\partial }{\partial \lambda}\left\{\eta V -\dot{\sigma}\frac{\partial U}{\partial \sigma}
- \frac{RT' }{a}\frac{\partial \tilde{\pi}}{\partial \lambda}   \right\} \\
&\qquad\qquad- \frac{1}{a}\frac{\partial }{\partial \mu}\left\{\eta U +\dot{\sigma}\frac{\partial V}{\partial \sigma}
+ RT' \frac{(1-\mu^2) }{a}\frac{\partial \tilde{\pi}}{\partial \mu} \right\} \\
&\qquad\qquad- \nabla^2 \left\{\frac{U^2+V^2}{2(1-\mu^2)}+R\bar{T} \tilde{\pi} +\phi \right\}\\
\end{align}
```
と書ける。ここで、
```math
\begin{align}
a F_U &\equiv \eta V -\dot{\sigma}\frac{\partial U}{\partial \sigma}
- \frac{RT' }{a}\frac{\partial \tilde{\pi}}{\partial \lambda}   \\
a F_V &\equiv -\eta U -\dot{\sigma}\frac{\partial V}{\partial \sigma}
- RT' \frac{(1-\mu^2) }{a}\frac{\partial \tilde{\pi}}{\partial \mu}\\
\end{align}
```
とおけば、$\alpha$-演算子を用いて渦度・発散方程式は、
```math
\begin{align}
\frac{\partial \eta }{\partial t} &= \alpha(F_V,-F_U)　　\\
\frac{\partial \delta}{\partial t} &=\alpha(F_U,F_V) - \nabla^2 \left\{\frac{U^2+V^2}{2(1-\mu^2)}+R\bar{T} \tilde{\pi} +\phi \right\}\\ 　
\end{align}
```
そして、熱力学方程式は、
 ```math
\begin{align}
\frac{\partial T'}{\partial t} &=- \frac{1}{a}\alpha(UT',VT') + T'\delta - \dot{\sigma}\frac{\partial T}{\partial \sigma}+\kappa T\frac{\omega}{P} + \frac{Q}{C_p}\\
\end{align}
```
となる。

# まとめ

　球極座標におけるプリミティブ方程式は以下となる。
### 予報方程式
 ```math
\begin{align}
\frac{\partial \eta }{\partial t} &= \alpha(F_V,-F_U)　　\\
\frac{\partial \delta}{\partial t} &=\alpha(F_U,F_V) - \nabla^2 \left\{\frac{U^2+V^2}{2(1-\mu^2)}+R\bar{T} \tilde{\pi} +\phi \right\}\\ 　

\frac{\partial \tilde{\pi}}{\partial t} &= -\int_0^{1} \mathrm{d}\hat{\sigma}A \\
\frac{\partial T'}{\partial t} &=- \frac{1}{a}\alpha(UT',VT') + T'\delta - \dot{\sigma}\frac{\partial T}{\partial \sigma}+\kappa T\frac{\omega}{P} + \frac{Q}{C_p}\\
\end{align}
```

### 診断方程式

```math
\begin{align}
A &= \delta + \frac{U}{a (1-\mu^2)} \frac{\partial \tilde{\pi}}{ \partial \lambda}+ \frac{V}{a}\frac{\partial \tilde{\pi}}{ \partial \mu} \\
 \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+{\bf u}\cdot\nabla \tilde{\pi} \\
\dot{\sigma}(\sigma) &= \sigma\int_0^{1} \mathrm{d}\hat{\sigma}A -\int_0^{\sigma} \mathrm{d}\hat{\sigma}A
\end{align}
```

### その他
```math
\begin{align}
 \alpha(A,B) &\equiv\frac{1}{1-\mu^2}\frac{\partial A}{\partial \lambda} + \frac{\partial B}{\partial \mu}\\

U &\equiv  a^{-1}\left\{(\mu^2-1) \frac{\partial \psi}{ \partial \mu}+\frac{\partial \chi}{ \partial \lambda} \right\} \\
V &\equiv a^{-1} \left\{\frac{\partial \psi}{ \partial \lambda}+ (1-\mu^2)\frac{\partial \chi}{ \partial \mu}\right\} \\

\zeta &= \frac{1}{a}\alpha(V,-U) ={\bf \nabla }^2 \psi \\
\delta &= \frac{1}{a}\alpha(U,V) ={\bf \nabla }^2 \chi \\

a F_U &\equiv \eta V -\dot{\sigma}\frac{\partial U}{\partial \sigma}
- \frac{RT' }{a}\frac{\partial \tilde{\pi}}{\partial \lambda}   \\
a F_V &\equiv -\eta U -\dot{\sigma}\frac{\partial V}{\partial \sigma}
- RT' \frac{(1-\mu^2) }{a}\frac{\partial \tilde{\pi}}{\partial \mu}\\

\eta &= \zeta +f \\
f &= 2\Omega \mu \\
T(\lambda,\mu,\sigma,t) &=  \bar{T}(\sigma) +T'(\lambda,\mu,\sigma,t) 
\end{align}
```

　解くべき方程式は準備できたので、次回は方程式の離散化方法について紹介したい。

# 参考文献

以下の書籍を参考にしています。

https://www.amazon.co.jp/-/en/Martin-Ehrendorfer/dp/1611971985

-