% 表題 中島健介対流モデル解説  運動方程式の計算(TSTEP)
% 
% 履歴 1998/07/12 小高正嗣
%      1999/06/14  小高正嗣
% 
\subsection{運動方程式の計算}

運動方程式の計算を行なう.

\subsubsection{浮力項の計算(CLBUOY)}

サブルーチン {\tt CLBUOY} を呼んで浮力項の計算を行なう.

\subsubsection*{CLBUOY}

\begin{enumerate}
  \item 機能
    \begin{quote}
      運動方程式における浮力項を計算する.
    \end{quote}
  \item パラメータ
    \begin{quote}
      \begin{tabular}{lllp{10cm}}
        {\tt TPOT} & (実数) & {\tt ( -2:NXB, -2:NZB, NROT )} & 温位 \\
        {\tt QVAP} & (実数) & {\tt ( -2:NXB, -2:NZB, NROT )} & 水蒸気混合比 \\
        {\tt QCLW} & (実数) & {\tt ( -2:NXB, -2:NZB, NROT )} & 雲水混合比 \\
        {\tt QRAI} & (実数) & {\tt ( -2:NXB, -2:NZB, NROT )} & 雨水混合比 \\
        {\tt IBY}  & (整数) & & 時間方向次元 \\
        {\tt BUOY} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 浮力項\\
        {\tt BUOYW} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 
        浮力項(湿潤量からの寄与) \\
        {\tt TPOT0} & (実数) & {\tt ( -2:NZB )} & 基本場の温位 \\
        {\tt GRAV}  & (実数) & & 重力加速度 \\
        {\tt RMVAP} & (実数) & & 湿潤空気分子量 \\
        {\tt FDZ}   & (実数) & {\tt ( -2:NZB )} & z 方向格子間隔 \\
        {\tt NX  } & (整数) & & 水平方向格子点数 \\
        {\tt NZ  } & (整数) & & 鉛直方向格子点数 \\
        {\tt NXB } & (整数) & & 水平方向次元数 {\tt (=nx+2)} \\
        {\tt NZB } & (整数) & & 鉛直方向次元数 {\tt (=nz+2)} \\
        {\tt NROT} & (整数) & & 時間方向次元数 \\
      \end{tabular}
    \end{quote}
  \item 解説

    運動方程式中の浮力項

    \[
      b \equiv g\left( \frac{\theta }{\Theta _{0}} + 
        0.608q_{v} - q_{c} - q_{r}\right)
    \]

    を計算する.

    \begin{enumerate}
       \item パラメータ {\tt IBY} には {\tt IT1} が渡される.
       \item 変数 {\tt BUOY} は $g\theta /\Theta _{0}$ であり, {\tt
         BUOYW} は残りの項に対応する.
         温位や混合比の格子点と z 方向の風の格子点が異なるので, 適当に補間して
         浮力項を計算する. ここでは 4 点平均して求める. 例えば

         \[
         \left( g\frac{\theta }{\Theta _{0}}\right)_{i,j+\frac{1}{2}}
         = \frac{g}{(\Theta _{0,j+1}+\Theta _{0,j})/2}
         \left( -\frac{1}{16}\theta _{i,j+2} +
           \frac{9}{16}\theta _{i,j+1} +
           \frac{9}{16}\theta _{i,j} -
           \frac{1}{16}\theta _{i,j-1} \right),
         \]

         などとする.
       \item 境界条件はサブルーチン {\tt BOUND} で与える.
     \end{enumerate}
\end{enumerate}

\subsubsection{浮力による運動エネルギー生成の計算}

浮力による運動エネルギーの生成と湿潤過程に伴うポテンシャルエネルギーの消滅の時空間積算量を計算する. 各時間ステップ毎に,

\begin{eqnarray*}
  && \rho _{0}g\frac{\theta }{\Theta _{0}}w \Dd x \Dd z \\
  && - \rho _{0}g(0.608q_{v} - q_{c} - q_{r})w \Dd x \Dd z 
\end{eqnarray*}

を全計算領域で積分し, それに時間刻($\Delta t$)をかけてそれぞれ変数
{\tt CTHKIN}, {\tt CPWKIN} に足し込んでいく.

\subsubsection{コリオリ力の計算(CLCORI)}

サブルーチン {\tt CLCOLI} を呼んでコリオリ力項の計算を行なう.

\subsubsection*{CLCORI}

\begin{enumerate}
  \item 機能
    \begin{quote}
      運動方程式中のコリオリ力項の計算を行なう.
    \end{quote}
  \item パラメータ
    \begin{tabular}{lllp{10cm}}
      {\tt U}    & (実数) & {\tt ( -2:NXB, -2:NZB, NROT )} & 風速 x 成分 \\
      {\tt V}    & (実数) & {\tt ( -2:NXB, -2:NZB, NROT )} & 風速 y 成分 \\
      {\tt DUCOLI} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 
      風速 x 成分コリオリ力 \\
      {\tt DVCOLI} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 
      風速 y 成分コリオリ力 \\
      {\tt ICL} & (整数) & & 時間方向次元 \\
      {\tt FCOLI} & (実数) & {\tt ( -2:NXB )} & コリオリ定数 \\
      {\tt ZGRID} & (実数) & {\tt ( -2:NZB )} & 鉛直方向格子点座標 \\
      {\tt NX  } & (整数) & & 水平方向格子点数 \\
      {\tt NZ  } & (整数) & & 鉛直方向格子点数 \\
      {\tt NXB } & (整数) & & 水平方向次元数 {\tt (=nx+2)} \\
      {\tt NZB } & (整数) & & 鉛直方向次元数 {\tt (=nz+2)} \\
      {\tt NROT} & (整数) & & 時間方向次元数 \\
    \end{tabular}
  \item 解説
    
    運動方程式中のコリオリ力項

    \[
      fv, -fu 
    \]
    
    を計算する. 風速の x, y 成分の格子点は同じなので簡単に計算できる.

    \begin{enumerate}
      \item パラメータ {\tt ICL} は {\tt IT1} が渡される.
      \item 境界条件は {\tt BOUND} で設定する.
    \end{enumerate}
\end{enumerate}

\subsubsection{地表摩擦の計算}

地表摩擦項

\begin{eqnarray*}
  \frac{1}{\rho _{0}}\DP{}{z}F_{x} &=& \frac{1}{\rho _{0(0)}}
  \frac{F_{x,j+\frac{1}{2}}}{\Delta z _{(0)}}, \\
  \frac{1}{\rho _{0}}\DP{}{z}F_{y} &=& \frac{1}{\rho _{0(0)}}
  \frac{F_{y,j+\frac{1}{2}}}{\Delta z _{(0)}}, 
\end{eqnarray*}

を計算する. それぞれ変数 {\tt DUFRIC}, {\tt DVFRIC} に代入される.


\subsubsection{圧力項以外の運動量変化の計算}

圧力項以外の運動量変化(移流, 乱流拡散, コリオリ, 人工粘性, 浮力, 地表
摩擦)を合計して, x, y, z 方向それぞれ {\tt DUTOTL}, {\tt DVTOTL}, {\tt
DWTOTL} に代入する. 境界条件はサブルーチン {\tt BOUND} で設定する.


\subsubsection{水平平均量の計算(CLZPRF)}

サブルーチン {\tt CLZPROF} を呼んで各種水平平均量を計算する.
計算するのは
\begin{quote}
\begin{tabular}{lllp{10cm}}
    {\tt U}    & (実数) & {\tt ( -2:NXB, -2:NZB, NROT )} & 風速 x 成分 \\
    {\tt V}    & (実数) & {\tt ( -2:NXB, -2:NZB, NROT )} & 風速 y 成分 \\
    {\tt DUADV} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 風速 x 成分変化量(移流)\\
    {\tt DUVIS} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 風速 x 成分変化量(拡散)\\
    {\tt DUNLV} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 
    風速 x 成分変化量(人工粘性)\\
    {\tt DVADV} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 風速 y 成分変化量(移流)\\
    {\tt DVVIS} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 風速 y 成分変化量(拡散)\\
    {\tt DVNLV} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 
    風速 y 成分変化量(人工粘性)\\ 
    {\tt DUCOLI} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 風速 x 成分コリオリ力 \\
    {\tt DVCOLI} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 風速 y 成分コリオリ力 \\
    {\tt DUFRIC} & (実数) & {\tt ( -2:NXB, 2 )} & 風速 x 成分摩擦力 \\
    {\tt DVFRIC} & (実数) & {\tt ( -2:NXB, 2 )} & 風速 y 成分摩擦力 
\end{tabular}
\end{quote}
の水平平均である.

\subsubsection*{CLZPROF}

\begin{enumerate}
  \item 機能
    \begin{quote}
      各種水平平均量を計算する.
    \end{quote}
  \item パラメータ
    \begin{quote}
      \begin{tabular}{lllp{10cm}}
        {\tt XDATA} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 変数 \\
        {\tt ICL} & (整数) & & 時間方向次元 \\
        {\tt IST} & (整数) & & 水平平均量番号 \\
        {\tt IAC} & (整数) & & 内部パラメータ \\
        {\tt FDT} & (実数) & & 時間刻 \\
        {\tt ZPROF}& (実数) & {\tt ( -2:NZB, NZPROF )} &  各種水平平均量 \\
        {\tt NX  } & (整数) & & 水平方向格子点数 \\
        {\tt NZ  } & (整数) & & 鉛直方向格子点数 \\
        {\tt NXB } & (整数) & & 水平方向次元数 {\tt (=nx+2)} \\
        {\tt NZB } & (整数) & & 鉛直方向次元数 {\tt (=nz+2)} \\
        {\tt NROT} & (整数) & & 時間方向次元数 \\
        {\tt NZPROF} & (整数) & & 水平平均量次元数
      \end{tabular}
    \end{quote}
  \item 解説
    \begin{enumerate}
      \item パラメータ {\tt ICL} には変数そのものの平均をとる場合 {\tt IC1},
        収支を計算する場合は {\tt 1} を渡す.
      \item パラメータ {\tt IST} には計算する水平平均量に対応する 4.3
        節で定義した番号を渡す.
      \item パラメータ {\tt IAC} は内部の {\tt IF} 分岐に用いる. 変数
        そのものの水平平均をとる場合は {\tt IAC} に {\tt 0} を渡す.
        収支を計算する場合は {\tt 1} を渡す.
      \item パラメータ {\tt FDT} は {\tt IAC} に {\tt 0} を渡す場合は
        {\tt 0.0D0}, {\tt IAC} に {\tt 1} を渡す場合は {\tt DTBGT} を
        渡す.
      \item 収支を計算する場合には時間刻 {\tt DTBGT} を掛けた量を計算する.
        このときの各 {\tt  ZPROF} は時間方向に積算された値となる.
    \end{enumerate}
\end{enumerate}

\subsubsection{各格子点での質量収束の計算(CLCONM)}

スカラー量の格子点で質量収束を計算する.

\subsubsection*{CLCONM}

\begin{enumerate}
  \item 機能
    \begin{quote}
      スカラー量の格子点で質量収束を計算する.
    \end{quote}
  \item パラメータ
    \begin{quote}
      \begin{tabular}{lllp{10cm}}
        {\tt U}    & (実数) & {\tt ( -2:NXB, -2:NZB, NROT )} & 風速 x 成分 \\
        {\tt V}    & (実数) & {\tt ( -2:NXB, -2:NZB, NROT )} & 風速 y 成分 \\
        {\tt IC1} & (整数) & & 時間方向次元 \\
        {\tt DTI} & (実数) & & 時間刻 \\
        {\tt CONVM} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 
        スカラー量格子点での質量収束 \\
        {\tt DENS0} & (実数) & {\tt ( -2:NZB )} & 基本場の密度 \\
        {\tt FDZ}   & (実数) & {\tt ( -2:NZB )} & z 方向格子間隔 \\
        {\tt NX  } & (整数) & & 水平方向格子点数 \\
        {\tt NZ  } & (整数) & & 鉛直方向格子点数 \\
        {\tt NXB } & (整数) & & 水平方向次元数 {\tt (=nx+2)} \\
        {\tt NZB } & (整数) & & 鉛直方向次元数 {\tt (=nz+2)} \\
        {\tt NROT} & (整数) & & 時間方向次元数 
      \end{tabular}
    \end{quote}
  \item 解説
\end{enumerate}


\subsubsection{圧力方程式の計算(CLPRES)}

圧力の診断方程式をサブルーチン {\tt CLPRES} を呼んで計算する.
さらに圧力変化に伴う運動量変化を計算する.

\subsubsection*{CLPRES}

\begin{enumerate}
  \item 機能
    \begin{quote}
      圧力の診断方程式と圧力変化に伴う運動量変化を計算する.
    \end{quote}
  \item パラメータ
    \begin{quote}
      \begin{tabular}{lllp{10cm}}
        {\tt DUTOTL} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 
        風速 x 成分総変化量 \\
        {\tt DWTOTL} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 
        風速 z 成分総変化量 \\
        {\tt CONVM} & (実数) & {\tt ( -2:NXB, -2:NZB )} &
        スカラー量格子点での質量収束 \\
        {\tt PFUNC} & (実数) & {\tt ( -2:NXB, -2:NZB )} & \\
        {\tt PPAI} & (実数) & {\tt ( -2:NXB, -2:NZB, NROT )} & 
        無次元圧力関数 \\
        {\tt PDIV} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 計算用配列1\\ 
        {\tt PDIVM} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 計算用配列2\\ 
        {\tt PFUNCM} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 計算用配列3\\ 
        {\tt ALPHA} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 計算用配列4\\ 
        {\tt GAMMA} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 計算用配列5\\ 
        {\tt VMODE} & (実数) & {\tt ( -2:NZB, NZ )} & 鉛直方向固有ベクトル行列\\
        {\tt TMODE} & (実数) & {\tt ( -2:NZB, NZ )} & {\tt VMODE} の逆行列\\
        {\tt FMODE} & (実数) & {\tt ( -2:NXB, NZ )} & 水平方向係数逆行列\\
        {\tt TRIGS} & (実数) & {\tt ( NX+15 )} & \\
        {\tt IFAC}  & (実数) & {\tt ( 15 )} & \\
        {\tt UUFT}  & (実数) & {\tt ( NX, NZ )} & フーリエ変換用配列\\
        {\tt ANFT}  & (実数) & {\tt ( NX, NZ )} & フーリエ逆変換用配列\\
        {\tt DX}    & (実数) & & x 方向格子間隔 \\
        {\tt CPDRY} & (実数) & & 乾燥空気定圧比熱 \\
        {\tt FDZM}  & (実数) & {\tt ( -2:NZB )} & z 方向半整数格子間隔 \\
        {\tt TPOT0} & (実数) & {\tt ( -2:NZB )} & 基本場の温位 \\
        {\tt NX  } & (整数) & & 水平方向格子点数 \\
        {\tt NZ  } & (整数) & & 鉛直方向格子点数 \\
        {\tt NXB } & (整数) & & 水平方向次元数 {\tt (=nx+2)} \\
        {\tt NZB } & (整数) & & 鉛直方向次元数 {\tt (=nz+2)} 
      \end{tabular}
    \end{quote}
  \item 解説

    解くべき圧力の診断式は

    \[
      \left(\DP[2]{}{x} + \frac{1}{\rho _{0}}\DP{}{z}\rho _{0}\DP{}{z}\right)
      c_{p}\Theta 
      \pi = \DP{}{x}F_{u} +\frac{1}{\rho _{0}}\DP{}{z}\rho _{0}F_{w},
    \]

    である. ただし $\Theta $ はの鉛直変化は非常に小さいとする. 右辺の
    $F_{u}, F_{w}$ は $u, w$ の式に現れる移流項, コリオリ項, 拡散項,
    浮力項の和である.
    
    壁での境界条件は $w=0$. 運動方程式から,

    \[
      -c_{p}\Theta \DP{\pi}{z} + F_{w} = 0,
    \]

    となる. これより $\pi $ についての境界条件は時間変化する. この時間
    依存性を含まないようにして解くために以下のような方法で解くことにする.

    まず,

    \[
      \gamma \equiv \int _{0}^{z} \frac{F_{w}} \Dd z,
    \]

    を定義する. 次に 

    \[
      P \equiv c_{p}\Theta \pi - \gamma ,
    \]

    を用意する. これを鉛直方向の運動方程式に代入すると

    \[
      \DP{w}{t} = -\DP{P}{z},
    \]

    となる. 壁で $w=0$ を考えると $P$ についての境界条件は時間依存しな
    いことになる. よって $\pi $ の代わりに $P$ を解くことにする.
    $\pi $ の診断式に $\pi = P +\gamma $ を代入すると,

    \[
      \left(\DP[2]{}{x} +\frac{1}{\rho _{0}}\DP{}{z}\rho _{0}\DP{}{z}\right)
      P = \DP{}{x}F_{u} - \DP[2]{}{x}\gamma ,
    \]

    となる. $\alpha \equiv F_{u} - \DP{}{x}\gamma $ を用意すれば,

    \[
      \left(\DP[2]{}{x} +\frac{1}{\rho _{0}}\DP{}{z}\rho _{0}\DP{}{z}\right)
      P = \DP{}{x}\alpha ,
    \]

    である. この $\alpha $ を用いると水平方向の運動方程式は

    \[
      \DP{u}{t} = -\DP{P}{x} + \alpha ,
    \]

    となる.

    \begin{enumerate}
      \item $\gamma $ を計算する.
      \item $\alpha $ を計算する. $\gamma $ の微分は 2 次中心差分で行なう.
      \item $P$ についてのポアソン方程式を解く.  サブルーチン {\tt
        PSETUP} において $z$ 微分項の係数行列の固有関数行列
        $\Dvect{V}$({\tt VMODE})及びその逆行列 $\Dvect{V}^{-1}$({\tt
        TMODE}), $x$ 方向係数逆行列 $\Dvect{F}$({\tt FMODE}) は既に計
        算されているので, $P$ の $x$ 依存部分 $\Dvect{H}$ を

        \[
        \Dvect{H} = \Dvect{F}\Dvect{V}^{-1}\Dvect{S}
        \]

        から計算する(4.10節 {\tt PSETUP} の解説参照).

        \begin{enumerate}
        \item $\Dvect{V}^{-1}\Dvect{S}$ を計算し, それをフーリエ変換す
          る. $\Dvect{S}$ を計算する際 $\alpha $ の $x$ 微分は 4 次中
          心差分で行なう.
          \item $P$ の $x$ 依存部分 $\Dvect{H}$ を計算する.
          \item $\Dvect{H}$ と $\Dvect{V}$ から $P$ を計算する.
        \end{enumerate}

      \item $\pi $ を $P$ から計算する.
   \end{enumerate}

\end{enumerate}

\subsubsection{時間積分}

圧力項を含めた運動量変化から 1 ステップ時間積分を行なう. 境界条件はサ
ブルーチン{\tt BOUND} で設定する.





