\chapter{序論}
% AGCM における静力学コアから非静力学コアへの移行
全球的な数値天気予報や気候予測において, 
大気大循環モデル (Atmospheric General Circulation Model: AGCM) は重要な道具である. 
継続的な計算機技術の発展は, 全球規模の気象シミュレーションにおいても水平解像度が数 km 以下の高解像度計算を可能にしつつある. 
伝統的な AGCM は静水圧平衡近似を用いるプリミティブ方程式を使用してきたが,  
数 km 以下の水平スケールの運動に対しては静水圧平衡の仮定が成り立たないため, 
近年開発される AGCM において, 静水圧近似を仮定しない非静力学方程式が使われる傾向にある. 
研究領域における非静力学モデルの歴史は 30 年ほどあるが, 
現業の数値予報分野で非静力学モデルが使われるようになったのは 10 年前ほどからであり, さらに領域モデルに限られていた. 
このような全球非静力学モデルとして, 
日本では, 東京大学気候システム研究センター(CCSR)と地球環境フロンティア研究センター(FRCGC)で共同開発された, 
非静力学正二十面格子大気モデル NICAM\citep{Tomita2004} がある. 
他の正二十面格子を使った全球非静力学モデルには, 
アメリカ海洋大気庁(NOAA) による, 非静力学正二十面格子モデル NIM や海洋 - 陸 - 大気モデル OLAM などが挙げられる. 
また, 立方球面格子を用いた全球モデルでは, 地球流体力学研究所(GFDL) と アメリカ航空宇宙局(NASA) のゴダード宇宙飛行センター (GSFC) による, 
GEOS\_FV\_CUBED\citep{PutmanAndLin2008}がある.   
一方, ヨーロッパ中期予報センター(ECMWF) は
セミラグランジュ・スペクトル変換法による全球非静力学モデル IFS\citep{ECMWF_IFS}を開発している.

% 非静力学モデルにおける数値解法の動向
伝統的にほとんどの AGCM の数値解法には, 格子点法に比べ高い数値精度を持つ, 球面調和関数展開に基づくスペクトル変換法が採用されてきた. 
しかし, 高性能計算(High Performance Computing: HPC)で実現されるような高解像度計算では,  
ルジャンドル変換の O($n^3$) ($n$は緯度経度格子における一次元的な格子数)の計算量や大規模並列計算機におけるノード間データ転送がネックになることが指摘されている.  
一方, アメリカ国立大気研究センター(NCAR)の WARF\citep{SkamarockAndKlemp2008}, 
イギリス気象庁(UK Met Office)の UM\citep{UKMet_UM}など現業の多くの領域非静力学モデルは, 
空間離散化に格子点法の一種である有限差分法を用いる. 
%, 時間離散化に時間分割法(HEVE法, HEVI法, HIVI 法)を用いている. 
有限差分法は効率的であるが, 位相誤差, 領域形状に対する柔軟性, 並列計算におけるスケーラビリティの問題が存在する. 

これらの欠点を克服するために, 近年開発されている領域および全球非静力学モデルでは, 
有限体積法(Finite Volume Method: FVM), 有限要素法(Finite Element Method: FEM), スペクトル要素法(Spectral Element Method: SEM), 
不連続 Galerkin 法(Discontinuous Galerkin Method: DGM)などがしばしば用いられる.  
これらはもともと数値流体力学(Computational Fluid Dynamics: CFD)分野で活発に研究されてきたものであり, 
気象あるいは地球流体力学分野での応用は 10 年前ほどからである.  
気象分野での FVM の応用例は, 正二十面格子全球非静力学モデル NICAM などが挙げられる.   
全球モデルの難しさの一つは極の特異性であり, 緯度経度格子の極における格子線の収束は時間積分の時間スッテプの大きさに厳しい制約を与える.  
NICAM では六角型の非構造格子を水平格子に採用し, 有限体積法を用いて空間離散化することで, 極の特異性を回避している. 

% 次世代非静力学モデルの数値解法に求められるもの  
次世代の非静力学モデル開発における数値解法と関連した問題は, 
%%
\begin{enumerate}
 \setlength{\itemsep}{-2pt}
 \item 空間高精度化, 
 \item 計算領域形状に対する柔軟性, 
 \item 優れた分散特性, 
 \item 並列計算機におけるスケーラビリティ
\end{enumerate}
%%
などである. 
1 の空間高精度化は, 計算コストを抑えながら, 大気の流れの詳細な構造を解像するためには不可欠である. 
また, 2 は急峻な地形を始めとする複雑な境界の取り扱いが必要な場合に重要である. 
1,2 と関連することとしては, 静的あるいは動的な格子の局所的な細分化が挙げられる. 
3 の分散特性については, 低解像度計算においても大気波動の分散関係を正しく再現できる性質が望まれる.  
4 と関係することとしては, 
(LINPACK ベンチマークによる)スーパーコンピュータTop500 が反映しているように, 1990 年代後半からの HPC の傾向として, 
ベクトル型よりもスカラ型のスーパーコンピュータ(例えばクラスタ)の開発に重点が置かれている事実である.  
ベクトル型と比較したスカラ型のスーパーコンピュータの特徴は, 
安価な汎用マルチコアプロセッサを大量に用いて並列計算を行うところである. 
例えば, 2012 年に完成した理化学研究所の京(けい)の総コア数は 50 万個以上に及ぶ. 
このような, 膨大な数のプロセッサコアを持つ計算機の利点を生かすには, 
一プロセッサに対する演算量が多く, ノード間でのデータ通信の少ない, より局所的な数値解法が望ましいであろう. 
したがって, このような HPC の動向も考慮すると, 
新しく作られる全球非静力学モデルに採用される数値解法には, 
(メモリ空間上で)局所的なデータに対して多くの演算を行うことにより高精度化する, 
コンパクトステンシルを用いる高精度 FVM や FEM, SEM, DGM といった手法がやはり適していると思われる. 
また, これらの手法は 1, 2, 3 の問題も得意とするところである. 

前述したように, 地球流体力学的問題に対して FVM, FEM, SEM, DGM を適用した研究は, 2000 年代に入ってから見られるが, 
当初は全球浅水モデルであったり, 全球非静力学モデルであっても古典的な FVM (空間 2 次精度)に限られていた.    
特に, 領域または全球非静力学モデルに対する適用は数年前からであり, まだ広く研究されていない.  
これらの研究例には, 高精度 FVM による全球非静力学モデル Mcore\citep{Ullrich2012}, SEM または DGM による領域非静力学モデル\citep{Giraldo2008}などが挙げられる. 
このため, 本研究では CENO 法を用いた高精度 FVM による全球非静力学モデル DECAMO(DEnnou Cubed sphere grid point Atmospheric MOdel)を提案する. 

ここからは, 本モデルの特徴を簡単に述べることにする. 
まず, 全球モデルに伴う極問題を解決するために, 格子系には立方球面格子を用いる. 
なお, \citet{Taylor2008}や\citet{Putman2007}によって, 
立方球面格子を用いた数値モデルは, 超並列計算機においてスケーラビリティが保たれることが示されている. 
また, 本モデルの支配方程式は, 保存形式の密度, 運動量, 全エネルギーの時間発展式で構成される完全圧縮方程式系であり,  
これらの空間離散化には有限体積法を用いる. 
これにより, 本質的には系の全質量, 全運動量, 全エネルギーは保存される. 
ただし, 非線形方程式の反復計算の残差や浮動小数演算の丸め誤差のために, 非常に小さな変動は存在する.  

有限体積法の高精度化において, セル界面での場の値を正確に与えるための再構築(reconstruction)の方法と, セル界面における数値流束の評価方法の選択が重要となる. 
再構築において, 検査体積内のサブグリッドスケールの分布は, 
近傍にある複数の検査体積(ステンシル)の体積積分平均値の情報を用いて生成された, 補間多項式で近似される. 
並列計算の効率性の観点から, 再構築においてできるだけコンパクトなステンシルが設定されることが望ましい. 
CFD 分野においてこの再構築の技術は古くからよく研究されており, 
近年では WENO 法(Weighted Essentially Non-Oscillatory Method)\citep{Liu1994}などが有名である. 
これらの手法には, 単調性を保証する工夫と一緒になっている場合が多い. 
本研究では, WENO 法よりもより効率的で安定な Central ENO 法\citep{IvanPHDTheisis, Ivan-et-al2010_ICCFD7}を用いる. 
これらの手法は, 原理的には, より高次の補間多項式を生成することによって任意空間精度の有限体積法を実現することができる. 

次に, 検査体積の界面における数値流束の評価方法の選択について述べる. 
数値流束の評価する方法としては, 近似リーマン解法(例えば, Roe スキーム\citep{Roe1981})や比較的新しいものには風上分離法がある. 
伝統的にはこれらの方法は, 大気のような低マッハ数の流れよりもむしろ, 遷音速や超音速のレジームにある流れを模倣するために設計されている. 
そのため, 低マッハ数の流れに対して近似リーマン解法を用いた場合に方程式系が stiff になるため, 
陰解法における反復計算が収束しなかったり, また陽解法でも数値精度が悪いなどの問題がある. 
しかし, 1990 年代後半から, 近似リーマン解法や風上分離法は低マッハ数流れに対しても適用できるように拡張されており,  
前処理を用いた近似リーマン解法\citep{WeissAndSmith1995}や低マッハ数流れ用の風上分離法\citep{Liou2006}などが挙げられる. 
本研究では, \citet{Ullrich2012}と同様に, 計算効率と数値精度の観点から, \citet{Liou2006}で提案されたAUSM$^{+}$-up スキームを数値流束の評価に用いる.  

全球モデルでは, 依然として鉛直格子幅は水平格子幅に比べ, $10^{-2} \sim 10^{-3}$倍ほど小さい. 
特に, 速い位相速度で鉛直方向に伝播する音波や重力波を含む非静力学モデルでは, 
この相対的に小さな鉛直方向の格子幅によりCFL条件(Courant-Friedrichs-Lewy Condition: CFL condtion)を満たす時間ステップ幅は, 
総観規模の運動がもつ時間スケールに比べて非常に小さくなる. 
したがって, 多くの全球モデルでなされるように, 本研究でも鉛直方向の運動に対する時間離散化は陰的に行う. 
時空間精度を保持しながら, 
水平方向の運動に対する陽的な積分と鉛直方向の運動に対する陰的な積分を組み合わすために, 
\citet{Ullrich2012}で採用された演算子分割 Runge-Kutta-Rosenbrock スキームを用いる. 
今回示す数値解法では, 水平方向に関しては空間精度は 3 次精度, 鉛直方向に関しては空間 2 次精度の再構築を行う.  
時間精度に関しては, 水平方向の陽的な時間積分には, 4 次の Runge-Kutta 法を用いる. 
演算子分割 Runge-Kutta-Rosenbrock スキームは, 時間積分一ステップあたりに陽解法と陰解法を複数回組み合わせることで, 
次元分割していても最高 3 次精度の時間積分を可能にする.  
今回は \citet{Ullrich2012} に習って, 全体としての空間・時間精度が 2 次精度となるような方法を適用した.  
そのため, 結果的に全体としては水平空間精度の 3 次収束は保証されないが, \citet{Ullrich2012}が示すように, 
惑星波から中速度の重力波の振動数を持つ水平伝播する波の分散関係の表現が改善される. 
なお, 鉛直方向に対してもより高次の再構築を用い, かつ全体で 2 次精度以上を保つような演算子分割 Runge-Kutta-Rosenbrock スキームを適用すれば, 
時空間 2 次精度以上を実現できる. 

開発したモデルの性能を検証するために, 非静力学モデルの力学コアに対する標準テストケースである, 音波や内部重力波のテストケース, 傾圧不安定のテストケースを行う. 
これらの数値計算の結果は, 高精度 FVM と低マッハ数流れ向けの数値流束の評価方法のおかげで, 
古典的な FVM による非静力学モデルに比べて場の急勾配をより捕捉できる特性を示している. 

最後に, この修士論文の構成について説明する. 
第2章は, 立方球面格子を導入する. 
第3章では, 立方球面格子座標系において本モデルの支配方程式の定式化を行う.  
第4章において, 支配方程式を解くために用いられる数値解法について述べる. 
第5章では, 開発したモデルの数値精度や特性を検証するために行った, 
非静力学コアのための標準テストケースの結果を示す.  
第6章では, 結論および今後の開発方針や本モデルを用いた数値シミュレーションの可能性を議論する. 
