%
% 表題   SPMODEL サンプルプログラムドキュメント
%        thermal_responce.tex
%
% 履歴   2003/05/21 山田 由貴子 
%        2004/06/10 小高 正嗣  
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\documentclass[a4j,12pt]{jarticle}
%
% プリアンブル
%
\usepackage{Dennou6}                      % dennou-sty-6 を使用


\Dtitle[{\small 水路領域赤道 $\beta$ 面浅水モデル}]  % ヘッダ部タイトル, モデル名と同じ
{ {\large SPMODEL サンプルプログラム} \\  % ここは変更しない 
   水路領域赤道 $\beta$ 面浅水モデル: \\
   局所熱源に対する非線形応答問題
  {\large thermal\_responce.f90}                   % プログラム名
}
%
\Dauthor[山田 由貴子, 小高 正嗣]{山田 由貴子, 小高 正嗣}            % 著者
\Ddate{2004 年 8 月 25 日}                % 日付
\Dnoparindent                             % 段落の字下げをしない
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 本文開始
%
% 基本的な章節構成は以下の通り. 
%
%   1. はじめに 
%   2. 支配方程式系
%   3. 離散化
%   4. 利用モジュールとその他の設定
%   5. 数値実験
%   6. 参考文献
%   -. 謝辞
% 
% 必要に応じて構成を変更すること.
%
\begin{document}

\maketitle                                % タイトルの作成
\tableofcontents                          % 目次の作成
\pagebreak                                

%----------------------------------------------------------------------
\Dparskip
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage                                  % 改ページ
\section{概要}

SPMODEL サンプルプログラム『thermal\_responce.f90』に用いられている基礎
方程式と境界条件, および, このプログラムを用いた数値実験の方法について解
説する. 基礎方程式は赤道 $\beta$ 面非線型浅水方程式である. 境界条件は,
$x$ 方向は周期境界, $y$ 方向は壁面境界としている. 計算はスペクトル法を用
いて行い, 展開関数は水平方向にはフーリエ級数, 鉛直方法には境界条件にあわ
せてフーリエ正弦級数またはフーリエ余弦級数を用いる. スペクトル変換と逆変
換および微分演算には, SPMODEL ライブラリ(spml)を用いている. 数値実験とし
て, Gill(1980) および Heckley and Gill(1984) と同様の質量源分布を与えた
場合の応答を計算する.


{\bf プログラム名} \\                  % プログラム名
{\footnotesize
thermal\_responce.f90
}

{\bf プログラム取得元}\\               % プログラム取得先
{\footnotesize
http:\slash \slash www.gfd-dennou.org\slash arch\slash spmodel\slash 2d-channel-esc\slash shallow-equator\slash thermal-responce\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} &=& 
  -u\DP{u}{x} -v\DP{u}{y}
  + \beta yv 
  -g\DP{h}{x} 
  -r_{M}u 
  + \nu_{M}\left(\DP[2]{}{x} +\DP[2]{}{y} \right) u 
  + F_x, \Deqlab{x方向運動方程式}  \\
 \DP{v}{t} &=& 
  -u\DP{v}{x} -v\DP{v}{y}
  -\beta yu 
  -g\DP{h}{y} 
  -r_{M}v 
  + \nu_{M}\left(\DP[2]{}{x} +\DP[2]{}{y} \right) v 
  + F_y, \Deqlab{y方向運動方程式}\\
 \DP{h}{t} &=& 
  -\DP{}{x}\left[ (h+H_{0})u \right]
  -\DP{}{y}\left[ (h+H_{0})v \right]
  -r_{H}h 
  + \nu_{H}\left(\DP[2]{}{x} +\DP[2]{}{y} \right) 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$   && ニュートン冷却係数\\
   $\nu_M$ && 粘性拡散係数    \\ 
   $\nu_H$ && 熱拡散係数      \\ 
   $F_x$   && 経度方向外力    \\ 
   $F_y$   && 緯度方向外力    \\ 
   $Q$     && 熱源(質量源)    \\  \hline
  \end{tabular} 
 \end{center}
 \caption{変数, 物理定数の定義}
% \Dtablab{}
\end{table}

\subsection{境界条件}

$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$ 方向の境界条件は剛体壁面で $v=0$, 応力なしとする.
すなわち $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} &=& 
                     -u_{i,j}\left(\DP{u}{x}\right)_{i,j}
                     -v_{i,j}\left(\DP{u}{y}\right)_{i,j}
                     +\beta y_j v_{i,j}\nonumber\\ 
 &&                  - g\left(\DP{h}{x}\right)_{i,j}
		     - r_{M}u_{i,j}
                     + \nu _{M}\left[
				\left(\DP[2]{u}{x}\right)_{i,j}
				+ \left(\DP[2]{u}{y}\right)_{i,j}
			     \right]
                     + F_{x_{i,j}},   
                    \Deqlab{離散化されたx方向運動方程式} \\
 \DP{v_{i,j}}{t} &=& 
                     -u_{i,j}\left(\DP{v}{x}\right)_{i,j}
                     -v_{i,j}\left(\DP{v}{y}\right)_{i,j}
                     - \beta y_j u_{i,j}\nonumber\\ 
 &&                  - g\left(\DP{h}{y}\right)_{i,j}
		     - r_{M}v_{i,j}
                     + \nu _{M}\left[
				\left(\DP[2]{v}{x}\right)_{i,j}
				+ \left(\DP[2]{v}{y}\right)_{i,j}
			     \right]
                     + F_{y_{i,j}}, 
                    \Deqlab{離散化されたy方向運動方程式} \\
 \DP{h_{i,j}}{t} &=& 
  -\left\{ \DP{}{x}\left[(h_{i,j}+H_{0})u_{i,j}\right] \right\}_{i,j}
  -\left\{ \DP{}{y}\left[(h_{i,j}+H_{0})v_{i,j}\right] \right\}_{i,j} 
  \nonumber\\
&&                   -r_{H}h_{i,j} 
                     + \nu _{H}\left[
				\left(\DP[2]{h}{x}\right)_{i,j}
				+ \left(\DP[2]{h}{y}\right)_{i,j}
			     \right]
                     - 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)
  u_{i,j}
  \left(
  \DP{u}{x}
 \right)_{i,j}
  \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)
  v_{i,j}
  \left(
   \DP{u}{y}
  \right)_{i,j}
  \nonumber\\
 & & 
  + \frac{\beta}{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\\
 & & 
  - g\left(\frac{2\pi ik}{x_m}\right) \hat{h}_{k,l} 
  -r_M \hat{u}_{k,l}
  - \nu_M\left[
	  \left(\frac{2\pi ik}{x_m}\right)^2\hat{u}_{k,l}
	  - \left(\frac{\pi l}{y_m}\right)^2\hat{u}_{k,l}
	\right] \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)
  F_{x_{i,j}},  \\
%-----------------------------------------------------------
 \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)
  u_{i,j}
  \left(
   \DP{v}{x}
  \right)_{i,j}
  \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)
  \sin\left(\frac{\pi ly_j}{y_m} \right)
  v_{i,j}
  \left(
   \DP{v}{y}
  \right)_{i,j}
  \nonumber\\
 && 
 - \frac{\beta}{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\\
 && 
  + g\left(\frac{\pi l}{y_m}\right)\hat{h}_{k,l} 
  -r_M \hat{v}_{k,l}
  - \nu_M\left[
	  \left(\frac{2\pi ik}{x_m}\right)^2\hat{v}_{k,l}
	  - \left(\frac{\pi l}{y_m}\right)^2\hat{v}_{k,l}
	\right]
  \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)
  \sin\left(\frac{\pi ly_j}{y_m} \right)
  F_{y_{i,j}}, 
  \\
%-----------------------------------------------------------
 \DP{\hat{h}_{k,l}}{t} 
  &=&- 
  \left(\frac{2\pi ik}{x_m}\right)
   \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}+H_{0}) u_{i,j}
  \nonumber\\
  && -  
\left(\frac{\pi l}{y_m}\right)
   \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)
  (h_{i,j}+H_{0}) v_{i,j}
  \nonumber\\
 && -r_H \hat{h}_{k,l}
  - \nu_H\left[
	  \left(\frac{2\pi ik}{x_m}\right)^2\hat{h}_{k,l}
	  - \left(\frac{\pi l}{y_m}\right)^2\hat{h}_{k,l}
	\right]
  \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}

ここで, 
\begin{eqnarray*}
 \left(\DP{u}{x}\right)_{i,j}
  &=&  
  \sum^{K}_{k=-K}
  \sum^{L}_{l=0}
  \left(\frac{2\pi ik}{x_m}\right)
  \exp\left(\frac{2\pi ikx_i}{x_m}\right)
  \cos\left(\frac{\pi ly_j}{y_m} \right)
  \hat{u}_{k,l}, 
  \\
 \left(\DP{u}{y}\right)_{i,j} 
  &=&   
  - \sum^{K}_{k=-K}
  \sum^{L}_{l=1}
  \left(\frac{\pi l}{y_m}\right)
  \exp\left(\frac{2\pi ikx_i}{x_m}\right)
  \sin\left(\frac{\pi ly_j}{y_m} \right)
  \hat{u}_{k,l}, 
  \\
 \left(\DP{v}{x}\right)_{i,j} 
  &=&  
  \sum^{K}_{k=-K}
  \sum^{L}_{l=1}
  \left(\frac{2\pi ik}{x_m}\right)
  \exp\left(\frac{2\pi ikx_i}{x_m}\right)
  \sin\left(\frac{\pi ly_j}{y_m} \right)
  \hat{v}_{k,l}, 
  \\
 \left(\DP{v}{y}\right)_{i,j} 
  &=&  
  \sum^{K}_{k=-K}
  \sum^{L}_{l=0}
  \left(\frac{\pi l}{y_m}\right)
  \exp\left(\frac{2\pi ikx_i}{x_m}\right)
  \cos\left(\frac{\pi ly_j}{y_m} \right)
  \hat{v}_{k,l}
\end{eqnarray*}
である. 

\subsection{時間積分}
ここでは時間積分法について記述し, プログラム内で実際に用いられている方
程式を記述する. 以下では $\Delta t$ を時間格子間隔, 時刻 $\tau \Delta t$ 
における $\hat u_{k,l}$ の値を $\hat u_{k,l}^{\tau}$ 等と表す.

時間積分は陽解法を用いて行う. 
拡散項と減衰項は Euler スキーム, 
その他の項に対しては Leapfrog スキームを適用する. 
\begin{eqnarray}
 \hat{u}_{k,l}^{\tau+1} 
  &=& 
 \hat{u}_{k,l}^{\tau-1} \nonumber\\
&&
 + 2\Delta t\left\{ 
  - \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}^{\tau}
  \left(
  \DP{u}{x}
 \right)_{i,j}^{\tau}
\right . 
  \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)
  v_{i,j}^{\tau}
  \left(
   \DP{u}{y}
  \right)_{i,j}^{\tau}
  \nonumber\\
 & & 
  + \frac{\beta}{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}^{\tau} 
  \nonumber\\
 & & 
  - g\left(\frac{2\pi ik}{x_m}\right) \hat{h}_{k,l}^{\tau} 
  -r_M \hat{u}_{k,l}^{\tau-1}
  - \nu_M\left[
	  \left(\frac{2\pi ik}{x_m}\right)^2\hat{u}_{k,l}^{\tau-1}
	  - \left(\frac{\pi l}{y_m}\right)^2\hat{u}_{k,l}^{\tau-1}
	\right] \nonumber \\
 & & 
\left .
  + \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)
  F_{x_{i,j}}^{\tau}
  \right\},
  \\
%-----------------------------------------------------------
 \hat{v}_{k,l}^{\tau+1}
  &=&
 \hat{v}_{k,l}^{\tau-1} \nonumber\\
&&
 + 2\Delta t\left\{ 
  - \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)
  u_{i,j}^{\tau}
  \left(
   \DP{v}{x}
  \right)_{i,j}^{\tau}
  \right .
  \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)
  \sin\left(\frac{\pi ly_j}{y_m} \right)
  v_{i,j}^{\tau}
  \left(
   \DP{v}{y}
  \right)_{i,j}^{\tau}
  \nonumber\\
 && 
 - \frac{\beta}{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}^{\tau}
  \nonumber\\
 && 
  + g\left(\frac{\pi l}{y_m}\right)\hat{h}_{k,l}^{\tau} 
  -r_M \hat{v}_{k,l}^{\tau-1}
  - \nu_M\left[
	  \left(\frac{2\pi ik}{x_m}\right)^2\hat{v}_{k,l}^{\tau-1}
	  - \left(\frac{\pi l}{y_m}\right)^2\hat{v}_{k,l}^{\tau-1}
	\right]
  \nonumber\\
 && 
  \left .
    + \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)
  F_{y_{i,j}}^{\tau}
	  \right\},
  \\
%-----------------------------------------------------------
 \hat{h}_{k,l}^{\tau+1}
  &=&
 \hat{h}_{k,l}^{\tau-1} \nonumber\\
&&
 + 2\Delta t\left\{
  - \left(\frac{2\pi ik}{x_m}\right)
   \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}^{\tau}+H_{0}) u_{i,j}^{\tau}
 \right .
  \nonumber\\
  && -  
\left(\frac{\pi l}{y_m}\right)
   \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)
  (h_{i,j}^{\tau}+H_{0}) v_{i,j}^{\tau}
  \nonumber\\
 && -r_H \hat{h}_{k,l}^{\tau-1}
  - \nu_H\left[
	  \left(\frac{2\pi ik}{x_m}\right)^2\hat{h}_{k,l}^{\tau-1}
	  - \left(\frac{\pi l}{y_m}\right)^2\hat{h}_{k,l}^{\tau-1}
	\right]
  \nonumber\\
 && 
\left .
- \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} \right\}. 
\end{eqnarray}


\subsection{時間フィルター}

計算モードの増幅を抑制するため Asselin (1972) の時間フィルターを 1 ステッ
プ毎に適用する. 
\begin{eqnarray}
 \bar{\hat{u}}_{k,l}^{\tau} &=& 
  (1-2\epsilon_f)\hat{u}_{k,l}^{\tau}
  + \epsilon_f (\bar{\hat{u}}_{k,l}^{\tau-1} + \hat{u}_{k,l}^{\tau+1}),
  \\
 \bar{\hat{v}}_{k,l}^{\tau} &=& 
  (1-2\epsilon_f)\hat{v}_{k,l}^{\tau} 
  + \epsilon_f (\bar{\hat{v}}_{k,l}^{\tau-1} + \hat{v}_{k,l}^{\tau+1}),
  \\
 \bar{\hat{h}}_{k,l}^{\tau} &=& 
  (1-2\epsilon_f)\hat{h}_{k,l}^{\tau} 
  + \epsilon_f (\bar{\hat{h}}_{k,l}^{\tau-1} + \hat{h}_{k,l}^{\tau+1}). 
\end{eqnarray}

ここで $\epsilon_f$ はフィルターの係数である. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\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.02$, 計算ステップ数は 1,000 ステップとする.
Asselin (1972) の時間フィルターの係数 $\epsilon_f$ には, 標準的な値
$0.05$ を用いる.  


\begin{table}[h]
 \begin{center}
  \begin{tabular}{c l l  } \hline
   パラメータ    &\qquad\qquad&  数値   \\ \hline 
   $g$   && 1.0  \\
   $H_{0}$   && 1.0  \\
   $\beta$   && 1.0  \\
   $r_M$   && 0.1  \\
   $r_H$   && 0.1  \\
   $\nu_M$ && 0.02  \\ 
   $\nu_H$ && 0.02  \\ 
   $F_x$   && 0.0  \\ 
   $F_y$   && 0.0  \\ 
   $Q_{0}$ && 0.01 \\  
   $a$     && 2.0  \\
   $b$     && 2.0 \\
   $x_{m}$ && 18.0 \\
   $y_{m}$ && 9.0  \\ \hline
  \end{tabular} 
 \end{center}
 \caption{使用したパラメータの値}
 \Dtablab{使用したパラメータの値}
\end{table}

\Dfigref{thermal_responce.f90 の結果}に厚さ分布の時間変化の様子を示す.
東向きにケルビン応答, 西向きにロスビー応答が伝播していることが分かる.
ロスビー応答の伝播距離はケルビン応答の約 1/3 である.

\begin{figure}[p]
\begin{center}
\Depsf[15cm][]{./figs/thermal_responce.eps}
\vspace*{-1.0cm}
\caption{厚さの分布の時間変化の様子. 上から順に $t=2, 4, 6, 8$ の結果.}
\Dfiglab{thermal_responce.f90 の結果}
\end{center}
\end{figure}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage                                  % 改ページ
\newpage
\section{参考文献}
\begin{description}
 \item Asselin, R. A., 1972: 
       Frequency filter for time integrations.
       {\it Mon. Wea. Rev.}, {\bf 100}, 487--490.

 \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}

