%
% 表題   SPMODEL サンプルプログラムドキュメント
%        spshallow_zd.tex
%
% 履歴   2004/06/16 小高 正嗣  
%        2005/06/24 小高 正嗣  
%        2005/09/21 竹広 真一 大幅改訂
%        2005/09/21 竹広 真一 拡散係数変更
%        2005/09/28 竹広 真一 拡散係数変更
%        2005/10/07 竹広 真一 標準実験解説を旧資源より移行
%        2015/07/02 竹広 真一 強制問題解説へと改訂
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\documentclass[a4j,12pt]{jarticle}
%
% プリアンブル
%
\usepackage{Dennou6}                      % dennou-sty-6 を使用

\Dtitle[回転球面浅水強制乱流]        % ヘッダ部タイトル, モデル名と同じ
{回転球面浅水強制乱流問題の定式化}
%
\Dauthor[竹広 真一]{竹広 真一}            % 著者
\Ddate{2015 年 7 月 2 日}                 % 日付
\Dnoparindent                             % 段落の字下げをしない
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 本文開始
%
% 基本的な章節構成は以下の通り. 
%
%   1. はじめに 
%   2. 支配方程式系
%   3. 離散化
%   4. 利用モジュールとその他の設定
%   5. 数値実験
%   6. 参考文献
%   -. 謝辞
% 
% 必要に応じて構成を変更すること.
%
\begin{document}

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

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

この文書では, 回転球面浅水強制乱流問題の定式化と
数値計算のための離散化について記述する. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\newpage                                  % 改ページ
\section{支配方程式系}
ここでは支配方程式系と境界条件を記す. 

  \subsection{支配方程式系}
    支配方程式系は回転系での球面浅水方程式系である. 



  \subsection{渦度・発散を用いた支配方程式系}

    球面浅水方程式系を渦度と発散を用いて表した式はつぎのとおりである. 
    \begin{eqnarray}
      \DP{\zeta }{t} &=& -\frac{1}{a(1-\mu ^{2})}
                          \DP{(f+\zeta)U}{\lambda}
                        - \frac{1}{a}\DP{(f+\zeta)V}{\mu} \nonumber \\
                      & & -d_\zeta + F(\lambda,\mu,t),
                        \Deqlab{渦度方程式} \\
      \DP{D}{t} &=& \frac{1}{a(1-\mu ^{2})}
                    \DP{(f+\zeta)V}{\lambda}
                    - \frac{1}{a}\DP{(f+\zeta)U}{\mu} 
                    - \Dlapla \left[g (h+h_{s}) + E\right] \nonumber \\
                & & -d_D,
                    \Deqlab{発散方程式} \\
      \DP{h}{t} &=& - \frac{1}{a(1-\mu ^{2})}\DP{(hU)}{\lambda}
                    - \frac{1}{a}\DP{(hV)}{\mu}  
                    -d_h,
                    \Deqlab{質量保存の式} 
    \end{eqnarray}
    %
    各記号の定義は\Dtabref{変数, 物理定数の定義}に示した.
    渦度 $\zeta $ と発散 $D$ , 単位質量あたりの運動エネルギー $E$, 
    演算子 $\nabla ^{2}$ は以下のように定義される.
    %
    \begin{eqnarray}
      \zeta &=& \frac{1}{a(1-\mu ^{2})}\DP{V}{\lambda} - \frac{1}{a}\DP{U}{\mu}, \\
      D     &=& \frac{1}{a(1-\mu ^{2})}\DP{U}{\lambda} + \frac{1}{a}\DP{V}{\mu}, \\
      E     &=& \frac{U^{2}+V^{2}}{2(1-\mu ^{2})}, \\ 
      \nabla^{2} A &=& \frac{1}{a^{2}(1-\mu ^{2})}\DP[2]{A}{\lambda}
                  + \frac{1}{a^{2}}\DP{}{\mu}\left[(1-\mu^{2})\DP{A}{\mu}\right].
    \end{eqnarray}
    %
    ここで $A$ は任意のスカラーである. 流線関数 $\psi $ と速度ポテンシャル
    $\chi$ を導入すると,
    %
    \begin{eqnarray}
      U &=& -\frac{1-\mu^ {2}}{a}\DP{\psi}{\mu}
            +\frac{1}{a}\DP{\chi}{\lambda}, \\
      V &=& \frac{1}{a}\DP{\psi}{\lambda}
            + \frac{1-\mu^{2}}{a}\DP{\chi}{\mu}, \\
      \zeta &=& \nabla ^{2} \psi, \\
      D     &=& \nabla ^{2} \chi
    \end{eqnarray}
    と表される. 

    \Deqref{渦度方程式}右辺の $F(\lambda,\mu,t)$ は渦度強制項である. 
    その具体的な形は実験設定の節で記述する. 
    $d_\zeta,d_D,d_h$ は角変数に対する散逸項を表している(下記参照). 

    \begin{table}[h]
     \begin{center}
      \begin{tabular}{c l l  } \hline
       記号    &\qquad\qquad&  変数/物理定数   \\ \hline 
       $\lambda$&& 経度            \\
       $\phi$   && 緯度            \\
       $\mu $   && $\sin$ 緯度 ($\equiv \sin \phi$) \\
       $t$      && 時間            \\ 
       $u$      && 経度方向流速    \\ 
       $v$      && 緯度方向流速    \\ 
       $U$      && $u\cos \phi $  \\ 
       $V$      && $v\cos \phi $  \\ 
       $h$      && 流体の全厚さ    \\ 
       $\zeta$  && 渦度            \\ 
       $D$      && 発散            \\
       $E$      && 運動エネルギー \\
       $\psi$   && 流線関数        \\
       $\chi$   && 速度ポテンシャル \\
       $f$      && $2\Omega \sin \phi$ \\
       $\Omega$ && 自転角速度      \\
       $g$      && 重力加速度     \\
       $a$      && 惑星半径       \\
       $F$      && 渦度強制       \\
       $\nu_{m}$  && 水平粘性係数   \\
       $\kappa_{h}$  && 水平拡散係数   \\
       $N_{m}$  && 超粘性の次数(水平ラプラシアンの階数) \\
       $N_{h}$  && 超拡散の次数(水平ラプラシアンの階数) \\ \hline
      \end{tabular} 
     \end{center}
     \caption{変数, 物理定数の定義}
     \Dtablab{変数, 物理定数の定義}
    \end{table}

    後の時間積分スキームの構成のために, 
    $h$ を平均成分 $\overline{h}$ とそれからの擾乱 $h'$ を用いて
    $h=\overline{h}+h'(\lambda,\mu,t)$ と表しておくと, 
    %
    \begin{eqnarray}
      \DP{\zeta }{t} &=& -\frac{1}{a(1-\mu ^{2})}
                          \DP{(f+\zeta)U}{\lambda}
                        - \frac{1}{a}\DP{(f+\zeta)V}{\mu} \nonumber \\
                      & & -d_\zeta  + F(\lambda,\mu,t),
                        \Deqlab{擾乱渦度方程式} \\
      \DP{D}{t} &=& \frac{1}{a(1-\mu ^{2})}
                    \DP{(f+\zeta)V}{\lambda}
                    - \frac{1}{a}\DP{(f+\zeta)U}{\mu} 
                    - \Dlapla \left[g (h'+h_{s}) + E\right] \nonumber \\
                & & -d_D
                    \Deqlab{擾乱発散方程式} \\
      \DP{h}{t} &=& - \frac{1}{a(1-\mu ^{2})}\DP{(h'U)}{\lambda}
                    - \frac{1}{a}\DP{(h'V)}{\mu}  
                    - \bar{h}D
                    - d_h
                    \Deqlab{擾乱質量保存の式} 
    \end{eqnarray}

    散逸過程としてレイリー摩擦, ニュートン冷却および超粘性拡散を
    考慮する. その具体的な形は
    %
    \begin{eqnarray}
      d_\zeta &=& - \zeta/\tau_R + (-1)^{N_m-1}\nu_{N_m}(\Dlapla)^{N_m}\zeta,\\
      d_D &=& - D/\tau_R + (-1)^{N_m-1}\nu_{N_m}(\Dlapla)^{N_m}D,\\
      d_h &=& - h'/\tau_N + (-1)^{N_h-1}\kappa_{N_h}(\Dlapla)^{N_h}h',
    \end{eqnarray}
    %
    ここで $\tau_R,\tau_N$ はそれぞれレイリー摩擦およびニュートン冷却の時定数, 
    $N_m$, $N_h$ は流れ場および高度場に対する超粘性・拡散の次数, 
    $\nu_{N_m}, \kappa_{N_m}$ は超粘性・拡散係数である. 

    \subsection{境界条件}

    境界条件は緯度, 経度方向ともに周期境界条件である. 
    \begin{eqnarray}
      \zeta (\lambda + 2\pi, \phi)  &=& \zeta (\lambda, \phi), \\
      \zeta (\lambda , \phi + 2\pi) &=& \zeta (\lambda, \phi).
    \end{eqnarray}
    などと表される.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\newpage                                  % 改ページ
\section{水平離散化}

  この節では方程式系の空間離散化について説明, 
  プログラム内で実際に用いられている方程式を記述する.

  格子点を経度方向には等間隔, 緯度方向にはガウス緯度上に置き, 
  物理量 A を格子点上の値で定義したものを $A_{i,j}$ と表すことにする. 
  添字 $i,j$ は格子点 $(\lambda _{i}, \mu _{j})$ 上の値であることを示す.

  $A_{i,j}$ に対する離散スペクトル変換を $\hat{A}_{m,n}$ であらわす.
  すなわち, 
  %
  \begin{eqnarray}
   A_{i,j} &=& \sum _{n=0}^{M}\sum _{m=-n}^{n}
                   P_{n}^{m}(\mu)e ^{im\lambda} \hat A_{m,n}, \\
   \hat A_{m,n} &=& \frac{1}{2I} \sum _{j=1}^{J} \sum _{i=1}^{I}
                    A _{i,j}e^{-i m \lambda _{i}}P_{n}^{m}(\mu_{j})w_{j},
  \end{eqnarray}
  %
  支配方程式\Deqref{擾乱渦度方程式},\Deqref{擾乱発散方程式},
  \Deqref{擾乱質量保存の式},
  をスペクトル変換すると
  %
  \begin{eqnarray}
    \DP{\hat\zeta _{m,n}}{t} 
         &=& -\widehat{\left[\frac{1}{a(1-\mu^{2})}
                        \DP{(f+\zeta_)U}{\lambda}\right]}_{m,n}
             - \widehat{\left[\frac{1}{a}\DP{(f+\zeta)V}{\mu}\right]}_{m,n} 
             - \widehat{d_\zeta}_{n,m}
             + \hat{F}_{n,m}(t)
          \Deqlab{スペクトル渦度方程式} \\
    \DP{\hat{D}_{m,n}}{t} 
         &=&   \widehat{\left[\frac{1}{a(1-\mu^2)}
                            \DP{(f+\zeta)V}{\lambda}\right]}_{m,n}
             - \widehat{\left[
                  \frac{1}{a}\DP{(f+\zeta)U}{\mu} \right]}_{m,n}\nonumber \\
         & & -g\Dlapla\hat{h}'_{m,n}
             - \widehat{\left[\Dlapla\left\{ gh_s+E \right\}\right]}_{m,n}
             - \widehat{d_D}_{,n,m}, 
         \Deqlab{スペクトル発散方程式} \\
    \DP{\hat{h}_{m,n}}{t} 
         &=& - \widehat{\left[\frac{1}{a(1-\mu^2)}
                              \DP{(h'U)}{\lambda}\right]}_{m,n}
             - \widehat{\left[\frac{1}{a}\DP{(h'V)}{\mu}\right]}_{m,n}
             -\bar{h}\hat{D}_{m,n}
             - \widehat{d_h}_{n,m}
         \Deqlab{スペクトル質量保存の式} 
  \end{eqnarray}
  %
  非線形項のスペクトル変換は, 格子点上での値を用いて非線形項を評価し, 
  それをスペクトル変換する(変換法)を適用することにする. 

  散逸過程のスペクトル表現は, 
  % 
  \begin{eqnarray}
    & &
    \widehat{d_\zeta}_{m,n} = -\nu_{m,n}\hat{\zeta}_{m,n}, \
    \widehat{d_D}_{m,n} = -\nu_{m,n}\hat{D}_{m,n}, \
    \widehat{d_h}_{m,n} = -\kappa_{m,n}\hat{h'}_{m,n}, \\
    & &
    \nu_{m,n}    =  1/\tau_R + \nu_{N_m}[n(n+1)]^{N_m}, \quad
    \kappa_{m,n} =  1/\tau_N + \kappa_{N_m}[n(n+1)]^{N_m}.
  \end{eqnarray}
  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\newpage                                  % 改ページ
\section{時間積分}

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

  % \subsection{ルンゲ・クッタ \& クランク・ニコルソン法}

    時間積分の数値スキームとして, 
    散逸項にはクランク・ニコルソン法, それ以外の項を
    4 次のルンゲ・クッタ法を組み合わせて用いる場合の定式化を行なう. 

    まず強制項を陽解法(オイラースキーム)により積分する. 
    %
    \begin{equation}
      \DD{\hat\zeta_{m,n}}{t} = \hat{F}_{n,m}(t) 
    \end{equation}
    %
    次に散逸項以外の項に関してルンゲ・クッタ法で積分する. 
    \Deqref{スペクトル渦度方程式},
    \Deqref{スペクトル発散方程式}, 
    \Deqref{スペクトル質量保存の式}から
    強制項と散逸項以外の時間変化を
    %
    \begin{eqnarray}
      & &\left(\DD{\hat\zeta_{m,n}}{t}\right)_{ND}
       [\hat{\zeta}_{m,n},\hat{D}_{m,n},\hat{h'}_{m,n}]
        = -\widehat{\left[\frac{1}{a(1-\mu^{2})}
                          \DP{(f+\zeta)U}{\lambda}\right]}_{m,n}
         - \widehat{\left[\frac{1}{a}\DP{(f+\zeta)V}{\mu}\right]}_{m,n} 
            \nonumber\\\Deqlab{スペクトル渦度時間変化散逸なし}\\
      & & \left(\DD{\hat{D}_{m,n}}{t}\right)_{ND}
       [\hat{\zeta}_{m,n},\hat{D}_{m,n},\hat{h'}_{m,n}] 
           = \widehat{\left[\frac{1}{a(1-\mu^2)}
                              \DP{(f+\zeta)V}{\lambda}\right]}_{m,n}
             - \widehat{\left[
                  \frac{1}{a}\DP{(f+\zeta)U}{\mu} \right]}_{m,n}
              \nonumber \\
           & & \hspace{60mm}
              -g\Dlapla\hat{h}^{'\tau}_{m,n}
               - \widehat{\left[\Dlapla\left\{ g h_s+E \right\}\right]}_{m,n}
            \Deqlab{スペクトル発散時間変化散逸なし} \\
      & & \left(\DD{\hat{h}_{m,n}}{t}\right)_{ND}
       [\hat{\zeta}_{m,n},\hat{D}_{m,n},\hat{h'}_{m,n}]
            = - \widehat{\left[\frac{1}{a(1-\mu^2)}
                                \DP{(h'U)}{\lambda}\right]}_{m,n}
               - \widehat{\left[\frac{1}{a}\DP{(h'V)}{\mu}\right]}_{m,n}
               -\bar{h}\hat{D}_{m,n}
        \nonumber\Deqlab{スペクトル質量時間変化散逸なし}\\
    \end{eqnarray}
    %
    と定義する. ここで非線形項で用いられている格子点で定義された変数は
    各スペクトル変数
    $\hat{\zeta}_{m,n},\hat{D}_{m,n},\hat{h'}_{m,n}$
    からスペクトル逆変換して求められる. 

    ルンゲ・クッタ法による積分の手順はつぎのようになる. 
    渦度方程式だけを例に手順を記す. 
    第 1 段は
    %
    \begin{eqnarray}
      \left(\DD{\hat\zeta_{m,n}}{t}\right)^{k1}
       &=& \left(\DD{\hat\zeta_{m,n}}{t}\right)_{ND}
          [\hat{\zeta}^{\tau}_{m,n},\hat{D}^{\tau}_{m,n},
           \hat{h'}^{\tau}_{m,n}],\\ 
      \hat\zeta^{\tau+1/2}_{m,n} &=& \hat\zeta^{\tau}_{m,n} + 
        \left(\DD{\hat\zeta_{m,n}}{t}\right)^{k1}\cdot \frac{(\Delta t)}{2},
    \end{eqnarray}
    %
    第 2 段は
    \begin{eqnarray}
      \left(\DD{\hat\zeta_{m,n}}{t}\right)^{k2}
       &=& \left(\DD{\hat\zeta_{m,n}}{t}\right)_{ND}
          [\hat{\zeta}^{\tau+1/2}_{m,n},\hat{D}^{\tau+1/2}_{m,n},
           \hat{h'}^{\tau+1/2}_{m,n}], \\
      \hat\zeta^{*\tau+1/2}_{m,n} &=& \hat\zeta^{\tau}_{m,n} + 
        \left(\DD{\hat\zeta_{m,n}}{t}\right)^{k2}\cdot \frac{(\Delta t)}{2},
    \end{eqnarray}
    %
    第 3 段は
    %
    \begin{eqnarray}
      \left(\DD{\hat\zeta_{m,n}}{t}\right)^{k3}
       &=& \left(\DD{\hat\zeta_{m,n}}{t}\right)_{ND}
          [\hat{\zeta}^{*\tau+1/2}_{m,n},\hat{D}^{*\tau+1/2}_{m,n},
           \hat{h'}^{*\tau+1/2}_{m,n}], \\
      \hat\zeta^{*\tau+1}_{m,n} &=& \hat\zeta^{\tau}_{m,n} + 
        \left(\DD{\hat\zeta_{m,n}}{t}\right)^{k3}\cdot (\Delta t),
    \end{eqnarray}
    %
    第 4 段は
    %
    \begin{eqnarray}
      \left(\DD{\hat\zeta_{m,n}}{t}\right)^{k4}
       = \left(\DD{\hat\zeta_{m,n}}{t}\right)_{ND}
          [\hat{\zeta}^{*\tau+1}_{m,n},\hat{D}^{*\tau+1}_{m,n},
           \hat{h'}^{*\tau+1}_{m,n}], 
    \end{eqnarray}
    %
    最終的に
    %
    \begin{eqnarray}
      \hat\zeta^{\tau+1,ND}_{m,n} &=& \hat\zeta^{\tau}_{m,n} + 
        (\Delta t)\left[
            \Dinv{6}\left(\DD{\hat\zeta_{m,n}}{t}\right)^{k1}
          + \Dinv{3}\left(\DD{\hat\zeta_{m,n}}{t}\right)^{k2}
          + \Dinv{3}\left(\DD{\hat\zeta_{m,n}}{t}\right)^{k3}
          + \Dinv{6}\left(\DD{\hat\zeta_{m,n}}{t}\right)^{k4}
          \right], \nonumber\\
    \end{eqnarray}
    %
    以上の手順を発散 $\hat{D}_{m,n}$, 変位擾乱 $\hat{h'}_{m,n}$ に
    ついても平行して計算していくことで, 非散逸項の積分が完了する. 

    ルンゲ・クッタ法による積分ののち, 散逸項からの修正を行なう. 
    渦度を例に取ると, 
    \begin{displaymath}
      \frac{\zeta^{\tau+1}_{m,n}-\zeta^{\tau+1,ND}_{m,n}}{\Delta t}
       = -\nu_{m,n}\frac{\zeta^{\tau+1}_{m,n}+\zeta^{\tau+1,ND}_{m,n}}{2}
    \end{displaymath}
    %
    したがって
    \begin{equation}
      \zeta^{\tau+1}_{m,n} = 
      \frac{1- \nu_{m_n}(\Delta t)/2}{1+ \nu_{m_n}(\Delta t)/2}
         \cdot\zeta^{\tau+1,ND}_{m,n},
    \end{equation}
    %
    と修正される. 発散 $\hat{D}_{m,n}$, 変位擾乱 $\hat{h'}_{m,n}$ に
    ついても同様に計算する. 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\newpage                                  % 改ページ
\section{使用モジュールとその他の設定}

  スペクトル変換と逆変換, 微分演算は SPMODEL ライブラリ (spml) の
  w\_module\_sjpack に含まれる関数を用いて行う. 
  スペクトル変換と逆変換の際の数値積分は台形公式を用いて行う. 
  spml が下位で使用する ISPACK の仕様から, 
  格子点数 $I$ は $2^n$ でなければならない. 
  非線形項の計算によって生じるエリアジングを防ぐため, 
  格子点数 $I,J$ は切断波数 $M$ に対し $I>3M+1,J>3M/2$ を満たすよう
  に与える.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\newpage                                  % 改ページ

\section{Scott and Polvani (2008) および Saito and Ishioka (2015) の実験}

  回転球面浅水強制乱流問題の設定として, 
  Scott and Polvani (2008) および Saito and Ishioka (2015)により行われた
  実験を参考にして設定する. 

  強制項の表現は
  %
  \begin{eqnarray}
    & &
    F(\lambda,\mu,t) 
    = 2\mbox{Re}\left[\sum_{n=n_f-\Delta n/2}^{n_f+\Delta n/2}
                      \sum_{m=1}^n \hat{f}_{m,n}(t) P_n^m(\mu)e^{im\lambda}
                \right], \\
    & &
    \hat{f}_{m,n}(t) 
       = \sqrt{\frac{2n(n+1)\varepsilon_0}{(2n+1)\Delta n}}\xi_{m,n}(t), \quad
    \xi_{m,n}(t) = \frac{e^{i\theta_{m,n}(t)}}{\sqrt{\Delta t}}.
  \end{eqnarray}
  %
  ここで  $\varepsilon_0$ はエネルギー注入率の大きさを表すパラメター,
  $P_n^m(\mu)$ は $2$ に規格化されたルジャンドル陪函数, 
  $\theta_{m,n}(t)$ は区間 $[0,2\pi]$ での一様乱数であり, 
  各時間ステップにおいて更新される
  \footnote{
    w\_module\_sjpack のスペクトル変数を用いれば
    l\_nm(n,m) (m>0) がスペクトル係数の実数部, 
    l\_nm(n,m) (m<0) がスペクトル係数の虚数部を
    格納しているので, $m \geq 1$ に対して
    %
    \begin{displaymath}
      w\_F(l\_(n,m)) = \sqrt{\frac{2n(n+1)\varepsilon_0}
                            {(2n+1)\Delta n\Delta t}}
                       \cos\theta_{m,n}(t),\quad
      w\_F(l\_(n,-m)) = \sqrt{\frac{2n(n+1)\varepsilon_0}
                             {(2n+1)\Delta n\Delta t}}
                       \sin\theta_{m,n}(t).
    \end{displaymath}
  }. 

  長さスケールと時間スケールを球半径 $a$ と自転周期 $2\pi/\Omega$ に
  選ぶのが簡単である. 
  %
  重要な力学パラメターは
  変形半径 $\Ddsty \frac{L_D}{a}=\frac{\sqrt{g\bar{h}}}{2\Omega}$. 
  %
  \begin{table}[h]
    \begin{center}
      \begin{tabular}{|l|c|} \hline
        実験パラメター                   & 設定値\\\hline
        切断波数 $N_T$                   & 170, (683) \\
        時間ステップ $\Delta t$          & 0.0125 \\
        球半径 $a$                       & 1\\
        回転角速度 $\Omega$              & $2\pi$ \\
        重力加速度 $g$                   & ? \\
        平均深さ   $\bar{h}$             & ? \\
        変形半径   $L_D$                 & $(1),0.1,(0.025)$ \\
        エネルギー注入率 $\varepsilon_0$ & $5.0\times 10^7$ \\
        強制の中心波数 $n_f$             & 42 \\
        強制の波数幅 $\Delta n$          & 4  \\
        レイリー摩擦の時定数 $\tau_R$    & $\infty$ \\
        ニュートン冷却の時定数 $\tau_N$  & $0.25 L_D^{-2}$\\
        超粘性次数 $N_m$                 & 4 ($\nabla^8$)\\
        超拡散次数 $N_m$                 & 4 ($\nabla^8$)\\
        超粘性係数 $\nu_4$               & $10/[N_T(N_T+1)]^4$\\
        超拡散係数 $\nu_4$               & $10/[N_T(N_T+1)]^4$\\\hline
      \end{tabular}
      \caption{実験パラメター一覧}
      \Dtablab{実験パラメター}
    \end{center}
  \end{table}

  \begin{table}[h]
    \begin{center}
      \begin{tabular}{|lc|cccc|} \hline
                  & OpenMP 数 & 1    & 2 & 4 & 6 \\\hline
                  & real time & ...s &   &   &   \\ 
        Intel ??? & usr time  & ...s &   &   &   \\\hline
      \end{tabular}
      \caption{10 日計算にかかる切断波数 $T170$ での計算時間 (sec)}
      \Dtablab{計算時間}
    \end{center}
  \end{table}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\newpage                                  % 改ページ

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

 \item 石岡圭一, 2004:
       スペクトル法による数値計算入門, 東京大学出版会, 232pp.

 \item Krishnamurti, T. N., H. S. Bedi, and V. M. Hardiker, 1998:
       An introduction to global spectral modeling. 
       Ocford Univ. Press, 253 pp.

%  \item Satoh, M., 2004: Atmospheric circulation dynamics and circulation 
%        models.
%        Springer-Verlag, 643 pp.

 \item 竹広 真一, 佐々木 洋平, 石岡 圭一, 小高 正嗣, 高橋 芳幸, 
   中島 健介, 石渡 正樹, 林 祥介, spmodel 開発グループ, 2013:
   階層的地球流体力学スペクトルモデル集 (SPMODEL),
   http://www.gfd-dennou.org/library/spmodel/, 地球流体電脳倶楽部.

 \item Scott, R. K., Polvani, L. M., 2007: 
   Forced-dissipative shallow-water turbulence
   on the sphere and the atmospheric circulation of the giant planets. 
   {\it J. Atmos. Sci.}, {\bf 64}, 3158--3176.

 \item Scott, R. K., Polvani, L. M., 2008: 
   Equatorial superrotation in shallow water atmospheres. 
   {\it Geophys. Res. Lett.}, {\bf 35}, L24202.

 \item Saito, I., Ishioka, K., 2015:
   Mechanism for the Formation of Equatorial Superrotation
   in Forced Shallow-Water Turbulence with Newtonian Cooling.
   {\it J. Atmos. Sci.}, {\bf 72}, 1466--1483.

\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 2014. ). 
% 本資源は, 著作者の諸権利に抵触しない(迷惑をかけない)限りにおいて自由に
% 利用していただいて構わない. なお, 利用する際には今一度自ら内容を確かめ
% ることをお願いする(無保証無責任原則).

% 本資源に含まれる元資源提供者(図等の版元等を含む)からは,  直接的な形で
% の WEB 上での著作権または使用許諾を得ていない場合があるが,  勝手ながら, 
% 「未来の教育」のための実験という学術目的であることをご理解いただけるも
% のと信じ,  学術標準の引用手順を守ることで諸手続きを略させていただいて
% いる.  本資源の利用者には, この点を理解の上,  注意して扱っていただける
% ようお願いする.  万一, 不都合のある場合には
% \begin{center}
%  spmodel@gfd-dennou.org
% \end{center}
% まで連絡していただければ幸いである. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{document}

