% 表題   DCPAM5  乱流過程
%
% 履歴
%\Drireki{2010/04/15 高橋芳幸}
%\Drireki{2011/08/12 高橋芳幸}
%
%  \Dchapterhead
\section{数理表現}


鉛直拡散による運動方程式, 熱力学の式, 成分の式における変化率は
下のようにそれぞれ下のように書くことができる.
%
\begin{eqnarray}
          \DP[][VD]{u}{t}
    &=& - \frac{1}{\rho} \DP{F_{m,x}}{z} \\
    &=&   g \DP{F_{m,x}}{p}, 
\\
          \DP[][VD]{v}{t}
    &=&   g \DP{F_{m,y}}{p}, 
\\
          \DP[][VD]{T}{t}
    &=&   \frac{g}{C_p} \DP{F_{h}}{p}, 
\\
          \DP[][VD]{q}{t}
    &=&   g \DP{F_{q}}{p} .
\end{eqnarray}
%
ここで, $F_{m,x}$, $F_{m,y}$, $F_{h}$, $F_{q}$ はそれぞれ東西方向, 南北方向の
運動量フラックス, 熱フラックス, 水蒸気(物質)のフラックスであり, 下のように表現される.
%
\begin{eqnarray}
  F_{m,x} &=& - \rho K_m \DP{u}{z}, 
\\
  F_{m,y} &=& - \rho K_m \DP{v}{z}, 
\\
  F_{h}   &=& - C_p P \rho K_h \DP{\theta}{z}, 
\\
  F_{q}   &=& - \rho K_q \DP{q}{z}.
\end{eqnarray} 
ただし, 
\begin{eqnarray}
  \theta  &=& \frac{ T }{ P }, 
\\
  P       &=& \left( \frac{ p_{00} }{ p } \right)^\kappa, 
\\
  \kappa  &=& \frac{R}{C_p}
\\
  \rho    &=& \frac{p}{R T_v}
\end{eqnarray}
%
ここで, $T_v$ は仮温度である.
%
上部境界では
%
\begin{eqnarray}
  F_{m,x} &=& 0, 
\\
  F_{m,y} &=& 0, 
\\
  F_{h}   &=& 0, 
\\
  F_{q}   &=& 0
\end{eqnarray}
%
とし, 下部境界では, バルク法を用いてフラックスを評価する場合には, 
%
\begin{eqnarray}
  F_{m,x} &=& - \rho C_d |\Dvect{v}| u, 
\\
  F_{m,y} &=& - \rho C_d |\Dvect{v}| v, 
\\
  F_{h}   &=& - C_p P \rho C_h |\Dvect{v}| 
                \left( \theta - \theta_s  \right), 
\\
  F_{q}   &=& - \epsilon \rho C_q |\Dvect{v}| \left( q - q_s^* \right)
   \Deqlab{バルク法Fq}
\end{eqnarray}
%
とし, 摩擦の時定数 (運動量フラックスに対して) や一定値 (熱や物質フラックス
に対して) を与える場合には, 
%
\begin{eqnarray}
  F_{m,x} &=& - \frac{1}{\tau_f} u, 
\\
  F_{m,y} &=& - \frac{1}{\tau_f} v, 
\\
  F_{h}   &=& F_{h,s}, 
\\
  F_{q}   &=& F_{q,s}
\end{eqnarray}
%
とする.
ただし, \Deqref{バルク法Fq}に示した下部境界における物質のフラックス $F_q$ 
は水蒸気の
フラックスである. 
ここで, $p_{00}$ は基準圧力である. 
$K_m$, $K_h$, $K_q$ はそれぞれ運動量, 熱, 物質の拡散係数である. 
$C_d$, $C_h$, $C_q$ はそれぞれ運動量, 熱, 水蒸気のバルク係数である. 
また, $\tau_f$ は下部境界における摩擦の時定数, $F_{h,s}$, $F_{q,s}$ は
固定する熱フラックス, 水蒸気フラックスである.
$\epsilon$ は地表面の湿潤度である.


$K_m$, $K_h$, $K_q$, はそれぞれ Mellor and Yamada (1982) レベル 2 の方法に従って評価する.
これら拡散係数の具体的な評価方法については \Dsecref{鉛直拡散係数連続系表現}で述べる.
$C_d$, $C_h$, $C_q$, は Monin-Obukhov の相似則に基づくパラメタリゼーションによって評価する.
これらバルク係数の具体的な評価方法については 
\Dsecref{vdiff-math:bulkcoef:L82}, \Dsecref{vdiff-math:bulkcoef:BH91} で述べる.


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

Mellor and Yamada (1982) レベル 2 の方法を用いる場合, 
鉛直拡散係数, $K_m$, $K_h$, $K_q$, および乱流運動エネルギー, 
$\frac{q^2}{2}$, は診断的に計算される. 

Mellor and Yamada (1974, 1982) レベル 2 の方法に従って
鉛直拡散係数, $K_m$, $K_h$, $K_q$, を 
評価する場合, 下のように表現される. 
%
\begin{eqnarray}
  K_m &=& l^2 \left| \DP{\Dvect{v}}{z} \right| S_M' , 
     \Deqlab{Km連続系表現}
\\
  K_h &=& l^2 \left| \DP{\Dvect{v}}{z} \right| S_H' , 
     \Deqlab{Kh連続系表現}
\\
  K_q &=& K_h .
     \Deqlab{Kq連続系表現}  
\end{eqnarray}
%
ここで, $l$ は混合距離であり, 
%
\begin{eqnarray}
  l = \frac{ k \left( z - z_s \right) }{ 1 + k \left( z - z_s \right) / l_0 }
  \Deqlab{混合距離}
\end{eqnarray}
%
の表式を用いる. 
ここで, $z_s$ は地表面高度, $l_0$ は支配混合距離
\footnote{(2011-8-17 石渡) この用語は一般的か?}, 
$k$ はカルマン定数 (K\'{a}rm\'{a}n's constant) である.
%
また, 
%
\begin{eqnarray}
  S_M' = B_1^{\frac{1}{2}} \left( 1 - R_f \right)^{\frac{1}{2}}
   S_M^{\frac{1}{2}} S_M, 
\\
  S_H' = B_1^{\frac{1}{2}} \left( 1 - R_f \right)^{\frac{1}{2}}
   S_M^{\frac{1}{2}} S_H
\end{eqnarray}
%
である.
$S_H$, $S_M$ は, 
%
\begin{eqnarray}
S_H 
&=& \left\{ 
      \begin{array}{cc}
       {\displaystyle \frac{ \alpha_1 - \alpha_2 R_f }{ 1 - R_f } }
           & (R_f <R_{f,critical})\\
       S_{H,min} 
           & (R_f \geq R_{f,critical})\\
      \end{array}
       \right. , 
\\
S_M 
&=& \left\{ 
      \begin{array}{cc}
       {\displaystyle \frac{ \beta_1 - \beta_2 R_f }{ \beta_3 - \beta_4	R_f} 
          S_H }
           & (R_f <R_{f,critical})\\
       S_{M,min} 
           & (R_f \geq R_{f,critical})\\
      \end{array}
       \right. .
\end{eqnarray}
%
である
\footnote{
  (2011-08-26 石渡) オリジナルの Mellor and Yamada (1974) では
  臨界リチャードソン数を使った議論はなされていない.
  この場合分けに関しては, 離散化の部分で記述するべき
  かもしれない.\\
  (2013-08-10 高橋) Mellor and Yamada (1974), p.1801, 左段, に, 
  フラックスリチャードソン数が臨界値を超えた場合に乱流混合が起
  こらないと述べられている. この記述に対応する処理. 
  記述に正確に対応させるためには $S_{M,min}=0$ となる.
}. 
ここで, $R_f$ はフラックスリチャードソン数 
(flux Richardson number)
%
\begin{eqnarray}
  R_f = \frac{1}{2\beta_2} 
        \left\{ \beta_1 + \beta_4 R_i - \sqrt{ \left( \beta_1 + \beta_4 R_i \right)^2 
               - 4 \beta_2 \beta_3 R_i } \right\}
\end{eqnarray}
%
であり, $R_i$ はバルクリチャードソン数 (bulk Richardson number)
%
\begin{eqnarray}
  R_{i} &\equiv& \frac{ {\displaystyle \frac{g}{\theta_v} \DP{\theta_v}{z} } }
                      { {\displaystyle \left| \DP{\Dvect{v}}{z} \right|^2  }  }
  \Deqlab{Ri定義}
\end{eqnarray}
%
である
\footnote{
  ここでは, 水蒸気の効果を考慮して仮温位 $\theta_v$ を用いている. 
}. 
$R_{f,critical}$ は臨界リチャードソン数であり, 
\begin{eqnarray}
 R_{f,critical} = \frac{\gamma_1}{\gamma_1 + \gamma_2}
\end{eqnarray}
である.

なお, 乱流運動エネルギー, $\frac{q^2}{2}$, は, 
%
\begin{eqnarray}
  \frac{q^2}{2} 
    &=& \frac{1}{2} B_1 l^2 ( 1 - R_f ) 
          \left| \DP{\Dvect{v}}{z} \right|^2 S_M
\end{eqnarray}
%
のように表される
\footnote{
  (2013-08-13 高橋)
  Mellor and Yamada level 2.5 の表現との対応を考えた場合, 
  \Deqref{Km連続系表現}, \Deqref{Kh連続系表現} の拡散係数は, 
  この $q$ を使って下のように表現する方が良いと思う.
  %
  \begin{eqnarray}
    K_m &=& l q S_M, 
  \\
    K_h &=& l q S_H, 
  \end{eqnarray}
  %
  このとき, \Deqref{vdiff-math:my2.5:Km}, \Deqref{vdiff-math:my2.5:Kh} と
  同じ形となる. 
  いずれ統一したい.
}. 


また, 
%
\begin{eqnarray}
  \alpha_1 &=& 3 A_2 \gamma_1, 
\\
  \alpha_2 &=& 3 A_2 \left( \gamma_1 + \gamma_2 \right), 
\\
  \beta_1  &=& A_1 B_1 \left( \gamma_1 - C_1 \right), 
\\
  \beta_2  &=& A_1 \left[ B_1 \left( \gamma_1 - C_1 \right) + 6 A_1 + 3
		    A_2 \right], 
\\
  \beta_3  &=& A_2 B_1 \gamma_1, 
\\
  \beta_4  &=& A_2 \left[ B_1 \left( \gamma_1 + \gamma_2 \right) - 3 A_1
		   \right], 
\\
  \gamma_1 &=& \frac{1}{3} - \frac{2 A_1}{B_1}, 
\\
  \gamma_2 &=& \frac{B_2}{B_1} + \frac{6 A_1}{B_1}
\end{eqnarray}
%
であり, 
\begin{eqnarray}
  (A_1, B_1, A_2, B_2, C_1) &=& (0.92, 16.6, 0.74, 10.1, 0.08)
  \Deqlab{vdiff-math:my2.5:constants}
\end{eqnarray}
%
である (Mellor and Yamada, 1982). 


\subsection{乱流運動エネルギー, 鉛直拡散係数 2 (Mellor and Yamada level 2.5)}
\Dseclab{vdiff-math:my2.5}

Mellor and Yamada (1982) レベル 2.5 の方法を用いる場合, 
鉛直拡散係数, $K_m$, $K_h$, $K_q$, は, 予報変数である乱流運動エネルギー, 
$\frac{q^2}{2}$, から計算される. 
乱流運動エネルギーの支配方程式は下のように表される
\footnote{
  $P_s$, $P_b$ の元の形は以下.
  \begin{eqnarray}
    P_s &=& - \overline{u'w'}\DP{u}{z} - \overline{v'w'}\DP{v}{z},
\\
    P_b &=& \frac{g}{\theta} \overline{w'\theta}.
\\
    \epsilon_{TKE} &=& \frac{q^3}{\Lambda_1}
  \end{eqnarray}
}.
%
\begin{eqnarray}
  \DD{}{t}\left( \frac{q^2}{2} \right) 
    &=& - \frac{1}{\rho} \DP{F_{TKE}}{z} 
        + P_s + P_b - \epsilon_{TKE}
\\
    &=& g \DP{F_{TKE}}{p} 
        + P_s + P_b - \epsilon_{TKE}
\\
  F_{TKE} &=& - \rho K_{TKE} \DP{}{z}\left( \frac{q^2}{2} \right)
\\
  P_s &=& K_M \left( \DP{\Dvect{u}}{z} \right)^2
\\
      &=& K_M \left\{ \left( \DP{u}{z} \right)^2 + \left( \DP{v}{z} \right)^2 \right\}
  \Deqlab{vdiff-math:my2.5:Ps}
\\
  P_b &=& - K_H \frac{g}{\theta_v} \DP{\theta_v}{z}
  \Deqlab{vdiff-math:my2.5:Pb}
\\
  \epsilon_{TKE} &=& \frac{q^3}{B_1 l}
\end{eqnarray}
%
ただし, $P_s$, $P_b$ は, \Deqref{vdiff-math:my2.5:Ps using Sm Gm'}, 
\Deqref{vdiff-math:my2.5:Pb using Sm Gm'} において, 別の形で表現される.

なお, $\frac{q^2}{2}$ の境界条件は下のようになる.
%
\begin{eqnarray}
  \frac{q^2}{2} &=& \frac{B_1^\frac{2}{3} u_\tau^2}{2} \hspace{2cm} {\rm at} \hspace{5mm} z = 0
\\
  \frac{q^2}{2} &=& 0 \hspace{2cm} {\rm at} \hspace{5mm} z = \infty
\end{eqnarray}
%
ここで, $u_\tau$ は下部境界における摩擦速度である.

鉛直拡散係数, $K_m$, $K_h$, $K_q$, $K_{TKE}$ は
下のように表現される. 
%
\begin{eqnarray}
  K_m &=& l q S_M, 
     \Deqlab{vdiff-math:my2.5:Km}
\\
  K_h &=& l q S_H, 
     \Deqlab{vdiff-math:my2.5:Kh}
\\
  K_q &=& K_h .
     \Deqlab{vdiff-math:my2.5:Kq}
\\
  K_{TKE} &=& l q S_{TKE}, 
\end{eqnarray}
%
ここで, $l$ は混合距離であり, 
\Deqref{混合距離} により評価する.
ただし, $l$ は (経験的に?) 以下のような制限が必要となる
\footnote{
  この制限は, 安定成層下での乱流渦のサイズの制限を反映している 
  (Galperin et al., 1988).
}.
%
\begin{eqnarray}
  l &\le& \frac{0.53q}{N} = 0.53 \sqrt{\frac{2}{N^2}\left(\frac{q^2}{2}\right)}, 
  \Deqlab{vdiff-math:my2.5:limit of l}
\\
  N^2 &=& \frac{g}{\theta_v}\DP{\theta_v}{z}.
\end{eqnarray}


また, $S_M$, $S_H$ は, Galperin et al. (1988) による修正を考慮し, 
下のように与えられる. 
%
\begin{eqnarray}
  S_M 
    &=& \frac{ A_1 \left( 1 - 3C_1 - \frac{6A_1}{B_1} \right) + 9A_1 ( A_2 + 2A_1 ) G_H S_H }
             { 1 - 9 A_1 A_2 G_H },
\\
  S_H
    &=& \frac{ A_2 \left( 1 - \frac{6A_1}{B_1} \right) }
             { 1 - 3A_2 ( 6A_1 + B_2 ) G_H },
\end{eqnarray}
%
また, $S_{TKE} = 0.2$ である
\footnote{
  Mellor and Yamada (1982), p.862, 左段参照.
}.
%
ここで, $G_H$ は
%
\begin{eqnarray}
  G_H = - \frac{l^2}{q^2} \frac{g}{\theta_v} \DP{\theta_v}{z}
\end{eqnarray}
%
であり, 上記の式には表れないが, 
%
\begin{eqnarray}
  G_M = \frac{l^2}{q^2} \left\{ \left( \DP{u}{z} \right)^2 + \left( \DP{v}{z} \right)^2 \right\}
\end{eqnarray}
%
である
\footnote{
  $S_M$, $S_H$ が $G_H$ のみに依存し, $G_M$ に依存しないのは, 
  Galperin et al. (1988) による修正の結果である. これにより計算が安定する. 
}. 
%
また, これら $S_M$, $S_H$ を用いると, 
\Deqref{vdiff-math:my2.5:Ps}, 
\Deqref{vdiff-math:my2.5:Ps} の $P_{s}$, $P_{b}$ は下のように
書ける.
%
\begin{eqnarray}
  P_s &=& K_M \left( \DP{\Dvect{u}}{z} \right)^2
       = l q S_M \left( \DP{\Dvect{u}}{z} \right)^2
       = 2^\frac{1}{2} l S_M \left( \DP{\Dvect{u}}{z} \right)^2 \left( \frac{q^2}{2} \right)^\frac{1}{2}
  \Deqlab{vdiff-math:my2.5:Ps using Sm Gm'}
\\
  P_b &=& - K_H \left( \frac{g}{\theta_v} \DP{\theta_v}{z} \right)
       =  - l q S_H \left( \frac{g}{\theta_v} \DP{\theta_v}{z} \right)
       = - 2^\frac{1}{2} l S_H 
            \left( \frac{g}{\theta_v} \DP{\theta_v}{z} \right)
            \left( \frac{q^2}{2} \right)^\frac{1}{2}
  \Deqlab{vdiff-math:my2.5:Pb using Sm Gm'}
\end{eqnarray}

また, 計算の安定性を考慮し, 実際には $G_H$ に下のような制限を与える
\footnote{
  $G_H$ の下限は $l$ の制限 (\Deqref{vdiff-math:my2.5:limit of l}) に起因しており, 
  上限はレベル 2 (の条件) において $G_M \ge 0$ となる条件に起因している 
  (Galperin et al., 1988).
}.
%
\begin{eqnarray}
  - 0.53^2 \le G_H &\le& \frac{1}{A_2(12A_1 + B_1 + 3B_2)}.
\end{eqnarray}

定数 $(A_1, B_1, A_2, B_2, C_1)$ は 
\Deqref{vdiff-math:my2.5:constants} で与えられる. 



\subsection{バルク係数 1 (Louis et al., 1982)}
\Dseclab{vdiff-math:bulkcoef:L82}

Louis et al. (1982) の方法によると, バルク係数は下のように評価される.

\subsubsection{中立もしくは安定 ($R_i \ge 0$) な場合}

中立, もしくは安定 $R_i \ge 0$ な場合には, バルク係数は下のように評価する
%
\footnote{
なお, 元論文 (Louis et al., 1982) では下のように表記されている. 

\begin{eqnarray}
  C_d &=& a^2 \frac{1}{ 1 + 2 b R_i \frac{1}{ \sqrt{ 1 + d R_i } } }
\\
  C_h &=& a^2 \frac{1}{ 1 + 3 b R_i \sqrt{ 1 + d R_i } }
\\
  a   &=& \frac{k}{ \log\left( \frac{z+z_0}{z_0} \right) }
\end{eqnarray}
%
ここで, $b = 5$, $d = 5$ である. 
}. 
%
\begin{eqnarray}
  C_d &=& a_m^2 \frac{1}{ 1 + 10 R_i \frac{1}{ \sqrt{ 1 + 5 R_i } } }, 
      \Deqlab{Cd中立・安定連続系表現}
\\
  C_h &=& a_h a_m \frac{1}{ 1 + 15 R_i \sqrt{ 1 + 5 R_i } }, 
      \Deqlab{Ch中立・安定連続系表現}
\\
  a_m &=& \frac{k}{ \log\left( \frac{z + z_{0,m}}{z_{0,m}} \right) }.
\\
  a_h &=& \frac{k}{ \log\left( \frac{z + z_{0,h}}{z_{0,h}} \right) }.
\end{eqnarray}
%
ここで, $z$ は地面からの距離, $k$ はカルマン定数, 
$z_{0,m}$, $z_{0,h}$ はそれぞれ風速, 温度に対する粗度長である
\footnote{
  Louis et al. (1982) では, 風速と温度に対する粗度長には同じ値を用いて
  いるようである. ここでは, それぞれ別の値を用いることを許す形で記述し
  ておく.
}.
なお, $z$ が基準等ポテンシャル面 (地球の場合はジオイド) からの距離では
ないことに注意.


\subsubsection{不安定 ($R_i < 0$) な場合}

不安定 $R_i < 0$ な場合には, バルク係数は下のように評価する
%
\footnote{
なお, 元論文 (Louis et al., 1982) では下のように表記されている. 

\begin{eqnarray}
  C_d &=& a^2 \left( 1 - \frac{ 2 b R_i }{ 1 + 3 a^2 b c \sqrt{
	      \frac{z+z_0}{z_0} |R_i| } } \right), 
\\
  C_h &=& a^2 \left( 1 - \frac{ 3 b R_i }{ 1 + 3 a^2 b c \sqrt{
	      \frac{z+z_0}{z_0} |R_i| } } \right) . 
\end{eqnarray}
%
ここで, $b=5$, $c = 5$ である. 
}$^,$\footnote{
  ここでは, 仮温位を使うことで水蒸気の効果を考慮している.
}. 
%
\begin{eqnarray}
  C_d &=& a_m^2 \left( 1 - \frac{ 10 R_i }{ 1 + 75 a_m^2 \sqrt{ \frac{z + z_{0,m}}{z_{0,m}} |R_i| } } 
              \right), 
      \Deqlab{Cd不安定連続系表現}
\\
  C_h &=& a_h a_m \left( 1 - \frac{ 15 R_i }{ 1 + 75 a_m a_h \sqrt{
	       \frac{z + z_{0,h}}{z_{0,h}} |R_i| } } 
               \right) .
       \Deqlab{Ch不安定連続系表現}
\end{eqnarray}


\subsection{バルク係数 2 (Beljaars and Holtslag, 1991; Beljaars, 1994)}
\Dseclab{vdiff-math:bulkcoef:BH91}

まとまっていない (2013/08/20, yot).

Beljaars and Holtslag (1991) の方法によると, バルク係数は下のように評価される
\footnote{
  ここで書いている $C_h$ は, 粗度長部分が Beljaars (1994) の $C_h$ とは
  異なっている. 
  これは Beljaars (1994) の誤植か, 自分がわかっていないのか? たぶん誤植.
  (yot, 2013/09/08)
}.
%
\begin{eqnarray}
  C_d &=& \left\{
            \frac{k}
                 {   \log\left(\frac{z+z_{0,m}}{z_{0,m}}\right)
                   - \Psi_M\left(\frac{z+z_{0,m}}{L}\right)
                   + \Psi_M\left(\frac{z_{0,m}}{L}\right) }
          \right\}^2 
%  \Deqlab{vdiff-math:Cd:BH91}
\\
  C_h &=& 
            \frac{k}
                 {   \log\left(\frac{z+z_{0,m}}{z_{0,m}}\right)
                   - \Psi_M\left(\frac{z+z_{0,m}}{L}\right)
                   + \Psi_M\left(\frac{z_{0,m}}{L}\right) }
            \frac{k}
                 {   \log\left(\frac{z+z_{0,h}}{z_{0,h}}\right)
                   - \Psi_H\left(\frac{z+z_{0,h}}{L}\right)
                   + \Psi_H\left(\frac{z_{0,h}}{L}\right) }
%  \Deqlab{vdiff-math:Ch:BH91}
\end{eqnarray}
%
ここで, 
ここで, $z$ は地面からの距離, $k$ はカルマン定数, 
$z_{0,m}$, $z_{0,h}$ はそれぞれ風速と温度に対する粗度長である.
$L$ は Monin-Obukhov 長さであり, 
%
\begin{eqnarray}
  L &=& - \frac{1}{k} \frac{ \left| \overline{u'w'} \right|^\frac{3}{2} }{ \displaystyle \frac{g}{\theta_v} \left( \overline{w'\theta_v'} \right) }
\end{eqnarray}
%
である
\footnote{
  Beljaars and Holtslag (1991) では水蒸気の効果を考慮しておらず, 
  $\theta$ で定式化されてきた.
  ここでは, $\theta$ を用いた式と $\theta_v$ を用いた式の関係を考えてみる.
  水蒸気の効果を無視すると, 下のように書くことができる.
%
  \begin{eqnarray}
    \frac{\tau}{\rho} &=& - u_*^2 = - C_m U^2
  \\
    Q_0 &=& - u_* \theta_* = - C_h |U| ( \theta - \theta_s )
  \\
    \frac{u_* \theta_*}{u_*^2} &=& \frac{ C_h |U| (\theta-\theta_s) }{ C_m U^2 }
  \\
    \frac{\theta-\theta_s}{U^2}
    &=& \frac{ u_* \theta_* C_m |U| }{ u_*^2 C_h }
  \\
    &=& \frac{ \theta_* C_m |U| }{ u_* C_h }
  \\
    R_b &=& \frac{g}{\theta} z \frac{ \theta - \theta_s }{ U^2 }
  \\
    L &=& - \frac{1}{k} \frac{ \left( \overline{u'w'} \right)^{1.5} }{ \frac{g}{\theta} \left( \overline{w'\theta'} \right) }
%  \\
%      &=& - \frac{1}{k} \frac{ \left( \frac{\tau}{\rho} \right)^{1.5} }{ \frac{g}{\theta} Q_0 }
%  \\
%      &=& \frac{1}{k} \frac{ \left(C_m U^2\right)^{1.5} }{ \frac{g}{\theta} C_h |U| (\theta-\theta_s) }
%  \\
%      &=& \frac{1}{k} U^2 \frac{\theta}{g} \frac{1}{\theta-\theta_s} \frac{C_m^{1.5}}{C_h}
%  \\
%      &=& \frac{1}{k} \left( \frac{g}{\theta} \frac{\theta-\theta_s}{U^2} \right)^{-1} \frac{C_m^{1.5}}{C_h}
%  \\
%      &=& \frac{1}{k} z \left( \frac{g}{\theta} \frac{\theta-\theta_s}{U^2} z \right)^{-1} \frac{C_m^{1.5}}{C_h}
  \\
      &=& \frac{1}{k} z R_b^{-1} \frac{C_m^{1.5}}{C_h}
  \end{eqnarray}
%
  ここで, $\theta$ を $\theta_v$ で置き換えても, 
%
  \begin{eqnarray}
    L &=& \frac{1}{k} z R_b^{-1} \frac{C_m^{1.5}}{C_h}
  \end{eqnarray}
%
  のように同じ形で書ける. 
  ただし, このとき, 
%
  \begin{eqnarray}
    L &=& - \frac{1}{k} \frac{ \left( \overline{u'w'} \right)^{1.5} }{ \frac{g}{\theta_v} \left( \overline{w'\theta_v'} \right) }
  \\
    Q_0 &=& - u_* \theta_* = - C_h |U| ( \theta_v - \theta_{v,s} )
  \\
    R_b &=& \frac{g}{\theta_v} z \frac{ \theta_v - \theta_{v,s} }{ U^2 }
  \end{eqnarray}
%
  である.

  ただし, ここでは, $\overline{w' \theta'}$ の $\theta$ を単純に 
  $\theta_v$ で置き換えることで, 
  $\overline{w' \theta_v'}$ を計算しているが, 正しくは下のように
  書けるのだろう.
%
  \begin{eqnarray}
    \theta_v &=& \theta ( 1 + \alpha q_v )
  \\
    \theta_v' &=&   \overline{\theta} \alpha q_v' 
                  + \theta' ( 1 + \alpha \overline{q_v} ) + \theta' \alpha q_v'
  \\
    \overline{w' \theta_v'} 
      &=&   \overline{w' \overline{\theta} \alpha q_v'}
          + \overline{w' \theta' (1 + \alpha \overline{q_v})}
          + \overline{w' \theta' \alpha q_v'}
%  \\
%      &=&   \alpha \overline{\theta} \overline{w'q_v'}
%          + \overline{w'\theta'} (1 + \alpha \overline{q_v})
%          + \alpha \overline{w'\theta'q_v'}
%  \\
%      &=&   \overline{w'\theta'} (1 + \alpha \overline{q_v})
%          + \alpha \overline{\theta}\overline{w'q_v'}
%          + \alpha \overline{w'\theta'q_v'}
  \\
      &=&   \overline{w'\theta'}
          + \alpha \overline{q_v}\overline{w'\theta'}
          + \alpha \overline{\theta}\overline{w'q_v'}
          + \alpha \overline{w'\theta'q_v'}
  \\
      &\sim&   \overline{w'\theta'}
             + \alpha \overline{q_v}\overline{w'\theta'}
             + \alpha \overline{\theta}\overline{w'q_v'}
  \end{eqnarray}
  %
  しかし, Beljaars (1994) では, 
  %
  \begin{eqnarray}
    \overline{w' \theta_v'} 
      &=&   \overline{w'\theta'}
          + \alpha \overline{\theta}\overline{w'q_v'}
  \end{eqnarray}
  %
  としている. 右辺第二項は落としているのだろう. 

  これをさらにバルク式を使って変形すると, 
  %
  \begin{eqnarray}
    \overline{w' \theta_v'} 
      &=& - C_h |U| ( \theta - \theta_s ) - \alpha \overline{\theta} \epsilon C_q |U| ( q - q_s )
  \\
      &=& - C_h |U| ( \theta - \theta_s ) - \alpha \theta \epsilon C_q |U| ( q - q_s )
  \end{eqnarray}
  %
  となる. ここで, $\overline{\theta}$ の実態が $\theta$ であることに注意.
  ここで $C_h = C_q$ とし, $\overline{\theta} = \theta$ と書くことにすると, 
  %
  \begin{eqnarray}
    \overline{w' \theta_v'} 
      &=& - C_h |U| \left\{ \theta   \left( 1 + \alpha \epsilon q \right)
                          - \theta_s - \theta \alpha \epsilon q_s \right\}
  \end{eqnarray}
  %
  さらに, $\theta_s = \theta$ とすると, 
  %
  \begin{eqnarray}
    \overline{w' \theta_v'} 
      &=& - C_h |U| \left\{ \theta   \left( 1 + \alpha \epsilon q \right)
                          - \theta_s \left( 1 + \alpha \epsilon q_s \right) \right\}
  \end{eqnarray}
  %
  となり, $\epsilon = 1$ ならば, 
  %
  \begin{eqnarray}
    \overline{w' \theta_v'} 
      &=& - C_h |U| \left\{ \theta   \left( 1 + \alpha q \right)
                          - \theta_s \left( 1 + \alpha q_s \right) \right\}
  \\
      &=& - C_h |U| \left( \theta_v - \theta_{s,v,sat} \right)
  \end{eqnarray}
  %
  ここで, $\theta_{s,v,sat}$ は惑星表面温度での飽和仮温位 (?) である. 

  したがって, $\overline{w' \theta'}$ の $\theta$ を単純に 
  $\theta_v$ で置き換えたものと, 正しい (?) $\overline{w' \theta_v'}$ とは
  いくつかの点で違いがある. 
}.
%
$\Psi_M$, $\Psi_H$ の詳細は下に示す. 
なお, $z$ が基準等ポテンシャル面 (地球の場合はジオイド) からの距離では
ないことに注意.

また, Beljaars (1994) によると, 
%
\begin{eqnarray}
  |\Dvect{v}| &=& \left\{ u^2 + v^2 + \left( \beta w_* \right)^2 \right\}^\frac{1}{2}
\end{eqnarray}
%
として, 特に自由対流レジームにおける対流を想定し, それに伴う風速を考慮する. 
ここで, $\beta$ は定数であり, $w_*$ は下のように与える
\footnote{
  Beljaars (1994) では下のように書かれている.
  %
  \begin{eqnarray}
    w_* &=& \left\{ z_{BL} \frac{g}{T} \overline{w'\theta_v'} \right\}^\frac{1}{3}
  \end{eqnarray}
  %
  ここでは, 地球以外の条件への適応も想定し, $T$ は $\theta$ としておく.
}.
%
\begin{eqnarray}
  w_* &=& \left\{ z_{BL} \frac{g}{\theta} \overline{w'\theta_v'} \right\}^\frac{1}{3}
\end{eqnarray}
%
ここで, $z_{BL}$ は境界層高度である.


\subsubsection{中立もしくは安定 ($R_i \ge 0$) な場合の $\Psi_M$, $\Psi_H$}

中立, もしくは安定 $R_i \ge 0$ な場合には, $\Psi_M$, $\Psi_H$ は
下のように評価する
%
\begin{eqnarray}
  \Psi_M(x) &=& \log\left\{ \frac{ (1+x)^2 (1+x^2) }{ 8 } \right\}
                - 2 \tan^{-1} x + \frac{\pi}{2}
\\
  \Psi_H(x) &=& \log\left\{ \frac{ (1+x^2)^2 }{ 4 } \right\}
\\
  x &=& \left( 1 - 16 \zeta \right)^\frac{1}{4}
\\
  \zeta &=& \frac{z}{L}
\end{eqnarray}

\subsubsection{不安定 ($R_i < 0$) な場合の $\Psi_M$, $\Psi_H$}

不安定 $R_i < 0$ な場合には, $\Psi_M$, $\Psi_H$ は下のように評価する
%
\begin{eqnarray}
  \Psi_M(x) &=& - a \zeta 
                - b \left( \zeta - \frac{c}{d} \right) \exp\left( -d \zeta \right)
                - \frac{bc}{d}
\\
  \Psi_H(x) &=& - \left( 1 + \frac{2}{3}a\zeta \right)^\frac{3}{2}
                - b \left( \zeta - \frac{c}{d} \right) \exp\left( -d \zeta \right)
                - \frac{bd}{d}
                + 1
\end{eqnarray}
%
ここで, $a$, $b$, $c$, $d$ は定数であり, 
$a = 1$, $b = 0.667$, $c = 5$, $d = 0.35$ である.


\subsection{乱流過程で用いられるパラメータの値}

\begin{itemize}
 \item 支配混合距離 $l_0$ : 
       \Dmodel の現在 (2010/02/20) のデフォルト値は $l_0 = 300$ m である.

 \item カルマン定数 $k$ : \Dmodel のデフォルト値は $k = 0.4$ .

       カルマン定数の値は
       1950 年代から 1980 年代までに行われたいくつかの観測
       に基づいて決定されている (近藤, 1982).
%       (\Dtabref{カルマン定数観測値一覧}). 
       Businger et al (1971) 
       はカンザス洲の麦畑における観測で 0.35
       という値を得た. 
       これにより, 0.35 という値が広く使われたのだが, 
       Businger et al (1971) では, 
       風速計付近に設置した観測用電子機器を収めた箱の影響があることや
       三杯風速計の補正が十分ではないことなどがわかった.
       その後, 近藤純正氏のグループが乾燥した水田で観測を行い, 
       0.39 という値を得た (Kondo and Sato, 1982).

       木田・柳瀬 (1999) には, Dean (1978)
       \footnote{
          Dean, R. B., 1978:  
          Reynolds Number Dependence of Skin Friction and Other Bulk Flow
          Variables in Two-Dimensional Rectangular Duct Flow
          {\it J. Fluids Eng.}, {\bf 100}, 215--213.
       }
      が実験的にカルマン定数の値として 0.41 を得た, 
      と記載されている.

%      \begin{table}
%	\begin{center}
%	\begin{tabular}{cc} \hline \hline
%	 研究者 (年) & $k$ \\ \hline \hline
%	 ライダー (1954) &  0.41 \\
%	 ディーコン (1959) &  0.42 \\
%	 ダイヤー・ヒックス (1970) &  0.41 \\
%	 ゴダード (1970) &  0.41 \\
%	 ブシンガーら (1971) &  0.35 \\
%	 プルイトら (1973) &  0.42 \\
%	 近藤・佐藤 (1982) &  0.39 \\ \hline
%	\end{tabular}
%	\end{center}
%	\caption{カルマン定数の観測値一覧. 近藤 (1982) の表2.3 を転載したもの.
%		 }
%	\Dtablab{カルマン定数観測値一覧}
%      \end{table}

 \item 粗度長 $z_0$

       \Dmodel において陸面と海洋の差のみ考慮する場合には, 
       陸面では $z_0=0.1$ m, 海洋上では $z_0=10^{-4}$ m 
       というデフォルト値が設定されている.

       粗度長は地表面の状態によって変化する. 
       Sheppard (1947) は屋外観測によってに各種の表面における粗度長を決定した.
       そこでは, 滑らかな表面上 (氷や滑らかな雪など) において $z_0=9 \times 10^{-6}$ m, 
       成長した根菜類が存在する地面において $z_0=0.139$ m
       という数値が示されている (Shpeppard, 1947 の Table 2).
       Sheppard (1947) が決定した粗度長の値は
       Sutton (1953)
         \footnote{
                    Sutton, O.G., 1953: Micrometeorology.
                    McGrawHill. 333pp.
                  }
       および島貫 (1982) 
         \footnote{
                    島貫 陸, 1982:
                    乱流と気象 -変動とうずの種々相- (気象学のプロムナード 6), 
                    東京堂出版.
                  }
       にも引用されている.
\end{itemize}