% 表題   DCPAM5  コード解説 -- 力学過程 -- セミインプリシット時間積分(サブルーチン TimeIntegration)
%
% 履歴
%\Drireki{2009/10/05  高橋 芳幸}
%\Drireki{2009/07/12  森川 靖大}
%\Drireki{2009/02/24  森川 靖大}

\section[セミインプリシット時間積分]{セミインプリシット時間積分 \\ (サブルーチン {\tt TimeIntegration})}

注: 本節は書き換え中であるため, 支配方程式系とその離散化の文章とは対応していない. 

\ref{sec:セミインプリシット時間積分}節で
述べた時間積分の2つの手続きは
{\tt TimeIntegration} で行われる. ここでは {\tt TimeIntegration} 内の計算の流れを, 物理量とモデ
ル変数の対応を示しながら説明する.

時間積分の第1段階は
\Deqref{code-smiimp:発散項のみに関する単一のセミインプリシット方程式}
を解くことである.
時間に依存しない行列 $\tilde{\underline{M}}^{m}_{n}$ はサブルーチン
{\tt SemiImplMatrix} で計算され, LU分解されて 
$\Carray{wzz\_siMtxLU}$として {\tt TimeIntegration} で使用される.
本来{\tt SemiImplMatrix}は, サブルーチン{\tt Dynamics}の冒頭で呼ばれるが,
ここに$\tilde{\underline{M}}^{m}_{n}$
を表す式とモデル内の変数との対応を示すこととする.
%
%\\-------------------------------------------ここから消す予定
%
\begin{align}
 \underbrace{\tilde{\underline{M}}^{m}_{n}}_{\Carray{wzz\_siMtxM}}
   & \equiv
      ( 1-2 \!\!\!\!\!\! \overbrace{\Delta t}^{\Carray{DelTime}}
        \!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!
        \underbrace{\tilde{\cal D}_{H,n}^{m}}_{\Carray{wz\_HDifCoefH}}
        \!\!\!\!\!\!\!\!\!\! )
%
      ( \underline{I}-2\Delta t
        \!\!\!\!\!\!\!\!\!
        \overbrace{ \underline{ \tilde{\cal D}_{M,n}^{m} } }^{\Carray{wz\_DisCoefM}}
        \!\!\!\!\!\!\!\!\!\!
      ) \nonumber \\
   & \ \ \ \ 
      - ( \Delta t )^{2}
        [
              \underbrace{
                 \overbrace{\underline{W}}^{\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\Carray{zz\_siMtxW}} \
                 \overbrace{\underline{h}}^{\Carray{zz\_siMtxH}\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!}
              }_{\Carray{zz\_siMtxWH}}
           + ( 1-2\Delta t \tilde{\cal D}_{H,n}^{m} )
             \underbrace{
                 \overbrace{\Dvect{G}}^{\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\Carray{z\_siMtxG}} \
                 \overbrace{\Dvect{C}^{T}}^{\Carray{z\_DelSigma}\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!}
             }_{\Carray{zz\_siMtxGCt}}
        ]
        \biggl(\underbrace{ - \frac{n(n+1)}{a^2} }_{
        \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\Carray{wz\_LaplaEigVal}\quad}\!
        \biggr). 
\end{align}
%
%\\-------------------------------------------ここまで消す予定
%%
%\begin{align}
% \underbrace{\tilde{\underline{M}}^{m}_{n}}_{\Carray{wzz\_siMtxM}}
%   & \equiv
%      ( \underline{I}-2 \!\!\!\!\!\! \overbrace{\Delta t}^{\Carray{DelTime}}
%        \!\!\!\!\!\!\!\!\!\!\!\!\!\!
%        \underbrace{ \underline{ \tilde{\cal D}_{M,n}^{m} } }_{\Carray{wz\_DisCoefM}}
%        \!\!\!\!\!\!\!\!\!\!
%      ) \nonumber \\
%   & \ \ \ \ 
%      - ( \Delta t )^{2}
%        \biggl(\underbrace{ - \frac{n(n+1)}{a^2} }_{
%        \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\Carray{wz\_LaplaEigVal}\quad}\!
%        \biggr) 
%        [
%              \underbrace{
%                 \overbrace{\underline{W}}^{\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\Carray{zz\_siMtxW}} \
%      ( \underline{I} - 2 \Delta t
%        \!\!\!\!\!\!\!\!\!\!
%        \underbrace{\underline{ \tilde{\cal D}_{H} }_{n}^{m}}_{\Carray{wz\_HDifCoefH}}
%        \!\!\!\!\!\!\!\!\!\! )^{-1}
%                 \overbrace{\underline{h}}^{\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\Carray{zz\_siMtxH}\!\!\!\!\!\!\!\!\!\!\!\!\!\!}
%              }_{\Carray{zz\_siMtxWH}}
%           + \underbrace{
%                 \overbrace{\Dvect{G}}^{\!\!\!\!\!\!\!\!\!\!\!\Carray{z\_siMtxG}} \
%                 \overbrace{\Dvect{C}^{T}}^{\Carray{z\_DelSigma}\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!}
%             }_{\Carray{zz\_siMtxGCt}}
%        ].
%\end{align}
%
この$\tilde{\underline{M}}^{m}_{n}$を用い,
$\tilde{\Dvect{f}}^{m}_{n}$ を求める.
\Deqref{code-smiimp:ベクトルfの式} を再掲し,
各項の下にモデル内の変数名を記す.

%
%----------------------------------ここから消す予定
%
\renewcommand{\arraystretch}{0.5}
\begin{align}
% \Deqlab{code-tintgr:ベクトルf}
 \begin{split}
 \underbrace{\tilde{\Dvect{f}}^{m}_{n}}_{\Warray{wz\_siVectF}} \!\!\!\!\!\!\!
  &\equiv   ( 1-2 \!\!\!\!\!\overbrace{\Delta t}^{\Carray{DelTime}} \!\!\!\!\!
                \!\!\!\!\!\!\!
                \underbrace{\tilde{{\cal D}}_{H,n}^{m}}_{
                \!\!\!\!\!\!\Carray{wz\_HDifCoefH}\!\!\!\!\!\!} \!\!\!\!\!\!\!)
            (
              \underline{I} - \Delta t
                \!\!\!\!\!\!\!\!\!
                \overbrace{ \underline{ \tilde{{\cal D}}_{M} }_{n}^{m}}^{\!\!\!\!\Warray{\Carray{wz\_DisCoefM}}\!\!}
                \!\!\!\!\!\!\!\!
            )
%
     \underbrace{\tilde{\Dvect{D}}^{m,t-\Delta t}_{n}}_{\!\!\!\!\!\Iarray{wz\_DivB}\!\!\!\!\!\!\!\!\!\!\!}
 \\
 & \quad
  + ( 1-2\Delta t \tilde{{\cal D}}_{H,n}^{m} ) \Delta t 
        \underbrace{ \left( \DP{\tilde{\Dvect{D}}^{m}_{n}}{t} \right)^{\rm NG} }_{
                  \Iarray{wz\_DDivDtNG}}
 \\
 & \quad
     -  \Delta t \biggl(\underbrace{ - \frac{n(n+1)}{a^2} }_{
      \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\Carray{wz\_LaplaEigVal}\quad} \! \biggr)
                   \Biggl\{  ( 1-2\Delta t \tilde{\cal D}_{H,n}^{m} ) 
                            \underbrace{ \tilde{\Dvect{\Phi}}_{s,n}^{m} }_{
      \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\Iarray{w\_SurfGeoPot}}
  \\
  &                  \hspace*{20mm} 
     + \underbrace{ \underline{W} 
             \overbrace{ \Biggl[ ( 1-\Delta t \tilde{\cal D}_{H,n}^{m} ) 
                                  \underbrace{ \tilde{\Dvect{T}}^{m,t-\Delta t}_{n} }_{ 
          \!\!\!\!\!\!\Iarray{wz\_TempB}\!\!\!\!\!\!}
                                  + \Delta t 
                                      \underbrace{ \left( \DP{\tilde{\Dvect{T}}^{m}_{n}}{t}
                                      \right)^{\rm NG} }_{ 
      \!\!\!\!\!\!\!\!\!\!\!\!\Iarray{wz\_DTempDtNG}}
                         \Biggr] }^{\Warray{wz\_siTemp}}
       }_{\Warray{wz\_siPhi}}
 \\
  &                  \hspace*{20mm} 
                          + ( 1-2\Delta t \tilde{\cal D}_{H,n}^{m} )  \!\!\!\!\!\!
                            \underbrace{ \Dvect{G} }_{
                  \!\!\!\!\!\!\!\!\Carray{z\_siMtxG}\quad} \!\!\!\!\!\!
                            \underbrace{ \mbox{\huge \lower.2ex\hbox{[}} 
                                 \overbrace{ \tilde{\pi}^{m,t-\Delta t}_{n} }^{
              \!\!\Iarray{w\_PiB}\!\!\!\!\!\!\!\!}
                                  + \Delta t
                            \overbrace{ \left( \DP{\tilde{\pi}^{m}_{n}}{t} \right)^{\rm NG} }^{ 
                          \!\!\Iarray{w\_DPiDtNG}}
                            \mbox{\huge \lower.2ex\hbox{]}} }_{\Warray{w\_siPi}}
                   \Biggr\} . 
 \end{split}
\end{align}
%
%\\----------------------------------ここまで消す予定
%%
%\renewcommand{\arraystretch}{0.5}
%\begin{align}
% \Deqlab{code-tintgr:ベクトルf}
% \begin{split}
% \underbrace{\tilde{\Dvect{f}}^{m}_{n}}_{\Warray{wz\_siVectF}} \!\!\!\!\!\!\!
%  &\equiv
%%            ( 1-2 \!\!\!\!\!\overbrace{\Delta t}^{\Carray{DelTime}} \!\!\!\!\!
%%                \!\!\!\!\!\!\!
%%                \underbrace{\tilde{{\cal D}}_{H,n}^{m}}_{
%%                \!\!\!\!\!\!\Carray{wz\_HDifCoefH}\!\!\!\!\!\!} \!\!\!\!\!\!\!)
%            (
%              \underline{I} - \!\!\!\!\!\! \underbrace{\Delta t}_{\Carray{DelTime}}
%                \!\!\!\!\!\!\!\!\!\!\!\!\!
%                \overbrace{ \underline{ \tilde{{\cal D}}_{M} }_{n}^{m}}^{\!\!\!\!\Warray{\Carray{wz\_DisCoefM}}\!\!}
%                \!\!\!\!\!\!\!\!
%            )
%%
%     \underbrace{\tilde{\Dvect{D}}^{m,t-\Delta t}_{n}}_{\!\!\!\!\!\Iarray{wz\_DivB}\!\!\!\!\!\!\!\!\!\!\!}
%% \\
%% & \quad
%%  + ( 1-2\Delta t \tilde{{\cal D}}_{H,n}^{m} ) 
%  + \Delta t 
%        \underbrace{ \left( \DP{\tilde{\Dvect{D}}^{m}_{n}}{t} \right)^{\rm NG} }_{
%                  \Iarray{wz\_DDivDtN}}
% \\
% & \quad
%     -  \Delta t \biggl( \!\!\! \overbrace{ - \frac{n(n+1)}{a^2} }^{\Carray{wz\_LaplaEigVal}} \!\!\! \biggr)
%         \Biggl\{  
%%                            ( 1-2\Delta t \tilde{\cal D}_{H,n}^{m} ) 
%                  \!\!\!\!\!\!\!\!\!\!
%                  \underbrace{ \tilde{\Dvect{\Phi}}_{s,n}^{m} }_{\Iarray{w\_SurfGeoPot}}
%  \\
%  &               \hspace*{20mm} 
%                + \underbrace{ 
%                     \underline{W} 
%                     ( \underline{I}-2\Delta t \underline{ \tilde{\cal D}_{H} }_{n}^{m} )^{-1}
%                  \overbrace{ \Biggl[ ( \underline{I}-\Delta t \underline{ \tilde{\cal D}_{H} }_{n}^{m} ) 
%                                  \underbrace{ \tilde{\Dvect{T}}^{m,t-\Delta t}_{n} }_{ 
%          \!\!\!\!\!\!\Iarray{wz\_TempB}\!\!\!\!\!\!}
%                                  + \Delta t 
%                                      \underbrace{ \left( \DP{\tilde{\Dvect{T}}^{m}_{n}}{t}
%                                      \right)^{\rm NG} }_{ 
%      \!\!\!\!\!\!\!\!\!\!\!\!\Iarray{wz\_DTempDtN}}
%                         \Biggr] }^{\Warray{wz\_siTemp}}
%       }_{\Warray{wz\_siPhi}}
% \\
%  &                  \hspace*{20mm} 
%                          + \!\!\!\!\!\! \overbrace{ \Dvect{G} }^{\!\!\!\!\!\!\!\!\Carray{z\_siMtxG}\quad} 
%                            \!\!\!\!\!\!
%                            \underbrace{ \mbox{\huge \lower.2ex\hbox{[}} 
%                                 \overbrace{ \tilde{\pi}^{m,t-\Delta t}_{n} }^{
%              \!\!\Iarray{w\_PiB}\!\!\!\!\!\!\!\!}
%                                  + \Delta t
%                            \overbrace{ \left( \DP{\tilde{\pi}^{m}_{n}}{t} \right)^{\rm NG} }^{ 
%                          \!\!\Iarray{w\_DPiDtN}}
%                            \mbox{\huge \lower.2ex\hbox{]}} }_{\Warray{w\_siPi}}
%                   \Biggr\} . 
% \end{split}
%\end{align}
%%

\renewcommand{\arraystretch}{1}

\Deqref{code-smiimp:発散項のみに関する単一のセミインプリシット方程式}は
サブルーチン {\tt LUSolve} により解かれる. 結果は
$\overline{\tilde{\Dvect{D}}^{m}_{n}}^{t}$として
変数$\Warray{wz\_siDivAvrTime}$に格納される. 
{\tt LUSolve} の引数は, $\tilde{\Dvect{f}}^{m}_{n}$および
LU 分解された $\tilde{\underline{M}}^{m}_{n}$ とそのピボット
($\Warray{wz\_siVectF}$, $\Carray{wzz\_siMtxLU}$,
$\Carray{wz\_siMtxPiv}$)
である. 

\subsection{$t+\Delta t$の値の算出}

\Deqref{code-smiimp:Dの時間積分} と \Deqref{code-smiimp:D以外の時間積分}
を解く. 以下に式とモデル変数との対応を記す. 
比湿は非重力波成分しかないのでコードの変数名の末尾は
{\rm NG} ではなく現在時刻を意味する {\rm N} を
用いていることに注意されたい.
%
\begin{align}
\Deqlab{code-tintgr:地表面気圧のｔ＋Δｔの計算}
 \underbrace{\tilde{\pi}^{m,t+\Delta t}_{n}}_{\Oarray{w\_PiA}}
   &= 
        \underbrace{\tilde{\pi}^{m,t-\Delta t}_{n}}_{\Iarray{w\_PiB}}
      + 2\Delta t
        \underbrace{
          \Bigg[
                      \overbrace{\DP[][]{\tilde{\pi}^{m}_{n}}{t}^{\rm NG}}^{\Iarray{w\_DPiDtNG}}
                    - \underbrace{\Dvect{C}^{T}}_{\!\!\!\!\!\!\!\!\!\!\Carray{z\_siMtxC}}
                      \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!
                      \overbrace{\overline{\tilde{\Dvect{D}}^{m}_{n}}^{t}}^{\Warray{wz\_siDivAvrTime}}
         \Bigg]
        }_{\Warray{w\_siDPiDt}},
\end{align}
%
\begin{align}
\Deqlab{code-tintgr:渦度のｔ＋Δｔの計算}
 \underbrace{\tilde{\Dvect{\zeta}}^{m,t+\Delta t}_{n}}_{\Oarray{wz\_VorA}}
   &= \left( \underline{I}-2\Delta t \underline{ \tilde{\cal D}_{M} }_{n}^{m} \right)^{-1}
      \Biggl\{
        \underbrace{\tilde{\Dvect{\zeta}}^{m,t-\Delta t}_{n}}_{\Iarray{wz\_VorB}}
      + 2\Delta t
                      \underbrace{\DP[][]{\, \tilde{\Dvect{\zeta}}^{m}_{n}}{t}^{\rm NG}}_{\Iarray{wz\_DVorDtNG}}
      \Biggr\} \, ,
\end{align}
%
\begin{align}
\Deqlab{code-tintgr:発散のｔ＋Δｔの計算}
 \underbrace{\tilde{\Dvect{D}}^{m,t+\Delta t}_{n}}_{\Oarray{wz\_DivA}}
  &= 2   \overline{\tilde{\Dvect{D}}^{m}_{n}}^{t}
       - \underbrace{\tilde{\Dvect{D}}^{m,t-\Delta t}_{n}}_{\Iarray{wz\_DivB}},
\end{align}
%
\begin{align}
\Deqlab{code-tintgr:温度のｔ＋Δｔの計算}
 \underbrace{\tilde{\Dvect{T}}^{m,t+\Delta t}_{n}}_{\Oarray{wz\_TempA}}
   &= \Dinv{1-2\Delta t \tilde{\cal D}_{H,n}^{m}}
      \Biggl\{
        \underbrace{\tilde{\Dvect{T}}^{m,t-\Delta t}_{n}}_{\Iarray{wz\_TempB}\!\!\!\!\!\!\!\!\!}
      + 2\Delta t
         \underbrace{
           \Biggl[
                  \overbrace{\DP[][]{\, \tilde{\Dvect{T}}^{m}_{n}}{t}^{\rm NG}}^{\!\!\!\!\!\!\!\!\!\Iarray{wz\_DTempDtNG}\!\!\!\!\!\!\!\!\!}
                - \overbrace{\underline{h}}^{\!\!\!\!\!\Carray{zz\_siMtxH}\!\!\!\!\!\!\!\!\!} \overline{\tilde{\Dvect{D}}^{m}_{n}}^{t}
           \Biggr]
         }_{\Warray{wz\_siDTempDt}}
      \Biggr\} \, ,
\end{align}
%
\begin{align}
\Deqlab{code-tintgr:比湿のｔ＋Δｔの計算}
 \underbrace{\tilde{\Dvect{q}}^{m,t+\Delta t}_{n}}_{\Oarray{wzf\_QMixA}}
   &= \Dinv{1-2\Delta t \tilde{\cal D}_{q,n}^{m}}
      \Biggl\{
        \underbrace{\tilde{\Dvect{q}}^{m,t-\Delta t}_{n}}_{\Iarray{wzf\_QMixB}\!\!\!\!\!\!\!\!\!}
      + 2\Delta t
                  \underbrace{\DP[][]{\, \tilde{\Dvect{q}}^{m}_{n}}{t}^{\rm NG}}_{\!\!\!\!\!\!\!\!\!\Iarray{wzf\_DQMixDtN}\!\!\!\!\!\!\!\!\!}
      \Biggr\}.
\end{align}