% 表題   dcpam4  力学過程 -- 時間積分(差分)
%
% 履歴
%\Drireki{1993/03/18 沼口敦・保坂征宏}
%\Drireki{2008/06/15 森川靖大}


\section{時間積分}

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

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{leap frog による時間積分と時間フィルター}

移流項等の時間積分スキームとして leap frog を用いる.
水平拡散項には $2 \Delta t$ の後方差分を使用する.
物理過程の項(${\cal F}_\lambda, {\cal F}_\varphi, Q, S_q$)には 
$2 \Delta t$ の前方差分を使用する.
ただし, 鉛直拡散の時間変化項の計算に関しては後方差分的な取扱いをする.
また, 対流・凝結過程に関しては
一度それらを考慮せずに時間積分して,
その後で場を修正するという形式をとる(adjustment).
摩擦熱の項も補正として扱う.


各予報変数の総体を ${\cal A}$ と表すと,
%
\begin{align}
  \hat{\cal A}^{t+\Delta t} 
    =  \bar{\cal A}^{t-\Delta t}
    + 2 \Delta t 
      \dot{\cal A}_{adv}\left( {\cal A}^{t} \right)
    + 2 \Delta t 
      \dot{\cal A}_{dif}\left( \hat{\cal A}^{t+\Delta t} \right)
    + 2 \Delta t 
      \dot{\cal A}_{phy}\left( \bar{\cal A}^{t-\Delta t} \right) . 
\end{align}
%
$ \dot{\cal A}_{adv} $ は移流項等,
$ \dot{\cal A}_{dif} $ は水平拡散項,
$ \dot{\cal A}_{phy} $ は物理過程等による項である.

$ \hat{\cal A}^{t+\Delta t} $ には, 
摩擦熱($ \dot{\cal A}_{dis} $)および
対流・凝結過程($ \dot{\cal A}_{cnd} $)の補正が加えられ,
$ {\cal A}^{t+\Delta t} $ となる.
%
\begin{align}
  {\cal A}^{t+\Delta t} 
    =  \hat{\cal A}^{t+\Delta t}
    + 2 \Delta t 
      \dot{\cal A}_{dis}\left( \hat{\cal A}^{t+\Delta t} \right)
    + 2 \Delta t 
      \dot{\cal A}_{cnd}\left( \hat{\cal A}^{t+\Delta t} \right) .
\end{align}

leap frog における計算モードの除去のために 
Asselin(1972) の時間フィルターを毎ステップ適用する.
%
すなわち, 
%
\begin{align}
  \bar{\cal A}^{t}
    = ( 1-2 \epsilon_f ) {\cal A}^{t}
    +  \epsilon_f 
        \left( \bar{\cal A}^{t-\Delta t} + {\cal A}^{t+\Delta t} \right)
\end{align}
%
と$\bar{\cal A}$を求める.
$\epsilon_f$ はフィルターの係数であり, 標準的な値である 0.05 を使用する. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{semi-implicit 時間積分}

方程式系において, $T=\overline{T}_k$ であるような静止場を基本場とする
線型重力波項とそれ以外の項に分離する.
後者には非重力波項 (Non Gravity wave terms) として添字NGを付ける. 
鉛直方向のベクトル表現
$\Dvect{D}=\{ D_{k} \}$, $\Dvect{T}=\{ T_{k} \}$ を用いて,
%
\begin{align}
   \DP{\pi}{t} =
          \left( \frac{\partial \pi}{\partial t} \right)_{\rm NG}
     - \Dvect{C} \cdot \Dvect{D}  ,
\end{align}
%
\begin{align}
  \frac{\partial \Dvect{D}}{\partial t} = 
          \left( \DP{\Dvect{D}}{t} \right)_{\rm NG}  
          - \nabla^{2}_{\sigma} ( \Dvect{\Phi}_{S} 
                                  + \underline{W} \Dvect{T}
                                  + \Dvect{G} \pi )  
          + {\cal D}_M \Dvect{D} ,
\end{align}
%
\begin{align}
  \DP{\Dvect{T}}{t} 
     & =   \left( \DP{\Dvect{T}}{t}       \right)_{\rm NG}  
         - \underline{h} \Dvect{D}
         + {\cal D}_H \Dvect{T} .
\end{align}
%
ここで, 非重力波項は,
%
\begin{align}
  \Deqlab{Z項}
  \left( \DP{\pi}{t} \right)^{\rm NG}
   =   - \sum_{k=1}^{K} \Dvect{v}_{k} \cdot \nabla \pi  
       \Delta  \sigma_{k} 
   =   Z_{k} ,
\end{align}
%
\begin{align}
  \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} ,
\end{align}
%
\begin{align}
  \left( \DP{D}{t} \right)^{\rm NG}
   =
            \Dinv{a \cos \varphi}
            \left\{
                \DP{{u_A}_{,k}}{\lambda}
              + \DP{({v_A}_{,k} \cos \varphi)}{\varphi}
            \right\}
          - \nabla^{2}_{\sigma} (\mbox{\sl KE})_{k} 
          - \nabla^{2}_{\sigma} \sum_{k=1}^{K} W_{kl} ( T_{v,l}-T_{l} )
          + {\cal D}(D_{k}) , 
\end{align}
%
\begin{align}
  \left( \DP{T_{k}}{t} \right)^{\rm NG} 
      =   - \frac{1}{a(1-\mu^{2})} 
               \DP{U_k T'_k}{\lambda}
          - \frac{1}{a}
               \DP{V_k T'_k}{\mu}
          + H_{k} 
          + {\cal D}(T_{k}) 
          + {\cal D}^{\prime}(\Dvect{v}) , 
\end{align}
%
\begin{align}
 \begin{split}
  H_k & =     T_{k}^{\prime} D_{k}  \\
     & \quad
           - \frac{1}{\Delta \sigma_{k}} 
             [   \dot{\sigma}_{k-1/2} ( \hat{T^{\prime}}_{k-1/2} 
                                         - T^{\prime}_{k}   )
               + \dot{\sigma}_{k+1/2} ( T^{\prime}_{k}  
                                         - \hat{T^{\prime}}_{k+1/2} ) ]
                \\
     & \quad
           - \frac{1}{\Delta \sigma_{k}} 
             [   \dot{\sigma}^{\rm NG}_{k-1/2} ( \hat{\overline{T}}_{k-1/2} 
                                         - \overline{T}_{k}   )
               + \dot{\sigma}^{\rm NG}_{k+1/2} ( \overline{T}_{k}  
                                         - \hat{\overline{T}}_{k+1/2} ) ]
                \\
     & \quad
           + \hat{\kappa}_{k} T_{v,k} \Dvect{v}_{k} \cdot \nabla \pi
                \\
     & \quad
           - \frac{\alpha_{k}}{\Delta \sigma_{k} } T_{v,k}
             \sum_{l=k}^{K} \Dvect{v}_{l} \cdot \nabla \pi 
               \Delta \sigma_{l}
           - \frac{\beta_{k}}{\Delta \sigma_{k} } T_{v,k}
             \sum_{l=k+1}^{K} \Dvect{v}_{l} \cdot \nabla \pi 
               \Delta \sigma_{l}
                \\
     & \quad
           - \frac{\alpha_{k}}{\Delta \sigma_{k} } T'_{v,k}
             \sum_{l=k}^{K} D_l  \Delta \sigma_{l}
           - \frac{\beta_{k}}{\Delta \sigma_{k} } T'_{v,k}
             \sum_{l=k+1}^{K} D_l  \Delta \sigma_{l} .
 \end{split}
\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} .
\end{align}
%
ここで, 例えば $\delta_{k \leq l}$ は, 
$ k \leq l$ が成り立つとき 1, そうでないとき 0 となる関数である.

以下では, 次のような表現を使用する. 
%
\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}
%
方程式系に semi-implicit 法を適用すると,
%
\begin{align}
  \Deqlab{semi-imp pi}
  \delta_{t} \pi =
          \left( \DP{\pi}{t} \right)_{\rm NG}  
     - \Dvect{C} \cdot \overline{ \Dvect{D} }^{t} ,
\end{align}
%
\begin{align}
  \Deqlab{semi-imp D}
  \delta_{t} \Dvect{D} =
          \left( \DP{\Dvect{D}}{t} \right)_{\rm NG}  
          - \nabla^{2}_{\sigma} ( \Dvect{\Phi}_{S} 
                                  + \underline{W} 
                                     \overline{ \Dvect{T} }^{t}
                                  + \Dvect{G}
                                  \overline{\pi}^{t} ) 
          + {\cal D}_M ( \Dvect{D}^{t-\Delta t} 
                         + 2 \Delta t \delta_{t} \Dvect{D} ) ,
\end{align}
%
\begin{align}
  \Deqlab{semi-imp T}
  \delta_{t} \Dvect{T} =
        \left( \DP{\Dvect{T}}{t} \right)_{\rm NG}  
         - \underline{h} \overline{ \Dvect{D} }^{t} 
         + {\cal D}_H ( \Dvect{T}^{t-\Delta t}
                        + 2 \Delta t \delta_{t} \Dvect{T} ) . 
\end{align}
%
\Deqref{semi-imp pi}, \Deqref{semi-imp D}, \Deqref{semi-imp T}
より, 
%
\begin{align}
 \Deqlab{semi-imp barD}
 \begin{split}
 &   \left\{   ( 1-2\Delta t {\cal D}_H )( 1-2\Delta t {\cal D}_M )
                 \underline{I}  
             - ( \Delta t )^{2}
                 (   \underline{W} \ \underline{h} 
                   + \Dvect{G} \Dvect{C}^{T} ) \nabla^{2}_{\sigma}
     \right\}
      \overline{ \Dvect{D} }^{t} 
       \\
 & \qquad
     =   ( 1-2\Delta t {\cal D}_H )( 1-\Delta t {\cal D}_M ) 
           \Dvect{D}^{t-\Delta t}
       + ( 1-2\Delta t {\cal D}_H ) \Delta t 
           \left( \DP{\Dvect{D}}{t} \right)_{\rm NG}  
  \\
 & \qquad \quad
       -  \Delta t \nabla^{2}_{\sigma}     
                   \left\{ \quad
                            ( 1-2\Delta t {\cal D}_H ) \Dvect{\Phi}_{S} 
                          + \underline{W} 
                            \left[ ( 1-2\Delta t {\cal D}_H ) 
                                    \Dvect{T}^{t-\Delta t}
                                  + \Delta t 
                                      \left( \DP{\Dvect{T}}{t}     
                                      \right)_{\rm NG} \right]
                   \right.
  \\
  & \qquad
                   \left.  \hspace*{20mm} 
                          + ( 1-2\Delta t {\cal D}_H ) \Dvect{G} 
                            \left[ \pi^{t-\Delta t} 
                                  + \Delta t
                                     \left( \DP{\pi}{t} \right)_{\rm NG}  
                            \right]
                   \right\} . 
 \end{split}
\end{align}
%
球面調和関数展開を用いているので,
%
\begin{align}
     \nabla^{2}_{\sigma} = - \frac{n(n+1)}{a^{2}}
\end{align}
%
であり, \Deqref{semi-imp barD}
を$\overline{ \Dvect{D}_n^m }^{t}$ について解くことができる.
%
その後,
%
\begin{align}
   D^{t+\Delta t} = 2\overline{ \Dvect{D} }^{t} - D^{t-\Delta t}
\end{align}
%
および, \Deqref{semi-imp pi}, \Deqref{semi-imp T}
により$t+\Delta t$ における値$\hat{\cal A}^{t+\Delta t}$
が求められる.
