% 表題 DCPAM5 座標系・変換公式に関する解説 -- 波数切断 % % 履歴 %\Drireki{1991/11/12 保坂征宏} \section{波数切断} \Dseclab{波数切断} GCM では, 物理量を 球面調和函数 $P_n^m(\sin \varphi) \exp(im \lambda)$ で展開したり 波数空間で計算するときに, 計算資源の都合上, ある一定波数以下の波数のみを考慮して計算する. そのことを波数切断するという \footnote{ 後述するように, 現実的には波数切断を決めると同時に格子点数が決まる. すなわち, 以上の理由は 格子点数を大きくとれないことの理由でもある. } . 以下ではまず, 切断の基礎知識として切断の仕方・流儀を述べ, ついで, 切断における事情を述べた上で切断波数の決め方を記す. \subsection{波数切断の仕方} 波数切断の仕方については, 東西波数 ($m$) , 南北波数 ($n-m$) の それぞれの切断の方法にいくつかの流儀がある. 一般によく用いられるものは 三角形切断 (Triangle) , 平行四辺形切断 (Rhomboidal : 偏菱形) と 呼ばれるものである. 三角形切断の場合について計算する波数領域を波数平面上に書くと \Deqref{三角形切断}のようになる. 平方四辺形切断の場合は, \Deqref{平方四辺形切断}である. \begin{figure}[hbtp] \begin{center} \begin{picture}(300,100)(10,-10) \put(20 , 10){\vector(1,0){85}} \put(30 , 0){\vector(0,1){85}} \put(-5 , 80) {\shortstack{$n-m$}} \put(100 , 0) {\shortstack{$m$}} \put(30 , 60){\line(1,-1){50}} \put(10 , 55) {\shortstack{$N_{tr}$}} \put(70 , 0) {\shortstack{$N_{tr}$}} \put(150, 10){\vector(1,0){85}} \put(160 , 0){\vector(0,1){85}} \put(145 , 80){\shortstack{$n$}} \put(230, 0){\shortstack{$m$}} \put(160, 60){\line(1,0){50}} \put(160, 10){\line(1,1){50}} \put(140, 55) {\shortstack{$N_{tr}$}} \put(200, 0) {\shortstack{$N_{tr}$}} % \put(100,-20) {\shortstack{三角形切断}} \end{picture} \end{center} \caption{三角形切断の場合の波数領域} \Deqlab{三角形切断} \end{figure} \begin{figure}[hbtp] \begin{center} \begin{picture}(300,170)(10,-20) \put(20 , 10){\vector(1,0){85}} \put(30 , 0){\vector(0,1){85}} \put(-5 , 80) {\shortstack{$n-m$}} \put(100 , 0) {\shortstack{$m$}} \put(30 , 60){\line(1,0){50}} \put(80 , 10){\line(0,1){50}} \put(10 , 55) {\shortstack{$N_{tr}$}} \put(70 , 0) {\shortstack{$N_{tr}$}} \put(150, 10){\vector(1,0){85}} \put(160 , 0){\vector(0,1){135}} \put(145 ,130){\shortstack{$n$}} \put(230, 0){\shortstack{$m$}} \put(160, 60){\line(1,1){50}} \put(210, 60){\line(0,1){50}} \put(160, 10){\line(1,1){50}} \put(140, 55) {\shortstack{$N_{tr}$}} \put(135,105) {\shortstack{$2N_{tr}$}} \put(200, 0) {\shortstack{$N_{tr}$}} % \put(90,-20) {\shortstack{平方四辺形切断}} \end{picture} \end{center} \caption{平方四辺形切断の場合の波数領域} \Deqlab{平方四辺形切断} \end{figure} 三角形切断, 平行四辺形切断, という名称は 波数平面上 ($(n,m)$平面) での形状による \footnote{ 平方四辺形切断には, $n$ の最大値を $m$ の最大値の2倍に しないようなとり方もある. 詳しくは五角形切断に関する脚注参照.}. より一般的な切断方法は五角形切断 (\Deqref{五角形切断}) である. \begin{figure}[hbtp] \begin{center} \begin{picture}(300,170)(10,-20) \put(20 , 10){\vector(1,0){95}} \put(30 , 0){\vector(0,1){85}} \put(-5 , 80) {\shortstack{$n-m$}} \put(110 , 0) {\shortstack{$m$}} \put(30 , 60){\line(1, 0){40}} \put(70 , 60){\line(1,-1){20}} \put(90 , 10){\line(0, 1){30}} \put(10 , 55) {\shortstack{$J$}} \put(0 , 35) {\shortstack{{\footnotesize $K-M$ }}} \put(83 , 0) {\shortstack{$M$}} \put(50 , 0) {\shortstack{{\footnotesize $K-J$ }}} \put(150, 10){\vector(1,0){95}} \put(160 , 0){\vector(0,1){125}} \put(145 ,120){\shortstack{$n$}} \put(240, 0){\shortstack{$m$}} \put(160, 60){\line(1,1){40}} \put(200, 100){\line(1,0){20}} \put(220, 70){\line(0,1){30}} \put(160, 10){\line(1,1){60}} \put(145, 55) {\shortstack{$J$}} \put(145, 95) {\shortstack{$K$}} \put(213, 0) {\shortstack{$M$}} \put(180, 0) {\shortstack{ {\footnotesize $K-J$}}} % \put(100,-20) {\shortstack{五角形切断}} \end{picture} \end{center} \caption{五角形切断の場合の波数領域} \Deqlab{五角形切断} \end{figure} 三角形切断, 平行四辺形切断はそれぞれ, 五角形切断において \begin{itemize} \item 三角形切断 \ \ \ \ \ $ J=K=M =N_{tr} $ \item 平行四辺形切断 $ K=2N_{tr},\ J=M=N_{tr} $ \end{itemize} であるような特別な場合である \footnote{ 単に $K=J+M$ であるものも平方四辺形切断と呼ばれる. だが, 例えば R21 と呼ばれるものは, $K=42, J=M=21$ のものである. } . 三角形切断と平行四辺形切断の違いについて, 世の中では次のように言われている \footnote{ 気象庁予報部, 1982 の p.47 より. } . \begin{itemize} \item 三角形切断の水平分解能は, 経度方向のみならず緯度方向にも一定である \footnote{ 分解能が緯度方向に変化することについては, 平行四辺形切断に限らず, 三角形切断以外のどれでも起こる.} . 分解能を上げて スケールの細かい波を表現できるようになった場合を考える. 物理的に スケールの小さい波には指向性がないことと, 水平分解能に方向依存性がないこととは調和的である. \\ また, このことは, ある三角形波数切断した球面調和函数により表現される 球面上の分布は 極の位置を変えても 同じ三角形波数切断した球面調和函数により正確に表現される ことの言い替えでもある. \item 平行四辺形切断の場合, 各東西波数について同じだけの南北波数をとれる. \end{itemize} \subsection{切断波数の決め方} ここでは切断波数と南北格子点数の決め方について記す. これらは 切断の仕方を決めた後に, 使用する計算資源がネックになって決まる. その際, FFT の仕様, aliasing の回避, という２つの数値的な事情を 考慮した上で決める必要がある. FFT の仕様の事情というのは, 話は簡単で, 東西方向に 「格子 $\Leftrightarrow$ スペクトル」 変換する ために用いる FFT が 効率よく動くための格子点数・波数がある \footnote{ コード依存性がある. 通常, 2のべき乗が好ましいとされる. コードによっては, 2,3,5 のべき乗の積でもよいものもある.} ことである. 一方, aliasing に関する事情は複雑である. ここで扱っているスペクトルモデルでは, 格子点でのみ値を計算している. いわゆるスペクトルを使うのは, 単に格子点上での水平微分項の評価をする時のみである. その意味で, 「微分の評価にのみスペクトルを用いるグリッドモデル」 と言ってもよい. そのように受け止めると, 格子点値を"正しく"計算することを目指し, また, 考慮する波数は 厳密にスペクトルの係数と格子との変換を行なうことのできる波数, すなわち変換において情報の落ちないだけの波数を とらねばならないように思える. ところが実際には, スペクトルモデル的な配慮 --- ある波数以下についてのみ 正しく計算し, それ以上の波数については計算しない --- により 切断波数・格子点数が決められている. また, 後述する理由により 情報は (非線形 aliasing のことを考えずとも) 必ず落ちてしまうのである \footnote{ 実際の GCM では 格子点値からスペクトルに変換する際に情報は落ちている. したがって, 格子 - スペクトル - 格子という変換を 行なうと元にはもどらない. \\ 例えば T42 の場合, 自由度は $ 1+ (2 \times1 +1) + \cdots + (2 \times 42 +1) = 43^2 = 1849$ に対して 格子点数は $ 128 \times 64 = 8192 $ である. R21 の場合も, 自由度は $ (2 \times 21 +1) \times (21+1)= 946$ に対して, 格子点数は $ 64 \times 64 = 4096 $ である. すなわち, $3/4$ 以上の情報は 格子点値からスペクトルに変換するときに落ちている. \\ 工夫すれば 情報が落ちないうまい方法があるかも知れないが, 今のところ見つけていないし多分見つからない. \\ もちろん, スペクトル - 格子 - スペクトルという変換では 元にもどる (ように決めている) . } . さて, 以下では aliasing に関する事情を具体的に述べながら, 切断波数に対する格子点数の決め方を記そう. 球面上に連続分布している物理量を球面調和函数で展開する. ある波数 $M,N(m)$ 以下 (例えば, T42 ならば $M=42, N=42$ ) については 線形項・非線形項の両方について 厳密に計算できるように $I,J$ を決めることを目指す. $M,N$ を仮に固定したとして, まずは線形項について 切断波数以下のスペクトルの係数のわかっている物理量 $A$ を 格子点値に変換しさらにスペクトルの係数に正しくもどすことを考える. $A$ は $ -M \le m \le M,\ |m| \le n \le N(m)$ の $m,n$ については $\tilde{A}_n^m$ が わかっているとする. 格子点値は, $1 \le i \le I,1 \le j \le J$ について \begin{align} A_{ij} &\equiv \sum_{m=-M}^M \sum_{n=|m|}^N \tilde{A}_n^m P_n^{m}(\sin \varphi_j) \exp(im \lambda_i) \end{align} で与えられる. これらの格子点値から逆に $\tilde{A}_n^m ( -M \le m \le M, \ |m| \le n \le N) $ を計算する. 離散化した系での積分を Gauss の公式, Gauss-Legendre の公式で評価すれば, \begin{align} \tilde{A}_n^m &= \frac{1}{I} \sum_{i=1}^I \sum_{j=1}^J A_{ij} P_n^{m}(\sin \varphi_j) \exp(-im \lambda_i) w_j \end{align} である. ここで, $w_j$ は $\varphi_j$ における重みである. $A_{ij}$ の定義を代入すれば, \begin{align} \tilde{A}_n^m &= \frac{1}{I} \sum_{i=1}^I \sum_{j=1}^J \left( \sum_{m'=-M}^{M} \sum_{n'=|m'|}^N \tilde{A}_{n'}^{m'} P_{n'}^{m'}(\sin \varphi_j) \exp(im' \lambda) \right) P_n^{m}(\sin \varphi_j) \exp(-im \lambda_i) w_j \nonumber \\ &= \frac{1}{I} \sum_{m'=-M}^{M} \sum_{n'=|m'|}^N \tilde{A}_{n'}^{m'} \sum_{i=1}^I \exp(i(m'-m) \lambda) \sum_{j=1}^J P_n^{m}(\sin \varphi_j) P_{n'}^{m'}(\sin \varphi_j) w_j \end{align} となる. この計算が $\tilde{A}_n^m$ を 正しく評価している (すなわち元にもどる) ための $I,J$ の条件は, $ -M \le m \le M,\ |m| \le n \le N$ を満たす $m,n$ について \begin{align} & \sum_{i=1}^I \exp(i(m'-m) \lambda) = I \delta_{mm'} , \\ & \sum_{j=1}^J P_n^{m}(\sin \varphi_j) P_{n'}^{m}(\sin \varphi_j) w_j = \delta_{nn'} \end{align} が成り立つことである. 三角関数の和による評価が正しいための条件は, ここに登場する波数 $|m'-m|$ が最大で $2M$ の値をとるので, Gauss の公式の適用条件より, 格子点数 $I$ が $I \ge 2M+1$ を満たすことである. Legendre函数の積の和による評価が正しいための条件は, ここに登場する計算が $n+n'$ 次の多項式 \footnote{ ここで, 三角関数の和が $I \delta_{mm'}$ となることを用いた. 一般には ($m,m'$ の偶奇が一致しない場合には) $ P_n^m P_{n'}^{m'} $ は多項式にならない. } の評価であることから, Gauss - Legendre の公式の適用条件より, 格子点数 $J$ が $2J-1 \ge {\rm max}[n+n'] = 2 {\rm max}[N]$ を 満たすことである. ここで, ${\rm max}[n+n']$ は $n+n'$ の最大値を, ${\rm max}[N]$ は $N$ の最大値を表す. ちなみに, 格子点値からスペクトルの係数に変換し格子点値にもどす という立場からすれば, この Gauss-Legendre の公式の適用条件というのが 情報を落とさずには済まない理由である \footnote{ Gauss の公式の適用条件と情報欠落との関係について コメントしておく. 格子点数 $I$ が奇数の場合には, スペクトルで同じ情報量を持つためには 波数 $\frac{I-1}{2}$ までを考慮すればよいので, 情報は欠落しないことは明らかである. 一方, $I$ が偶数の場合には, 情報は欠落させないためには 波数 $\frac{I}{2}$ が必要であるが, この波数は Gauss の公式の適用条件を満たさない. しかしこの場合にも, (私は根拠を調べていないが, 少なくとも) 経験的には FFT および 逆FFT によって 格子 - スペクトル - 格子変換によって 情報が落ちないことが知られている. } . このことを以下に述べる. 情報を落とさずに 格子点値をスペクトルの係数に変換し格子点値にもどすには, あらゆる東西波数について 南北方向の格子点数 $J$ と同じだけの個数のLegendre函数が必要である. 東西波数 $m$ の場合, 登場するLegendre陪函数の $n$ は $n=|m|, |m|+1, \cdots, |m|+J-1$ である. $ P_n^m P_{n'}^{m} $ の次数は $n+n'$ であるから, 最大で $2J+2|m|-2$ である. これが $2J-1$ 以下になるのは $m=0$ の時のみである. $m \neq 0$ の場合は高次のLegendre函数は計算してはならない. つまり情報を落とさざるをえない \footnote{ この事情により, 非線形項の場合を考えて さらに著しく落とすことが必要になることが次節からわかる. } . 改めて $M,N$ を固定するという立場にもどって, 切断波数以下のスペクトルの係数のわかっている物理量 $B,C$ の積から それらの格子点値を用いて $B$ と $C$ との積 (非線形項) $A$ のスペクトルの係数を正しく求めるための $I,J$ の条件を考える. \begin{align} & A= BC , \\ & B= \sum_{m=-M}^{M} \sum_{n=|m|}^{N} \left( \tilde{B}_n^m \exp(im \lambda) \right) P_n^m(\sin \varphi) , \\ & C= \sum_{m=-M}^{M} \sum_{n=|m|}^{N} \left( \tilde{C}_n^m \exp(im \lambda) \right) P_n^m(\sin \varphi) \end{align} なる物理量 $A,B,C$ があるとする \footnote{ $A,B,C$ とも 実数である. すなわち, $\tilde{B}_n^m = \tilde{B}_n^{m*}, {\it etc.}$ となっている. } . $B$,$C$ の $-M \le m \le M,\ |m| \le n \le N$ における スペクトルの係数 $\tilde{B}_n^m,\tilde{C}_n^m$ を用いて $A$ の スペクトルの係数 $\tilde{A}_n^m$ を $0 \le m \le M,\ |m| \le n \le N$ については 正しく計算することを考える. \begin{align} \tilde{A}_n^m &\equiv \widetilde{(BC)_n^m} \nonumber \\ &= \frac{1}{I} \sum_{i=1}^I \sum_{j=1}^J B_{ij} C_{ij} P_n^{m}(\sin \varphi_j) \exp(-im \lambda_i) w_j \nonumber \\ &= \frac{1}{I} \sum_{i=1}^I \sum_{j=1}^J \left( \sum_{m'=-M}^{M} \sum_{n'=|m'|}^N \tilde{B}_{n'}^{m'} \exp(im' \lambda_i) P_{n'}^{m'}(\sin \varphi_j) \right) \nonumber \\ & \quad \times \left( \sum_{m''=-M}^{M} \sum_{n''=|m''|}^N \tilde{C}_{n''}^{m''} \exp(im'' \lambda_i) P_{n''}^{m''}(\sin \varphi_j) \right) P_n^{m}(\sin \varphi_j) \exp(-im \lambda_i) w_j \nonumber \\ &= \frac{1}{I} \sum_{m'=-M}^{M} \sum_{n'=|m'|}^N \sum_{m''=-M}^{M} \sum_{n''=|m''|}^N \tilde{B}_{n'}^{m'} \tilde{C}_{n''}^{m''} \nonumber \\ & \quad \times \sum_{i=1}^I \exp(i(m'+m''-m) \lambda_i) \sum_{j=1}^J P_{n'}^{m'}(\sin \varphi_j) P_{n''}^{m''}(\sin \varphi_j) P_n^{m}(\sin \varphi_j) w_j . \end{align} この計算が $\tilde{A}_n^m$ を $0 \le m \le M,\ |m| \le n \le N$ について 正しく評価しているための, $I,J$ の条件を 線形項の場合と同様に考えると, 格子点数 $I$ が $I \ge 3M+1$ を, 格子点数 $J$ が $2J-1 \ge {\rm max}[n+n'+n''] = 3 {\rm max}[N]$ を 満たすことである. ここで, ${\rm max}[n+n'+n'']$ は $n+n'+n''$ の最大値を, ${\rm max}[N]$ は $N$ の最大値を表す. 再び 格子点値からスペクトルの係数に変換し格子点値にもどす という立場からすれば, これらの $I,J$ に関する条件から, 南北成分のみならず, 東西成分についても 変換によって情報が落ちてしまうことがわかる. これまでに述べた $M,N$ を固定したときに 格子点数 $I,J$ がとらねばならない個数について, 線形項・非線形項の２つの場合のうち条件が厳しいのは, 明らかに非線形項の場合である. この条件以下の格子点数しかとらない場合には, aliasing をおこすことになる. 以上, FFT, aliasing という２つの事情を考えて 格子点数と切断波数とは同時に決められる. 具体的手順は以下のとおりである. \begin{enumerate} \item 波数切断の仕方を決める. \item FFT のかけやすい数を選ぶ. それを東西格子点数 $I$ とする. \item 東西方向の波数の最大値 $M$ を ${\displaystyle M= \left[ \frac{I-1}{3} \right] }$ にする. ただし $[ \ ]$ はそれを越えない最大の整数を表す記号である. \item 最大全波数 $N_{{\rm max}}$ を決める. 三角形切断ならば $N_{{\rm max}}=M$, 平行四辺形切断ならば $N_{{\rm max}}=2M$ である. \item 南北方向の格子点数 $J$ を $J \ge \frac{3N_{{\rm max}}+1}{2}$ を 満たす数に選ぶ. (\Dmodel では偶数でなくてはならない. ) \end{enumerate} 例えば, T42 の場合には $M=42, N=42$, 東西格子点数 $I$ が 128, 南北格子点数 $J$ が 64 である. R21 の場合には $M=21, N=42$, 東西格子点数 $I$ が 64, 南北格子点数 $J$ が 64 である. 参考までに, 線形モデルの場合について決め方を示しておく. \begin{enumerate} \item 波数切断の仕方を決める. \item FFT のかけやすい数を選ぶ. それを東西格子点数 $I$ とする. \item 東西方向の波数の最大値 $M$ を ${\displaystyle M= \left[ \frac{I}{2} \right] }$ にする. ただし $[ \ ]$ は それを越えない最大の整数を表す記号である \footnote{ ここで, $I$ が偶数のときについては Gauss の公式の適用条件を越えて 最大波数 $\frac{I}{2}$ まで 計算できるという知識を用いた. } . \item 最大全波数 $N_{{\rm max}}$ を決める. 三角形切断ならば $N_{{\rm max}}=M$, 平行四辺形切断ならば $N_{{\rm max}}=2M$ である. \item 南北方向の格子点数 $J$ を $J \ge \frac{2N_{{\rm max}}+1}{2}$ を 満たす数に選ぶ. \end{enumerate} 例えば, 三角形切断の場合には, $I=128$ とすると, $M=64$, $N=64$, $J=65$ となる. つまり T64 では $I=128, J=65$ である. 平方四辺形切断の場合には, $I=64$ とすると, $M=32$, $N=64$, $J \ge 65$ となる. つまり R32 では $I=64, J=65$ でよい \footnote{ これらの場合でも, 南北方向の細かい情報は 格子 - スペクトル - 格子変換によって 落ちていることに注意せよ. } .