% 表題   dcpam5  乱流過程
%
% 履歴
%\Drireki{2010/09/01 竹広真一}
%\Drireki{2010/04/15 高橋芳幸}
%
%  \Dchapterhead
\section{離散表現}

dcpam では, 鉛直拡散は陰解法を用いて計算する. 
運動量, 熱の鉛直拡散方程式は下のように離散化する. 
%
\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}} }
\\
  \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}} }
\\
  \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}} }
\end{eqnarray}
%

拡散フラックスは下のように離散化される.
%
\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
\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 }
\end{eqnarray}
%
$k = 1$ のとき, 
%
\begin{eqnarray}
  F_{m,x,k-\frac{1}{2}} &=& - (TC)_{m,k-\frac{1}{2}} u_1
\\
  F_{m,y,k-\frac{1}{2}} &=& - (TC)_{m,k-\frac{1}{2}} v_1
\\
%  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)
%\\
%            &=& - 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
\\
\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|
\\
  \rho_s &=& \frac{ p_s }{ R T_0 }
\end{eqnarray}
%
であり
\footnote
{
最後は $T_0$ (大気の温度) なのかね?
$T_s$ ではなくて?
たぶん, 考え方の問題だけ. どちらが悪いとも言えないだろうけど.
}, 
%
$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
\end{eqnarray}
%
となる.




\subsection{鉛直拡散係数}


鉛直拡散係数, $K_m$, $K_h$, は ... に示した式で計算する. 
ここでは, リチャードソン数と風速の鉛直シアーの離散表現を示すのみとする. 

リチャードソン数
%
\begin{eqnarray}
  R_{i} &=& \frac{ \frac{g}{\theta}\frac{\partial \theta}{\partial z} }
                 {\left| \frac{ \partial \Dvect{v}}{\partial z} \right| }
\end{eqnarray}
%
は, 下のように離散化する. 
%
\begin{eqnarray}
  R_{i,k+\frac{1}{2}} 
    &=& \frac{g}{\theta_{k+\frac{1}{2}}}
        \frac{\theta_{k+1} - \theta_{k}}{z_{k+1} - z_{k}}
        \left| \frac{ \partial \Dvect{v}}{\partial z} \right|_{k+\frac{1}{2}}^{-1}
\\
  \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}


\subsection{バルク係数}

バルク係数は, ... に示した式で計算する. ここでは, 地表面のリチャードソン数の
離散表現を示すのみとする. 

リチャードソン数
%
\begin{eqnarray}
  R_{i} &=& \frac{ \frac{g}{\theta}\frac{\partial \theta}{\partial z} }
                 {\left| \frac{ \partial \Dvect{v}}{\partial z} \right| }
\end{eqnarray}
%
は, 地表面においては, 下のように離散化する. 
%
\begin{eqnarray}
  R_{i,\frac{1}{2}} 
    &=& \frac{g}{\theta_{s}}
        \frac{\theta_{1} - \theta_s}{z_{k+1} - z_{s}}
        \left| \frac{ \partial \Dvect{v}}{\partial z} \right|_{\frac{1}{2}}^{-1}
\\
  \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_s &=& \frac{ T_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 \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},
      ..., 
      u_{k_{max}}^{t+\Delta t} - u_{k_{max}}^{t-\Delta t} \right), 
\\
  \Dvect{G}_u &=& \left( g_{u,1}, g_{u,2}, ..., 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$ のとき, 
%
\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},
      ..., 
      v_{k_{max}}^{t+\Delta t} - v_{k_{max}}^{t-\Delta t} \right), 
\\
  \Dvect{G}_v &=& \left( g_{v,1}, g_{v,2}, ..., 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{熱拡散の差分方程式の整理}

熱の鉛直拡散の式 ... を整理すると, 
$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}
%
となり, $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}
%
となる. 

これらをまとめると, 
%
\begin{eqnarray}
    \Dvect{B}_a \Dvect{x}_a = \Dvect{G}_a
    \label{eq:disc:heatdiff}
\end{eqnarray}
%
\begin{eqnarray}
  \Dvect{x}_a &=& \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},
      ..., 
      T_{k_{max}}^{t+\Delta t} - T_{k_{max}}^{t-\Delta t} \right), 
\\
  \Dvect{G}_a &=& \left( g_{a,1}, g_{a,2}, ..., g_{a,k_{max}} \right), 
\\
  g_{a,k} &=& - \left( F_{a,k+\frac{1}{2}}^{t-\Delta t} - F_{a,k-\frac{1}{2}}^{t-\Delta t} \right)
\end{eqnarray}
%
と書くことができる. 
ここで, $2 \le k \le k_{max}-1$ のとき, 
%
\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}}
                     + 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$ のとき, 
%
\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$ 列の行列であり, この式だけでは
未知数が方程式数よりも多いために閉じない. 
方程式を閉じるために, 以下に述べる惑星表面での熱収支式を用いる. 





