% 表題  回転球面上での 2 次元順圧流体の自由減衰乱流問題の定式化
%
% 履歴  2005/03/09 竹広 真一 
%       2007/11/25 竹広 真一  参照文書追加
%       2008/01/03 竹広 真一  エネルギースペクトル計算追加
%       2008/02/12 竹広 真一  ロスビー波振動数バグfix
%       2008/02/25 竹広 真一  全運動エネルギー, 全エンストロフィー時間変化追加
%       2008/03/02 竹広 真一  全運動エネルギー, 全エンストロフィー時間変化修正
%
\documentclass[12pt,a4j,notitlepage]{jarticle}
\usepackage{Dennou6}
\usepackage{psfrag}
\usepackage{ascmac}
\usepackage{times}
\Dtitle[回転球面 2 次元順圧自由減衰乱流]
       {回転球面上での 2 次元順圧流体の\\自由減衰乱流問題の定式化}
\Dauthor[竹広 真一]{竹広 真一, SPMODEL 開発グループ}

\Ddate{平成 20 年 3 月 4 日} 

\Dparskip

\pagestyle{Dheadings}

\begin{document}

\maketitle

\section{はじめに}

\section{支配方程式}

  支配方程式は, 回転球面上での渦度方程式である
  (「2 次元回転球面上での非発散順圧流体の定式化」参照):
  %
  \begin{equation}
    \Deqlab{渦度方程式}
    \DP{\zeta}{t} + \Dinv{a^2}J(\psi,\zeta) +
       \frac{2\Omega}{a^2}\DP{\psi}{\lambda} 
    = (-1)^{p+1}\nu_{2p}\left(\Dlapla + \frac{2}{a^2}\right)^p \zeta,
  \end{equation}
  %
  ここで $\zeta$ は渦度の動径成分, $t$ は時間, $\psi$ は流線関数, 
  $a$ は 球の半径, $\Omega$ は球の回転角速度, $\lambda$ は経度, 
  $\nu_p$ は高階粘性の係数である\footnote{$p=1$ で通常の粘性項になる}. 
  $\Ddsty J(f,g)\equiv \DP{f}{\lambda}\DP{g}{\mu}-\DP{f}{\mu}\DP{g}{\lambda}
   = \Dinv{\cos\varphi}\DP{f}{\lambda}\DP{g}{\varphi}
     -\DP{f}{\varphi}\Dinv{\cos\varphi}\DP{g}{\lambda}$
  はヤコビアンであり, $\varphi$ は緯度, 
  $\mu=\sin\varphi$ は sin 緯度 である. 

  流線関数と速度の経度・緯度成分 $u_\lambda,u_\varphi$ との関係は
  %
  \begin{equation}
    \Deqlab{速度流線関数関係}
    u_\lambda = - \Dinv{a}\DP{\psi}{\varphi}, \quad
    u_\varphi =  \Dinv{a\cos\varphi}\DP{\psi}{\lambda},
  \end{equation}
  %
  であり, したがって渦度と流線関数の関係は
  %
  \begin{eqnarray}
    \Deqlab{渦度流線関数関係}
    \zeta &=& (\Drot\Dvect{u})_r 
         =   \Dinv{a\cos\varphi}\DP{u_\varphi}{\lambda}
          - \Dinv{a\cos\varphi}\DP{(u_\lambda\cos\varphi)}{\varphi} \nonumber\\
        &=&   \Dinv{a^2\cos^2\varphi}\DP[2]{\psi}{\lambda}
          + \Dinv{a^2\cos\varphi}\DP{}{\varphi}
               \left(\cos\varphi\DP{\psi}{\varphi}\right)  \nonumber\\
        &=& \Dlapla \psi.
  \end{eqnarray}
  %
  ここで $\Ddsty \Dlapla \equiv \Dinv{a^2\cos^2\varphi}\DP[2]{}{\lambda}
          + \Dinv{a^2\cos\varphi}\DP{}{\varphi}
               \left(\cos\varphi\DP{}{\varphi}\right)$
  は半径 $a$ の球面上での 2 次元水平ラプラシアンである. 

\section{運動エネルギー・エンストロフィー}

  \subsection{定義}

     球面上の各点$(\lambda,\varphi)$における局所的な
     運動エネルギー $e_k$ は次のように定義される. 
     %
     \begin{equation}
       \Deqlab{局所運動エネルギー}
       e_k(\lambda,\varphi,t)=\frac{1}{2}|\Dvect{u}|^2
         =\frac{1}{2}(u_\lambda^2+u_\varphi^2).
     \end{equation}
     %
     全運動エネルギー $E_k$ はこれを球面全体で積分して得られる. 
     \begin{eqnarray*}
       \Deqlab{全運動エネルギー}
       E_k(t) &=& \Ddint_S\Dd S e_k(\lambda,\varphi,t)
               =  \Dinv{2}\Ddint_S\Dd S |\Dvect{u}|^2 \\
              &=& \Dinv{2}\int_0^{2\pi}\Dd \lambda
                  \int_{-\pi/2}^{\pi/2}\Dd \varphi \ a^2\cos\varphi
                  (u_\lambda^2+u_\varphi^2).
     \end{eqnarray*}

     球面上の各点$(\lambda,\varphi)$における局所的な
     エンストロフィーは $\Ddsty \frac{1}{2}\zeta^2$ で定義される. 
     全エンストロフィー $\cal Q$ はこれを球面全体で積分して得られる. 
     %
     \begin{equation}
       \Deqlab{全エンストロフィー}
       {\cal Q}(t) = \Ddint_S \Dd S \left(\frac{1}{2}\zeta^2\right) 
            = \Dinv{2}\int_0^{2\pi}\Dd \lambda
                  \int_{-\pi/2}^{\pi/2}\Dd \varphi \ a^2\cos\varphi\ \zeta^2
     \end{equation}

  \subsection{エネルギー・エンストロフィーのスペクトル分解}

    \Deqref{速度流線関数関係}を用いて全エネルギーを表現すると
    %
    \begin{eqnarray*}
      E_k(t) &=& \Dinv{2}\int_0^{2\pi}\Dd \lambda
                  \int_{-\pi/2}^{\pi/2}\Dd \varphi \ a^2\cos\varphi
                  \left[
                      \left(\Dinv{a}\DP{\psi}{\varphi}\right)^2
                    + \left(\Dinv{a\cos\varphi}\DP{\psi}{\lambda}\right)^2
                  \right] \\
             &=& \Dinv{2}\int_0^{2\pi}\Dd \lambda
                  \int_{-\pi/2}^{\pi/2}\Dd \varphi \ \cos\varphi
                  \left[
                      \left(\DP{\psi}{\varphi}\right)^2
                    + \left(\Dinv{\cos\varphi}\DP{\psi}{\lambda}\right)^2
                  \right].
     \end{eqnarray*}
     %
     ここで流線関数を球面調和函数で展開する. すなわち, 
     %
     \begin{equation}
       \psi(\lambda,\varphi)
           =\sum_{n=0}^\infty\sum_{m=-n}^n
               \tilde{\psi}_n^m(t)Y_n^m(\lambda,\varphi)
     \end{equation}
     %
     これを上の式に代入すると
     %
     \begin{eqnarray*}
       E_k(t) &=& \Dinv{2}\int_0^{2\pi}\Dd \lambda
                  \int_{-\pi/2}^{\pi/2}\Dd \varphi \ \cos\varphi
                  \left[
                    \left(\sum_{n=0}^\infty\sum_{m=-n}^n
                          \sum_{n'=0}^\infty\sum_{m'=-n'}^{n'}
                          \tilde{\psi}_n^m\DP{Y_n^m}{\varphi}
                          \tilde{\psi}_{n'}^{m'}
                             \DP{Y_{n'}^{m'}}{\varphi}\right) \right.\\
              & &\left.\hspace{35mm}
                  + \left(\Dinv{\cos^2\varphi}
                            \sum_{n=0}^\infty\sum_{m=-n}^n
                            \sum_{n'=0}^\infty\sum_{m'=-n'}^{n'}
                            \tilde{\psi}_n^m\DP{Y_n^m}{\lambda}
                            \tilde{\psi}_{n'}^{m'}
                             \DP{Y_{n'}^{m'}}{\lambda}\right)
                  \right]\\
              &=& \Dinv{2}\sum_{n=0}^\infty\sum_{m=-n}^n
                  \sum_{n'=0}^\infty\sum_{m'=-n'}^{n'}
                     \tilde{\psi}_n^m\tilde{\psi}_{n'}^{m'}\\
              & & \hspace*{25mm}
                  \int_0^{2\pi}\Dd \lambda
                  \int_{-\pi/2}^{\pi/2}\Dd \varphi \ \cos\varphi
                    \left(  \DP{Y_n^m}{\varphi}\DP{Y_{n'}^{m'}}{\varphi}
                           +\Dinv{\cos^2\varphi}
                              \DP{Y_n^m}{\lambda}\DP{Y_{n'}^{m'}}{\lambda}
                     \right)\\
              &=& \Dinv{2}\sum_{n=0}^\infty\sum_{m=-n}^n
                  \sum_{n'=0}^\infty\sum_{m'=-n'}^{n'}
                     \tilde{\psi}_n^m\tilde{\psi}_{n'}^{m'}
                  \int_0^{2\pi}\Dd \lambda
                      \left[Y_n^m\DP{Y_{n'}^{m'}}{\varphi}\cos\varphi
                      \right]_{-\pi/2}^{\pi/2}\\
              & & \hspace*{45mm}
                  -\int_0^{2\pi}\Dd \lambda
                   \int_{-\pi/2}^{\pi/2}\Dd \varphi \ 
                          Y_n^m\DP{}{\varphi}\cos\varphi
                               \DP{Y_{n'}^{m'}}{\varphi}\\
              & & \hspace*{45mm}
                  +\int_{-\pi/2}^{\pi/2}\Dd \varphi \ 
                        \left[
                            \Dinv{\cos^2\varphi}
                              Y_n^m\DP{Y_{n'}^{m'}}{\lambda}
                        \right]_0^{2\pi}\\
              & & \hspace*{45mm}
                  -\int_0^{2\pi}\Dd \lambda
                   \int_{-\pi/2}^{\pi/2}\Dd \varphi \ 
                            \Dinv{\cos^2\varphi}
                              Y_n^m\DP[2]{Y_{n'}^{m'}}{\lambda}
                   \\
              &=& -\Dinv{2}\sum_{n=0}^\infty\sum_{m=-n}^n
                  \sum_{n'=0}^\infty\sum_{m'=-n'}^{n'}
                     \tilde{\psi}_n^m\tilde{\psi}_{n'}^{m'}\\
              & & \hspace*{10mm}
                   \int_0^{2\pi}\Dd \lambda
                   \int_{-\pi/2}^{\pi/2}\Dd \varphi \cos\varphi Y_n^m
                     \left[\Dinv{\cos\varphi}
                             \DP{}{\varphi}                             
                             \left( \cos\varphi \DP{Y_{n'}^{m'}}{\varphi}
                             \right)
                            + \Dinv{\cos^2\varphi}
                              \DP[2]{Y_{n'}^{m'}}{\lambda}
                     \right]\\
              &=& \Dinv{2}\sum_{n=0}^\infty\sum_{m=-n}^n
                  \sum_{n'=0}^\infty\sum_{m'=-n'}^{n'}
                     \tilde{\psi}_n^m\tilde{\psi}_{n'}^{m'}
                   \int_0^{2\pi}\Dd \lambda
                   \int_{-\pi/2}^{\pi/2}\Dd \varphi \cos\varphi 
                           \ n'(n'+1)Y_n^m Y_{n'}^{m'}\\
              &=& \Dinv{2}\sum_{n=0}^\infty\sum_{m=-n}^n
                  \sum_{n'=0}^\infty\sum_{m'=-n'}^{n'}
                     \tilde{\psi}_n^m\tilde{\psi}_{n'}^{m'}
                  \delta_{nn'}\delta_{mm'} n'(n'+1) N_l^m\\
              &=& \Dinv{2}\sum_{n=0}^\infty\sum_{m=-n}^n
                  N_l^m n(n+1) |\tilde{\psi}_n^m(t)|^2. 
    \end{eqnarray*}
    %
    ここで $N_n^m$ は
    %
    \begin{equation}
      N_n^m \equiv \Ddint_S \Dd S |Y_n^m(\lambda,\varphi)|^2 
         = \int_0^{2\pi}\Dd\lambda\int_{-\pi/2}^{\pi/2} \Dd\varphi 
              \cos\varphi|Y_n^m(\lambda,\varphi)|^2,
    \end{equation}
    %
    であり, 球面調和函数の規格化によって値が異なる
    \footnote{ISPACK では 2 に規格化されたルジャンドル函数を用いているので
      $N_n^m=4\pi$ である(要確認)}.
    したがって, エネルギースペクトル密度 ${\mathcal E}(n,t)$を
    %
    \begin{equation}
      \Deqlab{エネルギースペクトル密度}
      {\mathcal E}(n,t) \equiv
         \Dinv{2}\sum_{m=-n}^n N_n^m n(n+1) |\tilde{\psi}_n^m(t)|^2
    \end{equation}
    %
    と定義できて, 
    %
    \begin{equation}
      E_k(t)=\sum_{n=0}^\infty {\mathcal E}(n,t),
    \end{equation}
    %
    となる. 

    特に, 球面調和函数が $4\pi$ で規格化されている場合, すなわち
    $N_n^m=4\pi$ の場合には, この規格化係数を省略したエネルギー密度を
    %
    \begin{equation}
      \Deqlab{エネルギースペクトル密度係数省略}
      {\mathcal E}(n,t) \equiv
         \Dinv{2}n(n+1)\sum_{m=-n}^n |\tilde{\psi}_n^m(t)|^2
    \end{equation}
    %
    と定義することができる. その場合, 
    %
    \begin{equation}
      E_k(t)=4 \pi \sum_{n=0}^\infty {\mathcal E}(n,t),
    \end{equation}
    %
    となる. 

    同様に, \Deqref{渦度流線関数関係}を用いて全エンストロフィーを表し, 
    流線関数を球面調和函数展開すれば, 
    %
    \begin{eqnarray*}
       {\cal Q}(t) 
          &=& \frac{1}{2}\int_0^{2\pi}\Dd \lambda
                  \int_{-\pi/2}^{\pi/2}\Dd \varphi \ a^2\cos\varphi\ 
                  (\Dlapla \psi)^2\\
          &=& \frac{1}{2}\int_0^{2\pi}\Dd \lambda
                  \int_{-\pi/2}^{\pi/2}\Dd \varphi \ a^2\cos\varphi\ 
                    \left(\Dlapla\sum_{n=0}^\infty\sum_{m=-n}^n
                          \tilde{\psi}_n^m Y_n^m
                    \right)
                    \left(\Dlapla 
                          \sum_{n'=0}^\infty\sum_{m'=-n'}^{n'}
                          \tilde{\psi}_{n'}^{m'}Y_{n'}^{m'}
                    \right)\\
          &=& \frac{1}{2}
              \sum_{n=0}^\infty\sum_{m=-n}^n
              \sum_{n'=0}^\infty\sum_{m'=-n'}^{n'}
              \int_0^{2\pi}\Dd \lambda
                  \int_{-\pi/2}^{\pi/2}\Dd \varphi \ a^2\cos\varphi\ \\
          & & \hspace{50mm}
              \times\left(
                         -\frac{n(n+1)}{a^2}\tilde{\psi}_n^m Y_n^m
                    \right)
                    \left(-\frac{n'(n'+1)}{a^2}
                          \tilde{\psi}_{n'}^{m'}Y_{n'}^{m'}
                    \right)\\
          &=& \frac{1}{2}
              \sum_{n=0}^\infty\sum_{m=-n}^n
              \sum_{n'=0}^\infty\sum_{m'=-n'}^{n'}
              \frac{n(n+1)n'(n'+1)}{a^2}
                   \tilde{\psi}_n^m\tilde{\psi}_{n'}^{m'}
              \int_0^{2\pi}\Dd \lambda
                  \int_{-\pi/2}^{\pi/2}\Dd \varphi \ \cos\varphi\ 
                    Y_n^m Y_{n'}^{m'}\\
          &=& \frac{1}{2}\sum_{n=0}^\infty\sum_{m=-n}^n
              \sum_{n'=0}^\infty\sum_{m'=-n'}^{n'}
              \frac{n(n+1)n'(n'+1)}{a^2}
                   \tilde{\psi}_n^m\tilde{\psi}_{n'}^{m'}
                   N_l^m \delta_{n,n'} \delta_{m,m'}\\
          &=& \frac{1}{2}\sum_{n=0}^\infty\sum_{m=-n}^n
              N_l^m \frac{n^2(n+1)^2}{a^2}|\tilde{\psi}_n^m(t)|^2. 
    \end{eqnarray*}
    %
    したがって, エンストロフィースペクトル密度 ${\mathcal Q}(n,t)$を
    %
    \begin{equation}
      \Deqlab{エンストロフィースペクトル密度}
      {\mathcal Q}(n,t) \equiv
         \Dinv{2a^2}\sum_{m=-n}^n N_n^m n^2(n+1)^2 |\tilde{\psi}_n^m(t)|^2
    \end{equation}
    %
    と定義できて, 
    %
    \begin{equation}
      Q(t)=\sum_{n=0}^\infty {\mathcal Q}(n,t),
    \end{equation}
    %
    となる. 

    エンストロフィーについても, 
    球面調和函数が $4\pi$ で規格化されている場合, すなわち
    $N_l^m=4\pi$ の場合には, 規格化係数を省略したエンストロフィー密度を
    %
    \begin{equation}
      \Deqlab{エンストロフィースペクトル密度係数省略}
      {\mathcal Q}(n,t) \equiv
         \frac{n^2(n+1)^2}{2a^2}\sum_{m=-n}^n |\tilde{\psi}_n^m(t)|^2
    \end{equation}
    %
    と定義することができる. その場合, 
    %
    \begin{equation}
      Q(t)=4 \pi \sum_{n=0}^\infty {\mathcal Q}(n,t),
    \end{equation}
    %
    となる. 
             
  \subsection{全運動エネルギー, 全エンストロフィーの時間変化}

    全運動エネルギーの式は渦度方程式\Deqref{渦度方程式}に $-\psi$ を
    かけて球面全体で積分することにより得ることができる. 
    %
    \begin{displaymath}
      -\int\int dS \psi\DP{\zeta}{t} -\int\int dS \psi\Dinv{a^2}J(\psi,\zeta)
         -\int\int dS \psi\frac{2\Omega}{a^2}\DP{\psi}{\lambda} 
      = -\int\int dS \psi(-1)^{p+1}\nu_{2p}\left(\Dlapla + \frac{2}{a^2}\right)^p \zeta,
    \end{displaymath}
    %
    ただし
    %
    \begin{displaymath}
      \int\int dS 
      =a^2\int_{0}^{2\pi}d\lambda\ \int_{-\pi/2}^{\pi/2}d\varphi\cos\varphi
      =  a^2\int_{0}^{2\pi}d\lambda\int_{-1}^{1}d\mu
    \end{displaymath}
    %
    は球面全積分である. 
    第 1 項目は
    \begin{eqnarray*}
      & &
      -\int\int dS \psi\DP{\zeta}{t} 
       = -\int\int dS \psi\DP{}{t} \Dlapla \psi\\
      & &
       = -\int\int dS \psi\DP{}{t} \left[
            \Dinv{a^2\cos^2\varphi}\DP[2]{\psi}{\lambda}
          + \Dinv{a^2\cos\varphi}\DP{}{\varphi}
               \left(\cos\varphi\DP{\psi}{\varphi}\right)\right]\\
      & &
       = \int\int dS  
            \Dinv{a^2\cos^2\varphi}\DP{\psi}{\lambda}
            \Dinv{a^2\cos^2\varphi}\DP{}{\lambda} \DP{\psi}{t}
        + \int\int dS 
            \Dinv{a^2}\DP{\psi}{\varphi}\DP{}{\varphi}\DP{\psi}{t}\\
      & &
       = \int\int dS  
          \frac{1}{2}\DP{}{t}\Dinv{a^2\cos^2\varphi}\left(\DP{\psi}{\lambda}\right)^2
        + \int\int dS 
           \frac{1}{2}\DP{}{t} \Dinv{a^2}\left(\DP{\psi}{\varphi}\right)^2\\
      & &
       = \DD{}{t}\int\int dS \frac{1}{2}(u_\lambda^2+u_\varphi^2)
       = \DD{E_k}{t}
    \end{eqnarray*}
    %
    第 2 項目は
    \begin{eqnarray*}
      \int\int dS \psi\Dinv{a^2}J(\psi,\zeta)
      = \int\int dS \Dinv{a^2}J(\frac{1}{2}\psi^2,\zeta)
    \end{eqnarray*}
    %
    ここでヤコビアン $J(f,g)$ は
    %
    \begin{eqnarray*}
      J(f,g) = \DP{f}{\lambda}\DP{g}{\mu}-\DP{f}{\mu}\DP{g}{\lambda}
             = \DP{}{\lambda}\left(f\DP{g}{\mu}\right)
              -\DP{}{\mu}\left(f\DP{g}{\lambda}\right)
    \end{eqnarray*}
    %
    と変形できるので全球積分すると 0 となる. 
    %
    第 3 項目は
    \begin{eqnarray*}
      -\int\int dS \psi\frac{2\Omega}{a^2}\DP{\psi}{\lambda} 
      = -\int\int dS \frac{\Omega}{a^2}\DP{\psi^2}{\lambda} 
      = 0.
    \end{eqnarray*}
    %
    最後の粘性項は球面調和函数展開して評価する. 
    %
    \begin{eqnarray*}
      & &
       -\int\int dS \psi(-1)^{p+1}\nu_{2p}\left(\Dlapla + \frac{2}{a^2}\right)^p \zeta 
       = \int\int dS \psi (-1)^p\nu_{2p}\left(\Dlapla + \frac{2}{a^2}\right)^p \Dlapla\psi\\
      & &
       = \int\int dS 
           \sum_{n'=0}^\infty\sum_{m'=-n'}^{n'}
               \tilde{\psi}_{n'}^{m'}(t)Y_{n'}^{m'}(\lambda,\varphi)\\
      & &\qquad\qquad
           \sum_{n=0}^\infty\sum_{m=-n}^n (-1)^p\nu_{2p}
               \left(\frac{-n(n+1) + 2}{a^2}\right)^p 
               \left(\frac{-n(n+1)}{a^2}\right)
                \tilde{\psi}_n^m(t)Y_n^m(\lambda,\varphi) \\
      & &
       = -\sum_{n=0}^\infty\sum_{m=-n}^n \nu_{2p}
               \left(\frac{n(n+1) - 2}{a^2}\right)^p 
               \left(\frac{n(n+1)}{a^2}\right) |\tilde{\psi}_n^m(t)|^2
                 \int\int dS Y_n^m(\lambda,\varphi)^2 \\
      & &
       = -\sum_{n=0}^\infty\sum_{m=-n}^n \nu_{2p}
               \left(\frac{n(n+1) - 2}{a^2}\right)^p 
               \left(\frac{n(n+1)}{a^2}\right) a^2 N_n^m |\tilde{\psi}_n^m(t)|^2\\
      & &
       = -\sum_{n=0}^\infty\sum_{m=-n}^n \nu_{2p}
               \left(\frac{n(n+1) - 2}{a^2}\right)^p 
               n(n+1) N_n^m |\tilde{\psi}_n^m(t)|^2
    \end{eqnarray*}
    %
    ただし, $N_n^m$ は球面調和函数の規格化定数である. 
    よって全運動エネルギーの時間変化は
    %
    \begin{equation}
      \DD{E_k}{t} 
        = -\sum_{n=0}^\infty\sum_{m=-n}^n \nu_{2p}
               \left(\frac{n(n+1) - 2}{a^2}\right)^p 
               n(n+1) N_n^m |\tilde{\psi}_n^m(t)|^2. 
    \end{equation}
    %
    特に球面調和函数の規格化定数が $N_n^m=4\pi$ の場合には
    %
    \begin{eqnarray}
      \DD{E_k}{t} 
        &=& -4\pi\sum_{n=0}^\infty\sum_{m=-n}^n \nu_{2p}
               \left(\frac{n(n+1) - 2}{a^2}\right)^p 
               n(n+1) |\tilde{\psi}_n^m(t)|^2 \\
        &=& -4\pi\sum_{n=0}^\infty\sum_{m=-n}^n \nu_{2p}
               \left(\frac{n(n+1) - 2}{a^2}\right)^p 
               2 {\mathcal E}(n,t).
    \end{eqnarray}
    %
    ただし ${\mathcal E}(n,t)$ は\Deqref{エネルギースペクトル密度係数省略}で定義さ
    れるスペクトル密度である. 
    全運動エネルギーの代わりに球面平均したエネルギー $\bar{E_k}$ の時間変化は
    全体を $4\pi a^2$ でわって
    %
    \begin{eqnarray}
      \DD{\bar{E_k}}{t} 
        &=& -\sum_{n=0}^\infty\sum_{m=-n}^n \nu_{2p}
               \left(\frac{n(n+1) - 2}{a^2}\right)^p 
               \frac{n(n+1)}{a^2} |\tilde{\psi}_n^m(t)|^2\\
        &=& -\sum_{n=0}^\infty\sum_{m=-n}^n \nu_{2p}
               \left(\frac{n(n+1) - 2}{a^2}\right)^p 
               \frac{2}{a^2} {\mathcal E}(n,t).
    \end{eqnarray}
    %
    $E_k, \bar{E_k}$ も球面調和函数展開で表しておくと
    %
    \begin{eqnarray}
      E_k &=& 4\pi\sum_{n=0}^\infty\sum_{m=-n}^n
                  \frac{n(n+1)}{2} |\tilde{\psi}_n^m(t)|^2
           = 4\pi\sum_{n=0}^\infty  {\mathcal E}(n,t), \\
      \bar{E_k} &=& \sum_{n=0}^\infty\sum_{m=-n}^n
                  \frac{n(n+1)}{2a^2} |\tilde{\psi}_n^m(t)|^2
           = \Dinv{a^2}\sum_{n=0}^\infty  {\mathcal E}(n,t).
    \end{eqnarray}

    全エンストロフィーの時間変化の式は渦度方程式\Deqref{渦度方程式}に 
    $\Ddsty \zeta=\Dinv{a^2}\Dlapla\psi$ をかけて
    球面全体で積分することにより得ることができる. 
    %
    \begin{displaymath}
      \int\int dS \zeta\DP{\zeta}{t} +\int\int dS \zeta\Dinv{a^2}J(\psi,\zeta)
         \int\int dS \zeta\frac{2\Omega}{a^2}\DP{\psi}{\lambda} 
      = \int\int dS \zeta(-1)^{p+1}\nu_{2p}
          \left(\Dlapla + \frac{2}{a^2}\right)^p \zeta,
    \end{displaymath}
    %
    第 1 項目は
    %
    \begin{displaymath}
      \int\int dS \zeta\DP{\zeta}{t}
      = \int\int dS \DP{}{t}\left(\frac{1}{2}\zeta^2\right)
      = \DD{Q}{t}.
    \end{displaymath}
    %
    第 2 項目は
    \begin{displaymath}
      \int\int dS \zeta\Dinv{a^2}J(\psi,\zeta)
      \int\int dS \Dinv{2a^2}J(\psi,\zeta^2) = 0.
    \end{displaymath}
    %
    第 3 項目は
    %
    \begin{eqnarray*}
      \int\int dS \zeta\frac{2\Omega}{a^2}\DP{\psi}{\lambda} 
      = \frac{2\Omega}{a^2}\int\int dS
         \left[
            \Dinv{a^2\cos^2\varphi}\DP[2]{\psi}{\lambda}
          + \Dinv{a^2\cos\varphi}\DP{}{\varphi}
               \left(\cos\varphi\DP{\psi}{\varphi}\right) \right]
         \DP{\psi}{\lambda} 
    \end{eqnarray*}
    %
    1 つ目の被積分関数は
    %
    \begin{eqnarray*}
      \Dinv{a^2\cos^2\varphi}\DP[2]{\psi}{\lambda} \DP{\psi}{\lambda} 
     = \Dinv{2a^2\cos^2\varphi}\DP[]{}{\lambda} \DP[][]{\psi}{\lambda}^2
    \end{eqnarray*}
    %
    と変形できて, この全球積分は 0 となる. 
    2 つ目の被積分関数は
    %
    \begin{eqnarray*}
      & &
              \Dinv{a^2\cos\varphi}\DP{}{\varphi}
               \left(\cos\varphi\DP{\psi}{\varphi}\right)\DP{\psi}{\lambda} \\
      & &
          =   \DP{}{\lambda}\left[\psi
                \Dinv{a^2\cos\varphi}\DP{}{\varphi}
                   \left(\cos\varphi\DP{\psi}{\varphi}\right)\right]
             - \psi\Dinv{a^2\cos\varphi}\DP{}{\varphi}
                   \left(\cos\varphi\frac{\partial^2 \psi}
                                         {\partial\varphi\partial\lambda}\right),
    \end{eqnarray*}
    %
    一方で
    \begin{eqnarray*}
      & &
              \Dinv{a^2\cos\varphi}\DP{}{\varphi}
               \left(\cos\varphi\DP{\psi}{\varphi}\right)\DP{\psi}{\lambda}
          =   \Dinv{a^2\cos\varphi}\DP{}{\varphi}\left[
                 \DP{\psi}{\lambda}
                   \left(\cos\varphi\DP{\psi}{\varphi}\right)\right]
             - \frac{\cos\varphi}{a^2\cos\varphi}\DP{\psi}{\varphi}
                    \frac{\partial^2 \psi}
                         {\partial\varphi\partial\lambda}\\
      & &
          =   \Dinv{a^2\cos\varphi}\DP{}{\varphi}\left[
                 \DP{\psi}{\lambda}
                   \left(\cos\varphi\DP{\psi}{\varphi}\right)\right]
             - \Dinv{a^2\cos\varphi}\DP{}{\varphi}\left[
                   \cos\varphi\psi
                    \frac{\partial^2 \psi}
                         {\partial\varphi\partial\lambda}\right]\\
      & & \qquad
             + \psi\Dinv{a^2\cos\varphi}\DP{}{\varphi}\left[
                   \cos\varphi
                    \frac{\partial^2 \psi}
                         {\partial\varphi\partial\lambda}\right].
    \end{eqnarray*}
    %
    これらを1/2 倍して足し合わせると最後の項がキャンセルして
    %
    \begin{eqnarray*}
      & &
      \Dinv{a^2\cos\varphi}\DP{}{\varphi}
          \left(\cos\varphi\DP{\psi}{\varphi}\right)\DP{\psi}{\lambda}\\
      & &
      = \Dinv{2}\DP{}{\lambda}\left[\psi
                \Dinv{a^2\cos\varphi}\DP{}{\varphi}
                   \left(\cos\varphi\DP{\psi}{\varphi}\right)\right]
      + \Dinv{2a^2\cos\varphi}\DP{}{\varphi}\left[
                 \DP{\psi}{\lambda}
                   \left(\cos\varphi\DP{\psi}{\varphi}\right)\right]\\
      & & \qquad
      - \Dinv{2a^2\cos\varphi}\DP{}{\varphi}\left[
                  \cos\varphi\psi
                    \frac{\partial^2 \psi}
                         {\partial\varphi\partial\lambda}\right].
    \end{eqnarray*}
    %
    これを全球積分すると 0 になる. 
    最後に粘性項を球面調和関数展開で表現すると
    %
    \begin{eqnarray*}
      & &
            \int\int dS \zeta (-1)^{p+1}\nu_{2p}
               \left(\Dlapla + \frac{2}{a^2}\right)^p \zeta
          = \int\int dS \Dlapla\psi (-1)^{p+1}\nu_{2p}
               \left(\Dlapla + \frac{2}{a^2}\right)^p \Dlapla\psi\\
      & &
          = \int\int dS (-1)^{p+1}\nu_{2p}
               \sum_{n'=0}^\infty\sum_{m'=-n'}^{n'}
               \left(\frac{-n'(n'+1)}{a^2}\right)
               \tilde{\psi}_{n'}^{m'}(t)Y_{n'}^{m'}(\lambda,\varphi)\\
      & & \qquad
            \sum_{n=0}^\infty\sum_{m=-n}^n (-1)^{p+1}\nu_{2p} 
               \left(\frac{-n(n+1) + 2}{a^2}\right)^p 
               \left(\frac{-n(n+1)}{a^2}\right)
                \tilde{\psi}_n^m(t)Y_n^m(\lambda,\varphi) \\
      & &
          = \sum_{n=0}^\infty\sum_{m=-n}^n (-1)^{p+1}\nu_{2p}
               \left(\frac{-n(n+1) + 2}{a^2}\right)^p 
               \left(\frac{-n(n+1)}{a^2}\right)^2 |\tilde{\psi}_n^m(t)|^2
            \int\int dS  |Y_n^m(\lambda,\varphi)|^2 \\
      & &
          = -\sum_{n=0}^\infty\sum_{m=-n}^n \nu_{2p}
               \left(\frac{n(n+1) - 2}{a^2}\right)^p 
               \frac{n^2(n+1)^2}{a^2}N_n^m|\tilde{\psi}_n^m(t)|^2 \\
      & &
          = -\sum_{n=0}^\infty \nu_{2p}
               \left(\frac{n(n+1) - 2}{a^2}\right)^p 
               2 {\mathcal Q}(n,t)
    \end{eqnarray*}
    %
    したがって全エンストロフィーの時間変化の式は
    %
    \begin{equation}
      \DD{Q}{t}
         =   -\sum_{n=0}^\infty\sum_{m=-n}^n \nu_{2p}
               \left(\frac{n(n+1) - 2}{a^2}\right)^p 
               \frac{n^2(n+1)^2}{a^2}N_n^m|\tilde{\psi}_n^m(t)|^2\\
    \end{equation}
    %
    特に球面調和函数の規格化係数 $N_n^m=4\pi$ の場合には
    全エンストロフィーの時間変化の式は
    %
    \begin{eqnarray*}
      \DD{Q}{t}
         &=& -4\pi\sum_{n=0}^\infty\sum_{m=-n}^n \nu_{2p}
               \left(\frac{n(n+1) - 2}{a^2}\right)^p 
               \frac{n^2(n+1)^2}{a^2}|\tilde{\psi}_n^m(t)|^2\\
         &=& -4\pi\sum_{n=0}^\infty \nu_{2p}
               \left(\frac{n(n+1) - 2}{a^2}\right)^p 
               2 {\mathcal Q}(n,t).
    \end{eqnarray*}
    %
    ただし ${\mathcal Q}(n,t)$ は\Deqref{エンストロフィースペクトル密度係数省略}で
    定義されるエンストロフィースペクトル密度である. 
    全エンストロフィーの代わりに球面平均したエンストロフィー $\bar{Q}$ の
    時間変化は全体を $4\pi a^2$ でわって
    %
    \begin{eqnarray*}
      \DD{\bar{Q}}{t}
         &=& -\sum_{n=0}^\infty\sum_{m=-n}^n \nu_{2p}
               \left(\frac{n(n+1) - 2}{a^2}\right)^p 
               \frac{n^2(n+1)^2}{a^4}|\tilde{\psi}_n^m(t)|^2\\
         &=& -\sum_{n=0}^\infty \nu_{2p}
               \left(\frac{n(n+1) - 2}{a^2}\right)^p 
               \frac{2}{a^2} {\mathcal Q}(n,t).
    \end{eqnarray*}
    %
    $Q,\bar{Q}$ も球面調和関数展開表現で表しておくと
    \begin{eqnarray}
       Q  &=& 4\pi\frac{1}{2}\sum_{n=0}^\infty\sum_{m=-n}^n
                   \frac{n^2(n+1)^2}{a^2}|\tilde{\psi}_n^m(t)|^2
           = 4\pi\sum_{n=0}^\infty {\mathcal Q}(n,t), \\
       \bar{Q}  &=& \frac{1}{2}\sum_{n=0}^\infty\sum_{m=-n}^n
                   \frac{n^2(n+1)^2}{a^4}|\tilde{\psi}_n^m(t)|^2
           = \frac{1}{a^2}\sum_{n=0}^\infty {\mathcal Q}(n,t), \\
    \end{eqnarray}


\section{数値計算 : 球面調和函数展開と時間積分}

  \subsection{球面調和函数展開}

    スペクトル法による数値計算を行うために, 
    各物理量を球面調和函数 $Y_n^m(\lambda,\varphi)$で
    展開する. すなわち, 
    %
    \begin{eqnarray}
      \zeta(\lambda,\varphi)
          &=&\sum_{n=0}^\infty\sum_{m=-n}^n
              \tilde{\zeta}_n^m(t)Y_n^m(\lambda,\varphi), \\
      \psi(\lambda,\varphi)
          &=&\sum_{n=0}^\infty\sum_{m=-n}^n
              \tilde{\psi}_n^m(t)Y_n^m(\lambda,\varphi)
    \end{eqnarray}
    %
    ここで, 球面調和函数は数値計算のため実関数として
    つぎのように定義されたものを扱う. 
    %
    \begin{eqnarray}
      Y_n^m(\lambda,\varphi) 
          &=& P_n^{m}(\cos\varphi)\cos (m\lambda) \quad (m \ge 0),\\
          &=& P_n^{|m|}(\cos\varphi)\sin (|m|\lambda) \quad (m < 0),
    \end{eqnarray}
    %
    これを\Deqref{渦度方程式}に代入し, 
    $Y_{n'}^{m'}(\lambda,\varphi)$ をかけて全球面で積分することにより, 
    球面調和函数の直交関係から以下の式が得られる. 
    %
    \begin{displaymath}
      \DD{\tilde{\zeta}_n^m}{t} 
      = - \Dinv{a^2}[J(\psi,\zeta)]_n^m 
         - \frac{2\Omega}{a^2}\left[\DP{\psi}{\lambda}\right]_n^m
         + \left[(-1)^{p+1}\nu_{2p}\left(\Dlapla + \frac{2}{a^2}\right)^p
         \zeta \right]_n^m
    \end{displaymath}
    %
    ただし $[\ldots]_n^m$ は球面調和函数の $n,m$ 成分であることを表している. 
    %
    散逸項はラプラシアンが球面調和函数の固有演算子であることから
    \begin{eqnarray*}
      \left[(-1)^{p+1}\nu_{2p}\left(\Dlapla + \frac{2}{a^2}\right)^p
             \zeta(\lambda,\varphi,t) \right]_n^m 
          &=& (-1)^{p+1}\nu_{2p}\left[\frac{-n(n+1)+ 2}{a^2}\right]^p
             \tilde{\zeta}_n^m(t)\\
          &=& -\nu^*(n,p)\tilde{\zeta}_n^m(t), 
    \end{eqnarray*}
    % 
    となる. ここで, 波数毎の超粘性係数 $\nu^*(n,p)$ を
    \begin{equation}
        \nu^*(n,p)\equiv (-1)^p\nu_{2p}
                               \left[\frac{-n(n+1)+2}{a^2}\right]^p
                    = \nu_{2p}\left[\frac{n(n+1)-2}{a^2}\right]^p
    \end{equation}
    %
    と定義した($\nu^*(n,p)$ は 0 以上であることに注意). 
    したがって, 
    \begin{equation}
      \Deqlab{スペクトル渦度方程式}
      \DD{\tilde{\zeta}_n^m}{t} 
      = - \Dinv{a^2}[J(\psi,\zeta)]_n^m 
        - \frac{2\Omega}{a^2}\left[\DP{\psi}{\lambda}\right]_n^m
        - \nu^*(n,p)\tilde{\zeta}_n^m
    \end{equation}


  \subsection{時間積分}

    \subsubsection{euler スキーム(1 次精度)}

      euler スキームの場合は時間変化項をすべて現在時間で評価して, 
      %
      \begin{eqnarray*}
        \frac{\tilde{\zeta}(t+\Delta t)_n^m-\tilde{\zeta}(t)_n^m}{\Delta t}
         &=&
           - \Dinv{a^2}[J(\psi(\lambda,\varphi,t),\zeta(\lambda,\varphi,t))]_n^m
           - \frac{2\Omega}{a^2}\left[\DP{\psi(\lambda,\varphi,t)}{\lambda}\right]_n^m \\
        && - \nu^*(n,p)\tilde{\zeta}_n^m(t)
      \end{eqnarray*}
      %
      すなわち, 
      % 
      \begin{eqnarray}
        \tilde{\zeta}(t+\Delta t)_n^m 
            =\tilde{\zeta}(t)_n^m 
             &+&\Delta t \times \left\{
                - \Dinv{a^2}[J(\psi(\lambda,\varphi,t),\zeta(\lambda,\varphi,t))]_n^m
                - \frac{2\Omega}{a^2}\left[\DP{\psi(\lambda,\varphi,t)}{\lambda}\right]_n^m
                \right. \nonumber\\ & &\left. \hspace{20mm}
                - \nu^*(n,p)\tilde{\zeta}_n^m(t)
             \right\}.
      \end{eqnarray}

    \subsubsection{Crank-Nicolson and Adams-Bashforth スキーム(2 次精度)}

      散逸項に Clank-Nicolson スキーム, それ以外の項に
      2 次の Adams-Bashforth スキームを適用する場合には, 
      時間積分を 2 段階に分けて
      %
      \begin{eqnarray*}
        \frac{\tilde{\zeta}^*(t+\Delta t)_n^m-\tilde{\zeta}(t)_n^m}{\Delta t}
         &=&  \frac{3}{2}\left[\left(\DD{\zeta}{t}\right)_{nd}(t)\right]_n^m
             -\frac{1}{2}\left[\left(\DD{\zeta}{t}\right)_{nd}(t-\Delta t)\right]_n^m,\\
       \frac{\tilde{\zeta}(t+\Delta t)_n^m-\tilde{\zeta}^*(t+\Delta t)_n^m}{\Delta t}
         & & = \frac{1}{2}\left\{
                 \left[\left(\DD{\zeta}{t}\right)_{d}(t+\Delta t)\right]_n^m
                +\left[\left(\DD{\zeta}{t}\right)_{d}(t)\right]_n^m
               \right\}, \\
      \end{eqnarray*}
      %
      ただし, $\Ddsty\left(\DD{\zeta}{t}\right)_{nd}, 
      \left(\DD{\zeta}{t}\right)_{d}$ はそれぞれ時間変化の
      非散逸項と散逸項を表している:
      %
      \begin{eqnarray*}
        \left[\left(\DD{\zeta}{t}\right)_{nd}(t)\right]_n^m &\equiv&
           - \Dinv{a^2}[J(\psi(\lambda,\varphi,t),\zeta(\lambda,\varphi,t))]_n^m
           - \frac{2\Omega}{a^2}\left[\DP{\psi(\lambda,\varphi,t)}{\lambda}\right]_n^m,\\
        \left[\left(\DD{\zeta}{t}\right)_{d}(t)\right]_n^m
           &\equiv& - \nu^*(n,p)\tilde{\zeta}_n^m(t)
      \end{eqnarray*}
      %
      すると, 
      %
      \begin{eqnarray*}
        \frac{\tilde{\zeta}^*(t+\Delta t)_n^m-\tilde{\zeta}(t)_n^m}{\Delta t}
         &=&  \frac{3}{2}\left[\left(\DD{\zeta}{t}\right)_{nd}(t)\right]_n^m
             -\frac{1}{2}\left[\left(\DD{\zeta}{t}\right)_{nd}(t-\Delta t)\right]_n^m\\
        \frac{\tilde{\zeta}(t+\Delta t)_n^m-\tilde{\zeta}^*(t+\Delta t)_n^m}{\Delta t}
         &=& \frac{-\nu^*(n,p)\tilde{\zeta}_n^m(t+\Delta t)
                     -\nu^*(n,p)\tilde{\zeta}_n^m(t)}{2},
      \end{eqnarray*}
      %
      したがって
      %
      \begin{eqnarray}
        \tilde{\zeta}^*(t+\Delta t)_n^m &=& \tilde{\zeta}(t)_n^m
           + \Delta t \times \left\{
                  \frac{3}{2}\left[\left(\DD{\zeta}{t}\right)_{nd}(t)\right]_n^m
                  -\frac{1}{2}\left[\left(\DD{\zeta}{t}\right)_{nd}(t-\Delta t)
             \right]_n^m \right\}, \\
  % %
        \tilde{\zeta}(t+\Delta t)_n^m 
             &=& \frac{1 - \nu^*(n,p)\Delta t/2}{1 + \nu^*(n,p)\Delta t/ 2}
                   \times\tilde{\zeta}^*(t+\Delta t)_n^m.
      \end{eqnarray}

    \subsubsection{演算子分割処理法(粘性項)と 4 次精度 Runge-Kutta スキーム}

      変数変換することによって粘性項を解析的に評価し, 
      その他の項を高次の数値積分スキームを適用する
      演算子分割法をとる場合を以下に示す. 
      渦度のスペクトル成分について
      %
      \begin{equation}
        \tilde{\zeta}(t)_n^m = \hat{\zeta}(t)_n^m e^{-\nu^*(n,p)t}
      \end{equation}
      %
      と変数変換すれば, 
      %
      \begin{displaymath}
        \DD{\tilde{\zeta}_n^m}{t} 
         =   \DD{\hat{\zeta}(t)_n^m}{t} e^{-\nu^*(n,p)t} 
           - \nu^*(n,p)\hat{\zeta}(t)_n^m e^{\nu^*(n,p)t}
         =   \DD{\hat{\zeta}(t)_n^m}{t} e^{-\nu^*(n,p)t} 
           - \nu^*(n,p)\tilde{\zeta}(t)_n^m
      \end{displaymath}
      %
      であるから, \Deqref{スペクトル渦度方程式}は
      \begin{equation}
        \DD{\hat{\zeta}(t)_n^m}{t} e^{-\nu^*(n,p)t} 
        = - \Dinv{a^2}[J(\psi,\zeta)]_n^m 
          - \frac{2\Omega}{a^2}\left[\DP{\psi}{\lambda}\right]_n^m.
      \end{equation}
      %
      \begin{equation}
        \DD{\hat{\zeta}(t)_n^m}{t} 
        = e^{\nu^*(n,p)t}\left[ 
            - \Dinv{a^2}[J(\psi,\zeta)]_n^m 
            - \frac{2\Omega}{a^2}\left[\DP{\psi}{\lambda}\right]_n^m 
          \right].
      \end{equation}
      %
      この左辺の時間積分を Runge-Kutta スキームで計算する. 
      すなわち, 
      \begin{eqnarray*}
        \tilde{\zeta}(t+\Delta t)_n^m 
          &=& \hat{\zeta}(t+\Delta t)_n^m e^{-\nu^*(n,p)\Delta t},\\
        \hat{\zeta}(t+\Delta t)_n^m 
          &=& \tilde{\zeta}(t)_n^m
           + \Delta t \times \left(  \frac{k_1}{6} +\frac{k_2}{3} 
                                   +\frac{k_3}{3} +\frac{k_4}{6} \right), \\
        k_1 &=& f[\zeta(t),\psi(t),0]_n^m,\\
        k_2 &=& f[\zeta(t_1),\psi(t_1),\Delta t/2]_n^m\\
        k_3 &=& f[\zeta(t_2),\psi(t_2),\Delta t/2]_n^m\\
        k_4 &=& f[\zeta(t_3),\psi(t_3),\Delta t]_n^m\\
        f[\zeta,\psi,t]_n^m &=& 
          e^{\nu^*(n,p)t}\left[ 
            - \Dinv{a^2}[J(\psi,\zeta)]_n^m 
            - \frac{2\Omega}{a^2}\left[\DP{\psi}{\lambda}\right]_n^m 
          \right],\\
        \tilde{\zeta}_n^m(t_1) &=& e^{-\nu^*(n,p)\Delta t/2}
            [\tilde{\zeta}(t)_n^m + k_1 \Delta t/2],\\
        \tilde{\zeta}_n^m(t_2) &=& e^{-\nu^*(n,p)\Delta t/2}
            [\tilde{\zeta}(t)_n^m + k_2 \Delta t/2],\\
        \tilde{\zeta}_n^m(t_3) &=& e^{-\nu^*(n,p)\Delta t}
            [\tilde{\zeta}(t)_n^m + k_3 \Delta t].
      \end{eqnarray*}
      %
      この定式化において 
      $\tilde{\zeta}(t)_n^m = \hat{\zeta}(t)_n^m$ であることを
      用いていることに注意されたい. 


    \subsubsection{演算子分割処理法(線形項)と 4 次精度 Runge-Kutta スキーム}

      粘性項だけでなく, 全ての線形項を演算子分割処理方によって解析的に扱い, 
      残りの非線形項のみ数値積分して解くための定式化を記す. 
      
      今, $m$ を正の整数として, \Deqref{スペクトル渦度方程式}の
      $(n,m)$ 成分と $(n,-m)$ 成分の式を書き下すと
      %
      \begin{eqnarray*}
        \DD{\tilde{\zeta}_n^m}{t} 
        &=& - \Dinv{a^2}[J(\psi,\zeta)]_n^m 
         - \frac{2\Omega}{a^2}\left[\DP{\psi}{\lambda}\right]_n^m
         - \nu^*(n,p)\tilde{\zeta}_n^m, \\
        \DD{\tilde{\zeta}_n^{-m}}{t} 
        &=& - \Dinv{a^2}[J(\psi,\zeta)]_n^{-m}
         - \frac{2\Omega}{a^2}\left[\DP{\psi}{\lambda}\right]_n^{-m}
         - \nu^*(n,p)\tilde{\zeta}_n^{-m}
      \end{eqnarray*}
      %
      ここで $\Dlapla\psi=\zeta$ より
      $\tilde{\psi}_n^m = -a^2\tilde{\zeta}_n^m /n(n+1)$ であることを
      用いると, 
      \begin{displaymath}
        \frac{2\Omega}{a^2}\left[\DP{\psi}{\lambda}\right]_n^m
           = \frac{2\Omega m}{n(n+1)} \tilde{\zeta}_n^{-m},\quad
        \frac{2\Omega}{a^2}\left[\DP{\psi}{\lambda}\right]_n^{-m}
           = -\frac{2\Omega m}{ n(n+1)} \tilde{\zeta}_n^m,
      \end{displaymath}
      %
      となるので
      %
      \begin{eqnarray*}
        \DD{\tilde{\zeta}_n^m}{t} 
        &=& - \Dinv{a^2}[J(\psi,\zeta)]_n^m 
         - \frac{2\Omega m}{n(n+1)} \tilde{\zeta}_n^{-m}
         - \nu^*(n,p)\tilde{\zeta}_n^m, \\
        \DD{\tilde{\zeta}_n^{-m}}{t} 
        &=& - \Dinv{a^2}[J(\psi,\zeta)]_n^{-m}
         + \frac{2\Omega m}{n(n+1)} \tilde{\zeta}_n^m
         - \nu^*(n,p)\tilde{\zeta}_n^{-m}
      \end{eqnarray*}
      %
      ロスビー波の位相速度 
      $\Ddsty\omega_R(n,m)\equiv - \frac{2\Omega |m|}{n(n+1)}$
      を用いて書き直すと, 
      %
      \begin{eqnarray*}
        \DD{\tilde{\zeta}_n^m}{t} 
        &=& - \Dinv{a^2}[J(\psi,\zeta)]_n^m 
         + \omega_R(n,m) \tilde{\zeta}_n^{-m}
         - \nu^*(n,p)\tilde{\zeta}_n^m, \\
        \DD{\tilde{\zeta}_n^{-m}}{t} 
        &=& - \Dinv{a^2}[J(\psi,\zeta)]_n^{-m}
         - \omega_R(n,m) \tilde{\zeta}_n^m 
         - \nu^*(n,p)\tilde{\zeta}_n^{-m}
      \end{eqnarray*}
      %
      解きやすくするために $\tilde{\zeta}_n^m$ と
      $\tilde{\zeta}_n^{-m}$ を複素変数 $Z_n^m$ でまとめて扱う. 
      すなわち 
      $ Z_n^m \equiv \tilde{\zeta}_n^m + i \tilde{\zeta}_n^{-m} $
      と置き換えると, 
      \begin{displaymath}
        \DD{Z_n^m}{t} 
        = G(\psi,\zeta)_n^m  -i \omega_R(n,m) Z_n^m - \nu^*(n,p)Z_n^m,
      \end{displaymath}
      %
      ただし $G(\psi,\zeta)_n^m$ は非線形項を複素表示したものであり, \\
      $\Ddsty G(\psi,\zeta)_n^m \equiv 
      -\left[ \Dinv{a^2}[J(\psi,\zeta)]_n^m
             +i\Dinv{a^2}[J(\psi,\zeta)]_n^{-m}\right]$
      である.

      線形項を消去すべく 
      $Z_n^m(t) = \hat{Z}_n^m(t) e^{-[i\omega_R(n,m)+\nu^*(n,p)]t}$
      と変換すれば, 
      \begin{eqnarray*}
        \DD{\hat{Z}_n^m}{t} e^{-[i\omega_R(n,m)+\nu^*(n,p)]t}
        &-&[i\omega_R(n,m)+\nu^*(n,p)]\hat{Z}_n^m(t) e^{-[i\omega_R(n,m)+\nu^*(n,p)]t}\\
        &=& G(\psi,\zeta)_n^m  -[i\omega_R(n,m) Z_n^m + \nu^*(n,p)]Z_n^m,
      \end{eqnarray*}
      したがって
      \begin{displaymath}
        \DD{\hat{Z}_n^m}{t}
        = e^{[i\omega_R(n,m)+\nu^*(n,p)]t} G(\psi,\zeta)_n^m 
      \end{displaymath}
      %
      実部と虚部に分けて 
      $Re[\hat{Z}_n^m]=\xi_n^m, Im[\hat{Z}_n^m]=\xi_n^{-m}$ とおくと
      %
      \begin{eqnarray}
        \DD{\xi_n^m}{t}
        &=&   e^{\nu^*(n,p)]t}\cos(\omega_R(n,m) t) Re[G(\psi,\zeta)_n^m]
            - e^{\nu^*(n,p)]t}\sin(\omega_R(n,m) t)Im[G(\psi,\zeta)_n^m]\nonumber\\
        &=&    e^{\nu^*(n,p)]t}\cos(\omega_R(n,m) t) 
                 \left[-\Dinv{a^2}J(\psi,\zeta)\right]_n^m
            - e^{\nu^*(n,p)]t}\sin(\omega_R(n,m) t)
                \left[-\Dinv{a^2}J(\psi,\zeta)\right]_n^{-m},\\
        \DD{\xi_n^{-m}}{t}
        &=&   e^{\nu^*(n,p)]t}\cos(\omega_R(n,m) t) Im[G(\psi,\zeta)_n^m]
            + e^{\nu^*(n,p)]t}\sin(\omega_R(n,m) t) Re[G(\psi,\zeta)_n^m]\nonumber\\
        &=&   e^{\nu^*(n,p)]t}\cos(\omega_R(n,m) t) 
                 \left[-\Dinv{a^2}J(\psi,\zeta)\right]_n^{-m}
            + e^{\nu^*(n,p)]t}\sin(\omega_R(n,m) t)
                \left[-\Dinv{a^2}J(\psi,\zeta)\right]_n^m,\\
      \end{eqnarray}
      %
      時間積分のための $\xi_n^m$, $\xi_n^{-m}$ の各積分ステップでの初期値は
      $Z_n^m$ と $\hat{Z}_n^m$ の関形において $t=0$ を適用すればよく, 
      $\xi_n^m = \tilde{\zeta}_n^m$, $\xi_n^{-m} = \tilde{\zeta}_n^{-m}$
      となる. 上の式を用いて数値積分し, あたらしい時間の
      $\xi_n^m$, $\xi_n^{-m}$ が求まれば, 
      対応して $\tilde{\zeta}_n^m$ がつぎのように求められる. 
      \begin{eqnarray}
        \tilde{\zeta}_n^m &=& Re[Z_n^m]
             = Re[\hat{Z}_n^m e^{-[i\omega_R(n,m)+\nu^*(n,p)]t}]\nonumber\\
            &=&   Re[\hat{Z}_n^m] e^{-\nu^*(n,p)t}\cos(\omega_R(n,m) t)
                + Im[\hat{Z}_n^m] e^{-\nu^*(n,p)t}\sin(\omega_R(n,m) t)\nonumber\\
            &=&   \xi_n^m e^{-\nu^*(n,p)t}\cos(\omega_R(n,m) t)
                + \xi_n^{-m} e^{-\nu^*(n,p)t}\sin(\omega_R(n,m) t),\\
        \tilde{\zeta}_n^{-m} &=& Im[Z_n^m]
             = Im[\hat{Z}_n^m e^{-[i\omega_R(n,m)+\nu^*(n,p)]t}]\nonumber\\
            &=&   Im[\hat{Z}_n^m] e^{-\nu^*(n,p)t}\cos(\omega_R(n,m) t)
                - Re[\hat{Z}_n^m] e^{-\nu^*(n,p)t}\sin(\omega_R(n,m) t)\nonumber\\
            &=&   \xi_n^{-m} e^{-\nu^*(n,p)t}\cos(\omega_R(n,m) t)
                - \xi_n^m e^{-\nu^*(n,p)t}\sin(\omega_R(n,m) t).
      \end{eqnarray}




\section{実験設定}

  \subsection{初期条件}

  \subsection{パラメター}

  \subsection{時間積分}

    数値計算が意味を持つための時間積分の時間ステップの目安を見積もる. 
    各時間変化項の最大値は
    %
    \begin{eqnarray*}
      \Dinv{a^2}J(\psi,\zeta) 
           & \sim & \frac{UL\cdot U/L}{L^2}
             \sim  \frac{U^2}{L^2},\\
      \frac{2\Omega}{a^2}\DP{\psi}{\lambda} 
           & \sim & \frac{\Omega}{a^2}\cdot Ua
             \sim   \frac{\Omega U}{a},\\
      (-1)^{p+1}\nu_{2p}\left(\Dlapla + \frac{2}{a^2}\right)^p \zeta
           & \sim & \nu_{2p}\left(\frac{N_{max}(N_{max}+1)}{a^2}\right)^p 
                     \cdot\frac{U}{L}
    \end{eqnarray*}
    %
    ここで, $U,L$ は現象の典型的な速度と水平スケールである. 
    一方, 時間変化項は 
    $\Ddsty \DP{\zeta}{t}\sim \frac{U}{(L)(\Delta t)}$
    と見積もられるので, 時間変化が十分に解像できるには
    $\Ddsty \DP{\zeta}{t}$ が各時間変化項の見積もりより大きくなるよう
    時間刻 $\Delta t$ をとらねばならない
    \footnote{どうもすっきりしない論理だ...}.
    したがって, 
    %
    \begin{eqnarray*}
      \frac{U}{L\Delta t}
        & \gg & \max\left[
                        \Dinv{a^2}J(\psi,\zeta), 
                        \frac{2\Omega}{a^2}\DP{\psi}{\lambda},
                        (-1)^{p+1}\nu_{2p}
                            \left(\Dlapla + \frac{2}{a^2}\right)^p \zeta
                     \right]\\
       &\sim &\max\left[
                     \frac{U^2}{L^2},
                     \frac{\Omega U}{a},
                     \nu_{2p}\left(\frac{N_{max}(N_{max}+1)}{a^2}\right)^p 
                         \cdot\frac{U}{L}
                   \right],\\
      \Delta t 
        &\ll& \min\left[
                \frac{L}{U},
                \frac{a}{\Omega L},
                \Dinv{\nu_{2p}\left(\frac{N_{max}(N_{max}+1)}{a^2}\right)^p}
              \right].
    \end{eqnarray*}
    %
    右辺の各項が最小となる水平スケールをそれぞれ選べば, 
    第 1 項目は $L \sim \Delta x$, 第 2 項目は $L \sim a$ 
    となる. したがって,
    \begin{equation}
      \Delta t 
        \ll \min\left[
                \frac{\Delta x}{U},
                \frac{1}{\Omega},
                \Dinv{\nu_{2p}\left(\frac{N_{max}(N_{max}+1)}{a^2}\right)^p}
              \right].
    \end{equation}

\begin{Dreference}
  \item SPMODEL 解説文書「2 次元回転球面上での非発散順圧流体の定式化」
    http://www.gfd-dennou.org/library/spmodel/2d-sphere-w/baro/pub/
\end{Dreference}

\end{document}

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: t
%%% End: 
