% 表題 中島健介対流モデル解説  乱流拡散係数の計算(TSTEP)
% 
% 履歴 1998/07/14 小高正嗣
%      1999/06/14  小高正嗣
% 
\subsection{乱流拡散係数の計算}

乱流拡散係数の予報方程式の計算を行なう.

\subsubsection{風から質量フラックスへの換算(風の格子点)(CLWFLX)}

風の場の格子点における質量フラックスをサブルーチン {\tt CLWFLX} を呼ん
で計算する.

\subsubsection*{CLWFLX}

\begin{enumerate}
  \item 機能
    \begin{quote}
      風の場の格子点における質量フラックスを計算する.
    \end{quote}
  \item パラメータ
    \begin{quote}
      \begin{tabular}{lllp{10cm}}
        {\tt U}    & (実数) & {\tt ( -2:NXB, -2:NZB, NROT )} & 風速 x 成分 \\
        {\tt W}    & (実数) & {\tt ( -2:NXB, -2:NZB, NROT )} & 風速 z 成分 \\
        {\tt ICL}  & (整数) & & 時刻パラメータ \\
        {\tt UFLX} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 
        x 方向質量フラックス \\
        {\tt WFLX} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 
        y 方向質量フラックス \\
        {\tt FDZ}   & (実数) & {\tt ( -2:NZB )} & z 方向格子間隔 \\
        {\tt DENS0} & (実数) & {\tt ( -2:NZB )} & 基本場の密度 \\
        {\tt NX  } & (整数) & & 水平方向格子点数 \\
        {\tt NZ  } & (整数) & & 鉛直方向格子点数 \\
        {\tt NXB } & (整数) & & 水平方向次元数 {\tt (=nx+2)} \\
        {\tt NZB } & (整数) & & 鉛直方向次元数 {\tt (=nz+2)} \\
        {\tt NROT} & (整数) & & 時間方向次元数 
      \end{tabular}
    \end{quote}
  \item 解説
    
    風の格子点上で $\rho _{0}u, \rho _{0}w$ を計算する.

    \begin{enumerate}
      \item 引数 {\tt ICL} には {\tt IT1} が渡される.
    \end{enumerate}
\end{enumerate}


\subsubsection{移流項と拡散項の計算(ADDIF1)}

拡散係数の予報式における移流項と拡散項を, サブルーチン {\tt ADDIF1} を
呼んで計算する.

\subsubsection*{ADDIF1}

\begin{enumerate}
  \item 機能
    \begin{quote}
      基本場からの偏差についてのスカラー量の移流と拡散を計算する.
    \end{quote}
  \item パラメータ
    \begin{quote}
      \begin{tabular}{lllp{10cm}}
        {\tt QADV}& (実数) &{\tt ( -2:NXB, -2:NZB, NROT )} & スカラー変数 \\
        {\tt UFLX} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 
        x 方向質量フラックス \\
        {\tt WFLX} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 
        y 方向質量フラックス \\
        {\tt CD1} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 乱流拡散係数(
        $(i+\frac{1}{2}), j+\frac{1}{2})$ 格子上) \\
        {\tt IAD}  & (整数) & & 時刻パラメータ(移流計算用) \\
        {\tt IDF}  & (整数) & & 時刻パラメータ(拡散計算用) \\
        {\tt FLQX} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 
        x 方向スカラー量フラックス \\
        {\tt FLQZ} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 
        z 方向スカラー量フラックス \\
        {\tt DTQADV} & (実数) & {\tt ( -2:NXB, -2:NZB )} 
        & スカラー変数変化量(移流) \\
        {\tt DTQDIF} & (実数) & {\tt ( -2:NXB, -2:NZB )} 
        & スカラー変数変化量(拡散)\\
        {\tt DENS0} & (実数) & {\tt ( -2:NZB )} & 基本場の密度 \\
        {\tt DX}   & (実数) & & x 方向格子間隔 \\
        {\tt FDZ}   & (実数) & {\tt ( -2:NZB )} & z 方向格子間隔 \\
        {\tt FDZM}  & (実数) & {\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 解説

    任意のスカラー量 $Q$ の移流と拡散

    \begin{eqnarray*}
      && - \frac{1}{\rho _{0}}\left[ \DP{}{x}\rho _{0}uQ + \DP{}{z}\rho _{0}wQ 
        \right] \\
      && \DP{}{x}\left(K_{m}\DP{Q}{x}\right) + 
         \frac{1}{\rho _{0}}\DP{}{z}\left(\rho _{0}K_{m}\DP{Q}{z}\right) 
    \end{eqnarray*}

    を計算する.

    風の格子点上でのスカラー量のフラックスを計算する際には, スカラー量
    の 4 点平均値を用いる. 例えば時刻 $(n\Delta t)$ の
    $((i+\frac{1}{2})\Delta x, j\Delta z)$ における x 方向のフラックス
    $FQ_{x(i+\frac{1}{2},j)}^{n}$ は,

    \[
    FQ_{x(i+\frac{1}{2},j)}^{n} = 
      \rho _{0,j}u_{i+\frac{1}{2},j}^{n}
      \left(-\frac{1}{16}Q_{i+2,j} + \frac{9}{16}Q_{i+1,j}
        + \frac{9}{16}Q_{i,j}-\frac{1}{16}Q_{i-1,j}\right)
    \]

    とする. このようにして求めたフラックスを用いて 4 次精度の移流を計
    算する. 例えば

    \[
    \left( \frac{1}{\rho _{0}}\DP{}{x}\rho _{0}uQ \right)_{i,j}^{n}
    = \frac{1}{\rho _{0,j}\Delta x}\left[
      - \frac{1}{24}FQ_{x(i+\frac{3}{2},j)}^{n}
      + \frac{9}{8}FQ_{x(i+\frac{1}{2},j)}^{n}
      - \frac{9}{8}FQ_{x(i-\frac{1}{2},j)}^{n}  
      + \frac{1}{24}FQ_{x(i-\frac{3}{2},j)}^{n}\right]
    \]

    とする.

    拡散項は 2 次精度で計算する. 例えば

    \begin{eqnarray*}
      \left[ \DP{}{x}\left( K_{m}\DP{Q}{x}\right) \right]_{i,j}^{n}
      &=& \frac{1}{(\Delta x)^{2}}\left[
        \frac{(K_{m(i+1,j)}^{n}+K_{m(i,j)}^{n})}{2}
              (Q_{i+1,j}^{n}-Q_{i,j}^{n}) - \right .\\
      && \left. \frac{(K_{m(i,j)}^{n}+K_{m(i-1,j)}^{n})}{2}
              (Q_{i,j}^{n}-Q_{i-1,j}^{n}) \right]
    \end{eqnarray*}

    とする.

    \begin{enumerate}
      \item 引数 {\tt IAD} には {\tt IT1} が, {\tt IDF} には {\tt IT2}
        が渡される.
      \item 境界条件はサブルーチン {\tt BOUND} で設定する.
    \end{enumerate}
\end{enumerate}


\subsubsection{人工非線形拡散項の計算(NLDIFF)}

人工非線形拡散項をサブルーチン {\tt NLDIFF} を呼んで計算する.

\subsubsection*{NLDIFF}

\begin{enumerate}
  \item 機能
    \begin{quote}
     スカラー量の人工非線形拡散を計算する.
    \end{quote}
  \item パラメータ
    \begin{quote}
      \begin{tabular}{lllp{10cm}}
        {\tt QADV}& (実数) &{\tt ( -2:NXB, -2:NZB, NROT )} & スカラー変数 \\
%        {\tt SCTURB} & (実数) & & \\
        {\tt SCALE} & (実数) & & \\
        {\tt IDF}  & (整数) & & 時刻パラメータ \\
        {\tt DTQNLV} & (実数) & {\tt ( -2:NXB, -2:NZB )}  
        & スカラー変数変化量(人工粘性) \\
        {\tt QLAP} & (実数) & {\tt ( -2:NXB, -2:NZB )} & \\ 
        {\tt CDQNLX} & (実数) & {\tt ( -2:NXB, -2:NZB )} & \\ 
        {\tt CDQNLZ} & (実数) & {\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 解説

    \begin{enumerate}
      \item 引数 {\tt IDF} には {\tt IT2} が渡される.
      \item {\tt SCALE} には {\tt SCTURB=2.0D0} が渡される.
        この引数 {\tt SCALE} に渡される値は計算されるスカラー量により異なる.
    \end{enumerate}

\end{enumerate}

\subsubsection{浮力とシアーによる乱流生成の計算(CLTURB)}

浮力とシアーによる乱流生成項をサブルーチン {\tt CLTURB} を呼んで計算す
る.

\subsubsection*{CLTURB}

\begin{enumerate}
  \item 機能
    \begin{quote}
      浮力とシアーによる乱流生成を計算する.
    \end{quote}
  \item パラメータ
    \begin{quote}
      \begin{tabular}{lllp{10cm}}
        {\tt SSTURB} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 
        浮力とシアーによる乱流生成項 \\
        {\tt U}    & (実数) & {\tt ( -2:NXB, -2:NZB, NROT )} & 風速 x 成分 \\
        {\tt W}    & (実数) & {\tt ( -2:NXB, -2:NZB, NROT )} & 風速 z 成分 \\
        {\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 CDTURB}& (実数) &{\tt ( -2:NXB, -2:NZB, NROT )} & 乱流拡散係数 \\
        {\tt IT1}  & (整数) & & 時刻パラメータ(シアー生成項用) \\
        {\tt IT2}  & (整数) & & 時刻パラメータ(安定度計算用) \\
        {\tt TPOT0} & (実数) & {\tt ( -2:NZB )} & 基本場の温位 \\
        {\tt QVAP0} & (実数) & {\tt ( -2:NZB )} & 基本場の水蒸気混合比 \\
        {\tt PPAI0} & (実数) & {\tt ( -2:NZB )} & 基本場の無次元圧力関数 \\
        {\tt ZGRID} & (実数) & {\tt ( -2:NZB )} & 鉛直方向格子点座標 \\
        {\tt FDZ}   & (実数) & {\tt ( -2:NZB )} & z 方向格子間隔 \\
        {\tt FDZM}  & (実数) & {\tt ( -2:NZB )} & z 方向半整数格子間隔 \\
        {\tt DX}   & (実数) & & x 方向格子間隔 \\
        {\tt GRAV}  & (実数) & & 重力加速度 \\
        {\tt CPDRY} & (実数) & & 乾燥空気定圧比熱 \\
        {\tt ALATNT}& (実数) & & 単位質量あたりの潜熱 \\
        {\tt RMVAP} & (実数) & & 湿潤空気分子量 \\
        {\tt STABZ} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 
        安定度(乱流拡散計算用)\\
        {\tt EDFM1} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 
        シアーによる乱流生成項1\\
        {\tt EDFM2} & (実数) & {\tt ( -2:NXB, -2:NZB )} & 
        シアーによる乱流生成項2\\
        {\tt NX  } & (整数) & & 水平方向格子点数 \\
        {\tt NZ  } & (整数) & & 鉛直方向格子点数 \\
        {\tt NXB } & (整数) & & 水平方向次元数 {\tt (=nx+2)} \\
        {\tt NZB } & (整数) & & 鉛直方向次元数 {\tt (=nz+2)} \\
        {\tt NROT} & (整数) & & 時間方向次元数 
      \end{tabular}
    \end{quote}
  \item 解説
    \begin{enumerate}
      \item 引数 {\tt IT1} には {\tt IT1}, {\tt IT2} には {\tt IT2} が渡
        される.
    \end{enumerate}
\end{enumerate}


\subsubsection{乱流拡散係数の時間変化計算}

これまでに求めた移流項, 拡散項, 人工粘性項, 生成項( {\tt DTQADV}, {\tt
DTQDIF}, {\tt DTQDNL}, {\tt SSTURB})に時間刻を掛けて 1 ステップ後の拡
散係数 {\tt CDTURB} を計算する. 但し拡散係数の大きさに下限(0
m${}^{2}$/sec)と上限(800 m${}^{2}$/sec)を与えておく.

