% 表題   DCPAM5  力学過程 -- 時間積分(差分)
%
% 履歴
%\Drireki{2010/04/15 高橋芳幸}
%\Drireki{2009/10/11 高橋芳幸}
%\Drireki{2009/07/11 森川靖大}
%\Drireki{2009/02/18 森川靖大}
%\Drireki{2008/06/15 森川靖大}
%\Drireki{1993/03/18 沼口敦・保坂征宏}


\section{離散表現: 時間離散化}


%ここでは時間積分スキームについて記す. 
%
%時間差分スキームは基本的に leap frog である.
%ただし, 拡散項および物理過程の項は後方差分もしくは前方差分とする.
%計算モードを抑えるために時間フィルター(Asselin, 1972)を用いる.
%さらに$\Delta t$ を大きくとるために,
%重力波の項に semi-implicit の手法を適用する(Bourke, 1988).
%
%-----------------------------------------------------


ここでは時間積分スキームについて記す. 

時間差分には, 複数の方法を組み合わせて用いる. 用いる方法の
概要を以下に示す. 
%
\begin{itemize}
\item 力学過程
  \begin{itemize}
    \item 水平拡散およびスポンジ層における減衰項には, 後方差分を用いる. 
    \item その他の項には, leap frog 法と Crank-Nicolson 法を組み合わせた 
          semi-implicit 法
          (Bourke, 1988) を用いる. 
  \end{itemize}
\item 物理過程
  \begin{itemize}
    \item 予報型の物理過程には, 前方差分を用いる. 
    \item 調節型の物理過程は, semi-implicit 法での力学過程積分後に計算された値を
          用いて計算する. 
  \end{itemize}
\item 時間フィルタ
  \begin{itemize}
    \item 力学過程, 物理過程のすべての計算後に, 力学過程で用いている
          leap frog 法を起源とする計算モード抑制のための, 
          Asselin (1972) による時間フィルター, または Williams (2009) に
          よる時間フィルタを適応する. 
  \end{itemize}
\end{itemize}

この方法は, 予報変数を ${\cal A}$ と表すと, 以下の式で表現される. 
%
\begin{align}
  \Deqlab{symbolic_timeintegration}
  \frac{ \hat{\cal A}^{t+\Delta t} - \bar{\cal A}^{t-\Delta t} }{ 2 \Delta t }
%    =   \dot{\cal A}_{dyn,G}
%          \left( \frac{ \bar{\cal A}^{t-\Delta t} + \hat{\cal A}^{t+\Delta t} }{2} \right)
    =   \frac{1}{2}
          \left\{
            \dot{\cal A}_{dyn,G} \left( \bar{\cal A}^{t-\Delta t} \right)
          + \dot{\cal A}_{dyn,G} \left( \hat{\cal A}^{t+\Delta t} \right)
          \right\}
      + \dot{\cal A}_{dyn,NG} 
          \left( {\cal A}^{t} \right)
  \nonumber \\
      + \dot{\cal A}_{dyn,dis }\left( \hat{\cal A}^{t+\Delta t} \right)
      + \dot{\cal A}_{phy,pred}\left( \bar{\cal A}^{t-\Delta t} \right) , 
\end{align}
%
\begin{align}
  {\cal A}^{t+\Delta t} 
    =  \hat{\cal A}^{t+\Delta t}
    + 2 \Delta t 
      \dot{\cal A}_{fric}\left( \hat{\cal A}^{t+\Delta t} \right)
    + 2 \Delta t 
      \dot{\cal A}_{phy,adj}\left( \hat{\cal A}^{t+\Delta t} \right) , 
\end{align}
%
これに続いて適用する Asselin (1972) による時間フィルタは下のように
表される:
%
\begin{align}
  \bar{\cal A}^{t}
    = {\cal A}^{t}
    +  \epsilon_f 
        \left( \bar{\cal A}^{t-\Delta t} - 2 {\cal A}^{t} + {\cal A}^{t+\Delta t} \right). 
\end{align}
%
また, Williams (2009) による時間フィルタは下のように
表される:
\footnote{未だ書いていない (yot, 2012/12/24).}
%
%\begin{align}
%  \bar{\cal A}^{t}
%    = {\cal A}^{t}
%    +  \epsilon_f 
%        \left( \bar{\cal A}^{t-\Delta t} - 2 {\cal A}^{t} + {\cal A}^{t+\Delta t} \right). 
%\end{align}
\\
%
ここで, $ \dot{\cal A}_{dyn,G} $, $ \dot{\cal A}_{dyn,NG} $ はそれぞれ, 
力学過程において semi-implicit 法で分離された重力波項 (線型項) と非重力波項 (非線型項), 
$ \dot{\cal A}_{dyn,dis} $ は水平拡散とスポンジ層における減衰項,
$ \dot{\cal A}_{phy,pred} $ は予報型の物理過程項である.
$ \dot{\cal A}_{fric} $, $ \dot{\cal A}_{phy,adj} $ は, それぞれ
摩擦熱による加熱項および調節型の物理過程項である. 
$\epsilon_f$ は時間フィルタの係数であり, \Dmodel での標準値は 0.05 としている. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\subsection{semi-implicit 時間積分}
\subsection{力学過程の方程式系の時間差分式}

まず, semi-implicit 法を用いるために, 方程式系を $T=\overline{T}_k$ である
静止場に基づいて線形重力波項とそれ以外の項に分離する. 
鉛直方向のベクトル表現 $\Dvect{A}=\{ A_{k} \}$,
および行列表現$\underline{A}=\{ A_{kl} \}$ を用いると, 連続の式, 発散方程式, 
熱力学の式は, 
%
\begin{align}
   \DP{\tilde{\pi}^{m}_{n}}{t} =
          \left( \DP{\tilde{\pi}^{m}_{n}}{t} \right)^{\rm NG}
     - \Dvect{C} \cdot \tilde{\Dvect{D}}^{m}_{n}  ,
\end{align}
%
\begin{align}
  \DP{\tilde{\Dvect{D}}^{m}_{n}}{t} =
          \left( \DP{\tilde{\Dvect{D}}^{m}_{n}}{t} \right)^{\rm NG}  
          - \left( - \frac{n(n+1)}{a^2} \right)
            (   \tilde{\Dvect{\Phi}}^{m}_{s,n} 
              + \underline{W} \tilde{\Dvect{T}}^{m}_{n}
              + \Dvect{G} \tilde{\pi}^{m}_{n} )  
          + \underline{ \tilde{\cal D}_M }_{n}^{m} \tilde{\Dvect{D}}^{m}_{n} ,
\end{align}
%
\begin{align}
  \DP{\tilde{\Dvect{T}}^{m}_{n}}{t} 
     & =   \left( \DP{\tilde{\Dvect{T}}^{m}_{n}}{t} \right)^{\rm NG}  
         - \underline{h} \tilde{\Dvect{D}}^{m}_{n}
         + \underline{ \tilde{\cal D}_{H} }_{n}^{m} \tilde{\Dvect{T}}^{m}_{n}
\end{align}
%
となる
    \footnote{
       念のため注記しておくと,
       $\tilde{\Dvect{\Phi}}^{m}_{s,n} =
         \left( \tilde{\Phi}^{m}_{s,n}, \tilde{\Phi}^{m}_{s,n}, \cdots, \tilde{\Phi}^{m}_{s,n}\right)$ である. 
    }. 
$\widetilde{( \hspace{0.5cm} )}^{m}_{n}$や
$\widetilde{[ \hspace{0.5cm} ]}^{m}_{n}$
といった表記については
%$\widetilde{( \partial / \partial \lambda )}^{m}_{n}$,
%$\widetilde{( \partial / \partial \mu )}^{m}_{n}$
\ref{sec:crd:水平スペクトル}節の
\Deqref{crd:スペクトル係数への変換}, 
\Deqref{crd:λ微分のスペクトル表現}, \Deqref{crd:μ微分のスペクトル表現}
を参照のこと. 
ここで, 
添字 NG の付いた項は, 非重力波項であり, 
以下のように表される.
%
\begin{align}
  \Deqlab{Z項}
  \left( \DP{\tilde{\pi}^{m}_{n}}{t} \right)^{\rm NG}
   =   \tilde{Z}^{m}_{n},
\end{align}
%
\begin{align}
  \Deqlab{発散変化の非重力波項}
 \begin{split}
   \left( \DP{\tilde{D}^{m}_{k,n}}{t} \right)^{\rm NG}
   &=
            \Dinv{a}
            \left(
                \widetilde{
                  \left[
                    \Dinv{1 - \mu^2} \DP{\, {U_A}_{,ijk}}{\lambda}
                  \right]^{m}_{n}
                }
              + \widetilde{
                  \left[
                    \DP{\, {V_A}_{,ijk}}{\mu}
                  \right]^{m}_{n}
                }
            \right) \\
   & \qquad
          - \left( - \frac{n(n+1)}{a^2} \right)
              \widetilde{
                \left[   (\mbox{\sl KE})_{k} 
                       + \sum_{l=1}^{K} W_{kl} ( T_{v,l}-T_{l} )
                \right]
              }^{m}_{n},
 \end{split}
\end{align}
%
\begin{align}
  \Deqlab{温度変化の非重力波項}
 \begin{split}
   \left( \DP{\tilde{T_{k}}^{m}_{,n}}{t} \right)^{\rm NG} 
     & =
          - \Dinv{a}
            \left(
                \widetilde{
                  \left[
                    \Dinv{1 - \mu^2} \DP{{U}_{ijk} {T'}_{ijk}}{\lambda}
                  \right]^{m}_{n}
                } 
              + \widetilde{
                  \left[
                    \DP{{V}_{ijk} {T'}_{ijk}}{\mu}
                  \right]^{m}_{n}
                } 
            \right)
  \\ & \qquad
          + \widetilde{
              \left[
                  H_{ijk}
%                + {\cal D}^{\prime}_{ijk}(\Dvect{v})
              \right]^{m}_{n}
            }.
 \end{split}
\end{align}
%
各項は以下の通りである. 簡単化のため経度, 緯度方向添字 $i,j$ の
表記を省略する. 
%
\begin{align}
 Z & = - \sum_{k=1}^{K} \Dvect{v}_{k} \cdot \nabla \pi  
         \Delta  \sigma_{k}, \\
 \begin{split}
  H_k & =     T_{k}^{\prime} D_{k}  \\
     & \quad
           - \frac{1}{\Delta \sigma_{k}} 
             \left[
                 \dot{\sigma}_{k-1/2} \left(  \hat{T^{\prime}}_{k-1/2} 
                                             - T^{\prime}_{k}   \right)
               + \dot{\sigma}_{k+1/2} \left(   T^{\prime}_{k}  
                                             - \hat{T^{\prime}}_{k+1/2} \right) \right]
                \\
     & \quad
           - \frac{1}{\Delta \sigma_{k}} 
             \left[
                 \dot{\sigma}^{\rm NG}_{k-1/2} \left( \hat{\overline{T}}_{k-1/2} 
                                         - \overline{T}_{k}   \right)
               + \dot{\sigma}^{\rm NG}_{k+1/2} \left( \overline{T}_{k}  
                                         - \hat{\overline{T}}_{k+1/2} \right) \right]
                \\
     & \quad
           + \hat{\kappa}_{k} T_{v,k} \Dvect{v}_{k} \cdot \nabla \pi
                \\
     & \quad
           - \frac{\alpha_{k}}{\Delta \sigma_{k} }
               \left[
                   T_{v,k}
                     \sum_{l=k}^{K} \Dvect{v}_{l} \cdot \nabla \pi 
                     \Delta \sigma_{l}
                 + T'_{v,k} \sum_{l=k}^{K} D_l  \Delta \sigma_{l}
               \right]
                \\
     & \quad
           - \frac{\beta_{k}}{\Delta \sigma_{k} } 
               \left[
                   T_{v,k}
                     \sum_{l=k+1}^{K} \Dvect{v}_{l} \cdot \nabla \pi 
                     \Delta \sigma_{l}
                 + T'_{v,k} \sum_{l=k+1}^{K} D_l  \Delta \sigma_{l}
               \right]
       \qquad (k=1, \cdots, K-1), \\
%
  H_K & =     T_{K}^{\prime} D_{K}  \\
     & \quad
           - \frac{1}{\Delta \sigma_{K}} 
             \left[
                 \dot{\sigma}_{K-1/2} \left(   \hat{T^{\prime}}_{K-1/2} 
                                             - T^{\prime}_{K}   \right)
               + \dot{\sigma}_{K+1/2} \left(   T^{\prime}_{K}  
                                             - \hat{T^{\prime}}_{K+1/2} \right) \right]
                \\
     & \quad
           - \frac{1}{\Delta \sigma_{K}} 
             \left[
                 \dot{\sigma}^{\rm NG}_{K-1/2} \left( \hat{\overline{T}}_{K-1/2} 
                                         - \overline{T}_{K}   \right)
               + \dot{\sigma}^{\rm NG}_{K+1/2} \left( \overline{T}_{K}  
                                         - \hat{\overline{T}}_{K+1/2} \right) \right]
                \\
     & \quad
           + \hat{\kappa}_{K} T_{v,K} \Dvect{v}_{K} \cdot \nabla \pi
                \\
     & \quad
           - \frac{\alpha_{K}}{\Delta \sigma_{K} }
               \left[
                   T_{v,K} \Dvect{v}_{K} \cdot \nabla \pi \Delta \sigma_{K}
                 + T'_{v,K} D_K  \Delta \sigma_{K}
               \right],
 \end{split}
\end{align}
%
\begin{align}
 \begin{split}
  \dot{\sigma}^{\rm NG}_{k-1/2}
  &= - \sigma_{k-1/2} \left( \DP{\pi}{t} \right)^{\rm NG}
     - \sum_{l=k}^{K} \Dvect{v}_{l} \cdot \nabla \pi \Delta  \sigma_{l} \\
  &=   \sigma_{k-1/2} \sum_{k=1}^{K} \Dvect{v}_{k} \cdot \nabla \pi \Delta  \sigma_{k}
     - \sum_{l=k}^{K} \Dvect{v}_{l} \cdot \nabla \pi \Delta  \sigma_{l}, 
 \end{split}
\end{align}
%
\begin{align}
 \hat{T}_{k-1/2}' = \left\{
 \begin{array}{ll}
   0 , & \text{($k = 1$)} \\
   \hat{T}_{k-1/2} - \hat{\overline{T}}_{k-1/2}, & \text{($k = 2, \cdots, K$)} \\
   0 , & \text{($k = K+1$)}
 \end{array} \right. 
\end{align}
%
\begin{align}
 \hat{\overline{T}}_{k-1/2} = \left\{
 \begin{array}{ll}
   0 , & \text{($k = 1$)} \\
   a_k \overline{T}_k + b_{k-1} \overline{T}_{k-1}, & \text{($k = 2, \cdots, K$)} \\
   0 . & \text{($k = K+1$)}
 \end{array} \right. 
\end{align}
%
また, 重力波項のベクトルおよび行列は以下のとおりである. 
%
\begin{align}
  \Deqlab{係数C}
  C_{k} &= \Delta \sigma_{k} ,
    \\
  W_{kl} &= C_{p} \alpha_{l} \delta_{k \geq l}
          + C_{p} \beta_{l} \delta_{k-1 \geq l} ,
    \\
  G_{k} &= \hat{\kappa}_{k} C_{p} \overline{T}_{k} ,
    \\
  \underline{h} &= \underline{Q}\underline{S} - \underline{R} ,
    \\
  Q_{kl} &=
           \frac{1}{\Delta \sigma_{k}} 
             ( \hat{\overline{T}}_{k-1/2} - \overline{T}_{k} ) \delta_{k=l} 
         + \frac{1}{\Delta \sigma_{k}} 
             ( \overline{T}_{k} - \hat{\overline{T}}_{k+1/2}  )
 \delta_{k+1=l}  ,
    \\
  S_{kl} &=  \sigma_{k-1/2} \Delta \sigma_{l} 
           - \Delta \sigma_{l} \delta_{k \leq l }  ,
    \\
  \Deqlab{係数R}
  R_{kl} &= - \left(  \frac{ \alpha_{k} }{ \Delta \sigma_{k} } 
                      \Delta \sigma_{l} \delta_{k \leq l} 
                    + \frac{ \beta_{k} }{ \Delta \sigma_{k} } 
                      \Delta \sigma_{l} \delta_{k+1 \leq l}  
              \right) \overline{T}_{k} ,
    \\
  (\tilde{ {\cal D}_{M,kl} } )_n^m &= 
                         - K_{HD} \left[ 
                             \left( \frac{-n(n+1)}{a^{2}} \right)^{N_D/2}
                             - \left( \frac{2}{a^2} \right)^{N_D/2}
                             \right]  
                             \delta_{k=l} \nonumber \\
                         & \hspace{1cm} 
                           - \gamma_{M,0,n}^m \left( \frac{\sigma_k}{\sigma_K} \right)^{N_{SL}} \delta_{k=l} \delta_{k \ge k_{SLlim}} .
    \\
  (\tilde{ {\cal D}_{H,kl} } )_n^m &= 
                         - K_{HD} 
                             \left( \frac{-n(n+1)}{a^{2}} \right)^{N_D/2}
                             \delta_{k=l} \nonumber \\
                         & \hspace{1cm} 
                           - \gamma_{H,0,n}^m \left( \frac{\sigma_k}{\sigma_K} \right)^{N_{SL}} \delta_{k=l} \delta_{k \ge k_{SLlim}} .
\end{align}
%
$\delta_{k \leq l}$ は, 
$ k \leq l$ が成り立つとき 1, そうでないとき 0 となる関数である.

なお, 渦度方程式には線型重力波項がないため, ここでは示さない. 
\footnote{
  ここは本当は方程式を書くべきだろう. 後で書く. (YOT, 2009/10/11)
}

これらの方程式に, 
%
\begin{itemize}
  \item 水平拡散とスポンジ層における減衰項には後退差分
  \item その他の項には, leap frog 法と中心差分を組み合わせた semi-implicit 法
\end{itemize}
%
を適応すると, 
%
\begin{align}
  \Deqlab{semi-imp pi}
  \delta_{t} \tilde{\pi}^{m}_{n} = 
          \left( \DP{\tilde{\pi}^{m}_{n}}{t} \right)^{\rm NG}  
     - \Dvect{C} \cdot \overline{ \tilde{\Dvect{D}}^{m}_{n} }^{t} ,
\end{align}
%
\begin{align}
  \Deqlab{semi-imp D}
  \delta_{t} \tilde{\Dvect{D}}^{m}_{n} & =
          \left( \DP{\tilde{\Dvect{D}}^{m}_{n}}{t} \right)^{\rm NG}  
          - \left( - \frac{n(n+1)}{a^2} \right)
            (    \tilde{\Dvect{\Phi}}^{m}_{s,n} 
               + \underline{W} 
                 \overline{ \tilde{\Dvect{T}}^{m}_{n} }^{t}
               + \Dvect{G} \overline{\tilde{\pi}^{m}_{n}}^{t}
            ) 
%%%%%%%%%%%%%%%%%%%%%% to be deleted
%            ) \nonumber
%  \\ & \hspace{6.5em} 
%          + \underline{ \tilde{\cal D}_{M} }_{n}^{m} ( \tilde{\Dvect{D}}^{m,t-\Delta t}_{n}
%                         + 2 \Delta t \delta_{t} \tilde{\Dvect{D}}^{m}_{n} ) , 
%%%%%%%%%%%%%%%%%%%%%% to be deleted
          + \underline{ \tilde{\cal D}_{M} }_{n}^{m} \tilde{\Dvect{D}}^{m,t+\Delta t}_{n} ,
\end{align}
%
\begin{align}
  \Deqlab{semi-imp T}
  \delta_{t} \tilde{\Dvect{T}}^{m}_{n} =
        \left( \DP{\tilde{\Dvect{T}}^{m}_{n}}{t} \right)^{\rm NG}  
         - \underline{h} \overline{ \tilde{\Dvect{D}}^{m}_{n} }^{t} 
%%%%%%%%%%%%%%%%%%%%%% to be deleted
%         + \underline{ \tilde{\cal D}_{H} }_{n}^{m} 
%             ( 
%                 \tilde{\Dvect{T}}^{m,t-\Delta t}_{n} 
%               + 2 \Delta t \delta_{t} \tilde{\Dvect{T}}^{m}_{n} 
%             ) . 
%%%%%%%%%%%%%%%%%%%%%% to be deleted
         + \underline{ \tilde{\cal D}_{H} }_{n}^{m} \tilde{\Dvect{T}}^{m,t+\Delta t}_{n}.
\end{align}
%
となる. ただし, 
%
\begin{align}
  \Deqlab{δπ}
  \delta_{t} {\cal A}
    & \equiv
        \frac{1}{2 \Delta t} 
          \left( {\cal A}^{t+\Delta t} - {\cal A}^{t-\Delta t} \right) ,
    \\
  \overline{\cal A}^{t}
    & \equiv \frac{1}{2} \left(   {\cal A}^{t+\Delta t} 
                                 + {\cal A}^{t-\Delta t} \right)
      =      {\cal A}^{t-\Delta t} + \delta_{t} {\cal A} \Delta t .
\end{align}
である. 

\Deqref{semi-imp pi}, \Deqref{semi-imp D}, \Deqref{semi-imp T}
より, $\overline{\tilde{\Dvect{D}}^{m}_{n}}^{t}$ について整理すると, 
%
\begin{align}
 \Deqlab{semi-imp barD}
 \begin{split}
 &   \Biggl[    
                ( \underline{I}-2\Delta t \underline{ \tilde{\cal D}_{M} }_{n}^{m} )
             - ( \Delta t )^{2}
                 \left( - \frac{n(n+1)}{a^2} \right)
                 \left\{
                    \underline{W} 
                    ( \underline{I} -2\Delta t \underline{ \tilde{\cal D}_{H} }_{n}^{m} )^{-1}
                    \underline{h} 
                  + \Dvect{G} \Dvect{C}^{T} 
                 \right\}
     \Biggr]
      \overline{\tilde{\Dvect{D}}^{m}_{n}}^{t} 
       \\
 & \qquad
     =  ( \underline{I}- \Delta t \underline{ \tilde{\cal D}_{M} }_{n}^{m} )
           \tilde{\Dvect{D}}^{m,t-\Delta t}_{n}
       + \Delta t \left( \DP{\tilde{\Dvect{D}}^{m}_{n}}{t} \right)^{\rm NG}  
  \\
 & \qquad \quad
       -  \Delta t \left( - \frac{n(n+1)}{a^2} \right) 
             \Biggl[
                \tilde{\Dvect{\Phi}}_{s,n}^{m}
             \Biggr.
  \\
 & \qquad  \hspace{20mm}
              + \underline{W} 
                ( \underline{I}-2\Delta t \underline{ \tilde{\cal D}_{H} }_{n}^{m} )^{-1}
                \left\{
                       ( \underline{I}-\Delta t \underline{ \tilde{\cal D}_{H} }_{n}^{m} ) 
                       \tilde{\Dvect{T}}^{m,t-\Delta t}_{n}
                     + \Delta t \left( \DP{\tilde{\Dvect{T}}^{m}_{n}}{t} \right)^{\rm NG}
                \right\}
  \\
  & \qquad
                \left.  \hspace*{20mm} 
                       + \Dvect{G} 
                         \left\{ \tilde{\pi}^{m,t-\Delta t}_{n}
                               + \Delta t
                                  \left( \DP{\tilde{\pi}^{m}_{n}}{t} \right)^{\rm NG}  
                         \right\}
                \right]
 \end{split}
\end{align}
%

となる. ここで$\underline{I}$は単位行列, $\Dvect{C}^{T}$は$\Dvect{C}$の
転置ベクトルである. 
\Deqref{semi-imp barD}
を $\overline{\tilde{\Dvect{D}}^{m}_{n}}^{t}$ について解き, 
%
\begin{align}
   \tilde{\Dvect{D}}^{m,t+\Delta t}_{n} =
      2\overline{\tilde{\Dvect{D}}^{m}_{n}}^{t}
    - \tilde{\Dvect{D}}^{m,t-\Delta t}_{n}
\end{align}
%
および, \Deqref{semi-imp pi}, \Deqref{semi-imp T}
により $\hat{\cal A}^{t+\Delta t}$ が求められる.