% 表題  中島健介対流モデル解説  圧力方程式解法ルーチンの初期設定
%
% 履歴  1998/07/03  小高正嗣
%       1998/08/28  小高正嗣
%       1999/06/14  小高正嗣
%
\subsection{圧力方程式解法ルーチンの初期設定}

サブルーチン {\tt PSETUP} を呼んで圧力方程式解法ルーチンの初期設定を行
なう.

\subsubsection{PSETUP}

\begin{enumerate}
  \item 機能
    \begin{quote}
      圧力方程式解法ルーチンの初期設定を行なう. 
    \end{quote}
  \item パラメータ
    \begin{quote}
      \begin{tabular}{lllp{10cm}}
        {\tt DENS0} & (実数) & {\tt ( -2:NZB )} & 基本場の密度 \\
        {\tt FDZ}   & (実数) & {\tt ( -2:NZB )} & z 方向格子間隔 \\
        {\tt FDZM}  & (実数) & {\tt ( -2:NZB )} & z 方向半整数格子間隔 \\
        {\tt VMODE} & (実数) & {\tt ( -2:NZB, NZ )} &鉛直方向固有ベクトル行列\\
        {\tt TMODE} & (実数) & {\tt ( -2:NZB, NZ )} & {\tt VMODE} の逆行列\\
        {\tt VEIGEN}& (実数) & {\tt ( NZ )}     & 鉛直方向固有値\\
        {\tt FMODE} & (実数) & {\tt ( -2:NXB, NZ )} & 水平方向係数逆行列\\
        {\tt TRIGS} & (実数) & {\tt ( NX+15 )} & \\
        {\tt IFAC}  & (実数) & {\tt ( 15 )} & 
      \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$ の式に現れる移流項, コリオリ項, 拡散項,
    浮力項の和である.
    
    ここでは Dimension Reduction 法を用いてこの式を解く. これは変数分
    離法の一種である. 
    圧力の式を適当に差分化すると, 

    \[
      \Dvect{D_{x}\pi} + \Dvect{D_{z}\pi} = \Dvect{S}
    \]

    と書くことができる. ただし $\Dvect{D_{x}}, \Dvect{D_{z}},
    \Dvect{S}$ は行列である.  行列 $\Dvect{D_{z}}$ の固有値を $\lambda
    _{i}$, 固有ベクトルを $V_{i}(z)$ とし, 固有列ベクトルを並べた行列を
    $\Dvect{V}$, 対角成分が $\lambda _{i}$ である行列を $
    \Dvect{\Lambda }$ とすると, $\Dvect{D_{z}V} = \Dvect{V\Lambda }$
    である.  $\Dvect{\pi } = \Dvect{V}\cdot \Dvect{H}$ と展開すると

    \[
      \Dvect{V}\Dvect{D_{x}H} + \Dvect{V\Lambda H} = \Dvect{S}
    \]
    
    よって

    \[
      \left(\Dvect{D_{x}} + \Dvect{\Lambda }\right) \Dvect{H}
      = \Dvect{V}^{-1}\Dvect{S}
    \]
    
    となる. 次は $x$ 方向に適当な固有関数に展開して各モードに対する展
    開係数を求めればよい. ここでは三角関数で展開するのでこの式をフーリエ
    変換して解く.


    \begin{enumerate}
      \item $z$ 方向の係数行列を求める. 
        
        \[
        \frac{1}{\rho _{0}}\DP{}{z}\rho _{0}\DP{}{z} \pi 
        \]
    
        を差分化するのだが, 最初の微分は 2 次中央差分, 次の微分は 4 次
        中央差分で解くことに注意する. これは連続の式を 4 次中央差分し,
        圧力の微分は 2 次中央差分であることによる. よって係数行列は 5
        重対角行列になる.
        係数行列の成分を $A_{ij}$ とすると,

        \begin{eqnarray*}
          A_{i,i+2} &=& -\frac{1}{\rho _{0}\Delta z_{i}}\left(
            \frac{1}{24\Delta zm_{i+2}}\frac{(\rho _
            {0,i+2}+\rho _{0,i+1})}{2}\right) \\
          A_{i,i+1} &=& \frac{1}{\rho _{0}\Delta z_{i}}\left(
            \frac{1}{24\Delta zm_{i+2}}\frac{(\rho _
            {0,i+2}+\rho _{0,i+1})}{2} + 
            \frac{9}{8\Delta zm_{i+1}}\frac{(\rho _
            {0,i+1}+\rho _{0,i})}{2}\right) \\
          A_{i,i} &=& - \frac{1}{\rho _{0}\Delta z_{i}}\left(
            \frac{9}{8\Delta zm_{i+1}}\frac{(\rho _
            {0,i+1}+\rho _{0,i})}{2} 
            + \frac{9}{8\Delta zm_{i}}\frac{(\rho _
            {0,i}+\rho _{0,i-1})}{2}\right) \\
          A_{i+1,i} &=& \frac{1}{\rho _{0}\Delta z_{i}}\left(
            \frac{1}{24\Delta zm_{i-1}}\frac{(\rho _
            {0,i-1}+\rho _{0,i-2})}{2} + 
            \frac{9}{8\Delta zm_{i}}\frac{(\rho _
            {0,i}+\rho _{0,i-1})}{2}\right) \\
          A_{i+2,i} &=& -\frac{1}{\rho _{0}\Delta z_{i}}\left(
            \frac{1}{24\Delta zm_{i-1}}\frac{(\rho _
          {0,i-1}+\rho _{0,i-2})}{2} \right)
        \end{eqnarray*}
        
        となる. ただし $\Delta z_{i}$ は風の鉛直成分の鉛直格子間隔, 
        $\Delta zm_{i}$ はスカラー量の鉛直格子間隔である.
        境界条件は壁で $\DP{}{z}=0$ となるようにする. 

      \item 係数行列 $A_{ii}$ の固有値 $\lambda _{i}$ と固有ベクトル
        $V_{i}(z)$ を求める. 計算は lapack のサブルーチン {\tt DGEEVX}
        で行なう. 
        
      \item 固有列ベクトルを並べた行列 $\Dvect{V}$ ({\tt VMODE}) の逆
        行列 ({\tt TMODE}) を求める. 計算は lapack のサブルーチン
        {\tt DGESVX}で行なう.

      \item 水平方向係数行列 $\Dvect{D_{x}}+\Dvect{\Lambda}$ をフーリ
        エ変換した行列の逆行列 ({\tt FMODE}) を計算する.


    \end{enumerate}

\end{enumerate}








