% 表題   DCPAM5  乱流過程
%
% 履歴
%\Drireki{2010/04/15 高橋芳幸}
%
%  \Dchapterhead
\section{離散表現}
\Dseclab{乱流過程離散表現}


\Dmodel では, 鉛直拡散は陰解法を用いて計算する. 
運動量, 熱の鉛直拡散方程式は下のように離散化する. 
%
\begin{eqnarray}
  \frac{ u^{t+\Delta t}_k - u^{t-\Delta t}_k }{ 2 \Delta t } 
    &=& g \frac{ F_{m,x,k+\frac{1}{2}}^{t+\Delta t} -
    F_{m,x,k-\frac{1}{2}}^{t+\Delta t} }{ p_{k+\frac{1}{2}} -
    p_{k-\frac{1}{2}} } , 
   \Deqlab{U鉛直拡散方程式離散表現}
\\
  \frac{ v^{t+\Delta t}_k - v^{t-\Delta t}_k }{ 2 \Delta t } 
    &=& g \frac{ F_{m,y,k+\frac{1}{2}}^{t+\Delta t} -
    F_{m,y,k-\frac{1}{2}}^{t+\Delta t} }{ p_{k+\frac{1}{2}} -
    p_{k-\frac{1}{2}} } , 
   \Deqlab{V鉛直拡散方程式離散表現}
\\
  \frac{ T^{t+\Delta t}_k - T^{t-\Delta t}_k }{ 2 \Delta t } 
    &=& \frac{1}{C_p} g \frac{ F_{h,k+\frac{1}{2}}^{t+\Delta t} -
    F_{h,k-\frac{1}{2}}^{t+\Delta t} }{ p_{k+\frac{1}{2}} -
    p_{k-\frac{1}{2}} } . 
   \Deqlab{T鉛直拡散方程式離散表現}
\end{eqnarray}
%
一方, 水蒸気の鉛直拡散に関しては, 最下層以外 ($k \ge 2$) では下のように離散化される. 
%
\begin{eqnarray}
  \frac{ q^{t+\Delta t}_k - q^{t-\Delta t}_k }{ 2 \Delta t } 
    &=& g \frac{ F_{q,k+\frac{1}{2}}^{t+\Delta t} -
    F_{q,k-\frac{1}{2}}^{t+\Delta t} }{ p_{k+\frac{1}{2}} -
    p_{k-\frac{1}{2}} } .
   \Deqlab{q鉛直拡散方程式離散表現(k>2)}
\end{eqnarray}
%
一方, 最下層 ($k = 1$) においては, 陰解法を用いて計算する場合の効率性を考慮し, 
2 つの離散化方法を用意している. 
%
1 つは,
%
\begin{eqnarray}
  \frac{ q^{t+\Delta t}_k - q^{t-\Delta t}_k }{ 2 \Delta t } 
    &=& g \frac{ F_{q,k+\frac{1}{2}}^{t+\Delta t} -
    F_{q,k-\frac{1}{2}}^{t+\Delta t} }{ p_{k+\frac{1}{2}} -
    p_{k-\frac{1}{2}} } \ \ \ (k=1)
\end{eqnarray}
%
であり, 1 つは, 
%
\begin{eqnarray}
  \frac{ q^{t+\Delta t}_k - q^{t-\Delta t}_k }{ 2 \Delta t } 
    &=& g \frac{ F_{q,k+\frac{1}{2}}^{t+\Delta t} -
    F_{q,k-\frac{1}{2}}^{t-\Delta t} }{ p_{k+\frac{1}{2}} -
    p_{k-\frac{1}{2}} } \ \ \ (k=1)
  \Deqlab{最下層水蒸気拡散離散表現 t-Δt版}
\end{eqnarray}
%
である. 前者の場合, 最下層の離散化方法は最下層以外の層 ($k \ge 2$) と同じように
離散化される. 
後者の場合, 惑星表面のフラックスのみ $t - \Delta t$ の時刻の値が使われる
%
\footnote{
  後者の方法を利用しなければいけないのは, 陰解法で離散化した結果を整理して
  得られる連立一次方程式の行列を三重対角行列にするため, そして, 有限の土壌
  水分を扱うためである. 

  地表面における上向き熱フラックスは, 大気側から見れば, 下部境界において
  大気に入る熱フラックスであり, この意味で, 大気中の熱収支は地表面および
  地下の土壌の熱収支と関係している. 
  さらに, 水蒸気が存在する系では, 地表面の熱収支は, 惑星表面における水蒸気
  の蒸発と凝結を介して水蒸気の収支とも関係している. 
  このため, 本来は, 熱の鉛直拡散, 惑星表面の熱収支, 地下の土壌の熱拡散, 
  水蒸気の鉛直拡散を陰解法で計算するためには, すべての方程式を連立して
  計算しなければならない. 
  素直に定式化すると, これらすべてを含む連立一次方程式の行列は三重対角
  行列にならず, 計算量が多くなってしまう. 
  三重対角行列にするためには, 熱の鉛直拡散, 地下の土壌の熱拡散 (惑星表面の
  熱収支を含む), 水蒸気の鉛直拡散のうちの一つを分離して解く必要があり, 
  現在の \Dmodel の定式化では, 水蒸気の鉛直拡散を分離して解くことにしている
  ($t - \Delta t$ の時刻の惑星表面の水蒸気フラックスを用いることで, 水蒸気
  の鉛直拡散は分離される). 

  また, 上では触れていないが, 本来は土壌水分量の収支も関係している. 
  しかし, 有限の土壌水分量を考える場合, 土壌が含む以上の量の水蒸気が蒸発する
  ことはないが, そのような条件を連立一次方程式に課すことは難しく, 現実的には
  それを連立して解くことはできない. このことも, 上で書いたように水蒸気の鉛直
  拡散を分離して解く理由である. 

  一方, 地下の土壌の熱拡散を計算しないモデルにおいては, 熱の鉛直拡散, 惑星
  表面の熱収支, 水蒸気の鉛直拡散を連立して得られる行列は三重対角行列になる
  ため, 問題は起こらない. これが前者の式が用いられる場合である. 
}.
%
なお, 水蒸気以外の熱収支に関わらない物質の鉛直拡散は, 
\Deqref{最下層水蒸気拡散離散表現 t-Δt版}
と同様に離散化する. 


拡散フラックスは下のように離散化される.
%
\begin{eqnarray}
  F_{m,x,k+\frac{1}{2}} 
    &=& - (TC)_{m,k+\frac{1}{2}} \left( u_{k+1} - u_{k} \right), 
%\\
%    &=& - (TC)_{m,k+\frac{1}{2}} u_{k+1} + (TC)_{m,k+\frac{1}{2}} u_k
\\
  F_{m,y,k+\frac{1}{2}} 
    &=& - (TC)_{m,k+\frac{1}{2}} \left( v_{k+1} - v_{k} \right), 
%\\
%    &=& - (TC)_{m,k+\frac{1}{2}} v_{k+1} + (TC)_{m,k+\frac{1}{2}} v_k
\\
  F_{h,k+\frac{1}{2}} 
    &=& - C_p P_{k+\frac{1}{2}} (TC)_{h,k+\frac{1}{2}}
          \left( \frac{ T_{k+1} }{ P_{k+1} } - \frac{ T_{k} }{ P_{k} } \right), 
%\\
%    &=& - C_p \frac{ P_{k+\frac{1}{2}} }{ P_{k+1} } (TC)_{h,k+\frac{1}{2}} T_{k+1} + C_p \frac{ P_{k+\frac{1}{2}} }{ P_{k} } (TC)_{h,k+\frac{1}{2}} T_k
\\
  F_{q,k+\frac{1}{2}} 
    &=& - (TC)_{q,k+\frac{1}{2}} \left( q_{k+1} - q_{k} \right).
%\\
%    &=& - (TC)_{q,k+\frac{1}{2}} q_{k+1} + (TC)_{q,k+\frac{1}{2}} q_k
\end{eqnarray}
ここで, $(TC)_{m,k+\frac{1}{2}}$, $(TC)_{h,k+\frac{1}{2}}$,
$(TC)_{q,k+\frac{1}{2}}$ 
は以下のように表現される
\footnote{コードのコメントでは, $(TC)$ に「輸送係数」と名前が付けられている.}.

上部境界では, 
\begin{eqnarray}
  (TC)_{m,k_{max}+\frac{1}{2}} &=& 0, 
\\
  (TC)_{h,k_{max}+\frac{1}{2}} &=& 0,
\\
  (TC)_{q,k_{max}+\frac{1}{2}} &=& 0.
\end{eqnarray}


$k = k_{max}$ のとき, 
%
\begin{eqnarray}
  F_{m,x,k_{max}+\frac{1}{2}} &=& 0, 
\\
  F_{m,y,k_{max}+\frac{1}{2}} &=& 0, 
\\
  F_{h,k_{max}+\frac{1}{2}} &=& 0, 
\\
  F_{q,kamx+\frac{1}{2}} &=& 0
\end{eqnarray}
%
となる.


$2 \le k \le k_{max}-1$ のとき, 
%
\begin{eqnarray}
  (TC)_{m,k+\frac{1}{2}} &=& \rho_{k+\frac{1}{2}} K_{m,k+\frac{1}{2}} \frac{1}{ z_{k+1} - z_k },
\\
  (TC)_{h,k+\frac{1}{2}} &=& \rho_{k+\frac{1}{2}} K_{h,k+\frac{1}{2}} \frac{1}{ z_{k+1} - z_k },
\\
  (TC)_{q,k+\frac{1}{2}} &=& \rho_{k+\frac{1}{2}} K_{q,k+\frac{1}{2}} \frac{1}{ z_{k+1} - z_k }.
\end{eqnarray}
$\rho_{k+\frac{1}{2}}$ は次式を用いて評価する.
\begin{eqnarray}
  \rho_{k+\frac{1}{2}} &=& \frac{ p_{k+\frac{1}{2}} }{ R T_{v,k+\frac{1}{2}} }
\end{eqnarray}
%
ここで $T_v$ は仮温度である.
%
$k = 1$ のとき, バルク法を用いてフラックスを評価する場合には, 
%
\begin{eqnarray}
  F_{m,x,k-\frac{1}{2}} &=& - (TC)_{m,k-\frac{1}{2}} u_1, 
    \Deqlab{vdiff-disc:sfcmomfluxx}
\\
  F_{m,y,k-\frac{1}{2}} &=& - (TC)_{m,k-\frac{1}{2}} v_1, 
    \Deqlab{vdiff-disc:sfcmomfluxy}
\\
%  F_{h,k-\frac{1}{2}} &=& - C_p (TC)_{h,k-\frac{1}{2}} \left( \frac{T_k}{P_{k}} - T_s \right)
%\\
%            &=& - C_p (TC)_{h,k-\frac{1}{2}} \frac{1}{P_{k}} T_k + C_p (TC)_{h,k-\frac{1}{2}} T_s
%\\
  F_{h,k-\frac{1}{2}} &=& - C_p P_{k-\frac{1}{2}} (TC)_{h,k-\frac{1}{2}}
   \left( \frac{T_k}{P_{k}} - \frac{T_s}{P_{k-\frac{1}{2}}} \right), 
    \Deqlab{vdiff-disc:sfcheatflux}
%\\
%            &=& - C_p (TC)_{h,k-\frac{1}{2}} \frac{P_{k-\frac{1}{2}}}{P_{k}} T_k + C_p (TC)_{h,k-\frac{1}{2}} T_s
\\
  F_{q,k-\frac{1}{2}} &=& - \epsilon (TC)_{q,k-\frac{1}{2}} \left( q_k - q_s^* \right).
%\\
%            &=& - \epsilon (TC)_{q,k-\frac{1}{2}} q_k + \epsilon (TC)_{q,k-\frac{1}{2}} q_s^*
\end{eqnarray}
%
\begin{eqnarray}
  (TC)_{m,k-\frac{1}{2}} &=& \rho_s C_d \left| \Dvect{v}_k \right|, 
\\
  (TC)_{h,k-\frac{1}{2}} &=& \rho_s C_h \left| \Dvect{v}_k \right|, 
\\
  (TC)_{q,k-\frac{1}{2}} &=& \rho_s C_q \left| \Dvect{v}_k \right|, 
\\
  \rho_s &=& \frac{ p_s }{ R T_{v,0} }.
\end{eqnarray}
%
である
\footnote
{
最後は $T_{v,0}$ (大気の温度) なのかね?
$T_s$ ($T_{s,v}$?) ではなくて?
たぶん, 考え方の問題だけ. どちらが悪いとも言えないだろうけど.
}.
%

また, 下部境界で摩擦の時定数, $\tau_f$, を与える場合には, 
\Deqref{vdiff-disc:sfcmomfluxx}, \Deqref{vdiff-disc:sfcmomfluxy} においては
%
\begin{eqnarray}
  (TC)_{m,k-\frac{1}{2}} &=& 
        - \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}} }{g} \frac{1}{\tau_f}, 
\end{eqnarray}
%
となる
\footnote{
  このとき, 下から二層目以上の(乱流)混合がなければ, 
%
  \begin{eqnarray}
    \left(\DP{u}{t}\right) &=& - \frac{1}{\tau_f} u
  \end{eqnarray}
%
  となる.
}. 
下部境界で温度を規定する場合には, 
\Deqref{vdiff-disc:sfcheatflux} において
%
\begin{eqnarray}
  (TC)_{h,k+\frac{1}{2}} &=& \rho_{k+\frac{1}{2}} K_{h,k+\frac{1}{2}} \frac{1}{ z_{k+1} - z_s },
\end{eqnarray}
%
とし, 下部境界で混合比値を規定する場合には, 
%
\begin{eqnarray}
  F_{q,k-\frac{1}{2}} &=& - (TC)_{q,k-\frac{1}{2}} \left( q_k - q_s \right).
\\
  (TC)_{q,k+\frac{1}{2}} &=& \rho_{k+\frac{1}{2}} K_{q,k+\frac{1}{2}} \frac{1}{ z_{k+1} - z_s }.
\end{eqnarray}
%
となる\footnote{ $z_s$ の $s$ は記号としては良くないかも. }.

また, 一定の熱フラックス, 物質フラックスを与える場合には, 
拡散フラックスは
%
\begin{eqnarray}
  F_{h,k-\frac{1}{2}} &=& F_{h,s}, 
\\
  F_{q,k-\frac{1}{2}} &=& F_{q,s}, 
\end{eqnarray}
%
となる.


\subsection{乱流運動エネルギー, 鉛直拡散係数 1 (Mellor and Yamada level 2) の離散表現}


% (2011-8-26 石渡) 以下の記述も必要
%  Ri の最小値が設定されていること.
% (d\Dvect{v}/dz)_{k+1/2} の最大値が設定さていること


鉛直拡散係数, $K_m$, $K_h$, $K_q$ は, それぞれ
\Deqref{Km連続系表現}, \Deqref{Kh連続系表現}, \Deqref{Kq連続系表現}
に示した式で計算する. 
そのために, リチャードソン数, 風速の鉛直シアー, 混合距離の離散表現が
必要となる. 
それらの表式は以下の通りである.

\Deqref{Ri定義} で定義したリチャードソン数
%
%\begin{eqnarray}
%  R_{i} &=& \frac{ \frac{g}{\theta}\frac{\partial \theta}{\partial z} }
%                 {\left| \frac{ \partial \Dvect{v}}{\partial z} \right|^2
%                 }
% \nonumber
%\end{eqnarray}
%
は, 地表面以外では下のように離散化する. 
%
\begin{eqnarray}
  R_{i,k+\frac{1}{2}} 
    &=& \frac{g}{\theta_{v,k+\frac{1}{2}}}
        \frac{\theta_{v,k+1} - \theta_{v,k}}{z_{k+1} - z_{k}}
        \left| \frac{ \partial \Dvect{v}}{\partial z}
	\right|_{k+\frac{1}{2}}^{-2}, 
\\
  \left| \frac{ \partial \Dvect{v}}{\partial z} \right|_{k+\frac{1}{2}}
    &=& \sqrt{ \left( \frac{ u_{k+1} - u_k }{ z_{k+1} - z_k } \right)^2
             + \left( \frac{ v_{k+1} - v_k }{ z_{k+1} - z_k } \right)^2
	     } .
\end{eqnarray}

混合距離 \Deqref{混合距離} は以下のように離散化する.
\begin{eqnarray}
  l_{k+\frac{1}{2}} = 
    \frac{k (z_{k+\frac{1}{2}} - z_{surf})}{1 + k (z_{k+\frac{1}{2}} - z_{surf})/l_0 } .
\end{eqnarray}
ここで, $z_{surf}$ は地表面高度である.


\subsection{乱流運動エネルギー, 鉛直拡散係数 2 (Mellor and Yamada level 2.5) の離散表現}

\iffalse

\begin{eqnarray}
  \DP{}{t}\left(\frac{q^2}{2}\right) 
    &=& \left\{ \DP{}{t} \left(\frac{q^2}{2}\right) \right\}_{adv}
      + \left\{ \DP{}{t} \left(\frac{q^2}{2}\right) \right\}_{vdiff}
      + \left\{ \DP{}{t} \left(\frac{q^2}{2}\right) \right\}_{src}
\\
  \left\{ \DP{}{t} \left(\frac{q^2}{2}\right) \right\}_{vdiff}
    &=& - \frac{1}{\rho} \DP{F_{TKE}}{z} 
\\
  \left\{ \DP{}{t} \left(\frac{q^2}{2}\right) \right\}_{src}
    &=& P_{s} + P_{b} - \epsilon_{TKE}
\\
  F_{TKE} &=& - \rho K_{TKE} \DP{}{z}\left( \frac{q^2}{2} \right)
\end{eqnarray}
%
まず最初に鉛直混合項のみによる時間変化率, 
$\left\{ \DP{}{t} \left(\frac{q^2}{2}\right) \right\}_{vdiff}$, 
を求める. 
次にこれを強制として, 生成項, 消滅項を含めた時間変化率, 
  $\left\{ \DP{}{t} \left(\frac{q^2}{2}\right) \right\}_{vdiff}
    + \left\{ \DP{}{t} \left(\frac{q^2}{2}\right) \right\}_{src}$, 
を求める. 
最後に, これらを強制として, 
移流項, $\left\{ \DP{}{t} \left(\frac{q^2}{2}\right) \right\}_{adv}$, を含めて
全時間変化率, $\DP{}{t}\left(\frac{q^2}{2}\right)$, を求める. 

鉛直混合による時間変化率, 
$\left\{ \DP{}{t} \left(\frac{q^2}{2}\right) \right\}_{vdiff}$, 
は, 時間に対して陰解法により求める. つまり, 
%
\begin{eqnarray}
  \left\{ \DP{}{t} \left(\frac{q^2}{2}\right) \right\}_{vdiff}
  &=&
  \frac{
      \left( \frac{q^2}{2} \right)^{t+\Delta t}_k - \left( \frac{q^2}{2} \right)^{t-\Delta t}_k 
  }{ 2 \Delta t } 
     =  g \frac{   F_{TKE,k+\frac{1}{2}}^{t+\Delta t}
                 - F_{TKE,k-\frac{1}{2}}^{t+\Delta t} }
               { p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}} }
%   \Deqlab{vdiff-disc:my2.5:TKEeq}
\end{eqnarray}


----------------------
\fi

乱流運動エネルギーの支配方程式のうち, 移流項を除いた部分は下のように
離散化される
\footnote{
  $P_s$, $P_b$ に対しては, どの時刻の値を使うかは選択による.
  $\epsilon_{TKE}$ に対しては, $t+\Delta t$ の値を使わないと上手く行かないと思われる.
}.
%
\begin{eqnarray}
  \frac{
      \left( \frac{q^2}{2} \right)^{t+\Delta t}_k - \left( \frac{q^2}{2} \right)^{t-\Delta t}_k 
  }{ 2 \Delta t } 
    &=& g \frac{   F_{TKE,k+\frac{1}{2}}^{t+\Delta t}
                 - F_{TKE,k-\frac{1}{2}}^{t+\Delta t} }
               { p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}} } \\ \nonumber
    & & + P_{s,k}^{t+\Delta t} + P_{b,k}^{t+\Delta t} - \epsilon_{TKE,k}^{t+\Delta t}.
   \Deqlab{vdiff-disc:my2.5:TKEeq}
\end{eqnarray}
%
ここで,
%
\begin{eqnarray}
  F_{TKE,k+\frac{1}{2}}
    &=& - (TC)_{TKE,k+\frac{1}{2}} 
            \left\{ \left( \frac{q^2}{2} \right)_{k+1} - \left( \frac{q^2}{2} \right)_k \right\}.
\\
  (TC)_{TKE,k+\frac{1}{2}} &=& \rho_{k+\frac{1}{2}} K_{TKE,k+\frac{1}{2}} \frac{1}{z_{k+1} - z_{k}},
\\
  P_{s,k}
    &=& K_{M,k} \left\{ \left( \DP{\Dvect{u}}{z} \right)^2 \right\}_k
\\
    &=& 2^\frac{1}{2} l_k S_{M,k} \left\{ \left( \DP{\Dvect{u}}{z} \right)^2 \right\}_k \left\{ \left( \frac{q^2}{2} \right)_k \right\}^\frac{1}{2}
\\
  P_{b,k}
    &=& - K_{H,k} \left( \frac{g}{\theta_v} \DP{\theta_v}{z} \right)_k 
\\
    &=& - 2^\frac{1}{2} l_k S_{H,k}
            \left( \frac{g}{\theta_v} \DP{\theta_v}{z} \right)_k
            \left\{ \left( \frac{q^2}{2} \right)_k \right\}^\frac{1}{2}
\\
  \epsilon_{TKE} &=& \frac{2^\frac{3}{2}}{B_1 l_k} \left( \frac{q^2}{2} \right)^\frac{3}{2}
\\
  K_{M,k} &=& 2^\frac{1}{2} l_k \left( \frac{q^2}{2} \right)^\frac{1}{2} S_{M,k}
\\
  K_{H,k} &=& 2^\frac{1}{2} l_k \left( \frac{q^2}{2} \right)^\frac{1}{2} S_{H,k}
\\
  K_{M,k+\frac{1}{2}} &=& \frac{1}{2} ( K_{M,k} + K_{M,k+1} )
\\
  K_{H,k+\frac{1}{2}} &=& \frac{1}{2} ( K_{H,k} + K_{H,k+1} )
\\
  K_{q,k+\frac{1}{2}} &=& K_{H,k+\frac{1}{2}}
\\
  K_{TKE,k+\frac{1}{2}} &=& K_{H,k+\frac{1}{2}}
\end{eqnarray}
%
\begin{eqnarray}
  \left\{ \left( \DP{\Dvect{u}}{z} \right)^2 \right\}_k
    &=&   \left( \frac{ u_{k+1} - u_{k  } }{ z_{k+1} - z_{k  } } \right)^2
        + \left( \frac{ v_{k+1} - v_{k  } }{ z_{k+1} - z_{k  } } \right)^2,
        \hspace{10mm}  k = 1
\\
  \left\{ \left( \DP{\Dvect{u}}{z} \right)^2 \right\}_k
    &=&   \left( \frac{ u_{k+1} - u_{k-1} }{ z_{k+1} - z_{k-1} } \right)^2
        + \left( \frac{ v_{k+1} - v_{k-1} }{ z_{k+1} - z_{k-1} } \right)^2,
        \hspace{10mm}  2 \le k \le k_{max}-1
\\
  \left\{ \left( \DP{\Dvect{u}}{z} \right)^2 \right\}_k
    &=&   \left( \frac{ u_{k  } - u_{k-1} }{ z_{k  } - z_{k-1} } \right)^2
        + \left( \frac{ v_{k  } - v_{k-1} }{ z_{k  } - z_{k-1} } \right)^2,
        \hspace{10mm}  k = k_{max}
\\
  \left( \frac{g}{\theta_v} \DP{\theta_v}{z} \right)_k
    &=& \frac{g}{\theta_{v,k}}
        \frac{\theta_{v,k+1} - \theta_{v,k  }}{z_{k+1} - z_{k  }}, \hspace{10mm}  k = 1
\\
  \left( \frac{g}{\theta_v} \DP{\theta_v}{z} \right)_k
    &=& \frac{g}{\theta_{v,k}}
        \frac{\theta_{v,k+1} - \theta_{v,k-1}}{z_{k+1} - z_{k-1}}, \hspace{10mm}  2 \le k \le k_{max}-1
\\
  \left( \frac{g}{\theta_v} \DP{\theta_v}{z} \right)_k
    &=& \frac{g}{\theta_{v,k}}
        \frac{\theta_{v,k  } - \theta_{v,k-1}}{z_{k  } - z_{k-1}}, \hspace{10mm}  k = k_{max}
\end{eqnarray}
%
となる
\footnote{
  実際には, 閉じるためには, 境界の値など他にも変数が必要になる.
}$^,$\footnote{
  (2013-08-13 高橋) 
  ここで示している離散表現は, 鉛直の格子配置の点で Mellor and Yamada level 2 
  の離散表現と異なる. 
  level 2, 2.5 の乱流モデルはともに, 鉛直層の境界における拡散係数を
  求めるために用いられる. 
  しかし, Mellor and Yamada level 2.5 では, 拡散係数を計算するために乱流
  運動エネルギーを計算する必要があり, 移流計算をすることを考えると, 現実
  的には乱流運動エネルギーの値は鉛直層の中心に配置せざるを得ない. 
  その結果, 何らかの内挿処理により鉛直層境界の拡散係数を計算する. 
  一方, Mellor and Yamada level 2 では, 乱流運動エネルギーの計算に移流過程
  が含まれず, 鉛直層の境界における乱流運動エネルギーを無理なく診断するこ
  とができる. 
  level 2 の場合も 2.5 の場合も同じように離散化する方が良いと思うが, 
  統一するためにわざわざ難しいことをやるのが良いかどうかを含めて検討が必要. 
}.

混合距離 \Deqref{混合距離} は以下のように離散化する.
\begin{eqnarray}
  l_k = 
    \frac{k (z_{k} - z_{surf})}{1 + k (z_{k} - z_{surf})/l_0 } .
  \Deqlab{vdiff-disc:MY2:l}
\end{eqnarray}
ここで, $z_{surf}$ は地表面高度である.


\Deqref{vdiff-disc:my2.5:TKEeq} は, さらに下のように整理される.
$2 \le k \le k_{max}-1$ のとき,
%
%
%%----- ここから消す
%ここから消す
%\begin{eqnarray}
%  & & \hspace{-1cm}
%       - (TC)_{TKE,k-\frac{1}{2}} \left\{ \left( \frac{q^2}{2} \right)_{k-1}^{t+\Delta t} - \left( \frac{q^2}{2} \right)_{k-1}^{t-\Delta t} \right\} \nonumber 
%\\
%    & & + \left\{
%            - \frac{1}{2 \Delta t} \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}}}{g}
%            + (TC)_{TKE,k+\frac{1}{2}} + (TC)_{TKE,k-\frac{1}{2}}
%          \right\}
%          \left\{ \left( \frac{q^2}{2} \right)_{k  }^{t+\Delta t}
%        - \left( \frac{q^2}{2} \right)_{k  }^{t-\Delta t} \right\} \nonumber
%\\
%    & & - (TC)_{TKE,k+\frac{1}{2}} 
%            \left\{ \left( \frac{q^2}{2} \right)_{k+1}^{t+\Delta t}
%          - \left( \frac{q^2}{2} \right)_{k+1}^{t-\Delta t} \right\} \nonumber
%\\
%    &=& - \left( F_{TKE,k+\frac{1}{2}}^{t-\Delta t} 
%               - F_{TKE,k-\frac{1}{2}}^{t-\Delta t} \right) 
%        - \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}}}{g}
%            \left( P_{s,k}^{t-\Delta t} + P_{b,k}^{t-\Delta t} \right)
%        + \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}}}{g}
%            \epsilon_{TKE,k}
%\end{eqnarray}
%ここまで消す
%%----- ここまで消す

%%----- ここから消す
%ここから消す
%\begin{eqnarray}
%  & & \hspace{-1cm}
%       - (TC)_{TKE,k-\frac{1}{2}} \left\{ \left( \frac{q^2}{2} \right)_{k-1}^{t+\Delta t} - \left( \frac{q^2}{2} \right)_{k-1}^{t-\Delta t} \right\} \nonumber 
%\\
%    & & + \left\{
%            - \frac{1}{2 \Delta t} \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}}}{g}
%            + (TC)_{TKE,k+\frac{1}{2}} + (TC)_{TKE,k-\frac{1}{2}}
%            - \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}}}{g} C_{d,2,k}
%          \right\} \nonumber
%\\
%    & & \hspace{9cm}
%        \times \left\{ \left( \frac{q^2}{2} \right)_{k  }^{t+\Delta t}
%                     - \left( \frac{q^2}{2} \right)_{k  }^{t-\Delta t} \right\} \nonumber
%\\
%    & & - (TC)_{TKE,k+\frac{1}{2}} 
%            \left\{ \left( \frac{q^2}{2} \right)_{k+1}^{t+\Delta t}
%          - \left( \frac{q^2}{2} \right)_{k+1}^{t-\Delta t} \right\} \nonumber
%\\
%    &=& - \left( F_{TKE,k+\frac{1}{2}}^{t-\Delta t} 
%               - F_{TKE,k-\frac{1}{2}}^{t-\Delta t} \right) \nonumber
%\\
%    & & - \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}}}{g}
%            \left( P_{s,k}^{t-\Delta t} + P_{b,k}^{t-\Delta t} \right)
%        + \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}}}{g}
%            C_{d,1,k}
%\end{eqnarray}
%ここまで消す
%%----- ここまで消す

\begin{eqnarray}
  & & \hspace{-1cm}
       - (TC)_{TKE,k-\frac{1}{2}} \left\{ \left( \frac{q^2}{2} \right)_{k-1}^{t+\Delta t} - \left( \frac{q^2}{2} \right)_{k-1}^{t-\Delta t} \right\} \nonumber 
\\
  & & + \left\{
          - \frac{1}{2 \Delta t} \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}}}{g}
          + (TC)_{TKE,k+\frac{1}{2}} + (TC)_{TKE,k-\frac{1}{2}} \right. \nonumber
\\
    & & \hspace{3cm} \left.
          + \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}}}{g}
             ( C_{s,2,k} + C_{b,2,k} - C_{d,2,k} )
        \right\} \nonumber
\\
    & & \hspace{6cm}
        \times \left\{ \left( \frac{q^2}{2} \right)_{k  }^{t+\Delta t}
                     - \left( \frac{q^2}{2} \right)_{k  }^{t-\Delta t} \right\} \nonumber
\\
    & & - (TC)_{TKE,k+\frac{1}{2}} 
            \left\{ \left( \frac{q^2}{2} \right)_{k+1}^{t+\Delta t}
          - \left( \frac{q^2}{2} \right)_{k+1}^{t-\Delta t} \right\} \nonumber
\\
    &=& - \left( F_{TKE,k+\frac{1}{2}}^{t-\Delta t} 
               - F_{TKE,k-\frac{1}{2}}^{t-\Delta t} \right) \nonumber
\\
    & & - \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}}}{g}
            \left( C_{s,1,k} + C_{b,1,k} - C_{d,1,k} \right)
\end{eqnarray}
%
ここで, 
%
\begin{eqnarray}
  C_{s,1,k}
%    &=&    2^\frac{3}{2} l_k \hat{S}_{M,k}
%             \left\{ \left( \DP{\Dvect{u}}{z} \right)^2 \right\}_k 
%             \left\{ \left( \frac{q^2}{2} \right)_k^{t-\Delta t} \right\}^\frac{3}{2} old
%\\
%    &=& 2^\frac{3}{2} l_k 
%             \left\{ \left( \DP{\Dvect{u}}{z} \right)^2 \right\}_k 
%             \hat{S}_{M,k}^{t-\Delta t}
%             \left\{ \left( \frac{q^2}{2} \right)_k^{t-\Delta t} \right\}^\frac{3}{2}
%\\
    &=& 2^\frac{1}{2} l_k 
             \left\{ \left( \DP{\Dvect{u}}{z} \right)^2 \right\}_k 
             S_{M,k}^{t-\Delta t}
             \left\{ \left( \frac{q^2}{2} \right)_k^{t-\Delta t} \right\}^\frac{1}{2}
\\
  C_{s,2,k}
%    &=&    \frac{3}{2} \cdot 2^\frac{3}{2} l_k \hat{S}_{M,k} 
%             \left\{ \left( \DP{\Dvect{u}}{z} \right)^2 \right\}_k 
%             \left\{ \left( \frac{q^2}{2} \right)_k^{t-\Delta t} \right\}^\frac{1}{2} old
%\\
%    &=&   2^\frac{3}{2} l_k 
%          \left\{ \left( \DP{\Dvect{u}}{z} \right)^2 \right\}_k 
%          \left[
%              \left\{ \DP{\hat{S}_{M}}{\left(\frac{q^2}{2}\right)} \right\}_k^{t-\Delta t}
%              \left\{ \left( \frac{q^2}{2} \right)_k^{t-\Delta t} \right\}^\frac{3}{2}
%            + \frac{3}{2} \hat{S}_{M,k}^{t-\Delta t}
%              \left\{ \left( \frac{q^2}{2} \right)_k^{t-\Delta t} \right\}^\frac{1}{2}
%          \right]
%\\
    &=& \frac{1}{2} \cdot 2^\frac{1}{2} l_k 
             \left\{ \left( \DP{\Dvect{u}}{z} \right)^2 \right\}_k 
             S_{M,k}^{t-\Delta t}
             \left\{ \left( \frac{q^2}{2} \right)_k^{t-\Delta t} \right\}^{-\frac{1}{2}}
\\
  C_{b,1,k}
%    &=& - 2^\frac{3}{2} l_k \hat{S}_{H,k}
%            \left( \frac{g}{\theta} \DP{\theta}{z} \right)_k
%            \left\{ \left( \frac{q^2}{2} \right)_k^{t-\Delta t} \right\}^\frac{3}{2} old
%\\
%    &=& - 2^\frac{3}{2} l_k 
%            \left( \frac{g}{\theta} \DP{\theta}{z} \right)_k
%            \hat{S}_{H,k}^{t-\Delta t}
%            \left\{ \left( \frac{q^2}{2} \right)_k^{t-\Delta t} \right\}^\frac{3}{2}
%\\
    &=& - 2^\frac{1}{2} l_k 
            \left( \frac{g}{\theta_v} \DP{\theta_v}{z} \right)_k
            S_{H,k}^{t-\Delta t}
            \left\{ \left( \frac{q^2}{2} \right)_k^{t-\Delta t} \right\}^\frac{1}{2}
\\
  C_{b,2,k}
%    &=& - \frac{3}{2} \cdot 2^\frac{3}{2} l_k \hat{S}_{H,k}
%            \left( \frac{g}{\theta} \DP{\theta}{z} \right)_k
%            \left\{ \left( \frac{q^2}{2} \right)_k^{t-\Delta t} \right\}^\frac{1}{2} old
%\\
%    &=& - 2^\frac{3}{2} l_k 
%            \left( \frac{g}{\theta} \DP{\theta}{z} \right)_k
%            \left[
%              \left\{ \DP{\hat{S}_{H}}{\left(\frac{q^2}{2}\right)} \right\}_k^{t-\Delta t}
%              \left\{ \left( \frac{q^2}{2} \right)_k^{t-\Delta t} \right\}^\frac{3}{2}
%              +  \frac{3}{2} 
%                 \hat{S}_{H,k}^{t-\Delta t}
%                 \left\{ \left( \frac{q^2}{2} \right)_k^{t-\Delta t} \right\}^\frac{1}{2}
%            \right]
%\\
    &=& - \frac{1}{2} \cdot 2^\frac{1}{2} l_k 
            \left( \frac{g}{\theta_v} \DP{\theta_v}{z} \right)_k
            S_{H,k}^{t-\Delta t}
            \left\{ \left( \frac{q^2}{2} \right)_k^{t-\Delta t} \right\}^{-\frac{1}{2}}
\\
  C_{d,1,k} &=& \frac{ 2^\frac{3}{2} }{B_1 l_k}
                \left\{ \left(\frac{q^2}{2}\right)^{t-\Delta t}_k \right\}^\frac{3}{2}
\\
  C_{d,2,k} &=& \frac{3}{2} \frac{ 2^\frac{3}{2} }{B_1 l_k}
                \left\{ \left(\frac{q^2}{2}\right)^{t-\Delta t}_k \right\}^\frac{1}{2}
\end{eqnarray}
%
である
\footnote{
  ここでは $C_{s,1}$, $C_{s,2}$, $C_{b,1}$, $C_{b,2}$ を
  示しているが, 少なくとも現状では
  $C_{s,2,k} = C_{b,2,k} = 0$ としている.
  これは $P_s$, $P_b$ として, $t-\Delta t$ の値を使っていること, 
  つまり, これらの項を陽的に扱うことに対応する.
  現状このようにする理由は, 線形化することで, この定式化の下で
  必ず正になる $P_s$ が負になることがあるためである.
}$^,$\footnote{
  $\epsilon_{TKE}$ を線形化した結果, $\epsilon_{TKE}$ 
  が負になることがある.
  この時は, $C_{d,*,k} = 0$ として解きなおす.
}$^,$\footnote{
$P_s$, $P_b$, $\epsilon_{TKE}$ は, 下のように時間に対して線形化
する (テイラー展開して一次の項までとる).
%
\begin{eqnarray}
  P_{s,k}^{t+\Delta t} 
    &=& 2^\frac{3}{2} l_k \hat{S}_{M,k}^{t+\Delta t} \left\{ \left( \DP{\Dvect{u}}{z} \right)^2 \right\}_k \left\{ \left( \frac{q^2}{2} \right)_k^{t+\Delta t} \right\}^\frac{3}{2}
\\
    &\sim& 2^\frac{3}{2} l_k \hat{S}_{M,k}^{t-\Delta t}
             \left\{ \left( \DP{\Dvect{u}}{z} \right)^2 \right\}_k 
             \left\{ \left( \frac{q^2}{2} \right)_k^{t-\Delta t} \right\}^\frac{3}{2} \nonumber
\\
    & & + 2^\frac{3}{2} l_k 
          \left\{ \left( \DP{\Dvect{u}}{z} \right)^2 \right\}_k 
          \left[
              \left\{ \DP{\hat{S}_{M}}{\left(\frac{q^2}{2}\right)} \right\}_k^{t-\Delta t}
              \left\{ \left( \frac{q^2}{2} \right)_k^{t-\Delta t} \right\}^\frac{3}{2}
            + \frac{3}{2} \hat{S}_{M,k}^{t-\Delta t}
              \left\{ \left( \frac{q^2}{2} \right)_k^{t-\Delta t} \right\}^\frac{1}{2}
          \right] \nonumber
\\
    & & \hspace{5cm} \times
          \left\{   \left( \frac{q^2}{2} \right)_k^{t+\Delta t} 
                  - \left( \frac{q^2}{2} \right)_k^{t-\Delta t} \right\}
\\
   &=& C_{s,1}
     + C_{s,2} \left\{ \left(\frac{q^2}{2}\right)^{t+\Delta t}_k 
                     - \left(\frac{q^2}{2}\right)^{t-\Delta t}_k \right\}
\\
%  old\ version ---
%\\
%  P_{s,k}^{t+\Delta t} 
%    &=& 2^\frac{3}{2} l_k \hat{S}_{M,k} \left\{ \left( \DP{\Dvect{u}}{z} \right)^2 \right\}_k \left\{ \left( \frac{q^2}{2} \right)_k \right\}^\frac{3}{2}
%\\
%    &\sim& 2^\frac{3}{2} l_k \hat{S}_{M,k}
%             \left\{ \left( \DP{\Dvect{u}}{z} \right)^2 \right\}_k 
%             \left\{ \left( \frac{q^2}{2} \right)_k^{t-\Delta t} \right\}^\frac{3}{2} \nonumber
%\\
%    & &  + \frac{3}{2} \cdot 2^\frac{3}{2} l_k \hat{S}_{M,k} 
%             \left\{ \left( \DP{\Dvect{u}}{z} \right)^2 \right\}_k 
%             \left\{ \left( \frac{q^2}{2} \right)_k^{t-\Delta t} \right\}^\frac{1}{2}
%             \left\{   \left( \frac{q^2}{2} \right)_k^{t+\Delta t} 
%                     - \left( \frac{q^2}{2} \right)_k^{t-\Delta t} \right\}
%\\
%   &=& C_{s,1}
%     + C_{s,2} \left\{ \left(\frac{q^2}{2}\right)^{t+\Delta t}_k 
%                     - \left(\frac{q^2}{2}\right)^{t-\Delta t}_k \right\}
%\\
%  --- old\ version
%\\
  P_{b,k}^{t+\Delta t} 
    &=& - 2^\frac{3}{2} l_k \hat{S}_{H,k}^{t+\Delta t}
            \left( \frac{g}{\theta_v} \DP{\theta_v}{z} \right)_k
            \left\{ \left( \frac{q^2}{2} \right)_k^{t+\Delta t} \right\}^\frac{3}{2}
\\
    &\sim& - 2^\frac{3}{2} l_k \hat{S}_{H,k}^{t-\Delta t}
               \left( \frac{g}{\theta_v} \DP{\theta_v}{z} \right)_k
               \left\{ \left( \frac{q^2}{2} \right)_k^{t-\Delta t} \right\}^\frac{3}{2} \nonumber
\\
    & &    - 2^\frac{3}{2} l_k 
               \left( \frac{g}{\theta_v} \DP{\theta_v}{z} \right)_k
               \left[
                 \left\{ \DP{\hat{S}_{H}}{\left(\frac{q^2}{2}\right)} \right\}_k^{t-\Delta t}
                 \left\{ \left( \frac{q^2}{2} \right)_k^{t-\Delta t} \right\}^\frac{3}{2}
                 +  \frac{3}{2} 
                    \hat{S}_{H,k}^{t-\Delta t}
                    \left\{ \left( \frac{q^2}{2} \right)_k^{t-\Delta t} \right\}^\frac{1}{2}
               \right] \nonumber
\\
   & & \hspace{5cm} \times
               \left\{   \left( \frac{q^2}{2} \right)_k^{t+\Delta t} 
                       - \left( \frac{q^2}{2} \right)_k^{t-\Delta t} \right\}
\\
   &=& C_{b,1}
     + C_{b,2} \left\{ \left(\frac{q^2}{2}\right)^{t+\Delta t}_k 
                     - \left(\frac{q^2}{2}\right)^{t-\Delta t}_k \right\}
%\\
%  old\ version ---
%\\
%  P_{b,k}^{t+\Delta t} 
%    &=& - 2^\frac{3}{2} l_k \hat{S}_{H,k}
%            \left( \frac{g}{\theta} \DP{\theta}{z} \right)_k
%            \left\{ \left( \frac{q^2}{2} \right)_k \right\}^\frac{3}{2}
%\\
%    &\sim& - 2^\frac{3}{2} l_k \hat{S}_{H,k}
%               \left( \frac{g}{\theta} \DP{\theta}{z} \right)_k
%               \left\{ \left( \frac{q^2}{2} \right)_k^{t-\Delta t} \right\}^\frac{3}{2} \nonumber
%\\
%    & &    - \frac{3}{2} \cdot 2^\frac{3}{2} l_k \hat{S}_{H,k}
%               \left( \frac{g}{\theta} \DP{\theta}{z} \right)_k
%               \left\{ \left( \frac{q^2}{2} \right)_k^{t-\Delta t} \right\}^\frac{1}{2}
%               \left\{   \left( \frac{q^2}{2} \right)_k^{t+\Delta t} 
%                       - \left( \frac{q^2}{2} \right)_k^{t-\Delta t} \right\}
%\\
%   &=& C_{b,1}
%     + C_{b,2} \left\{ \left(\frac{q^2}{2}\right)^{t+\Delta t}_k 
%                     - \left(\frac{q^2}{2}\right)^{t-\Delta t}_k \right\}
%\\
%  --- old\ version
\\
  \epsilon_{TKE,k}^{t+\Delta t}
    &=& \frac{ 2^\frac{3}{2} }{B_1 l_k}
          \left\{ \left(\frac{q^2}{2}\right)^{t+\Delta t}_k \right\}^\frac{3}{2}
\\
    &\sim& \frac{ 2^\frac{3}{2} }{B_1 l_k}
             \left\{ \left(\frac{q^2}{2}\right)^{t-\Delta t}_k \right\}^\frac{3}{2} \nonumber
\\
    & &  + \frac{3}{2} \frac{ 2^\frac{3}{2} }{B_1 l_k}
             \left\{ \left(\frac{q^2}{2}\right)^{t-\Delta t}_k \right\}^\frac{1}{2}
             \left\{ \left(\frac{q^2}{2}\right)^{t+\Delta t}_k 
                   - \left(\frac{q^2}{2}\right)^{t-\Delta t}_k \right\}
\\
    &=& C_{d,1}
      + C_{d,2} \left\{ \left(\frac{q^2}{2}\right)^{t+\Delta t}_k 
                      - \left(\frac{q^2}{2}\right)^{t-\Delta t}_k \right\}
\end{eqnarray}
}.

境界条件より上下境界では乱流運動エネルギーを固定するため, 下のようになる
\footnote{
  実際には, 下部境界では固定されない.
  下部境界における乱流運動エネルギーは摩擦速度に依存するが, 
  当然 $t+\Delta t$ における摩擦速度は $t-\Delta t$ における摩擦速度と異なる.
}$^,$\footnote{
  特に上部境界は扱いに不整合な点がある. 
  上部境界条件として乱流運動エネルギーの値を決めた上で, 上部境界での
  拡散係数をゼロと仮定すると, 境界条件の乱流運動エネルギーの値には意
  味がない. 
  別の言い方では, 境界条件として, 値を与えるのか微分値を与えるのかの
  問題. 両者を与えるとおかしなことになる.
  現状では, このような意味での整合性は取れていない. 
}.
%
$k = 1$ のとき,
%
\begin{eqnarray}
  & & \hspace{-1cm}
          \left\{
            - \frac{1}{2 \Delta t} \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}}}{g}
            + (TC)_{TKE,k+\frac{1}{2}} + (TC)_{TKE,k-\frac{1}{2}} \right. \nonumber
\\
    & & \hspace{3cm} \left.
            + \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}}}{g}
               ( C_{s,2,k} + C_{b,2,k} - C_{d,2,k} )
          \right\} \nonumber
\\
    & & \hspace{6cm}
        \times \left\{ \left( \frac{q^2}{2} \right)_{k  }^{t+\Delta t}
                     - \left( \frac{q^2}{2} \right)_{k  }^{t-\Delta t} \right\} \nonumber
\\
    & & - (TC)_{TKE,k+\frac{1}{2}} 
            \left\{ \left( \frac{q^2}{2} \right)_{k+1}^{t+\Delta t}
          - \left( \frac{q^2}{2} \right)_{k+1}^{t-\Delta t} \right\} \nonumber
\\
    &=& - \left( F_{TKE,k+\frac{1}{2}}^{t-\Delta t} 
               - F_{TKE,k-\frac{1}{2}}^{t-\Delta t} \right) \nonumber
\\
    & & - \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}}}{g}
            \left( C_{s,1,k} + C_{b,1,k} - C_{d,1,k} \right)
\end{eqnarray}
%
$k = k_{max}$ のとき,
%
\begin{eqnarray}
  & & \hspace{-1cm}
       - (TC)_{TKE,k-\frac{1}{2}} \left\{ \left( \frac{q^2}{2} \right)_{k-1}^{t+\Delta t} - \left( \frac{q^2}{2} \right)_{k-1}^{t-\Delta t} \right\} \nonumber 
\\
    & & + \left\{
            - \frac{1}{2 \Delta t} \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}}}{g}
            + (TC)_{TKE,k+\frac{1}{2}} + (TC)_{TKE,k-\frac{1}{2}} \right. \nonumber
\\
    & & \hspace{3cm} \left.
            + \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}}}{g}
               ( C_{s,2,k} + C_{b,2,k} - C_{d,2,k} )
          \right\} \nonumber
\\
    & & \hspace{6cm}
        \times \left\{ \left( \frac{q^2}{2} \right)_{k  }^{t+\Delta t}
                     - \left( \frac{q^2}{2} \right)_{k  }^{t-\Delta t} \right\} \nonumber
\\
    &=& - \left( F_{TKE,k+\frac{1}{2}}^{t-\Delta t} 
               - F_{TKE,k-\frac{1}{2}}^{t-\Delta t} \right) \nonumber
\\
    & & - \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}}}{g}
            \left( C_{s,1,k} + C_{b,1,k} - C_{d,1,k} \right)
\end{eqnarray}


これらをまとめると,
%
\begin{eqnarray}
    \Dvect{D} \Dvect{x}_{TKE} = \Dvect{G}_{TKE}
    \Deqlab{vdiff-disc:my2.5:tkediffeq}
\end{eqnarray}
%
と書くことができる.
ここで,
%
\begin{eqnarray}
  \Dvect{x}_{TKE} &=& \left(
      \left( \frac{q^2}{2} \right)_1^{t+\Delta t} - \left( \frac{q^2}{2} \right)_1^{t-\Delta t},
      \left( \frac{q^2}{2} \right)_2^{t+\Delta t} - \left( \frac{q^2}{2} \right)_2^{t-\Delta t}, 
      \cdots,
      \left( \frac{q^2}{2} \right)_{max}^{t+\Delta t} - \left( \frac{q^2}{2} \right)_{max}^{t-\Delta t} 
  \right)
\\
  \Dvect{G}_{TKE} &=& \left( g_{TKE,1}, g_{TKE,2}, \cdots, g_{TKE,k_{max}} \right),
\\
  g_{TKE,k}
    &=& - \left( F_{TKE,k+\frac{1}{2}}^{t-\Delta t} - F_{TKE,k-\frac{1}{2}}^{t-\Delta t} \right) \nonumber
\\
    & & - \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}}}{g}
            \left( C_{s,1,k} + C_{b,1,k} - C_{d,1,k} \right)
\end{eqnarray}
%
であり, $\Dvect{D}=(d_{m,n})$ の各成分は, $2 \le k \le k_{max}-1$ のとき, 
%
\begin{eqnarray}
  d_{k, k-1} &=& - (TC)_{TKE,k-\frac{1}{2}},
\\
  d_{k, k  } &=& - \frac{1}{2 \Delta t} \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}} }{g}
                 + (TC)_{TKE,k+\frac{1}{2}}
                 + (TC)_{TKE,k-\frac{1}{2}} \nonumber
\\
             & & + \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}}}{g}
                     ( C_{s,2,k} + C_{b,2,k} - C_{d,2,k} )
\\
  d_{k, k+1} &=& - (TC)_{TKE,k+\frac{1}{2}}.
\end{eqnarray}
%
$k = 1$ のとき,
%
\begin{eqnarray}
  d_{k, k  } &=&
            - \frac{1}{2 \Delta t} \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}}
}{g}
            + (TC)_{TKE,k+\frac{1}{2}}
            + (TC)_{TKE,k-\frac{1}{2}} \nonumber
\\           & &
            + \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}}}{g}
               ( C_{s,2,k} + C_{b,2,k} - C_{d,2,k} )
\\
  d_{k, k+1} &=& - (TC)_{TKE,k+\frac{1}{2}},
\end{eqnarray}
%
$k = k_{max}$ のとき
%
\begin{eqnarray}
  d_{k, k-1} &=& - (TC)_{TKE,k-\frac{1}{2}},
\\
  d_{k, k  } &=& - \frac{1}{2 \Delta t} \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}} }{g}
                 + (TC)_{TKE,k-\frac{1}{2}}
                 + (TC)_{TKE,k+\frac{1}{2}} \nonumber
\\
             & & + \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}}}{g}
                    ( C_{s,2,k} + C_{b,2,k} - C_{d,2,k} )
\end{eqnarray}
%
である.
%


\subsection{バルク係数 共通部分 (Louis et al., 1982; Beljaars and Holtslag, 1991) の離散表現}

バルク係数は, 
\Dsecref{vdiff-math:bulkcoef:L82}, \Dsecref{vdiff-math:bulkcoef:BH91} 
に示した式で計算する. 
そのために, 地表面のリチャードソン数の離散表現が必要となる.
その表式は以下の通りである.

\Deqref{Ri定義} で定義したリチャードソン数 
%
%\begin{eqnarray}
%  R_{i} &=& \frac{ \frac{g}{\theta}\frac{\partial \theta}{\partial z} }
%                 {\left| \frac{ \partial \Dvect{v}}{\partial z} \right|^2
%                 }
%  \nonumber
%\end{eqnarray}
%
は, 地表面においては, 下のように離散化する. 
%
\begin{eqnarray}
  R_{i,\frac{1}{2}} 
    &=& \frac{g}{\theta_{v,s}}
        \frac{\theta_{v,1} - \theta_{v,s}}{z_{k+1} - z_{s}}
        \left| \frac{ \partial \Dvect{v}}{\partial z}
	\right|_{\frac{1}{2}}^{-2}, 
\\
  \left| \frac{ \partial \Dvect{v}}{\partial z} \right|_{\frac{1}{2}}
    &=& \sqrt{ \left( \frac{ u_{k_1} - u_s }{ z_1 - z_s } \right)^2
             + \left( \frac{ v_{k_1} - v_s }{ z_1 - z_s } \right)^2 }, 
\\
  \theta_{v,s} &=& \frac{ T_{v,s} }{P_s}, 
\\
  P_s &=& \left( \frac{p_{00}}{p_s} \right)^\kappa .
\end{eqnarray}
%
ここで, $z_s$ は地表面の高度, $T_s$ は惑星表面温度, $p_s$ は惑星表面気圧
である
\footnote
{
  ここでは, $R_i$ の計算に惑星表面温度を用いているが, 惑星表面上の大気の温度
  を用いる方法もあるのかもしれない. 
  どちらが良いのかはよく分からない. 
}.


\subsection{バルク係数 2 (Beljaars and Holtslag, 1991; Beljaars, 1994) の離散表現}

\Dsecref{vdiff-math:bulkcoef:BH91} に示した Beljaars and Holtslag (1991), 
Beljaars (1994) の方法では, バルク係数は Monin-Obukhov 長さ, $L$, に依存する. 
しかし, $L$ は, 摩擦速度, 摩擦温度の関数であり, 
つまり, $L$ に依存する. 従って, 繰り返し法により $L$ を求める.
このとき, バルクリチャードソン数, $R_i$, を用いることにする.

$L$ と $R_i$ は, バルク係数, $C_d$, $C_h$ を用いて
%
\begin{eqnarray}
  L &=& \frac{1}{k R_i} \zeta \frac{C_m^\frac{3}{2}}{C_h}
\end{eqnarray}
%
のような関係にあるため, 下のように繰り返し法によって求める.
%
\begin{eqnarray}
  L^{n+1} &=& \frac{1}{k R_i} \zeta(L^n) \frac{C_m^\frac{3}{2}(L^n)}{C_h(L^n)}
\end{eqnarray}
%
ここで, 上付き添え字 $n$ は, $n$ 回の繰り返しによって得られた
値であることを示す. 繰り返し法で計算する際の初期値は $L^1 = 1$ としている.

また, 他に必要となる値は下のように離散化する
\footnote{
  \Deqref{vdiff-disc:beta w*} は, Beljaars (1994) の下の式, 
  %
  \begin{eqnarray}
    \overline{w'\theta'} 
      &=& \frac{k^2 \beta w_* (\theta_s - \theta_m)}
               {
                 \left\{ \log\left( -\frac{38.5L}{\gamma z_{0,m}} \right)
                         + \Psi_m\left(\frac{z_{0,m}}{L}\right) \right\}
                 \left\{ \log\left( -\frac{4L   }{\gamma z_{0,h}} \right)
                         + \Psi_m\left(\frac{z_{0,h}}{L}\right) \right\}
               }
  \\
    \overline{w'\theta'} 
      &=& b_h \left(\frac{g}{T}z_i\right)^\frac{1}{2}
            ( \theta_s - \theta_m )^\frac{3}{2}
  \\
    b_h
      &=& \frac{k^3 \beta^\frac{3}{2}}
               {
                 \left\{ \log\left( -\frac{38.5L}{\gamma z_{0,m}} \right)
                         + \Psi_m\left(\frac{z_{0,m}}{L}\right) \right\}^\frac{3}{2}
                 \left\{ \log\left( -\frac{4L   }{\gamma z_{0,h}} \right)
                         + \Psi_m\left(\frac{z_{0,h}}{L}\right) \right\}^\frac{3}{2}
               }
  \\
    z_i &=& - \frac{L}{k^2 \beta^3} 
                \left\{   \log\left(-\frac{38.5L}{\gamma z_{0,m}}\right)
                        + \Psi_m\left(\frac{z_{0,m}}{L}\right) \right\}^3
  \end{eqnarray}
  %
  から下のように導出した.
  %
  \begin{eqnarray}
    \beta w_*
      &=& \left[
            \frac{ \frac{g}{T} ( \theta_s - \theta_m ) k^2 \beta z_i }
                 {
                   \left\{ \log\left( -\frac{38.5L}{\gamma z_{0,m}} \right)
                           + \Psi_m\left(\frac{z_{0,m}}{L}\right) \right\}
                   \left\{ \log\left( -\frac{4L   }{\gamma z_{0,h}} \right)
                           + \Psi_m\left(\frac{z_{0,h}}{L}\right) \right\}
                 }
          \right]^\frac{1}{2}
  \end{eqnarray}
  %
  ただし, \Deqref{vdiff-disc:beta w*} では, $T$ を $\theta$ に置き換えている. 
  実質的には, 惑星表面付近の値なので差はほとんどないと考えられるが, 
  $( \theta_s - \theta_m )$ がかかっているため, $\theta$ にしておかないと, 
  参照気圧が惑星表面気圧から大きくずれているときに
  エクスナー関数分だけずれてしまう.
}.
%
\begin{eqnarray}
  |\Dvect{v_1}| &=& \left\{ u_1^2 + v_1^2 + \left( \beta w_* \right)^2 \right\}^\frac{1}{2}
\\
  \beta w_*
    &=& \left[
          \frac{ \frac{g}{\theta} ( \theta_s - \theta_m ) k^2 \beta z_{BL} }
               {
                 \left\{ \log\left( -\frac{38.5L}{\gamma z_{0,m}} \right)
                         + \Psi_m\left(\frac{z_{0,m}}{L}\right) \right\}
                 \left\{ \log\left( -\frac{4L   }{\gamma z_{0,h}} \right)
                         + \Psi_m\left(\frac{z_{0,h}}{L}\right) \right\}
               }
        \right]^\frac{1}{2}
  \Deqlab{vdiff-disc:beta w*}
\end{eqnarray}
%
ここで, $\beta = 1.2$, $z_{BL} = 1000 \ {\rm m}$ とする
\footnote{
  Beljaars (1994) の本文では下のように表現している.
%
  \begin{eqnarray}
    z_{BL} &=& - \frac{L}{k^2 \beta^3} 
                  \left\{   \log\left(-\frac{38.5L}{\gamma z_{0,m}}\right)
                          + \Psi_m\left(\frac{z_{0,m}}{L}\right) \right\}^3
  \end{eqnarray}
%
  しかし, $z_{BL}$ の値を固定してもほとんど影響がない 
  (e.g., ECMWF IFS documentation CY38r1)
  (該当箇所は, Part IV. Physical processes, 
    http://www.ecmwf.int/research/ifsdocs/CY38r1/IFSPart4.pdf
    の Section 3.2.1, p.36 である).
  また, 上の式を用いて計算すると, 計算が不安定になるようである.
}.



\subsection{運動量拡散の差分方程式の整理}

東西方向の運動量の鉛直拡散方程式 \Deqref{U鉛直拡散方程式離散表現} を整理すると, 
$2 \le k \le k_{max}-1$ のとき, 
%
\begin{eqnarray}
  & & \hspace{-1cm}
    - (TC)_{m,k-\frac{1}{2}} \left( u_{k-1}^{t+\Delta t} - u_{k-1}^{t-\Delta t} \right)
\nonumber 
\\
  & & + \left( 
        - \frac{1}{2 \Delta t} \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}} }{g}
        + (TC)_{m,k-\frac{1}{2}} 
        + (TC)_{m,k+\frac{1}{2}} 
      \right) \left( u_{k}^{t+\Delta t} - u_{k}^{t-\Delta t} \right)
\nonumber 
\\
  & & - (TC)_{m,k+\frac{1}{2}} \left( u_{k+1}^{t+\Delta t} - u_{k+1}^{t-\Delta t} \right) \nonumber
\\
  &=& - \left( F_{m,x,k+\frac{1}{2}}^{t-\Delta t} - F_{m,x,k-\frac{1}{2}}^{t-\Delta t} \right)
\end{eqnarray}
%
$k = 1$ のとき, 
%
\begin{eqnarray}
  & & \hspace{-1cm}
    \left( 
        - \frac{1}{2 \Delta t} \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}} }{g}
        + (TC)_{m,k-\frac{1}{2}} 
        + (TC)_{m,k+\frac{1}{2}} 
      \right) \left( u_{k}^{t+\Delta t} - u_{k}^{t-\Delta t} \right) \nonumber
\\
  & & - (TC)_{m,k+\frac{1}{2}} \left( u_{k+1}^{t+\Delta t} - u_{k+1}^{t-\Delta t} \right) \nonumber
\\
  &=& - \left( F_{m,x,k+\frac{1}{2}}^{t-\Delta t} - F_{m,x,k-\frac{1}{2}}^{t-\Delta t} \right)
\end{eqnarray}
%
$k = k_{max}$ のとき, 
%
\begin{eqnarray}
  & & \hspace{-1cm}
    - (TC)_{m,k-\frac{1}{2}} \left( u_{k-1}^{t+\Delta t} - u_{k-1}^{t-\Delta t} \right) \nonumber
\\
  & & + \left( 
        - \frac{1}{2 \Delta t} \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}} }{g}
        + (TC)_{m,k-\frac{1}{2}} 
      \right) \left( u_{k}^{t+\Delta t} - u_{k}^{t-\Delta t} \right) \nonumber
\\
  &=& - \left( F_{m,x,k+\frac{1}{2}}^{t-\Delta t} - F_{m,x,k-\frac{1}{2}}^{t-\Delta t} \right)
\end{eqnarray}
%
となる.

これらをまとめると, 
%
\begin{eqnarray}
    \Dvect{A} \Dvect{x}_u = \Dvect{G}_u
\end{eqnarray}
%
\begin{eqnarray}
  \Dvect{x}_u &=& \left( 
      u_1^{t+\Delta t} - u_1^{t-\Delta t}, 
      u_2^{t+\Delta t} - u_2^{t-\Delta t},
      \cdots, 
      u_{k_{max}}^{t+\Delta t} - u_{k_{max}}^{t-\Delta t} \right), 
\\
  \Dvect{G}_u &=& \left( g_{u,1}, g_{u,2}, \cdots, g_{u,k_{max}} \right), 
\\
  g_{u,k} &=& - \left( F_{m,x,k+\frac{1}{2}}^{t-\Delta t} - F_{m,x,k-\frac{1}{2}}^{t-\Delta t} \right)
\end{eqnarray}
%
ここで, $2 \le k \le k_{max}-1$ のとき, $\Dvect{A}=(a_{m,n})$ の各成分は, 
%
\begin{eqnarray}
  a_{k, k-1} &=& - (TC)_{m,k-\frac{1}{2}}, 
\\
  a_{k, k}   &=& - \frac{1}{2 \Delta t} 
                   \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}} }{g} 
                 + (TC)_{m,k-\frac{1}{2}} 
                 + (TC)_{m,k+\frac{1}{2}}, 
\\
  a_{k, k+1} &=& - (TC)_{m,k+\frac{1}{2}} .
\end{eqnarray}
%
$k = 1$ のとき, 
%
\begin{eqnarray}
  a_{k, k} &=& - \frac{1}{2 \Delta t} \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}} }{g} + (TC)_{m,k-\frac{1}{2}} + (TC)_{m,k+\frac{1}{2}},
\\
  a_{k, k+1} &=&                                                - (TC)_{m,k+\frac{1}{2}}.
\end{eqnarray}
%
$k = k_{max}$ のとき, 
%
\begin{eqnarray}
  a_{k, k-1} &=&                                                - (TC)_{m,k-\frac{1}{2}},
\\
  a_{k, k} &=& - \frac{1}{2 \Delta t} \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}} }{g} + (TC)_{m,k-\frac{1}{2}}
\end{eqnarray}
%
である. 

南北風に関しては, 東西風と同様に下のように書くことができる. 
%
\begin{eqnarray}
    \Dvect{A} \Dvect{x}_v = \Dvect{G}_v
\end{eqnarray}
%
\begin{eqnarray}
  \Dvect{x}_v &=& \left( 
      v_1^{t+\Delta t} - v_1^{t-\Delta t}, 
      v_2^{t+\Delta t} - v_2^{t-\Delta t},
      \cdots, 
      v_{k_{max}}^{t+\Delta t} - v_{k_{max}}^{t-\Delta t} \right), 
\\
  \Dvect{G}_v &=& \left( g_{v,1}, g_{v,2}, \cdots, g_{v,k_{max}} \right), 
\\
  g_{v,k} &=& - \left( F_{m,y,k+\frac{1}{2}}^{t-\Delta t} - F_{m,y,k-\frac{1}{2}}^{t-\Delta t} \right)
\end{eqnarray}
%
である. 


\subsection{熱拡散の差分方程式の整理}

熱の鉛直拡散の式 \Deqref{T鉛直拡散方程式離散表現} を整理すると, 
$2 \le k \le k_{max}-1$ のとき, 
%
\begin{eqnarray}
  & & \hspace{-2cm}
  - C_p \frac{ P_{k-\frac{1}{2}} }{ P_{k-1} } (TC)_{h,k-\frac{1}{2}} \left( T_{k-1}^{t+\Delta t} - T_{k-1}^{t-\Delta t} \right) \nonumber
\\
    & & + \left( 
            - C_p \frac{1}{2 \Delta t} \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}} }{g}
            + C_p \frac{ P_{k-\frac{1}{2}} }{ P_{k} } (TC)_{h,k-\frac{1}{2}} 
            + C_p \frac{ P_{k+\frac{1}{2}} }{ P_{k} } (TC)_{h,k+\frac{1}{2}} 
          \right) \left( T_{k}^{t+\Delta t} - T_{k}^{t-\Delta t} \right) \nonumber
\\
    & & - C_p \frac{ P_{k+\frac{1}{2}} }{ P_{k+1} } (TC)_{h,k+\frac{1}{2}} \left( T_{k+1}^{t+\Delta t} - T_{k+1}^{t-\Delta t} \right) \nonumber
\\
    &=& - \left( F_{h,k+\frac{1}{2}}^{t-\Delta t} - F_{h,k-\frac{1}{2}}^{t-\Delta t} \right)
\end{eqnarray}
%
のとき, $k = 1$ のとき, バルク法でフラックスを評価する場合には, 
%
\begin{eqnarray}
  & & \hspace{-2cm}
        - C_p (TC)_{h,k-\frac{1}{2}} \left( T_s^{t+\Delta t} - T_s^{t-\Delta t} \right) \nonumber
\\
  & &   + \left( 
            - C_p \frac{1}{2 \Delta t} \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}} }{g}
            + C_p \frac{ P_{k+\frac{1}{2}} }{ P_{k  } } (TC)_{h,k+\frac{1}{2}} 
            + C_p \frac{P_{k-\frac{1}{2}}}{P_{k}}  (TC)_{h,k-\frac{1}{2}}
          \right) \left( T_{k}^{t+\Delta t} - T_{k}^{t-\Delta t} \right) \nonumber
\\
    & & - C_p \frac{ P_{k+\frac{1}{2}} }{ P_{k+1} } (TC)_{h,k+\frac{1}{2}} \left( T_{k+1}^{t+\Delta t} - T_{k+1}^{t-\Delta t} \right) \nonumber
\\
    &=& - \left( F_{h,k+\frac{1}{2}}^{t-\Delta t} - F_{h,k-\frac{1}{2}}^{t-\Delta t} \right)
\end{eqnarray}
%
下部境界での温度を規定する場合には, 
%
\begin{eqnarray}
  & & \hspace{-2cm}
          \left( 
            - C_p \frac{1}{2 \Delta t} \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}} }{g}
            + C_p \frac{ P_{k+\frac{1}{2}} }{ P_{k  } } (TC)_{h,k+\frac{1}{2}} 
            + C_p \frac{P_{k-\frac{1}{2}}}{P_{k}}  (TC)_{h,k-\frac{1}{2}}
          \right) \left( T_{k}^{t+\Delta t} - T_{k}^{t-\Delta t} \right) \nonumber
\\
    & & - C_p \frac{ P_{k+\frac{1}{2}} }{ P_{k+1} } (TC)_{h,k+\frac{1}{2}} \left( T_{k+1}^{t+\Delta t} - T_{k+1}^{t-\Delta t} \right) \nonumber
\\
    &=& - \left( F_{h,k+\frac{1}{2}}^{t-\Delta t} - F_{h,k-\frac{1}{2}}^{t-\Delta t} \right)
\end{eqnarray}
%
一定値の熱フラックスを与える場合には, 
%
\begin{eqnarray}
  & & \hspace{-2cm}
          \left( 
            - C_p \frac{1}{2 \Delta t} \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}} }{g}
            + C_p \frac{ P_{k+\frac{1}{2}} }{ P_{k  } } (TC)_{h,k+\frac{1}{2}} 
          \right) \left( T_{k}^{t+\Delta t} - T_{k}^{t-\Delta t} \right) \nonumber
\\
    & & - C_p \frac{ P_{k+\frac{1}{2}} }{ P_{k+1} } (TC)_{h,k+\frac{1}{2}} \left( T_{k+1}^{t+\Delta t} - T_{k+1}^{t-\Delta t} \right) \nonumber
\\
    &=& - \left( F_{h,k+\frac{1}{2}}^{t-\Delta t} - F_{h,k-\frac{1}{2}}^{t-\Delta t} \right)
\end{eqnarray}
%
となる. また, $k = k_{max}$ のとき, 
%
\begin{eqnarray}
  & & \hspace{-2cm}
  - C_p \frac{ P_{k-\frac{1}{2}} }{ P_{k-1} } (TC)_{h,k-\frac{1}{2}} \left( T_{k-1}^{t+\Delta t} - T_{k-1}^{t-\Delta t} \right) \nonumber
\\
    & & + \left( 
            - C_p \frac{1}{2 \Delta t} \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}} }{g}
            + C_p \frac{ P_{k-\frac{1}{2}} }{ P_{k} } (TC)_{h,k-\frac{1}{2}} 
          \right) \left( T_{k}^{t+\Delta t} - T_{k}^{t-\Delta t} \right) \nonumber
\\
    &=& - \left( F_{h,k+\frac{1}{2}}^{t-\Delta t} - F_{h,k-\frac{1}{2}}^{t-\Delta t} \right)
\end{eqnarray}
%
となる. 

% (2011-9-2 石渡)
%   以下, \Dvect{B}_a \Dvect{x}_a = \Dvect{G}_a
%   を \Dvect{B} \Dvect{x}_h = \Dvect{G}_h
%   と変更した.  b_{a,k, k} だったものを b_{k, k} にした.

これらをまとめると, 惑星表面におけるフラックスをバルク法で評価する場合には, 
%
\begin{eqnarray}
    \Dvect{B}_a \Dvect{x}_a = \Dvect{G}_a
    \label{eq:disc:heatdiff}
\end{eqnarray}
%
\begin{eqnarray}
\Deqlab{vdiff-disc:explanation of array and vector in simul. lin. eq. for heat 1}
  \Dvect{x}_h &=& \left( 
      T_s^{t+\Delta t} - T_s^{t-\Delta t}, 
      T_1^{t+\Delta t} - T_1^{t-\Delta t}, 
      T_2^{t+\Delta t} - T_2^{t-\Delta t},
      \cdots, 
      T_{k_{max}}^{t+\Delta t} - T_{k_{max}}^{t-\Delta t} \right), 
\\
\Deqlab{vdiff-disc:explanation of array and vector in simul. lin. eq. for heat 2}
  \Dvect{G}_a &=& \left( g_{h,1}, g_{h,2}, \cdots, g_{h,k_{max}} \right), 
\\
\Deqlab{vdiff-disc:explanation of array and vector in simul. lin. eq. for heat 3}
  g_{h,k} &=& - \left( F_{a,k+\frac{1}{2}}^{t-\Delta t} - F_{a,k-\frac{1}{2}}^{t-\Delta t} \right)
\end{eqnarray}
%
と書くことができる. 
\footnote{$\Dvect{B}_a$ の下つき添字の $a$ は「大気」を表すラベルである.
          \Dchapref{惑星表面・地下の熱収支} では地表面の熱収支を
          扱い, そこでは $\Dvect{B}_s$ を用いる.}
ここで, $2 \le k \le k_{max}-1$ のとき, $\Dvect{B}_a=(b_{a,m,n})$ の各成分は, 
%
\begin{eqnarray}
  b_{a, k, k-1} &=& - C_p \frac{ P_{k-\frac{1}{2}} }{ P_{k-1} } (TC)_{h,k-\frac{1}{2}},
  \Deqlab{b_{k,k-1}一般形}
\\
  b_{a, k, k} &=& - C_p \frac{1}{2 \Delta t} \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}} }{ g } 
                     + C_p \frac{ P_{k-\frac{1}{2}} }{ P_k } (TC)_{h,k-\frac{1}{2}}
                     + C_p \frac{ P_{k+\frac{1}{2}} }{ P_k } (TC)_{h,k+\frac{1}{2}},
\\
  b_{a, k, k+1} &=& - C_p \frac{ P_{k+\frac{1}{2}} }{ P_{k+1} } (TC)_{h,k+\frac{1}{2}}
  \Deqlab{b_{k,k+1}一般形}
\end{eqnarray}
%
であり, $k = 1$ のとき, 
%
\begin{eqnarray}
  b_{a, k, k-1} &=& - C_p (TC)_{h,k-\frac{1}{2}},
\\
  b_{a, k, k} &=& - C_p \frac{1}{2 \Delta t} \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}} }{ g } 
               + C_p \frac{ P_{k+\frac{1}{2}} }{ P_k } (TC)_{h,k+\frac{1}{2}}
               + C_p \frac{P_{k-\frac{1}{2}}}{P_{k}} (TC)_{h,k-\frac{1}{2}},
\\
  b_{a, k, k+1} &=& - C_p \frac{ P_{k+\frac{1}{2}} }{ P_{k+1} } (TC)_{h,k+\frac{1}{2}}
\end{eqnarray}
%
であり, $k = k_{max}$ のとき, 
%
\begin{eqnarray}
  b_{a, k, k-1} &=& - C_p \frac{ P_{k-\frac{1}{2}} }{ P_{k-1} } (TC)_{h,k-\frac{1}{2}},
\\
  b_{a, k, k} &=& - C_p \frac{1}{2 \Delta t} \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}} }{ g } 
               + C_p \frac{ P_{k-\frac{1}{2}} }{ P_k } (TC)_{h,k-\frac{1}{2}}
\end{eqnarray}
%
である. 

ここで, $\Dvect{B}_a$ は $k_{max}$ 行 $k_{max}+1$ 列の行列であり, この式だけでは
未知数が方程式数よりも多いために閉じない. 
方程式を閉じるために, 以下に述べる惑星表面での熱収支式や地下の熱収支式, 
もしくは水蒸気の式を用いる. 
これらの式とあわせて同時に解く際に用いる行列の形式に関しては, 
\Dchapref{熱収支を統合した連立方程式の構成}
を参照せよ.



また, 惑星表面におけるフラックスに一定値を与える場合には, 
同じように式を変形して整理すると, $k = 1$ のとき, 
%
\begin{eqnarray}
  b_{a, k, k} &=& - C_p \frac{1}{2 \Delta t} \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}} }{ g } 
               + C_p \frac{ P_{k+\frac{1}{2}} }{ P_k } (TC)_{h,k+\frac{1}{2}},
\\
  b_{a, k, k+1} &=& - C_p \frac{ P_{k+\frac{1}{2}} }{ P_{k+1} } (TC)_{h,k+\frac{1}{2}}
\end{eqnarray}
となる. $k > 1$ の場合には, 
\Deqref{b_{k,k-1}一般形}$\sim$\Deqref{b_{k,k+1}一般形}
と同じである. 
この場合には, $\Dvect{B}_a$ は $k_{max}$ 行 $k_{max}$ 列の行列であり, 未知数が
方程式数と等しいため, この式のみで解くことができる. 



\subsection{水蒸気 (物質) 拡散の差分方程式の整理}

ここでは, 水蒸気の鉛直拡散の式の離散化方程式を整理する. 

\Dsecref{乱流過程離散表現} の最初の部分で
述べたように, 水蒸気の鉛直拡散は, 用いる惑星表面の水蒸気フラックスの時刻に
よって 2 通りの離散化方法を用いる. 

惑星表面の水蒸気フラックスとして $t + \Delta t$ の時刻の値を用いる場合, 
水蒸気の鉛直拡散の式 \Deqref{q鉛直拡散方程式離散表現(k>2)} を整理すると, 
$2 \le k \le k_{max}-1$ のとき, 
%
\begin{eqnarray}
  & & \hspace{-1cm}
       - (TC)_{q,k-\frac{1}{2}} \left( q_{k-1}^{t+\Delta t} - q_{k-1}^{t-\Delta t} \right) \nonumber
\\
    & & + \left( 
            - \frac{1}{2 \Delta t} \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}} }{g}
            + (TC)_{q,k+\frac{1}{2}} + (TC)_{q,k-\frac{1}{2}} 
          \right) 
          \left( q_{k  }^{t+\Delta t} - q_{k  }^{t-\Delta t} \right) \nonumber
\\
    & & - (TC)_{q,k+\frac{1}{2}} \left( q_{k+1}^{t+\Delta t} - q_{k+1}^{t-\Delta t} \right) \nonumber
\\
    &=& - \left( F_{q,k+\frac{1}{2}}^{t-\Delta t} - F_{q,k-\frac{1}{2}}^{t-\Delta t} \right)
\end{eqnarray}
%
となり, $k = 1$ のとき, 
%
\begin{eqnarray}
  & & \hspace{-1cm}
        - \epsilon (TC)_{q,k-\frac{1}{2}}
          \frac{\partial q_s^*}{\partial T_s} \left( T_s^{t+\Delta t} - T_s^{t-\Delta t} \right) \nonumber
\\
    & & + \left(
            - \frac{1}{2 \Delta t} \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}} }{g}
            + (TC)_{q,k+\frac{1}{2}} 
            + \epsilon (TC)_{q,k-\frac{1}{2}} 
          \right)
          \left( q_{k  }^{t+\Delta t} - q_{k  }^{t-\Delta t} \right) \\
    & & - (TC)_{q,k+\frac{1}{2}} \left( q_{k+1}^{t+\Delta t} - q_{k+1}^{t-\Delta t} \right) \nonumber
\\
    &=& - \left(  F_{q,k+\frac{1}{2}}^{t-\Delta t} - F_{q,k-\frac{1}{2}}^{t-\Delta t} \right)
\end{eqnarray}
%
下部境界の混合比を規定する場合には, 
%
\begin{eqnarray}
  & & \hspace{-1cm}
          \left(
            - \frac{1}{2 \Delta t} \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}} }{g}
            + (TC)_{q,k+\frac{1}{2}} 
            + (TC)_{q,k-\frac{1}{2}} 
          \right)
          \left( q_{k  }^{t+\Delta t} - q_{k  }^{t-\Delta t} \right) \\
    & & - (TC)_{q,k+\frac{1}{2}} \left( q_{k+1}^{t+\Delta t} - q_{k+1}^{t-\Delta t} \right) \nonumber
\\
    &=& - \left(  F_{q,k+\frac{1}{2}}^{t-\Delta t} - F_{q,k-\frac{1}{2}}^{t-\Delta t} \right)
\end{eqnarray}
%
となり, 
%
%メモ書き:
%惑星表面フラックスとして $t - \Delta t$ の時刻の値を用いる場合には, 
%左辺第一項の係数はゼロ, 左辺第二項内の第三項 ($\epsilon$ を含む項) はゼロである. 
%
$k = k_{max}$ のとき, 
%
\begin{eqnarray}
  & & \hspace{-1cm}
       - (TC)_{q,k-\frac{1}{2}} \left( q_{k-1}^{t+\Delta t} - q_{k-1}^{t-\Delta t} \right) \nonumber
\\
    & & + \left( 
            - \frac{1}{2 \Delta t} \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}} }{g}
            + (TC)_{q,k-\frac{1}{2}} 
          \right) 
          \left( q_{k  }^{t+\Delta t} - q_{k  }^{t-\Delta t} \right) \\
    &=& - \left( F_{q,k+\frac{1}{2}}^{t-\Delta t} - F_{q,k-\frac{1}{2}}^{t-\Delta t} \right)
\end{eqnarray}
%
となる. 

これらをまとめると, 
%
\begin{eqnarray}
    \Dvect{C} \Dvect{x}_q = \Dvect{G}_q
    \label{eq:disc:watervapordiff}
\end{eqnarray}
%
と書くことができる. 
ここで, 
%
\begin{eqnarray}
  \Dvect{x}_q &=& \left( 
      T_s^{t+\Delta t} - T_s^{t-\Delta t}, 
      q_1^{t+\Delta t} - q_1^{t-\Delta t}, 
      q_2^{t+\Delta t} - q_2^{t-\Delta t},
      \cdots, 
      q_{k_{max}}^{t+\Delta t} - q_{k_{max}}^{t-\Delta t} \right), 
\\
  \Dvect{G}_q &=& \left( g_{q,1}, g_{q,2}, \cdots, g_{q,k_{max}} \right), 
\\
  g_{q,k} &=& - \left( F_{q,k+\frac{1}{2}}^{t-\Delta t} - F_{q,k-\frac{1}{2}}^{t-\Delta t} \right)
\end{eqnarray}
%
であり, $2 \le k \le k_{max}-1$ のとき, $\Dvect{C}=(c_{m,n})$ の各成分は, 
%
\begin{eqnarray}
  c_{k, k-1} &=& - (TC)_{q,k-\frac{1}{2}}, 
  \Deqlab{c_{k,k-1}t+Δt版}
\\
  c_{k, k  } &=& - \frac{1}{2 \Delta t} \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}} }{g}
                 + (TC)_{q,k+\frac{1}{2}} 
                 + (TC)_{q,k-\frac{1}{2}}, 
\\
  c_{k, k+1} &=& - (TC)_{q,k+\frac{1}{2}}.
  \Deqlab{c_{k,k+1}t+Δt版}
\end{eqnarray}
%
$k = 1$ のとき, 
%
\begin{eqnarray}
  c_{k, k-1} &=& - \epsilon (TC)_{q,k-\frac{1}{2}} \frac{\partial q_s^*}{\partial T_s},
\\
  c_{k, k  } &=& 
            - \frac{1}{2 \Delta t} \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}} }{g}
            + (TC)_{q,k+\frac{1}{2}} 
            + \epsilon (TC)_{q,k-\frac{1}{2}}, 
\\
  c_{k, k+1} &=& - (TC)_{q,k+\frac{1}{2}}, 
\end{eqnarray}
%
$k = k_{max}$ のとき
%
\begin{eqnarray}
  c_{k, k-1} &=& - (TC)_{q,k-\frac{1}{2}}, 
\\
  c_{k, k  } &=& - \frac{1}{2 \Delta t} \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}} }{g} 
                 + (TC)_{q,k-\frac{1}{2}}
\end{eqnarray}
%
である. 

ここで, \Dvect{C} は $k_{max}$ 行 $k_{max}+1$ 列の行列であり, この式だけでは
未知数が方程式数よりも多いために閉じない. 
方程式を閉じるために, 熱の鉛直拡散の式や惑星表面での熱収支式や地下の熱収支式を
同時に解く. 
同時に解く際に用いる行列の形式に関しては, 
\Dchapref{熱収支を統合した連立方程式の構成}
を参照せよ.


なお, 惑星表面フラックスとして $t - \Delta t$ の時刻の値を用いる場合には, 
同じように式を変形して整理すると, 
$k = 1$ のとき, 
%
\begin{eqnarray}
  c_{k, k-1} &=& 0, 
\\
  c_{k, k  } &=& 
            - \frac{1}{2 \Delta t} \frac{ p_{k+\frac{1}{2}} - p_{k-\frac{1}{2}} }{g}
            + (TC)_{q,k+\frac{1}{2}} ,
\\
  c_{k, k+1} &=& - (TC)_{q,k+\frac{1}{2}}
\end{eqnarray}
%
となる. $k \ge 2$ においては,  
\Deqref{c_{k,k-1}t+Δt版} 〜\Deqref{c_{k,k+1}t+Δt版} と同様である. 
この場合には, \Dvect{C} は $k_{max}$ 行 $k_{max}$ 列の行列であり, この式だけで
閉じる. 

なお, 惑星表面フラックスとして一定値を用いる場合にも同様の方法で解くことができる.