\chapter{物質移流} \section{はじめに} 本章ではセミラグランジュ法による物質移流の計算手順について述べる． 本章で述べるのは, 物質混合比$q$の移流方程式を解く, 非保存型のセミラグランジュ法である． 一般に, 非保存型セミラグランジュ法では, \begin{enumerate} \item 時刻$t+\Delta t$で格子点上にくる流体粒子の, 時刻$t$での位置（上流点）を求め, \item その位置の時刻$t$での$q$の値を補間計算で求め, \end{enumerate} それを時刻$t+\Delta t$での格子点での値とする計算を行う． すなわち, 計算手順は{\bf 上流点探索}と{\bf 補間}に分けられる． セミラグランジュ法を３次元で一度に解こうとすると上流点探索も補間計算も, 非常に複雑になり 計算コストが増大する． このため, DCPAMでは方向分離という手法を用いている．すなわち水平方向の物質移流を計算し, その分布を用いて鉛直方向の移流計算を行う． なお, セミラグランジュ法の一般的な利点のひとつはCFL条件を超えて時間刻み幅を大きくとれることである． しかし, DCPAMでは他の力学や放射などの計算と合わせるために, CFL条件を超えた時間刻み幅を用いることはない(2013年1月現在)． \section{移流方程式の数理表現} 大気密度を$\rho$とするとき, 単位体積中に含まれる物質の質量の時間変化は \begin{equation} \frac{\partial \rho q}{\partial t} = -\nabla \cdot (\bm{v}\rho q) + S \label{eq:1} \end{equation} で表される．ただし, $t$は時間, $\bm{v}$は流速を表し, $S$は生成・消滅項である． 物質の生成・消滅を考えない場合, (\ref{eq:1})式と$\rho$の時間変化式 \begin{equation} \frac{\partial \rho}{\partial t} = -\nabla \cdot (\bm{v}\rho) \label{eq:2} \end{equation} から, $q$の時間変化式（移流方程式）が得られる． \begin{equation} \frac{\partial q}{\partial t} = -\bm{v} \cdot \nabla q \label{eq:3} \end{equation} \section{上流点探索} DCPAMの時間積分はリープフロッグ法であるため, 時刻$t+\Delta t$ ($A$; After) の混合比を求めるために, $t$ ($N$; Now)および$t-\Delta t$ ($B$; Before)の情報を用いることが出来る． このため, $B$から$A$まで間の平均流速場を$N$の流速場で近似して, 上流点を探索する(Williamson and Rasch, 1989)． すなわち, 時刻$A$に格子点$\bm{x}^G_A$に流れてくる流体粒子の時刻$N$での位置を$\bm{x}_N$（中間点）, 時刻$B$での位置を$\bm{x}_B$（上流点）とすると, これらは次式を満たす． \begin{align} \bm{x}^G_A &= \bm{x}_N+\Delta t\bm{v}_N(\bm{x}_N) \label{eq:4}\\ \bm{x}^G_A & =\bm{x}_B+2\Delta t\bm{v}_N(\bm{x}_N) \label{eq:5} \end{align} 一般に$\bm{x}_N$, $\bm{x}_B$は格子点上の位置ではないので, (\ref{eq:4})式を満たす $\bm{x}_N$と$\bm{v}_N(\bm{x}_N)$を補間の繰り返し計算（イタレーション）によって求める．$\bm{v}_N(\bm{x}_N)$が求まれば, (\ref{eq:5})式より上流点$\bm{x}_B$が求められる． \subsection{水平} 水平方向の上流点探索は$u_N$, $v_N$に対して{\bf ２次元ラグランジュ３次補間}を 繰り返し行う． すなわち, 中間点位置の現在推定値($\lambda_N$, $\phi_N$)での水平流速 $u_N(\lambda_N$, $\phi_N)$, $u_N(\lambda_N$, $\phi_N)$を ２次元ラグランジュ３次補間で求め, (\ref{eq:4})式を出来るだけ満たすように 推定値を修正していく． なお, ２次元ラグランジュ３次補間とは, 補間したい座標($\lambda_N$, $\phi_N$)を 取り囲む16グリッドでの値を用いて, 多項式 \begin{equation*} f(\lambda, \phi) = \sum^3_{i=0} \sum^3_{j=0}a_{i,j}\lambda^i \phi^j \end{equation*} で分布を近似し, 係数$a_{i,j}$を求めて, 補間値$f(\lambda_N, \phi_N)$を求める 方法である． \subsection{鉛直} 鉛直方向の上流点探索はより単純で, $\dot{\sigma}_N$に対して {\bf １次元ラグランジュ３次補間}を繰り返す． \section{混合比補間} \subsection{水平} 球面上の高精度な補間方法として, スペクトル双３次補間(Enomoto, 2008)が知られている．これは, 球面調和関数変換を利用して, 混合比分布のグリッド上の空間微分値 $q_\lambda, q_\phi, q_{\lambda\phi}$ %\begin{equation*} %\frac{\partial q}{\partial \lambda}, \quad %\frac{\partial q}{\partial \phi}, \quad %\frac{\partial^2 q}{\partial \lambda\phi} %\end{equation*} を求め, 補間したい座標$(\lambda_B, \phi_B)$の周囲４点での$q, q_\lambda, q_\phi, q_{\lambda\phi}$の情報から多項式 \begin{equation} f(\lambda, \phi) = \sum^3_{i=0} \sum^3_{j=0}a_{i,j}\lambda^i \phi^j \label{eq:6} \end{equation} の係数$a_{i,j}$を求めて, 補間値$f(\lambda_B, \phi_B)$を求める方法である． この方法の特徴は, 微分値をスペクトル変換を利用して求めている点であり, このために精度の高い補間が可能になっている． DCPAMでは, 上述の方法を発展させた {\bf 方向分離型スペクトル変則エルミート５次補間}で混合比の水平方向の補間計算を行う． \subsubsection{方向分離} 双３次補間では補間したい座標の周囲の混合比分布を多項式(\ref{eq:6})で近似して, 周囲４グリッド上の計16個の情報から16個の係数$a_{i,j}$を求めていた． これは, 16元の連立方程式を解くことに相当する． ところが, 以下ように経度方向と 緯度方向の補間を分離することで, より簡便に補間することができる． ここで, 補間したい座標を$(\lambda_B, \phi_B)$, その周囲の４グリッドの座標を $(\lambda_i, \phi_j)$, $(\lambda_{i+1}, \phi_j)$, $(\lambda_i, \phi_{j+1})$, $(\lambda_{i+1}, \phi_{j+1})$とする． \begin{enumerate} \item $(\lambda_i, \phi_j)$と$(\lambda_{i+1}, \phi_j)$上の $q, q_\lambda$を用いて, 経度方向の１次元エルミート３次補間により座標$(\lambda_B, \phi_j)$での$q$を求める． \item $(\lambda_i, \phi_j)$と$(\lambda_{i+1}, \phi_j)$上の $q_\phi, q_{\lambda\phi}$を用いて, 補間により座標$(\lambda_B, \phi_j)$での$q_\phi$を求める． \item $(\lambda_i, \phi_{j+1})$と$(\lambda_{i+1}, \phi_{j+1})$上の $q, q_\lambda$を用いて, 補間により座標$(\lambda_B, \phi_{j+1})$での$q$を求める． \item $(\lambda_i, \phi_{j+1})$と$(\lambda_{i+1}, \phi_{j+1})$上の $q_\phi, q_{\lambda\phi}$を用いて, 補間により座標$(\lambda_B, \phi_{j+1})$での$q_\phi$を求める． \item 上の1--4で求めた$q(\lambda_B, \phi_j)$, $q_\phi(\lambda_B, \phi_j)$, $q(\lambda_B, \phi_{j+1})$, $q_\phi(\lambda_B, \phi_{j+1})$ を用いて, 緯度方向にエルミート３次補間を行えば, $q(\lambda_B, \phi_B)$が求められる． \end{enumerate} この方法も双３次補間と同じく周囲４グリッド上の計16個の情報を用いており, 補間の精度は同じになる．一方で, 方向分離することで, 計算式やコードが簡単になるため 高速化しやすくなる．さらに各補間ごとに任意のフィルタ処理を施すことが可能となる(後述)． \begin{figure}[h] \begin{center} \includegraphics[width=8.5cm]{sltt/sltt-images/Fig1.eps} \caption{方向分離の模式図．Tは補間したい座標を表す．} \end{center} \end{figure} なお, エルミート３次補間とは１次元の補間方法であり, 補間したい座標の周囲２点の 値と空間微分値の計４個の情報から補間多項式 \begin{equation*} f(x) = a_3x^3 + a_2x^2 +a_1x +a_0 \end{equation*} の係数$a_i$を求めて, 補間値を求める方法である． ラグランジュ補間が補間多項式の次数を上げるために, より広い範囲のグリッド値を用いるのに対して, エルミート補間は補間多項式の次数を上げるために, より高階の空間微分値を利用する． \subsubsection{変則エルミート５次補間} DCPAMでは混合比補間では補間多項式を５次の多項式として, 補間精度を上げる． 通常の１次元エルミート５次補間は補間したい座標の 周囲２点の$q$, $q_x$, $q_{xx}$と, 空間２階微分まで必要となる． そして, これを上述の方向分離の方法で２次元の補間を行うためには, $q$, $q_\lambda$, $q_{\phi}$, $q_{\lambda\lambda}$, $q_{\lambda\phi}$, $q_{\phi\phi}$, $q_{\lambda\lambda\phi}$, $q_{\lambda\phi\phi}$, $q_{\lambda\lambda\phi\phi}$と計算すべき微分値が大幅に増えてしまう． そこで, 次のような変則エルミート５次補間を用いる．つまり, 補間したい座標の周囲２点の$q$, $q_x$に加え, それぞれもうひとつ外側のグリッドの$q$の値を用いる． \begin{figure}[h] \begin{center} \includegraphics[width=10cm]{sltt/sltt-images/Fig2.eps} \caption{一般的なエルミート５次補間と変則エルミート５次補間．×印が補間したい座標を表す．} \end{center} \end{figure} %こうすることで, 補間のコンパクトさは減じるものの, 新たに微分値を計算する必要はなくなる． 以下に具体的な計算手順を, 格子点が等間隔の場合(経度)と不等間隔(緯度, 鉛直) の場合に分けて述べる． \paragraph{等間隔格子の場合} 補間したい座標のすぐ左側(座標が減少する側)の点の座標を原点として, 座標を取り直すと, 周囲４点の座標は$x=-\Delta x$, $0$, $\Delta x$, $2\Delta x$である． いま, 補間多項式$f(x)$を \begin{equation*} f(x) = a_5x^5 + a_4x^4 + a_3x^3+a_2x^2+a_1x +a_0 \end{equation*} とすると, その導関数$g(x)$は \begin{equation*} g(x) = 5a_5x^4 + 4a_4x^3 + 3a_3x^2+ 2a_2x+a_1 \end{equation*} となる．既知の値は$f(-\Delta x)$, $f(0)$, $f(\Delta x)$, $f(2\Delta x)$, $g(0)$, $g(\Delta x)$なので, 以下の連立方程式が成り立つ． \begin{align} f(-\Delta x) &= -a_5\Delta x^5 + a_4 \Delta x^4 - a_3 \Delta x^3+a_2\Delta x^2-a_1\Delta x +a_0 \label{eq:7}\\ f(0) &= a_0\label{eq:8}\\ f(\Delta x) &= a_5\Delta x^5 + a_4 \Delta x^4 + a_3 \Delta x^3+a_2\Delta x^2+a_1\Delta x +a_0\label{eq:9}\\ f(2\Delta x) &= 32a_5\Delta x^5 + 16a_4 \Delta x^4 + 8a_3 \Delta x^3+4a_2\Delta x^2+2a_1\Delta x +a_0\label{eq:10}\\ g(0) &= a_1\label{eq:11}\\ g(\Delta x) &= 5a_5\Delta x^4 + 4a_4\Delta x^3 + 3 a_3\Delta x^2 + 2a_2\Delta x+a_1\label{eq:12} \end{align} 式(\ref{eq:8})と(\ref{eq:11})より, $a_0$, $a_1$は自明なので, 上記の方程式系は実質的には４元連立方程式である． 演算(\ref{eq:9})$-$(\ref{eq:7})から \begin{equation} f(\Delta x) -f(-\Delta x) = 2a_5\Delta x^5 + 2a_3 \Delta x^3 + 2a_1x \label{eq:13} \end{equation} 演算(\ref{eq:10})$-4\times$(\ref{eq:9})から \begin{equation} f(2\Delta x) -4f(\Delta x) = 28a_5\Delta x^5 +12a_4\Delta x^4 + 4a_3 \Delta x^3 -2a_1\Delta x -3a_0\label{eq:14} \end{equation} 演算$\Delta x \times $(\ref{eq:12})$-2\times$(\ref{eq:9})から \begin{equation} \Delta x g(\Delta x) -2f(\Delta x) = 3a_5\Delta x^5 +2a_4\Delta x^4 + a_3 \Delta x^3 -a_1\Delta x -2a_0\label{eq:15} \end{equation} 演算(\ref{eq:14})$-6\times$(\ref{eq:15})から \begin{equation} f(2\Delta x) + 8f(\Delta x) - 6\Delta x g(\Delta x) = 10a_5\Delta x^5 - 2a_3 \Delta x^3 +4a_1\Delta x +9a_0\label{eq:16} \end{equation} 演算(\ref{eq:13})$+$(\ref{eq:16})から \begin{equation} f(2\Delta x) + 9f(\Delta x) - f(-\Delta x) -6\Delta x g(\Delta x) = 12a_5\Delta x^5 +6a_1\Delta x +9a_0\label{eq:17} \end{equation} これより \begin{equation} a_5 = \frac{1}{12\Delta x^5}\left[ f(2\Delta x) + 9f(\Delta x) - f(-\Delta x) -6\Delta x g(\Delta x) -6a_1\Delta x -9a_0\right]\label{eq:18} \end{equation} と, $a_5$が求まる．数値計算なので, $a_5$を計算した後は, これを既知の値として, 他の係数を求めるために使うことができる．すなわち$a_3$は式(\ref{eq:13})より, \begin{equation} a_3 = \frac{1}{2\Delta x^3}\left[ f(\Delta x) - f(-\Delta x) - 2a_5\Delta x^5 - 2a_1\Delta x \right] \end{equation} 同様にして, 式(\ref{eq:15})より \begin{equation} a_4 = \frac{1}{2\Delta x^4}\left[ \Delta x g(\Delta x) - 2f(\Delta x) - 3a_5\Delta x^5 - a_3\Delta x^3 +a_1\Delta x + 2a_0 \right] \end{equation} そして, 演算(\ref{eq:9})$+$(\ref{eq:7})で得られる式より \begin{equation} a_2 = \frac{1}{2\Delta x^2}\left[ f(\Delta x) + f(-\Delta x) - 2a_4\Delta x^4 - 2a_0\right] \end{equation} が得られる．以上の計算により係数$a_i$を求めたのち, 次式で補間値を求める． \begin{equation} f(\xi) = a_5\xi^5 + a_4\xi^4 + a_3\xi^3+a_2\xi^2+a_1\xi +a_0 \end{equation} ここで$x = \xi$が補間したい場所の座標である． \paragraph{不等間隔格子の場合} 補間したい座標のすぐ左側(座標が減少する側)の点の座標を原点として, 座標を取り直し, 周囲４点の座標を$x=\Delta x_{21}$, $0$, $\Delta x_{23}$, $\Delta x_{24}$と表す．ここで$\Delta x_{21}<0$であることに注意したい．不等間隔格子の 場合も原理的には, 等間隔格子のときと同様に４元連立方程式を解けばよいのだが, そのまま解こうとすると式が煩雑になる． そこで, $\Delta x_{23}$が大きさ1になるように, 座標変換を施す．変換後の 座標を$X$で表すと, $X=r$, $0$, $1$, $t$である．ただし, $r \equiv \Delta x_{21}/\Delta x_{23}$, $t \equiv \Delta x_{24}/\Delta x_{23}$とする． いま, 補間多項式$F(X)$を \begin{equation*} F(X) = a_5X^5 + a_4X^4 + a_3X^3+a_2X^2+a_1X +a_0 \end{equation*} とすると, その導関数$G(X)$は \begin{equation*} G(X) = 5a_5X^4 + 4a_4X^3 + 3a_3X^2+ 2a_2X+a_1 \end{equation*} となる．既知の値は$F(r) = f(\Delta x_{21})$, $F(0)= f(0)$, $F(1)=f(\Delta x_{23})$, $F(t)=f(\Delta x_{24})$および, $G(0) = g(0)\Delta x_{23}$, $G(1) = g(\Delta x_{23})\Delta x_{23}$である． 変換により座標が伸縮したために, 元々の空間微分値に$\Delta x_{23}$をかけている． (\ref{eq:8})と(\ref{eq:11})に対応する自明な式は次のように書かれる． \begin{align*} a_0 &= f(0)\\ a_1 &= g(0)\Delta x_{23} \end{align*} さらに \begin{align*} F_1 &\equiv F(r) - a_0 -a_1r \\ F_3 &\equiv F(1) - a_0 -a_1 \\ F_4 &\equiv F(t) - a_0 -a_1t \\ G_3 &\equiv G(1) - a_1 \end{align*} と置けば, 解くべき連立方程式は \begin{align} a_5 + a_4 + a_3 + a_2 &= F_3\\ a_5r^5 + a_4r^4 + a_3r^3 + a_2r^2 &= F_1\\ a_5t^5 + a_4t^4 + a_3t^3 + a_2t^2 &= F_4\\ 5a_5 + 4a_4 + 3a_3 + 2a_2 &= G_3 \end{align} となる．この連理方程式の数値解は \begin{align} a_5 &= \frac{F_1}{(r-1)^2r^2(r-t)} - \frac{G_3}{(t-1)^2t^2(r-t)} - \frac{[4 + 2rt -3(r+t)]F_4}{(r-1)^2(t-1)^2} + \frac{G_3}{ (r-1)(t-1)}\\ % a_4 &= -\frac{(t+2)F_3}{((r-1)^2r^2(r-t) ) } + \frac{(r+2)F_4}{(t-1)^2t^2(r-t)}\nonumber \\ & \quad \quad \quad \quad \quad \quad + \frac{[5 - 3(r^2 + rt + t^2) +2rt(r+t)]F_3}{(r-1)^2(t-1)^2} - \frac{(r+t+1)G_3}{(r-1)(t-1)}\\ % a_3 &= -2F_3 + 1 -3a_5 - 2a_4 \\ a_2 &= F_3 - a_5 - a_4 - a_3 \end{align} で計算できる．補間値は次式で求められる． \begin{equation} F(\Xi) = a_5\Xi^5 + a_4\Xi^4 + a_3\Xi^3+a_2\Xi^2+a_1\Xi +a_0 \end{equation} ここで$X = \Xi = \xi/\Delta x_{23}$が補間したい場所の座標である． ２次元補間をする場合は, 補間したい座標の周囲16グリッドでの混合比とその 空間微分値から, 方向分離の方法を用いて 図 \ref{fig:3}のように１次元補間を７回行うことで求められる(図\ref{fig:3})．この際, 緯度方向の補間を最後に行うことで, より計算コストの高い不等間隔格子での補間回数を最少限に抑えられる． \begin{figure}[h] \begin{center} \includegraphics[width=13cm]{sltt/sltt-images/Fig3.eps} \caption{変則エルミート５次補間の方向分離による２次元補間の模式図．Tは補間したい座標を表す．緯度方向は不等間隔格子になるので注意が必要．\label{fig:3}} \end{center} \end{figure} \subsection{鉛直} 鉛直方向の混合比補間にも変則エルミート５次補間を用いる． ただし鉛直の場合, スペクトル変換を利用した微分計算が出来ないので, ４次精度の中央差分法で微分値を計算する． この際, 鉛直グリッドは一般に不等間隔格子なので, 注意が必要である． \subsubsection{不等間隔格子の４次精度中央差分} 不等間隔格子の４次精度中央差分の計算には, 微分値を計算したい点と その上下２点, 計５点の情報は必要となる．これら５点の座標をそれぞれ $\sigma_{k-2}, \sigma_{k-1}, \sigma_{k}, \sigma_{k+1}, \sigma_{k+2}$ とする．さらに以下のように格子点間隔を定義する． \begin{align*} s_1 &\equiv \sigma_k - \sigma_{k-1}\\ t_1 &\equiv \sigma_{k+1} - \sigma_{k}\\ s_2 &\equiv \sigma_k - \sigma_{k-2}\\ t_2 &\equiv \sigma_{k+2} - \sigma_{k} \end{align*} 次の４次のオーダーまでのテイラー展開を考える．（簡単のため$f(\sigma_k) = f_k$とする．） \begin{align} f_{k+1} &= f_k + t_1f^\prime_k + \frac{t_1^2}{2}f^{\prime\prime}_k +\frac{t_1^3}{6}f^{\prime\prime\prime}_k + \frac{t_1^4}{24}f^{\prime\prime\prime\prime}_k \label{eq:32}\\ f_{k-1} &= f_k - s_1f^\prime_k + \frac{s_1^2}{2}f^{\prime\prime}_k -\frac{s_1^3}{6}f^{\prime\prime\prime}_k + \frac{s_1^4}{24}f^{\prime\prime\prime\prime}_k\label{eq:33}\\ f_{k+2} &= f_k + t_2f^\prime_k + \frac{t_2^2}{2}f^{\prime\prime}_k +\frac{t_2^3}{6}f^{\prime\prime\prime}_k + \frac{t_2^4}{24}f^{\prime\prime\prime\prime}_k\label{eq:34}\\ f_{k-2} &= f_k - s_2f^\prime_k + \frac{s_2^2}{2}f^{\prime\prime}_k -\frac{s_2^3}{6}f^{\prime\prime\prime}_k + \frac{s_2^4}{24}f^{\prime\prime\prime\prime}_k\label{eq:35} \end{align} 演算$s_1^2\times$(\ref{eq:32})$-t_1^2\times$(\ref{eq:33})より \begin{equation*} f_{k+1}s_1^2 - f_k(s_1^2 -t_1^2) - f_{k-1}t_1^2 = f^{\prime}_k s_1t_1(s_1+t_1) + \frac{s_1^2t_1^2(s_1+t_1)}{6}f^{\prime\prime\prime}_k +\frac{s_1^2t_1^2(t_1^2-s_1^2)}{24}f^{\prime\prime\prime\prime}_k \end{equation*} 変形して \begin{equation} f^{\prime}_k = F_{11} - \frac{s_1t_1}{6}f^{\prime\prime\prime}_k -\frac{s_1t_1(t_1-s_1)}{24}f^{\prime\prime\prime\prime}_k \label{eq:36} \end{equation} ただし, \begin{equation} F_{11} \equiv \frac{f_{k+1}s_1^2 - f_k(s_1^2 -t_1^2) - f_{k-1}t_1^2}{s_1t_1(s_1+t_1)} \end{equation} となる．同様にして, 演算$s_2^2\times$(\ref{eq:34})$-t_2^2\times$(\ref{eq:35})より \begin{equation} f^{\prime}_k = F_{22} - \frac{s_2t_2}{6}f^{\prime\prime\prime}_k -\frac{s_2t_2(t_2-s_2)}{24}f^{\prime\prime\prime\prime}_k \end{equation} ただし, \begin{equation} F_{22} \equiv \frac{f_{k+2}s_2^2 - f_k(s_2^2 -t_2^2) - f_{k-2}t_2^2}{s_2t_2(s_2+t_2)} \label{eq:39} \end{equation} が得られる．演算$s_2t_2\times$(\ref{eq:36})$-s_1t_1\times$(\ref{eq:39}) により \begin{equation} f_k^\prime(s_2t_2- s_1t_1) = F_{11}s_2t_2 - F_{22}s_1t_1 -\frac{s_1t_1s_2t_2(t_1-s_1-t_2+s_2)}{24}f^{\prime\prime\prime\prime}_k\label{eq:40} \end{equation} が得られる． また, 演算$s_2^2\times$(\ref{eq:32})$-t_1^2\times$(\ref{eq:35})より \begin{equation} f^{\prime}_k = F_{21} - \frac{s_2t_1}{6}f^{\prime\prime\prime}_k -\frac{s_2t_1(t_1-s_2)}{24}f^{\prime\prime\prime\prime}_k \end{equation} ただし, \begin{equation} F_{21} \equiv \frac{f_{k+1}s_2^2 - f_k(s_2^2 -t_1^2) - f_{k-2}t_1^2}{s_2t_1(s_2+t_1)} \label{eq:42} \end{equation} さらに, 演算$s_1^2\times$(\ref{eq:33})$-t_2^2\times$(\ref{eq:34})より \begin{equation} f^{\prime}_k = F_{12} - \frac{s_1t_2}{6}f^{\prime\prime\prime}_k -\frac{s_1t_2(t_2-s_1)}{24}f^{\prime\prime\prime\prime}_k \end{equation} ただし, \begin{equation} F_{12} \equiv \frac{f_{k+2}s_1^2 - f_k(s_1^2 -t_2^2) - f_{k-1}t_2^2}{s_1t_2(s_1+t_2)} \label{eq:44} \end{equation} を得る．演算$s_1t_2\times$(\ref{eq:42})$-s_2t_1\times$(\ref{eq:44}) により \begin{equation} f_k^\prime(s_1t_2- s_2t_1) = F_{21}s_1t_2 - F_{12}s_2t_1 -\frac{s_1t_1s_2t_2(t_1-s_2-t_2+s_1)}{24}f^{\prime\prime\prime\prime}_k\label{eq:45} \end{equation} が得られる．ここで, \begin{equation*} r_1 \equiv t_1-s_1-t_2+s_2, \quad r_2 \equiv t_1-s_2-t_2+s_1 \end{equation*} (ただし, $r_1, r_2 \neq 0$)とすれば, 演算$r_2\times$(\ref{eq:40})$-r_1\times$(\ref{eq:45})より \begin{equation} f_k^\prime\left[(s_2t_2- s_1t_1)r_2 - (s_1t_2- s_2t_1)r_1\right] = (F_{11}s_2t_2 - F_{22}s_1t_1)r_2 - (F_{21}s_1t_2 - F_{12}s_2t_1)r_1 \end{equation} よって \begin{equation} f_k^{\prime} =\frac{(F_{11}s_2t_2 - F_{22}s_1t_1)r_2 - (F_{21}s_1t_2 - F_{12}s_2t_1)r_1}{(s_2t_2- s_1t_1)r_2 - (s_1t_2- s_2t_1)r_1} \label{eq:47} \end{equation} が不等間隔格子の４次精度の中央差分式になる． なお, 鉛直格子が等間隔の場合$2s_1 = 2t_1 = s_2 = t_1$なので, $r_1 = r_2 = 0$となり, 式(\ref{eq:47})は不定形になる．この場合, 式(\ref{eq:40})は \begin{equation*} 3s_1^2f_k^\prime = 4F_{11}s_1^2 - F_{22}s_1^2 \end{equation*} になるので, \begin{equation*} f_k^\prime = \frac{4F_{11} - F_{22}}{3} = \frac{-f_{k+2} +8f_{k+1} -8f_{k-1}+f_{k-2} }{12s_1} \end{equation*} が４次精度の中央差分式になる． \subsection{上下端境界} 上下端付近の計算には上下端の向こうに鏡像関係にある仮想点をおいて計算している． すなわち$q(\sigma_0) = q(\sigma_1)$, $q(\sigma_{-1}) = q(\sigma_2)$, $q(\sigma_{\rm kmax+1}) = q(\sigma_{\rm kmax})$, $q(\sigma_{\rm kmax+2}) = q(\sigma_{\rm kmax-1})$となるように 配列を拡張して計算している． ただし, 上流点が最上層($q$は半整数グリッド上で定義されている)よりも上に なった場合には, 補間値としては最上層の値$q(\sigma_{\rm kmax})$を そのまま使用する． 同様に, 上流点が最下層よりも下になった場合には, 補間値として最下層の値$q(\sigma_{1})$をそのまま使用する． \section{単調フィルタ} 変則エルミート５次補間によって急峻な混合比分布でも精度よく補間できる一方で, 補間多項式の次数が高いために, 分布が平坦な箇所では 人工的な短波(振動)が生じやすくなる．また, 微分値の計算にスペクトル変換を利用しているために, ギブス現象に由来する人工的な短波も発生してしまう． このため, 混合比分布に人工的な極大・極小が生じて, 移流スキームに求められる「単調性」が損なわれてしまう． このような人工的な短波を除去するために, Sun et al.(1996)の単調フィルタを 一部修正して用いる． \subsection{Sun et al. (1996) フィルタ} １次元で考え, 補間した値を$q_T$, 補間した場所の両側４点の値を$q_1$, $q_2$, $q_3$, $q_4$とする．また, $q$分布の全領域の最大値／最小値を大域的最大値($q_{\rm GMax}$)／最小値($q_{\rm GMin}$)と呼び, $q_2$と$q_3$のうち大きい方を局所的最大値($q_{\rm LMax}$), 小さい方を局所的最小値($q_{\rm LMin}$)と呼ぶ． このとき, \begin{enumerate} \item $(q_2 - q_1)(q_4 -q_3) \geq 0 $の場合 \begin{itemize} \item もし $q_T < q_{\rm LMin}$ ならば, $q_T$を$q_{\rm LMin}$で置換える \item もし $q_T > q_{\rm LMax}$ ならば, $q_T$を$q_{\rm LMax}$で置換える \end{itemize} \item $(q_2 - q_1)(q_4 -q_3) < 0 $の場合 \begin{itemize} \item もし $q_T < q_{\rm LGin}$ ならば, $q_T$を$q_{\rm GMin}$で置換える \item もし $q_T > q_{\rm GMax}$ ならば, $q_T$を$q_{\rm GMax}$で置換える \end{itemize} \end{enumerate} とするのがSun et al. (1996)の単調フィルタである． このフィルタによって, 人工的な短波が除去できる． しかし, 上の2.の処理によって, 混合比分布の最大値が徐々に減少したり, 最小値が増加する副作用がある． そこでDCPAMでは, \underline{上の2.の処理は行わず, 1.のみ実行する．} こうすることで, 混合比分布の最大(小)値が減少(増加)せずにすむ． ところが, 今度は逆にオーバーシュート, アンダーシュートの効果で 最大値が増加したり, 最小値が減少しうる．特に, 最小値だアンダーシュートによって 負の値になることが大きな問題になる． なぜなら, 混合比は非負で定義される量であり, 負の混合比は物理的にあり得ないからである． \subsection{Arcsine変換フィルタ} そこで, 非負を保証するためのフィルタとして, Arcsine 変換を利用する(Kashimura et al, 2013)． $n+1$ステップ目の移流計算を行う前に, 混合比$q_{n}$を次式に従い変数変換する． \begin{equation} Q_{n} = \frac{1}{2}\arcsin\left( \frac{2q_{n}}{\alpha q^{\rm max}_{n} } - 1\right); \quad -\frac{\pi}{4}\leq Q_n \leq \frac{\pi}{4} \label{eq:48} \end{equation} ここで, $q_{n}^{\rm max}$は$q_{n}$の最大値を表し, $\alpha$は１以上の定数とする． $Q_{n}$の分布をもとに上流点で補間計算をして$Q_{n+1}$を求め, 以下の逆変換により$q_{n+1}$を求める． \begin{equation} q_{n+1} = \frac{\alpha q_{n}^{\rm max}}{2}\left[\sin(2Q_{n+1}) + 1 \right] \label{eq:49} \end{equation} 以上の手続きにより$q$の非負性が担保される． 補間計算の結果, $Q_{n+1}$が(\ref{eq:48})の不等式を満たさなくても, 逆変換(\ref{eq:49})により, $0\leq q_{n+1} \leq \alpha q_{n}^{\rm max}$になるからである(図\ref{fig:4})． 上述の変換を施すことで, 図\ref{fig:4}において傾きの変化率が大きい$q= 0$, $\alpha q^{\rm max}$付近では, 逆拡散の効果が働く(矢部他, 2003). $q^{\rm max}$付近の歪みは$\alpha=1.05$程度とすることで抑制できる（DCPAMでは$\alpha=1.05$としている）． $q=0$付近では通常, 補間操作に伴う数値拡散が強いため, 変換による逆拡散性は数値拡散を抑えるように働く． 上述のArcsine変換とSun et al.フィルタの一部を併用することで, 非負性を担保しながら, 人工的な短波を除去することが出来る． \begin{figure}[t] \begin{center} \includegraphics[width=6cm]{sltt/sltt-images/asin.eps} {\footnotesize \caption{Arcsine 変換(\ref{eq:4})．横軸が$q$で縦軸が$Q$．\label{fig:4}} } \end{center} \end{figure} \section{ネームリストによる制御} 本章で述べたセミラグランジュ法による物質移流スキームに関係するネームリスト変数を 以下に記す． \begin{center} \begin{tabular}{lcc} \verb|dynamics_hspl_vas83_nml| &&\\ 変数名& デフォルト値& 備考\\ \hline \vspace{3mm} FlagSLTT & .false. & セミラグスキーム使用フラグ\\ \verb|sltt_nml| & & \\ \hline FlagSLTTArcsine & .true. & Arcsine変換フィルタ使用フラグ\\ SLTTIntHor & ``HQ'' & 水平補間法$^*$\\ SLTTIntVer & ``HQ'' & 鉛直補間法$^*$ \end{tabular} \end{center} $^*$HQ: エルミート５次補間, HC: エルミート３次補間． \section*{参考文献} \begin{itemize} \item Enomoto, T.: Bicubic interpolation with spectral derivatives, SOLA, 4, pp.~5--8, 2008. \item Kashimura, H., T.~Enomoto, Y.~O.~Takahashi: Non-negative filter using arcsine transformation for tracer advection with semi-Lagrangian scheme, NCTAM, 62, 2013. \item Sun, W.-Y., K.-S.~Yeh, and R.-Y.~Sun: A simple semi-Lagrangian scheme for advection equations, QJRMS, 122, pp.~1211--1226, 1996. \item Williamson, D.~L., and Rasch, P.~J.: Two-dimensional semi-Lagrangian transport with shape-preserving interpolation, Mon.~Wea.~Rev., 117, pp.~102--129, 1989 \item 矢部孝, 内海隆行, 尾形陽一: CIP法---原子から宇宙スケールまでを解くマルチスケール解法, 森北出版, 2003. \end{itemize}