% % 表題 SPMODEL サンプルプログラムドキュメント % shallow_eqbeta1.tex % % 履歴 2003/05/21 山田 由貴子 % 2004/06/10 小高 正嗣 % 2004/08/27 山田 由貴子 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \documentclass[a4j,12pt]{jarticle} % % プリアンブル % \usepackage{Dennou6} % dennou-sty-6 を使用 \Dtitle[{\small 水路領域赤道 $\beta$ 面浅水モデル}] % ヘッダ部タイトル, モデル名と同じ { {\large SPMODEL サンプルプログラム} \\ % ここは変更しない 水路領域赤道 $\beta$ 面浅水モデル: \\ 局所熱源に対する線形応答問題 {\large shallow\_eqbeta1.f90} % プログラム名 } % \Dauthor[山田 由貴子, 小高 正嗣]{山田 由貴子, 小高 正嗣} % 著者 \Ddate{2004 年 8 月 27 日} % 日付 \Dnoparindent % 段落の字下げをしない % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % 本文開始 % % 基本的な章節構成は以下の通り. % % 1. はじめに % 2. 支配方程式系 % 3. 離散化 % 4. 利用モジュールとその他の設定 % 5. 数値実験 % 6. 参考文献 % -. 謝辞 % % 必要に応じて構成を変更すること. % \begin{document} \maketitle % タイトルの作成 \tableofcontents % 目次の作成 \pagebreak %---------------------------------------------------------------------- \Dparskip %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \newpage % 改ページ \section{概要} SPMODEL サンプルプログラム『shallow\_eqbeta1.f90』に用いられている基礎方 程式と境界条件, および, このプログラムを用いた数値実験の方法について解説 する. 基礎方程式は線形化された赤道 $\beta$ 面浅水方程式である. 境界条件 は, $x$ 方向は周期境界, $y$ 方向は壁面境界としている. 計算はスペクトル法 を用いて行い, 展開関数は水平方向にはフーリエ級数, 鉛直方法には境界条件に あわせてフーリエ正弦級数またはフーリエ余弦級数を用いる. スペクトル変換と 逆変換および微分演算には, SPMODEL ライブラリ(spml)を用いている. 数値実験 として, Gill(1980) および Heckley and Gill(1984) と同様の質量源分布を与 えた場合の応答を計算する. {\bf プログラム名} \\ % プログラム名 {\footnotesize shallow\_eqbeta1.f90 } {\bf プログラム取得元}\\ % プログラム取得先 {\footnotesize http:\slash \slash www.gfd-dennou.org\slash arch\slash spmodel\slash 2d-channel-esc\slash shallow-equator\slash sample1\slash SIGEN.htm } {\bf SPMODEL サンプルプログラム目次}\\ % サンプルプログラム目次(変更しない) {\footnotesize http:\slash \slash www.gfd-dennou.org\slash arch\slash spmodel\slash sample.htm } {\bf SPMODEL の使い方}\\ % SPMODEL の使い方(変更しない) {\footnotesize http:\slash \slash www.gfd-dennou.org\slash arch\slash spmodel\/ } %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \newpage % 改ページ \section{支配方程式系} ここでは, 支配方程式系と境界条件を記す. \subsection{支配方程式系} 支配方程式系は線形化された赤道 $\beta$ 面浅水方程式系である. \begin{eqnarray} \DP{u}{t} &=& \beta yv -g\DP{h}{x} -r_{M}u, \Deqlab{x方向運動方程式}\\ \DP{v}{t} &=& -\beta yu -g\DP{h}{y} -r_{M}v, \Deqlab{y方向運動方程式}\\ \DP{h}{t} &=& - H_{0}\left( \DP{u}{x} + \DP{v}{y}\right) -r_{H}h - Q. \Deqlab{質量保存の式} \end{eqnarray} 各記号は以下の量をあらわす. \begin{table}[h] \begin{center} \begin{tabular}{c l l } \hline 記号 &\qquad\qquad& 変数/物理定数 \\ \hline $x$ && 経度方向座標 \\ $y$ && 緯度方向座標 \\ $t$ && 時間 \\ $u$ && 経度方向流速 \\ $v$ && 緯度方向流速 \\ $H_{0}$ && 流体の平均厚さ \\ $h$ && 流体の厚さの偏差 \\ $\beta$ && 赤道 $\beta$ パラメタ \\ $g$ && 重力加速度 \\ $r_M$ && レイリー減衰係数 \\ $r_H$ && ニュートン冷却係数\\ $Q$ && 熱源(質量源) \\ \hline \end{tabular} \end{center} \caption{変数, 物理定数の定義} % \Dtablab{} \end{table} \subsection{無次元化} 支配方程式\Deqref{x方向運動方程式}, \Deqref{y方向運動方程式}, \Deqref{質量保存の式}を無次元化する. 長さスケールとして赤道変形半径 \begin{equation} L_{R} \equiv \left(\frac{\sqrt{gH_{0}}}{\beta}\right)^{\frac{1}{2}} \end{equation} を導入する. 速度スケールは $\sqrt{gH_{0}}$, 厚さスケールは $H_{0}$ とする. これらを用いると, 長さ, 時間, 速度, 厚さは以下のように無次元化される. \[ x_{*}=x/L_{R}, \quad y_{*}=y/L_{R}, \quad t_{*}=t/(L_{R}/\sqrt{gH_{0}}). \] \[ u_{*} = u/\sqrt{gH_{0}}, \quad v_{*} = v/\sqrt{gH_{0}}, \quad h_{*} = h/H_{0}. \] ここで添字 $*$ は無次元量を表す. 以上の無次元変数を用いて\Deqref{x方向運動方程式}, \Deqref{y方向運動方程式}, \Deqref{質量保存の式}を無次元化すると, 以下の式が得られる. \begin{eqnarray} \DP{u_{*}}{t_{*}} &=& y_{*}v_{*} -\DP{h_{*}}{x_{*}} -r_{M*}u_{*}, \Deqlab{無次元化されたx方向運動方程式} \\ \DP{v_{*}}{t_{*}} &=& -y_{*}u_{*} -\DP{h_{*}}{y_{*}} -r_{M*}v_{*}, \Deqlab{無次元化されたy方向運動方程式}\\ \DP{h_{*}}{t_{*}} &=& - \left( \DP{u_{*}}{x_{*}} + \DP{v_{*}}{y_{*}}\right) -r_{H*}h_{*} - Q_{*}. \Deqlab{無次元化された質量保存の式} \end{eqnarray} ここで, \[ r_{M*} = r_{M}/(L_{R}/\sqrt{gH_{0}}), \quad r_{H*} = r_{H}/(L_{R}/\sqrt{gH_{0}}), \quad Q_{*} = Q/(H_{0}\sqrt{gH_{0}}/L_{R}) \] である. \subsection{境界条件} $x$ 方向の境界条件は周期境界条件とする. $y$ 方向の境界条件は $y=\pm y_{m}/2$ に置いた剛体壁面で $v=0$, 応力なしとする. $x_{*}$方向の境界条件は, $x_{*}$ 方向の計算領域を $x_{m*}$ とすると \begin{eqnarray} u_{*}(x_{*}+x_{m*},y_{*}) &=& u_{*}(x_{*},y_{*}), \\ v_{*}(x_{*}+x_{m*},y_{*}) &=& v_{*}(x_{*},y_{*}), \\ h_{*}(x_{*}+x_{m*},y_{*}) &=& h_{*}(x_{*},y_{*}) \end{eqnarray} と表される. $y_{*}$ 方向の境界条件は $y_{*}=\pm y_{m_{*}}/2$ において, \begin{eqnarray} && v_{*}=0, \\ && \DP{u_{*}}{y_{*}} = 0, \\ && \DP{h_{*}}{y_{*}} = 0 \end{eqnarray} である\footnote{ただしこの境界条件はふさわしくない. 本来ならば $x, \, y$ 方向ともに周期境界とするのが正しい. しかし数値計算に SPMODEL ライブラリの esc\_module を用いているため, このような境界条件に決まる. }. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \newpage % 改ページ \section{離散化} この節では方程式系の空間離散化および使用した時間積分法について説明, プログラム内で実際に用いられている方程式を記述する. 以下では無次元量を示す添字 $*$ は省略する. \subsection{水平離散化} 支配方程式\Deqref{無次元化されたx方向運動方程式}, \Deqref{無次元化された y方向運動方程式}, \Deqref{無次元化された質量保存の式}の空間離散表現は以 下である. \begin{eqnarray} \DP{u_{i,j}}{t} &=& y_j v_{i,j} - \left(\DP{h}{x}\right)_{i,j} - r_{M}u_{i,j}, \Deqlab{離散化されたx方向運動方程式} \\ \DP{v_{i,j}}{t} &=& - y_j u_{i,j} - \left(\DP{h}{y}\right)_{i,j} - r_{M}v_{i,j}, \Deqlab{離散化されたy方向運動方程式} \\ \DP{h_{i,j}}{t} &=& - \left(\DP{u}{x}\right)_{i,j} - \left(\DP{v}{y}\right)_{i,j} -r_{H}h_{i,j} - Q_{i,j}. \Deqlab{離散化された質量保存の式} \end{eqnarray} ここで添字 $i,j$ は格子点 $(x_{i}, y_{j})$ 上の値であることを表す. \subsection{水平方向のスペクトル表現} 空間離散化した支配方程式\Deqref{離散化されたx方向運動方程式}, \Deqref{離 散化されたy方向運動方程式}, \Deqref{離散化された質量保存の式}をスペクト ル法を用いて表現する. スペクトル展開は東西方向にフーリエ級数, 南北方向に は境界条件にあわせてフーリエ正弦級数またはフーリエ余弦級数を用いて行う. 非線形項を扱う場合は, 先に格子点上での非線形項の値を計算し, その値のスペ クトルを求める方法(変換法) を用いる. 以下では $k, l$ をそれぞれ$x, y$ 方 向波数, $K, L$ を切断波数, $I, J$ を格子点数とする. $u_{i,j}, v_{i,j}, H_{i,j}$ はスペクトル逆変換によって以下のように展開される. \begin{eqnarray} u_{i,j} &=& \sum^{K}_{k=-K} \sum^{L}_{l=0} \exp\left(\frac{2\pi ikx_i}{x_m}\right) \cos\left(\frac{\pi ly_j}{y_m} \right) \hat{u}_{k,l}, \\ v_{i,j} &=& \sum^{K}_{k=-K} \sum^{L}_{l=1} \exp\left(\frac{2\pi ikx_i}{x_m}\right) \sin\left(\frac{\pi ly_j}{y_m} \right) \hat{v}_{k,l}, \\ h_{i,j} &=& \sum^{K}_{k=-K} \sum^{L}_{l=0} \exp\left(\frac{2\pi ikx_i}{x_m}\right) \cos\left(\frac{\pi ly_j}{y_m} \right) \hat{h}_{k,l}. \end{eqnarray} スペクトル係数 $\hat u _{i,j}, \hat v _{i,j}, \hat h _{i,j}$ は スペクトル変換によって以下のように与えられる. \begin{eqnarray} \hat{u}_{k,l} &=& \frac{1}{x_m y_m} \sum^{I-1}_{i=0}\sum^{J-1}_{j=0} \exp\left(-\frac{2\pi ikx_i}{x_m}\right) \cos\left(\frac{\pi ly_j}{y_m} \right) u_{i,j}, \\ \hat{v}_{k,l} &=& \frac{1}{x_m y_m} \sum^{I-1}_{i=0} \sum^{J-1}_{j=0} \exp\left(-\frac{2\pi ikx_i}{x_m}\right) \sin\left(\frac{\pi ly_j}{y_m} \right) v_{i,j}, \\ \hat{h}_{k,l} &=& \frac{1}{x_m y_m} \sum^{I-1}_{i=0} \sum^{J-1}_{j=0} \exp\left(-\frac{2\pi ikx_i}{x_m}\right) \cos\left(\frac{\pi ly_j}{y_m} \right) h_{i,j}. \end{eqnarray} 以上を用いると, \Deqref{離散化されたx方向運動方程式}, \Deqref{離散化され たy方向運動方程式}, \Deqref{離散化された質量保存の式}のスペクトル表現は 以下のようになる. \begin{eqnarray} \DP{\hat{u}_{k,l}}{t} &=& \frac{1}{x_m y_m} \sum^{I-1}_{i=0} \sum^{J-1}_{j=0} \exp\left(-\frac{2\pi ikx_i}{x_m}\right) \cos\left(\frac{\pi ly_j}{y_m} \right) y_j v_{i,j} \nonumber\\ & & - \left(\frac{2\pi ik}{x_m}\right) \hat{h}_{k,l} -r_M \hat{u}_{k,l}, \\ %----------------------------------------------------------- \DP{\hat{v}_{k,l}}{t} &=& - \frac{1}{x_m y_m} \sum^{I-1}_{i=0} \sum^{J-1}_{j=0} \exp\left(-\frac{2\pi ikx_i}{x_m}\right) \sin\left(\frac{\pi ly_j}{y_m} \right) y_j u_{i,j} \nonumber\\ && + \left(\frac{\pi l}{y_m}\right)\hat{h}_{k,l} -r_M \hat{v}_{k,l} \\ %----------------------------------------------------------- \DP{\hat{h}_{k,l}}{t} &=&- \left(\frac{2\pi ik}{x_m}\right) \hat{u}_{k,l} - \left(\frac{\pi l}{y_m}\right) \hat{v}_{k,l} -r_H \hat{h}_{k,l} \nonumber\\ && - \frac{1}{x_m y_m} \sum^{I-1}_{i=0} \sum^{J-1}_{j=0} \exp\left(-\frac{2\pi ikx_i}{x_m}\right) \cos\left(\frac{\pi ly_j}{y_m} \right) Q_{i,j}. \end{eqnarray} \subsection{時間積分} ここでは時間積分法について記述し, プログラム内で実際に用いられている方 程式を記述する. 以下では $\Delta t$ を時間格子間隔, 時刻 $\tau \Delta t$ における $\hat u_{k,l}$ の値を $\hat u_{k,l}^{\tau}$ 等と表す. 時間方向の離散化は Euler スキームを用いて行う\footnote{ プログラムは Euler スキームではなく, 半端な time split のようになってし まっており, 本ドキュメントの定式化と異なっている. }. 時空間方向に離散化され た方程式は以下のように表される. \begin{eqnarray} \hat u_{k,l}^{\tau +1} &=& \hat u _{k,l}^{\tau} + \Delta t A_{u}^{\tau}, \\ \hat v_{k,l}^{\tau +1} &=& \hat v _{k,l}^{\tau} + \Delta t A_{v}^{\tau}, \\ \hat h_{k,l}^{\tau +1} &=& \hat h _{k,l}^{\tau} + \Delta t A_{h}^{\tau} \end{eqnarray} \begin{eqnarray} A_{u}^{\tau} &=& \frac{1}{x_m y_m} \sum^{I-1}_{i=0} \sum^{J-1}_{j=0} \exp\left(-\frac{2\pi ikx_i}{x_m}\right) \cos\left(\frac{\pi ly_j}{y_m} \right) y_j \hat v_{i,j}^{\tau} \nonumber \\ & & - \left(\frac{2\pi ik}{x_m}\right) \hat{h}_{k,l}^{\tau} -r_M \hat{u}_{k,l}^{\tau}, \\ %----------------------------------------------------------- A_{v}^{\tau} &=& - \frac{1}{x_m y_m} \sum^{I-1}_{i=0} \sum^{J-1}_{j=0} \exp\left(-\frac{2\pi ikx_i}{x_m}\right) \sin\left(\frac{\pi ly_j}{y_m} \right) y_j \hat u_{i,j}^{\tau} \nonumber\\ && + \left(\frac{\pi l}{y_m}\right)\hat{h}_{k,l}^{\tau} -r_M \hat{v}_{k,l}^{\tau}, \\ %----------------------------------------------------------- A_{h}^{\tau} &=&- \left(\frac{2\pi ik}{x_m}\right) \hat{u}_{k,l}^{\tau} - \left(\frac{\pi l}{y_m}\right) \hat{v}_{k,l}^{\tau} -r_H \hat{h}_{k,l}^{\tau} \nonumber\\ && - \frac{1}{x_m y_m} \sum^{I-1}_{i=0} \sum^{J-1}_{j=0} \exp\left(-\frac{2\pi ikx_i}{x_m}\right) \cos\left(\frac{\pi ly_j}{y_m} \right) Q_{i,j}^{\tau}. \end{eqnarray} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \newpage % 改ページ \section{使用モジュールとその他の設定} スペクトル変換と逆変換, 微分演算は SPMODEL ライブラリ (spml) の esc\_module に含まれる関数を用いて行う. フーリエ正弦および余弦変換, それ らの逆変換の際の数値積分は台形公式を用いて行う. spml が下位で使用する ISPACK の仕様から, 格子点数 $I,J$ は偶数で, かつ $I/2, J/2=2^{a}3^{b}5^{c}$ ($a, b, c$ は 0 または整数) でなければならない. 非 線形項の計算によって生じるエリアジングを防ぐため, 格子点数 $I,J$ と切断 波数 $K,L$ は $I>3K,J>3K/2$ を満たすように与える. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \newpage % 改ページ \section{数値実験} 数値実験として, Gill (1980) および Heckley and Gill (1984) と同様の定常 な質量源分布を与えた場合の応答を計算する. 質量源分布は以下のように与える. \begin{equation} Q = \left\{ \begin{array}{lc} Q_{0}\exp (y^{2}/b) \cos \left(\frac{\pi x}{2a}\right) & - a \le x \le a \\ 0 & x < -a, x > a \end{array} \right. \end{equation} 初期条件は $u=v=h=0$ である. パラメータは\Dtabref{使用したパラメータの値} にまとめた値を用いる. 格子点数 $I, J$ と切断波数 $K, L$ はそれぞれ $I=64, J=32, K=L=16$ とする. 時間格子間隔 $\Delta t$ は $0.2$, 計算ステップ数は 1,000 ス テップである. \begin{table}[h] \begin{center} \begin{tabular}{c l l } \hline パラメータ &\qquad\qquad& 数値 \\ \hline $r_M$ && 0.0 \\ $r_H$ && 0.0 \\ $Q_{0}$ && 1.0 \\ $a$ && 20 \\ $b$ && 0.2 \\ $x_{m}$ && 400 \\ $y_{m}$ && 10 \\ \hline \end{tabular} \end{center} \caption{使用したパラメータの値} \Dtablab{使用したパラメータの値} \end{table} \Dfigref{shallow-eqbeta1.f90 の結果}に厚さ分布の時間変化の様子を示す. 東向きにケルビン応答, 西向きにロスビー応答が伝播していることが分かる. ロスビー応答の伝播距離はケルビン応答の約 1/3 である. \begin{figure}[p] \begin{center} \Depsf[15cm][]{./figs/shallow_eqbeta1.ps} \vspace*{-1.0cm} \caption{厚さの分布の時間変化の様子. 上から順に $t=20, 40, 60, 80$ の結果.} \Dfiglab{shallow-eqbeta1.f90 の結果} \end{center} \end{figure} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \newpage % 改ページ \newpage \section{参考文献} \begin{description} \item Gill, A. E., 1980: Some simple solution for heat-induced tropical circulation. {\it Qurt. J. R. Meteorol. Soc.}, {\bf 106}, 447--462. \item Heckley, W. A., and A. E. Gill, 1984: Some simple analytical solutions to the ploblems of forced equatorial long waves. {\it Quart. J. Roy. Met. Soc.}, {\bf 110}, 203--217. \item 竹広真一, 石岡圭一, 豊田英司, 石渡正樹, 林祥介, SPMODEL 開発グループ, 2002: 階層的地球流体力学スペクトルモデル集 (SPMODEL), http:\slash \slash www.gfd-dennou.org\slash arch\slash spmodel\slash , 地球流体電脳倶楽部. \end{description} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \newpage % 改ページ \section*{謝辞} % 謝辞. このまま利用する. \markright{謝辞} % ヘッダに表示 本資源は, 地球流体電脳倶楽部のインターネット上での学術知識の集積と活用 の実験の一環として \begin{center} http:\slash \slash www.gfd-dennou.org\slash arch\slash spmodel\slash \end{center} において公開されているものである (\copyright 地球流体電脳倶楽部スペクトルモデルプロジェクト spmodel@gfd-dennou.org 2002. ). 本資源は, 著作者の諸権利に抵触しない(迷惑をかけない)限りにおいて自由に 利用していただいて構わない. なお, 利用する際には今一度自ら内容を確かめ ることをお願いする(無保証無責任原則). 本資源に含まれる元資源提供者(図等の版元等を含む)からは, 直接的な形で の WEB 上での著作権または使用許諾を得ていない場合があるが, 勝手ながら, 「未来の教育」のための実験という学術目的であることをご理解いただけるも のと信じ, 学術標準の引用手順を守ることで諸手続きを略させていただいて いる. 本資源の利用者には, この点を理解の上, 注意して扱っていただける ようお願いする. 万一, 不都合のある場合には \begin{center} spmodel@gfd-dennou.org \end{center} まで連絡していただければ幸いである. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \end{document}