% 表題: 基礎方程式系
% 
% 履歴: 2005/11/15 杉山耕一朗
% 履歴: 2006/03/12 杉山耕一朗
% 履歴: 2007/06/22 小高正嗣
% 履歴: 2009/09/07 山下達也
% 履歴: 2009/11/24 山下達也
% 履歴: 2009/12/01 山下達也
% 履歴: 2009/12/08 山下達也
% 履歴: 2010/03/25 山下達也
% 履歴: 2010/04/13 山下達也
% 履歴: 2010/04/21 山下達也
% 履歴: 2010/04/27 山下達也

地球大気における湿潤対流の定式化同様, 大気の乾燥成分と湿潤成分の
分子量の差を密度の式には考慮するが, 熱の式には考慮しないような
系を考える. また気体は理想気体みなせるものとする. 
このような系では温位 $\theta$ が非凝結時の保存量として使える. 


\subsection{
%温位 $\theta$, エクスナー関数 $\Pi$, 風速 $u, w$ を予報変数とする場合
元となる方程式系
}

3 次元大気の状態を
気温 $T$, 圧力 $p$, 風速 $u, v, w$, 密度 $\rho$ で表現する場合, 
一般的な圧縮性流体の方程式系は以下のようになる. 
%
\begin{description}
%
 \item[運動方程式]~
%
\begin{eqnarray}
&& \DD{u}{t} = - \Dinv{\rho}\DP{p}{x} + D_u 
  \Deqlab{equations:u} \\
&& \DD{v}{t} = - \Dinv{\rho}\DP{p}{y} + D_v 
  \Deqlab{equations:v} \\
&& \DD{w}{t} = - \Dinv{\rho}\DP{p}{z} -g + D_w
  \Deqlab{equations:w} 
\end{eqnarray}
%
 \item[熱力学第一法則]~
%
%\begin{eqnarray}
% \langle c_{p} \rangle \DD{T}{t} - \Dinv{\rho_{g}} \DD{p}{t} 
% = Q_{cond} + Q_{rad} + Q_{dis}
%  \Deqlab{equations:T}
%\end{eqnarray}
%
\begin{eqnarray}
 \hat{c_p} \DD{T}{t} - \Dinv{\rho} \DD{p}{t}
= \hat{c_p} \left( Q_{cond} + Q_{rad} + Q_{dis} \right) 
  \Deqlab{theta1}
\end{eqnarray}
%
%
% \item[連続の式]~
%
%\begin{eqnarray}
% \DP{\rho}{t} + \DP{}{x_j} \left( \rho u_j \right) = M_{fall}
%  (\rho_r) 
%  \Deqlab{equations:continue}
%\end{eqnarray}
%
% \item[全密度の式]~
%
%\begin{eqnarray}
%\rho = \rho_d + \sum \rho_v + \sum \rho_c + \sum \rho_r
% \rho = \frac{p}{R_{d} T_{v}}
%  \Deqlab{equations:rho}
%\end{eqnarray}
%
 \item[状態方程式]~
  \begin{eqnarray}
   p_d &=& \rho_d R_d T, 
    \Deqlab{equations:eqsdry} \\
   p_v &=& \rho_v R_v T. 
    \Deqlab{equations:eqsvapor}
  \end{eqnarray}
%
% \begin{eqnarray}
%\rho
%= \frac{p}{\langle R \rangle \Pi \theta_v}
%= \frac{p_0 \Pi^{\langle c_{v} \rangle / \langle R \rangle}}{\langle R
%\rangle \theta_v}
% \Deqlab{equations:eqs}
%\end{eqnarray}
%
% \item[凝結成分の混合比保存式]~
 \item[密度の時間発展方程式]~
%
\begin{eqnarray}
\DP{\rho_d}{t} + \DP{}{x_j} (\rho_d u_j) &=& D_{\rho_d}, \Deqlab{renzoku1}
 \\
\DP{\rho_v}{t} + \DP{}{x_j} (\rho_v u_j) &=& M_{src}(\rho_v) +
 D_{\rho_v}, \Deqlab{renzoku2}
 \\
\DP{\rho_c}{t} + \DP{}{x_j} (\rho_c u_j) &=& M_{src}(\rho_c) +
 D_{\rho_c}, \Deqlab{renzoku3}
 \\
\DP{\rho_r}{t} + \DP{}{x_j} (\rho_r u_j) &=& M_{src}(\rho_r) +
 M_{fall}(\rho_r) +  D_{\rho_r} \Deqlab{renzoku4}
\end{eqnarray}
%
%\begin{eqnarray}
%\DD{q_{d}}{t}
% &=& - \frac{q_{d}}{\rho} \sum M_{fall}(\rho_r)
%     + D_{q_d}, 
% \\
%\DD{q_{v}}{t}
% &=& - \frac{q_{v}}{\rho} \sum M_{fall}(\rho_r)
%     + \Dinv{\rho} M_{src}(\rho_v)
%     + D_{q_v}, \Deqlab{equations:qv}
% \\
%\DD{q_{c}}{t}
% &=& - \frac{q_{c}}{\rho} \sum M_{fall}(\rho_r)
%     + \Dinv{\rho} M_{src}(\rho_c)
%     + D_{q_c}, \Deqlab{equations:qc}
% \\
%\DD{q_r}{t}
% &=& - \frac{q_r}{\rho} \sum M_{fall}(\rho_r)
%     + \Dinv{\rho} M_{fall}(\rho_r) + \Dinv{\rho} M_{src}(\rho_r)
%     + D_{q_r}. \Deqlab{equations:qr}
%\end{eqnarray}
%
%
\end{description}
%
ここで $\hat{c_p}$ は凝結物も含んだ単位質量の気塊の定圧比熱, 
$Q$ は非断熱加熱による温度変化率, $q_{v}$ は凝結性気体の比湿, 
$q_{c}$ は雲水比湿, 
$q_{r}$ は雨水比湿である. 
$q_{v}, q_{r}, q_{c}$ は, 凝結成分の数だけ存在する. 
$D_{\ast}$, $M_{src}$, $M_{fall}$ はそれぞれ
乱流拡散項, 生成消滅項, 落下項を意味する. 
%$T_{v}$ は仮温度であり, ここでは以下のように定義する. 
%%
%\begin{eqnarray}
%&& T_{v} \equiv \frac{T}{f(q_{v},q_{x})}  \\
%&&f(q_{v}, q_{x})  \equiv 
%% f \equiv  
%  \left( 1 - \frac{\sum q_{v}/M_{v}}{1/M_{d} + \sum q_{v}/M_{v}	}\right)
%  (1 + \sum q_{v} + \sum q_{c} + \sum q_{r} )
%  \Deqlab{equations:moist:f}
%\end{eqnarray}
%%
%但し M は乾燥成分と凝結成分の分子量である. 
以下では, 温位 $\theta$, エクスナー関数 $\Pi$, 風速 $u, v, w$ を予報変数と
する場合の基礎方程式系を導出する. 


\subsection{
密度の時間発展方程式の書き換え -- 比湿の時間発展方程式の導出 -- } 


$D_{\rho_d} + \sum D_{\rho_v} + \sum D_{\rho_c} + \sum D_{\rho_r} = 0$ 
となると仮定して 
\Deqref{renzoku1} -- \Deqref{renzoku4} の和をとると, 
\begin{equation}
\DP{\rho}{t} + \DP{}{x_j} (\rho u_j) = \sum M_{fall} (\rho_r)
 \Deqlab{renzoku5}
\end{equation}
が得られる. 
但し $M_{src} (\rho_v) + M_{src} (\rho_c) + M_{src} (\rho_r) = 0$ となる
ことを用いた. 
\Deqref{renzoku1} -- \Deqref{renzoku4} 及び \Deqref{renzoku5} より
\begin{eqnarray}
%
\DP{}{t} \left( \frac{\rho_d}{\rho} \right)
 &=& - \frac{\rho_d}{\rho^2} \DP{\rho}{t} + \Dinv{\rho} \DP{\rho_d}{t}
\nonumber \\
 &=& - \frac{\rho_d}{\rho^2} \left[ - \DP{}{x_j} (\rho u_j) + \sum
  M_{fall} (\rho_r) \right]
  + \Dinv{\rho} 
    \left[ - \DP{}{x_j} (\rho_d u_j) + D_{\rho_d} \right] \nonumber \\
 &=& - u_j \DP{}{x_j} \left( \frac{\rho_d}{\rho} \right)
    - \frac{\rho_d}{\rho^2} \sum M_{fall} (\rho_r)
    + \Dinv{\rho} D_{\rho_d}, 
  \Deqlab{sphm2} \\
%
\DP{}{t} \left( \frac{\rho_v}{\rho} \right)
 &=& 
  - \frac{\rho_v}{\rho^2} \DP{\rho}{t} + \Dinv{\rho}
  \DP{\rho_v}{t} \nonumber \\
 &=& 
  - \frac{\rho_v}{\rho^2} \left[ - \DP{}{x_j} (\rho u_j) + \sum
  M_{fall} (\rho_r) \right] \nonumber \\
 & &
  + \Dinv{\rho} \left[ - \DP{}{x_j} (\rho_v
  u_j) + M_{src} (\rho_v) + D_{\rho_v} \right] \nonumber \\
 &=& - u_j \DP{}{x_j} \left( \frac{\rho_v}{\rho} \right)
    - \frac{\rho_v}{\rho^2} \sum M_{fall} (\rho_r)
    + \Dinv{\rho} M_{src} (\rho_v)
    + \Dinv{\rho} D_{\rho_v}, 
    \Deqlab{sphm3} \\
%
\DP{}{t} \left( \frac{\rho_c}{\rho} \right)
 &=& 
  - \frac{\rho_c}{\rho^2} \DP{\rho}{t} + \Dinv{\rho}
  \DP{\rho_c}{t} 
   \nonumber \\
 &=& 
  - \frac{\rho_c}{\rho^2} \left[ - \DP{}{x_j} (\rho u_j) + \sum
  M_{fall} (\rho_r) \right] 
   \nonumber \\
 & &
  + \Dinv{\rho} \left[ - \DP{}{x_j} (\rho_c
  u_j) + M_{src} (\rho_c) + D_{\rho_c} \right] \nonumber \\
 &=& - u_j \DP{}{x_j} \left( \frac{\rho_c}{\rho} \right)
    - \frac{\rho_c}{\rho^2} \sum M_{fall} (\rho_r)
    + \Dinv{\rho} M_{src} (\rho_c)
    + \Dinv{\rho} D_{\rho_c}, 
   \Deqlab{sphm4} \\
%
\DP{}{t} \left( \frac{\rho_r}{\rho} \right)
 &=& - \frac{\rho_r}{\rho^2} \DP{\rho}{t} + \Dinv{\rho}
 \DP{\rho_r}{t} \nonumber \\
 &=& - \frac{\rho_r}{\rho^2} \left[ - \DP{}{x_j} (\rho u_j) + \sum
  M_{fall} (\rho_r) \right] \nonumber \\
 & & + \Dinv{\rho} \left[ - \DP{}{x_j} (\rho_r u_j) + M_{src} (\rho_r)
  + M_{fall} (\rho_r) + D_{\rho_r} \right] \nonumber \\
 &=& - u_j \DP{}{x_j} \left( \frac{\rho_r}{\rho} \right)
    - \frac{\rho_r}{\rho^2} \sum M_{fall} + \Dinv{\rho}
    M_{src} (\rho_r) + \Dinv{\rho} M_{fall} (\rho_r)
    + \Dinv{\rho} D_{\rho_r}. \nonumber \\
   \Deqlab{sphm5}
\end{eqnarray}
%ここで
%\begin{eqnarray}
% \Dinv{\rho} D_{\rho_d} &=& D_{q_d}, \\ 
% \Dinv{\rho} D_{\rho_v} &=& D_{q_v}, \\ 
% \Dinv{\rho} D_{\rho_c} &=& D_{q_c}, \\ 
% \Dinv{\rho} D_{\rho_r} &=& D_{q_r}
%\end{eqnarray}
%と置くと, 
\Deqref{sphm2} -- \Deqref{sphm5} を比湿で表現すると, 
\begin{eqnarray}
\DD{q_{d}}{t}
 &=& - \frac{q_{d}}{\rho} \sum M_{fall}(\rho_r)
     + \Dinv{\rho} D_{\rho_d}, 
     \Deqlab{sphm6} \\
\DD{q_{v}}{t}
 &=& - \frac{q_{v}}{\rho} \sum M_{fall}(\rho_r)
     + \Dinv{\rho} M_{src}(\rho_v)
     + \Dinv{\rho} D_{\rho_v}, 
     \Deqlab{sphm7} \\
\DD{q_{c}}{t}
 &=& - \frac{q_{c}}{\rho} \sum M_{fall}(\rho_r)
     + \Dinv{\rho} M_{src}(\rho_c)
     + \Dinv{\rho} D_{\rho_c}, 
     \Deqlab{sphm8} \\
\DD{q_r}{t}
 &=& - \frac{q_r}{\rho} \sum M_{fall}(\rho_r)
     + \Dinv{\rho} M_{fall}(\rho_r) + \Dinv{\rho} M_{src}(\rho_r)
     + \Dinv{\rho} D_{\rho_r}.
     \Deqlab{sphm9}
\end{eqnarray}
となる.
但し, $q_d + \sum q_v + \sum q_c + \sum q_r = 1$ の関係が成り立つので,
$q_d$ については時間発展方程式を解かずに, 診断的に求めることとする. 


\subsection{ 熱の式の導出 }

A.1.3 節では, Satoh(2004) に従い, Clausius-Clapeyron の式, 各カテゴリー
の比エンタルピーと比内部エネルギーの導出, 比内部エネルギーの時間発展方程
式の導出を行なう. 
更にそれらの結果を用いて, 熱の式の導出を行なう. 
A.1.3 節においては, 雨粒または氷粒の存在を無視する. 
凝結性気体と凝結物は相平衡状態にあるとする. 
また各気体と凝結物の定圧比熱の温度・圧力依存性を無視する. 
凝結性気体は非凝結性気体に比べて圧倒的に多いか, 或いは圧倒的に少ないもの
とする. 

\subsubsection{ Clausius-Clapeyron の式の導出 }

各カテゴリーの比エンタルピー, 比内部エネルギー, 比エントロピー, 比容, 化
学ポテンシャルをそれぞれ $H_i$, $U_i$, $S_i$, $v_i$, $\mu_i$とする. 
また飽和蒸気圧を $p^{*}$ と表すことにする. 
相平衡状態が実現されているとき, 
\begin{eqnarray}
\mu_c (p^{*}, T) &=& \mu_v (p^{*}, T), 
 \Deqlab{Cl-Cl1} \\
\mu_c (p^{*} + d p^{*}, T + d T) 
 &=& 
  \mu_v (p^{*} + d p^{*}, T + d T), 
 \Deqlab{Cl-Cl2}
\end{eqnarray}
が成り立つ. 
\Deqref{Cl-Cl1}, \Deqref{Cl-Cl2} の差をとると, 
\begin{eqnarray}
\mu_c (p^{*} + d p^{*}, T + d T) - \mu_c (p^{*}, T)
 &=&
\mu_v (p^{*} + d p^{*}, T + d T) - \mu_v (p^{*}, T)
 \Deqlab{Cl-Cl3}
  \nonumber \\
\end{eqnarray}
となる. 
$\mu_v (p^{*} + d p^{*}, T + d T)$ 及び $\mu_c (p^{*} + d p^{*}, T + d
T)$ に関する $(p^{*}, T)$ 近傍での 1 次の Taylor 展開を考えると, 
\begin{eqnarray}
\mu_c (p^{*} + d p^{*}, T + d T)
 &\simeq& 
  \mu_c (p^{*}, T)
  + \left( \DP{\mu_c}{p^{*}} \right)_T d p^{*}
  + \left( \DP{\mu_c}{T} \right)_{p^{*}} d T, 
   \Deqlab{Cl-Cl4}
    \nonumber \\
\mu_c (p^{*} + d p^{*}, T + d T)
 &\simeq& 
  \mu_v (p^{*}, T)
  + \left( \DP{\mu_v}{p^{*}} \right)_T d p^{*}
  + \left( \DP{\mu_v}{T} \right)_{p^{*}} d T 
   \Deqlab{Cl-Cl5}
\end{eqnarray}
となる. 
\Deqref{Cl-Cl4}, \Deqref{Cl-Cl5} を \Deqref{Cl-Cl3} に代入すると, 
\begin{eqnarray}
 \left( \DP{\mu_c}{p^{*}} \right)_T d p^{*}
 + \left( \DP{\mu_c}{T} \right)_{p^{*}} d T
 =
  \left( \DP{\mu_v}{p^{*}} \right)_T d p^{*}
  + \left( \DP{\mu_v}{T} \right)_{p^{*}} d T 
   \Deqlab{Cl-Cl6}
\end{eqnarray}
となる. 
ここで
\begin{eqnarray}
\left( \DP{\mu_i}{T} \right)_p &=& - S_i, 
   \Deqlab{Cl-Cl7} \\
\left( \DP{\mu_i}{p} \right)_T &=& v_i, 
   \Deqlab{Cl-Cl8}
\end{eqnarray}
となることに注意すると, 
\begin{eqnarray}
 v_c d p^{*}
 - s_c d T
 =
 v_v d p^{*}
 - s_v d T, 
  \nonumber
\end{eqnarray}
すなわち
\begin{eqnarray}
\DD{p^{*}}{T} &=& \frac{s_v - s_c}{v_v - v_c}
   \Deqlab{Cl-Cl9}
\end{eqnarray}
となる. 
潜熱 $L$ を
\begin{eqnarray}
L &\equiv& T (s_v - s_c)
   \Deqlab{Cl-Cl10}
\end{eqnarray}
と定義すると, 
\begin{eqnarray}
\DD{p^{*}}{T} &=& \frac{L}{T(v_v - v_c)}
   \Deqlab{Cl-Cl11}
\end{eqnarray}
となる. 
ここで $v_v \gg v_c$ であることを用い, さらに理想気体の状態方程式を用い
ると, 
\begin{eqnarray}
\DD{p^{*}}{T} 
 &\simeq& \frac{L}{T v_v}
   \nonumber \\
 &=& \frac{L p^{*}}{R_v T^2}
   \Deqlab{Cl-Cl12}
\end{eqnarray}
となり, Clausius-Clapeyron の式が得られる. 

\subsubsection{ 比エンタルピーの導出 }

定圧比熱の定義
\begin{eqnarray}
c_{pi} = T \left( \DP{S_i}{T} \right)_p
 \Deqlab{thermodyn1}
\end{eqnarray}
及び Maxwell の関係式
\begin{eqnarray}
- \left( \DP{S_i}{p} \right)_T
 = \left( \DP{v_i}{T} \right)_p
 \Deqlab{thermodyn2}
\end{eqnarray}
を用いると, 
\begin{eqnarray}
c_{pi} (p,T)
 &=& 
  c_{pi} (p_0,T) 
   + \int_{p_0}^{p} \left( \DP{c_{pi}}{p^{\prime}} \right)_T d
   p^{\prime}
   \nonumber \\
 &=&
  c_{pi} (p_0,T) 
   + \int_{p_0}^{p} 
    \left\{ \DP{}{p^{\prime}}
     \left[
      T \left( \DP{S_i}{T} \right)_{p^{\prime}} 
     \right]_T
    \right\} d p^{\prime}
   \nonumber \\
 &=&
  c_{pi} (p_0,T) 
   + T \int_{p_0}^{p} 
     \left[ \DP{}{T}
       \left( \DP{S_i}{p^{\prime}} \right)_{T} 
     \right]_{p^{\prime}}
    d p^{\prime}
   \nonumber \\
 &=&
  c_{pi} (p_0,T) 
   - T \int_{p_0}^{p} 
     \left( \DP[2]{v_i}{T} \right)_{p^{\prime}}
    d p^{\prime}
 \Deqlab{thermodyn3}
\end{eqnarray}
となる. 
但し $p_0$ は圧力の基準値である. 
特に理想気体の場合, 状態方程式より $( \partial^2 v_i / \partial
T^2 )_p = 0$ となるので, 
\begin{eqnarray}
c_{pi} (p,T)
 &=&
  c_{pi} (p_0,T) 
 \Deqlab{thermodyn4}
\end{eqnarray}
となる. 
\Deqref{thermodyn4} より理想気体の比熱は圧力に依存しないことが分かる. 

エンタルピーに関して
\begin{eqnarray}
d H_i 
 &=& T d S_i + v_i d p + \mu_i d q_i
  \nonumber \\
 &=& 
  T \left[ 
     \left( \DP{S_i}{T} \right)_p d T  
     + \left( \DP{S_i}{p} \right)_T d p  
    \right]
   + v_i d p + \mu_i d q_i
  \nonumber \\
 &=& 
  T  \left( \DP{S_i}{T} \right)_p d T  
   + \left[
      T \left( \DP{S_i}{p} \right)_T d p + v_i
     \right] d p
   + \mu_i d q_i
 \Deqlab{thermodyn5}
\end{eqnarray}
より
\begin{eqnarray}
\left( \DP{H_i}{T} \right)_p
 &=& T  \left( \DP{S_i}{T} \right)_p,  
  \Deqlab{thermodyn6} \\
\left( \DP{H_i}{p} \right)_T
 &=& T \left( \DP{S_i}{p} \right)_T d p + v_i
  \Deqlab{thermodyn7}
\end{eqnarray}
が成り立つ. 
$T_0$ を温度の基準値とすると, \Deqref{thermodyn1}, \Deqref{thermodyn2},
\Deqref{thermodyn6}, \Deqref{thermodyn7} より
\begin{eqnarray}
H_i (p,T)
 &=& 
  H_i (p_0, T_0) 
   + \int_{T_0}^{T} \left( \DP{H_i}{T^{\prime}} \right)_{p_0} 
   d T^{\prime}
   + \int_{p_0}^{p} \left( \DP{H_i}{p^{\prime}} \right)_{T} 
   d p^{\prime}
  \nonumber \\
 &=&
  H_i (p_0, T_0) 
   + \int_{T_0}^{T} T^{\prime} \left( \DP{S_i}{T^{\prime}} \right)_{p_0} 
   d T^{\prime}
   + \int_{p_0}^{p} 
    \left[
     T \left( \DP{S_i}{p^{\prime}} \right)_{T} + v_i 
    \right]
   d p^{\prime}
  \nonumber \\ 
 &=&
  H_i (p_0, T_0) 
   + \int_{T_0}^{T} c_{pi} (p_0, T^{\prime})
   d T^{\prime}
   + \int_{p_0}^{p} 
    \left[
     - T \left( \DP{v_i}{T} \right)_{p^{\prime}} + v_i 
    \right]
   d p^{\prime}
  \nonumber \\  
 &=&
  H_i (p_0, T_0) 
   + \int_{T_0}^{T} c_{pi} (p_0, T^{\prime})
     d T^{\prime}
   - T^2 \int_{p_0}^{p} 
    \left[
     \Dinv{T} \left( \DP{v_i}{T} \right)_{p^{\prime}}
     - \frac{v_i}{T^2} 
    \right]
     d p^{\prime}
   \nonumber \\
 &=&
  H_i (p_0, T_0) 
   + \int_{T_0}^{T} c_{pi} (p_0, T^{\prime})
     d T^{\prime}
   - T^2 \int_{p_0}^{p} 
    \left\{ \DP{}{T}
      \left[ \frac{v_i (p^{\prime}, T)}{T} \right]_{p^{\prime}}
    \right\}
   d p^{\prime}
   \nonumber \\
 \Deqlab{thermodyn8}
\end{eqnarray}
となる. 
特に理想気体の場合, 状態方程式より $\partial (v_i / T)_{p} / \partial T
= 0 $ となるので, 
\begin{eqnarray}
H_i (p,T)
 &=& 
  H_i (p_0, T_0) 
   + \int_{T_0}^{T} c_{pi} (p_0, T^{\prime})
     d T^{\prime}
 \Deqlab{thermodyn9}
\end{eqnarray}
となる. 

以下, 相平衡時における $H_c$ と $H_v$ の関係式を求める. 
$H_c$ はそれぞれ
\begin{eqnarray}
H_c (p, T) &=& \mu_c (p, T) + T S_c (p, T), 
 \Deqlab{thermodyn10} \\
H_v (p, T) &=& \mu_v (p, T) + T S_v (p, T)
 \Deqlab{thermodyn11}
\end{eqnarray}
と表される. 
\Deqref{thermodyn10}, \Deqref{thermodyn11} の差をとると, 
\begin{eqnarray}
H_c (p, T) - H_v (p, T)
 &=& 
  \mu_c (p, T) - \mu_v (p, T) 
  + T \left[ S_c (p, T) - S_v (p, T) \right]
 \Deqlab{thermodyn12}  
\end{eqnarray}
となる. 
\Deqref{thermodyn12} に
\begin{eqnarray}
S_i (p, T) &=& - \left( \DP{\mu_i}{T} \right)_p
 \Deqlab{thermodyn13}
\end{eqnarray}
を用いると, 
\begin{eqnarray}
H_c (p, T) - H_v (p, T)
 &=& 
  \mu_c (p, T) - \mu_v (p, T) 
  - T \DP{}{T} \left[ \mu_c (p, T) - \mu_v (p, T) \right]_p
  \nonumber \\
 \Deqlab{thermodyn14}  
\end{eqnarray}
となる. 
\Deqref{Cl-Cl7} において $i=v, c$ とし, $p^{*}$ から $p$ まで積分すると, 
\begin{eqnarray}
\mu_c (p, T) 
 &=& 
  \mu_c (p^{*}, T) 
  + \int_{p^{*}}^{p} v_c (p^{\prime}, T) d p^{\prime}, 
 \Deqlab{thermodyn15} \\
\mu_v (p, T) 
 &=& 
  \mu_v (p^{*}, T) 
  + \int_{p^{*}}^{p} v_v (p^{\prime}, T) d p^{\prime}
 \Deqlab{thermodyn16}
\end{eqnarray}
となる. 
\Deqref{thermodyn15}, \Deqref{thermodyn16} の差をとり, $\mu_c (p^{*}, T)
= \mu_v (p^{*}, T) $ となることを用いると, 
\begin{eqnarray}
\mu_c (p, T) 
 - \mu_v (p, T) 
 &=& 
   - \int_{p^{*}}^{p} 
    \left[ v_v (p^{\prime}, T) - v_c (p^{\prime}, T) \right]
    d p^{\prime}
 \Deqlab{thermodyn17}
\end{eqnarray}
となる. 
\Deqref{thermodyn17} を \Deqref{thermodyn14} に代入し, $\mu_v$, $\mu_c$
を消去すると, 
\begin{eqnarray}
H_c (p, T) - H_v (p, T)
 &=& 
   - \int_{p^{*}}^{p} 
    \left[ v_v (p^{\prime}, T) - v_c (p^{\prime}, T) \right]
    d p^{\prime}
  \nonumber \\
 & & 
  + T \DP{}{T} 
    \left[  
     \int_{p^{*}}^{p} 
      \left[ v_v (p^{\prime}, T) - v_c (p^{\prime}, T) \right]
      d p^{\prime}
    \right]_p
  \nonumber \\
 &\simeq&
   - \int_{p^{*}}^{p} 
      v_v (p^{\prime}, T) d p^{\prime}
  + T \DP{}{T} 
    \left[  
     \int_{p^{*}}^{p} 
       v_v (p^{\prime}, T) d p^{\prime}
    \right]_p
  \nonumber \\
 &=&
   - \int_{p^{*}}^{p} 
      v_v (p^{\prime}, T) d p^{\prime}
  + T \int_{p^{*}}^{p} 
      \DP{ v_v (p^{\prime}, T)}{T} d p^{\prime}
  - T \DD{p^{*}}{T} v_v (p^{*}, T)
  \nonumber \\
 &=&
   - \int_{p^{*}}^{p} 
      v_v (p^{\prime}, T) d p^{\prime}
   +  \int_{p^{*}}^{p} 
      v_v (p^{\prime}, T) d p^{\prime}
  - L(T)
  \nonumber \\
 &=&
  - L(T)
 \Deqlab{thermodyn18}  
\end{eqnarray}
が得られる. 

さらに相平衡時における $c_{pc}$ と $c_{pv}$ の関係式を求める. 
\Deqref{thermodyn1} に \Deqref{thermodyn9}, \Deqref{thermodyn10},
\Deqref{thermodyn13}, \Deqref{thermodyn18} を適用すると, 
\begin{eqnarray}
c_{pc} (p, T)
 &=& 
  T \left( \DP{S_c}{T} \right)_p
   \nonumber \\
 &=&
  T 
   \left\{ 
    \DP{}{T} \left[ \frac{H_c (p, T) - \mu_c (p, T)}{T} \right] 
   \right\}_p
   \nonumber \\
 &=&
  T 
   \left[
    - \frac{S_c}{T} 
    + \Dinv{T} \left( \DP{H_c}{T} \right)_p
    + \Dinv{T} \left( \DP{\mu_c}{T} \right)_p
   \right]
   \nonumber \\
 &=&
  \left( \DP{H_c}{T} \right)_p
   \nonumber \\
 &=&
  \DP{}{T} 
   \left[ H_v (p, T) - L(T) \right]_p
   \nonumber \\
 &=&
  \DP{}{T}
   \left[ 
    H_v (p_0, T_0) 
    + \int_{T_0}^{T} c_{pv} (p_0, T^{\prime}) d T^{\prime}
    - L(T)  
   \right]
   \nonumber \\
 &=&
  c_{pv} (p_0, T) - \DD{L}{T}
 \Deqlab{thermodyn19}  
\end{eqnarray}
が得られる. 
特に定圧比熱が温度に依存しない場合, \Deqref{thermodyn19} を $T$ について
積分することにより, 
\begin{eqnarray}
L(T) 
 &=&
  L(T_0) + \left( c_{pv} - c_{pc} \right) (T - T_0)
   \nonumber \\
 &=&
  L_{00} + \left( c_{pv} - c_{pc} \right) T
 \Deqlab{thermodyn20}  
\end{eqnarray}
が得られる. 
但し
\begin{eqnarray}
L_{00} 
 &=& L (T_0) - \left( c_{pv} - c_{pc} \right) T_0
 \Deqlab{thermodyn21}  
\end{eqnarray}
である. 

\Deqref{thermodyn9}, \Deqref{thermodyn18}, \Deqref{thermodyn19},
\Deqref{thermodyn20} を用いて定圧比熱が温度に依存しない場合の比エンタル
ピーを求めると, 
\begin{eqnarray}
H_d (p, T)
 &=&
  H_d (p_0, T_0)
  + \int_{T_0}^{T} c_{pd} (p_0, T_0) d T^{\prime}
   \nonumber \\
 &=&
  H_d (p_0, T_0) + c_{pd} (T - T_0)
   \nonumber \\
 &=& c_{pd} T, 
  \Deqlab{thermodyn22} \\
H_c (p, T)
 &=&
  H_v (p, T) - L(T)
   \nonumber \\
 &=&
  H_v (p_0, T_0)
  + \int_{T_0}^{T} c_{pv} (p_0, T^{\prime}) d T^{\prime}
  - L(T)
   \nonumber \\
 &=&
  H_v (p_0, T_0)
  + \int_{T_0}^{T} 
    \left[ 
     c_{pc} (p_0, T^{\prime}) + \DD{L}{T^{\prime}} 
    \right] d T^{\prime}
  - L(T)
   \nonumber \\
 &=&
  H_v (p_0, T_0) - L(T_0) + c_{pc} (T - T_0)
   \nonumber \\
 &=&
  H_c (p_0, T_0) - c_{pc} T_0 + c_{pc} T
   \nonumber \\
 &=& c_{pc} T
  \Deqlab{thermodyn23} \\
H_v (p, T)
 &=&
  H_c (p, T) + L(T)
   \nonumber \\
 &=&
  c_{pc} T + L_{00} + (c_{pv} - c_{pc}) T
   \nonumber \\
 &=&
  c_{pv} T + L_{00}
  \Deqlab{thermodyn24}
\end{eqnarray}
となる. 

\subsubsection{ 比内部エネルギーの導出 }

エンタルピーの定義より
\begin{eqnarray}
U_i (p, T)
 &=& H_i (p, T) - p v_i (p, T)
  \Deqlab{thermodyn25}
\end{eqnarray}
となる. 
従って, 定圧比熱が温度に依存しない場合の比内部エネルギーを求めると, 
\Deqref{thermodyn22}, \Deqref{thermodyn23}, \Deqref{thermodyn24} より
\begin{eqnarray}
U_d (p, T)
 &=& H_d (p, T) - p v_d (p, T)
   \nonumber \\
 &=& c_{pd} T - R_d T
   \nonumber \\
 &=& c_{vd} T, 
  \Deqlab{thermodyn26} \\
U_v (p, T)
 &=& H_v (p, T) - p v_v (p, T)
   \nonumber \\
 &=& c_{pv} T + L_{00}- R_v T
   \nonumber \\
 &=& c_{vv} T + L_{00}, 
  \Deqlab{thermodyn27} \\
U_c (p, T)
 &=& U_v (p, T) + H_c (p, T) - H_v (p, T) - p (v_c - v_v)
  \nonumber \\
 &\simeq& U_v (p, T) + H_c (p, T) - H_v (p, T) + p v_v
  \nonumber \\
 &=& H_c (p, T)
  \nonumber \\
 &=& c_{pc} T
  \Deqlab{thermodyn28}
\end{eqnarray}
となる. 

\subsubsection{ 熱の式の導出 }

Satoh(2004)と Landau and Lifshitz(1987) を参考にして熱の式
\Deqref{theta1} を導出する. 

比内部エネルギー $U$ の時間発展方程式は
\begin{equation}
 \rho \DD{U}{t}
  + \nabla \cdot \Dvect{F}^{rad}
  + \nabla \cdot \Dvect{F}^{therm}
 =
  - p \nabla \cdot \Dvect{u}
  + \varepsilon
 \Deqlab{netu1}
\end{equation}
と表される. 
ここで $F_i^{rad}$, $F_i^{therm}$, $\varepsilon$ はそれぞれ放射フラック
ス, 熱拡散フラックス, 消散率であり, 
\begin{eqnarray}
 \Dvect{F}^{therm}
  &=& 
   H_d \Dvect{i}_d + \sum H_v \Dvect{i}_v - \kappa \nabla T
    \nonumber \\
  &=&
   c_{pd} T \Dvect{i}_d 
   + \sum \left(c_{pv} T + L_{00} \right) \Dvect{i}_v 
   - \kappa \nabla T, 
 \Deqlab{netu1-A} \\
 \varepsilon
  &=&
   \Dinv{2} \eta 
    \left( 
     \DP{u_i}{x_j} + \DP{u_j}{x_i} 
      - \frac{2}{3} \delta_{ij} \nabla \cdot \Dvect{u} \right)^2
   + \zeta \left( \nabla \cdot \Dvect{u} \right)^2
 \Deqlab{netu1-B}
\end{eqnarray}
である. 
$H_d$, $H_v$ はそれぞれ非凝結性気体と凝結性気体の比エンタルピー, 
$\Dvect{i}_d$, $\Dvect{i}_v$ はそれぞれ非凝結性気体と凝結性気体の拡散流
束密度, 
$\kappa$ は熱伝導率, $\eta$, $\zeta$ は第一粘性係数, 第二粘性係数である. 
$L_{00}$ は潜熱 $L(T)$ の定数部分であり,
\begin{equation}
L_{00} = L(T) - (c_{pv} - c_{pc}) T
  \Deqlab{netu5}
\end{equation}
である.
$U$ は比エンタルピー $H$ を用いて
\begin{equation}
U = H - \frac{p}{\rho}
  \Deqlab{netu2}
\end{equation}
と表される.
\Deqref{netu2} を \Deqref{netu1} に代入し, $dH/dt$ について整理すると, 
\begin{eqnarray}
 & &
 \rho \DD{H}{t}
  + \frac{p}{\rho} \DD{\rho}{t}
  - \DD{p}{t}
  + \nabla \cdot \Dvect{F}^{rad}
  + \nabla \cdot \Dvect{F}^{therm}
  + p \DP{u_i}{x_i}
  - \varepsilon
 \nonumber \\
 &=&
 \rho \DD{H}{t}
  - p \DP{u_i}{x_i}
  - \DD{p}{t}
  + \nabla \cdot \Dvect{F}^{rad}
  + \nabla \cdot \Dvect{F}^{therm}
  + p \DP{u_i}{x_i}
  - \varepsilon
 \nonumber \\
 &=&
 \rho \DD{H}{t}
  - \DD{p}{t}
  + \nabla \cdot \Dvect{F}^{rad}
  + \nabla \cdot \Dvect{F}^{therm}
  - \varepsilon
 \nonumber \\
 &=&
 0, 
\end{eqnarray}
すなわち
\begin{eqnarray}
\DD{H}{t}
 = 
  \Dinv{\rho} 
   \left( 
    \DD{p}{t}
    - \nabla \cdot \Dvect{F}^{rad}
    - \nabla \cdot \Dvect{F}^{therm}
    + \varepsilon
   \right)
  \Deqlab{netu3}
\end{eqnarray}
となる. 
比エンタルピーを各カテゴリーの比エンタルピーの和で表現すると
\begin{eqnarray}
H
 &=&
  q_d H_d + \sum q_v H_v + \sum q_c H_c
  \nonumber \\
 &=&
  q_d c_{pd} T + \sum q_v \left( c_{pv} T + L_{00} \right)
  + \sum q_c c_{pc} T
  \Deqlab{netu4}
\end{eqnarray}
となる.
\Deqref{netu4} のラグランジュ微分をとると,
\begin{eqnarray}
\DD{H}{t}
 &=&
 \left(q_d c_{pd} + \sum q_v c_{pv} + \sum q_c c_{pc} \right) \DD{T}{t}
  \nonumber \\
 & &
    + c_{pd} T \DD{q_d}{t} 
    + \sum \left( c_{pv} T + L_{00} \right) \DD{q_v}{t} 
    + \sum c_{pc} T \DD{q_c}{t}
  \nonumber \\
 &=&
  \hat{c_p} \DD{T}{t}
    + c_{pd} T \DD{q_d}{t} 
    + \sum \left( c_{pv} T + L_{00} \right) \DD{q_v}{t} 
    + \sum c_{pc} T \DD{q_c}{t}
  \Deqlab{netu6}
\end{eqnarray}
となる.
但し
\begin{equation}
\hat{c_p}
 = q_d c_{pd} + \sum q_v c_{pv} + \sum q_c c_{pc}
  \Deqlab{netu7}
\end{equation}
と置いた.
\Deqref{netu3} を \Deqref{netu6} に代入し, $\hat{c_p} dT/dt$ について整
理すると, 
\begin{eqnarray}
  \hat{c_p} \DD{T}{t}
 &=&
  \Dinv{\rho}
  \left(
  \DD{p}{t}
  - \nabla \cdot \Dvect{F}^{rad}
  - \nabla \cdot \Dvect{F}^{therm}
  + \varepsilon
  \right)
  \nonumber \\
 & &
    - c_{pd} T \DD{q_d}{t} 
    - \sum \left( c_{pv} T + L_{00} \right) \DD{q_v}{t} 
    - \sum c_{pc} T \DD{q_c}{t}
  \Deqlab{netu8}
\end{eqnarray}
雨粒又は氷粒の存在を無視する場合, 比湿の時間発展方程式は
\begin{eqnarray}
\DD{q_{d}}{t}
 &=& \Dinv{\rho} D_{\rho_d}, 
     \Deqlab{sphm6-A} \\
\DD{q_{v}}{t}
 &=& \Dinv{\rho} M_{src}(\rho_v)
     + \Dinv{\rho} D_{\rho_v}, 
     \Deqlab{sphm7-A} \\
\DD{q_{c}}{t}
 &=& - \Dinv{\rho} M_{src}(\rho_v)
     + \Dinv{\rho} D_{\rho_c}, 
     \Deqlab{sphm8-A}
\end{eqnarray}
と表される. 
$D_{\rho_d}$, $D_{\rho_v}$, $D_{\rho_c}$ はそれぞれ
\begin{eqnarray}
D_{\rho_d} 
  &=& - \nabla \cdot \Dvect{i}_d, 
  \Deqlab{DiffFlux1} \\
D_{\rho_v} 
  &=& - \nabla \cdot \Dvect{i}_v, 
  \Deqlab{DiffFlux2} \\
D_{\rho_c} 
  &=& 0
  \Deqlab{DiffFlux3}
\end{eqnarray}
と表されるので, \Deqref{sphm6-A}, \Deqref{sphm7-A}, \Deqref{sphm8-A} は
\begin{eqnarray}
\DD{q_{d}}{t}
 &=& - \Dinv{\rho} \nabla \cdot \Dvect{i}_d, 
     \Deqlab{sphm6-B} \\
\DD{q_{v}}{t}
 &=& \Dinv{\rho} M_{src}(\rho_v)
     - \Dinv{\rho} \nabla \cdot \Dvect{i}_v,  
     \Deqlab{sphm7-B} \\
\DD{q_{c}}{t}
 &=& - \Dinv{\rho} M_{src}(\rho_v)
     \Deqlab{sphm8-B}
\end{eqnarray}
となる. 
 \Deqref{sphm6-B}, \Deqref{sphm7-B},
\Deqref{sphm8-B} を \Deqref{netu8} に代入すると,
\begin{eqnarray}
  \hat{c_p} \DD{T}{t}
 &=&
  \Dinv{\rho}
  \left(
  \DD{p}{t}
  - \nabla \cdot \Dvect{F}^{rad}
  - \nabla \cdot \Dvect{F}^{therm}
  + \varepsilon
  \right)
  \nonumber \\
 & &
    + \Dinv{\rho} c_{pd} T \nabla \cdot \Dvect{i}_d 
    - \Dinv{\rho} 
      \sum \left( c_{pv} T + L_{00} \right) 
      \left[ M_{src} (\rho_v) - \nabla \cdot \Dvect{i}_v \right]
  \nonumber \\
 & &
    + \Dinv{\rho} \sum c_{pc} T M_{src} (\rho_v)
  \Deqlab{netu9}
\end{eqnarray}
となる. 
$\Dvect{i}_d$, $\Dvect{i}_v$ に関して
\begin{eqnarray}
\Dvect{i}_d + \sum \Dvect{i}_v &=& 0
 \Deqlab{netu10} \\
\Dvect{i}_v
 &=& - \rho D_v 
     \left( 
      \nabla q_v + \frac{k_T}{T} \nabla T + \frac{k_p}{p} \nabla p
     \right)
 \Deqlab{netu11}
\end{eqnarray}
が成り立つ. 
但し $D_v$, $k_T$, $k_p$ はそれぞれ拡散係数, 熱拡散比, barodiffusion 係
数である. 
$q_v \rightarrow 0$ 又は $q_v \rightarrow 1$ の極限において $k_T$, $k_p$
はゼロに収束することが知られている. 
即ち, 
\begin{eqnarray}
\Dvect{i}_v
 &\simeq& - \rho D_v \nabla q_v 
 \Deqlab{netu11-1}
\end{eqnarray}
となる. 
\Deqref{netu1-A} 及び \Deqref{netu11-1} を用いて, \Deqref{netu9} から 
$\Dvect{F}^{therm}$, $\Dvect{i}_d$, $\Dvect{i}_v$ を消去すると, 
\begin{eqnarray}
  \hat{c_p} \DD{T}{t}
 &=&
  \Dinv{\rho} \DD{p}{t}
  - \Dinv{\rho} \nabla \cdot \Dvect{F}^{rad}
  - \Dinv{\rho} \nabla \cdot \Dvect{F}^{therm}
  + \frac{\varepsilon}{\rho} 
  \nonumber \\
 & &
    + \Dinv{\rho} c_{pd} T \nabla \cdot \Dvect{i}_d 
    + \Dinv{\rho} 
      \sum \left( c_{pv} T + L_{00} \right) 
      \nabla \cdot \Dvect{i}_v
  \nonumber \\
 & &
  - \Dinv{\rho} 
    \sum 
      \left[
       L_{00} + \left( c_{pv} - c_{pc} \right) T 
      \right] M_{src} (\rho_v)
  \nonumber \\
 &=&
  \Dinv{\rho} \DD{p}{t}
  - \Dinv{\rho} \nabla \cdot \Dvect{F}^{rad}
  + \frac{\varepsilon}{\rho} 
  \nonumber \\
 & &
  - \Dinv{\rho} \nabla \cdot 
    \left[
     c_{pd} T \Dvect{i}_d 
     + \sum \left(c_{pv} T + L_{00} \right) \Dvect{i}_v 
     - \kappa \nabla T
    \right]
  \nonumber \\
 & &
    + \Dinv{\rho} c_{pd} T \nabla \cdot \Dvect{i}_d 
    + \Dinv{\rho} 
      \sum \left( c_{pv} T + L_{00} \right) 
      \nabla \cdot \Dvect{i}_v
  \nonumber \\
 & &
  - \Dinv{\rho} 
    \sum 
       L M_{src} (\rho_v)
  \nonumber \\
 &=&
  \Dinv{\rho} \DD{p}{t}
  - \Dinv{\rho} \nabla \cdot \Dvect{F}^{rad}
  + \frac{\varepsilon}{\rho} 
  + \Dinv{\rho} \nabla \cdot \left( \kappa \nabla T \right)
  \nonumber \\
 & &
    - \Dinv{\rho} \Dvect{i}_d \cdot \nabla \left( c_{pd} T \right) 
    - \Dinv{\rho} \sum \Dvect{i}_v \cdot
      \nabla \left( c_{pv} T + L_{00} \right) 
   - \Dinv{\rho} 
    \sum 
       L M_{src} (\rho_v)
  \nonumber \\
 &=&
  \Dinv{\rho} \DD{p}{t}
  - \Dinv{\rho} \nabla \cdot \Dvect{F}^{rad}
  + \frac{\varepsilon}{\rho} 
  + \Dinv{\rho} \nabla \cdot \left( \kappa \nabla T \right)
  \nonumber \\
 & &
    + \Dinv{\rho} \sum \Dvect{i}_v \cdot \nabla \left( c_{pd} T \right) 
    - \Dinv{\rho} \sum \Dvect{i}_v \cdot
      \nabla \left( c_{pv} T \right) 
  - \Dinv{\rho} 
    \sum 
       L M_{src} (\rho_v)
  \nonumber \\
 &=&
  \Dinv{\rho} \DD{p}{t}
  - \Dinv{\rho} \nabla \cdot \Dvect{F}^{rad}
  + \frac{\varepsilon}{\rho} 
  + \Dinv{\rho} \nabla \cdot \left( \kappa \nabla T \right)
  \nonumber \\
 & &
    - \sum D_v \nabla q_v \cdot
      \nabla \left[ ( c_{pd} - c_{pv} ) T \right]
  - \Dinv{\rho} 
    \sum 
       L M_{src} (\rho_v)
  \nonumber \\
  \Deqlab{netu12}
\end{eqnarray}
となる. 
ここで
\begin{eqnarray}
\hat{c_p} Q_{rad}
  &=& 
    - \Dinv{\rho} \nabla \cdot \Dvect{F}^{rad}, 
      \Deqlab{netu13} \\
\hat{c_p} Q_{dis}
  &=& 
      \frac{\varepsilon}{\rho}, 
      \Deqlab{netu14} \\
\hat{c_p} Q_{cond}
  &=& 
    - \Dinv{\rho} L M_{src} (\rho_v), 
      \Deqlab{netu15} \\
\hat{c_p} D_{T}
 &=& 
    \Dinv{\rho} \nabla \cdot \left( \kappa \nabla T \right)
    - \sum D_v \nabla q_v \cdot
      \nabla \left[ ( c_{pd} - c_{pv} ) T \right]
  \Deqlab{netu16}
\end{eqnarray}
と置くと, 
\begin{equation}
\hat{c_p} \DD{T}{t} - \Dinv{\rho} \DD{p}{t}
 =
  \hat{c_p} \left( Q_{cond} + Q_{rad} + Q_{dis} \right)
  + \hat{c_p} D_{T}
  \Deqlab{netu17}
\end{equation}
が得られる.



\subsection{ 状態方程式の書き換え }


\Deqref{equations:eqsdry}, \Deqref{equations:eqsvapor} より
\begin{eqnarray}
p 
 &=& 
  p_d + \sum p_v 
  \nonumber \\
 &=& 
  \left( \rho_d R_d + \sum \rho_v R_v \right) T 
  \nonumber \\
 &=&
  \left( \rho_d + \sum \rho_v \right)
   \frac{q_d R_d + \sum q_v R_v}{q_d + \sum q_v} T. 
\end{eqnarray}
ここで
\begin{eqnarray}
\langle R \rangle 
 &\equiv& 
  \frac{ q_d R_d + \sum q_v R_v}{ q_d + \sum q_v} \\
\rho_g 
 &\equiv& 
  \rho_d + \sum \rho_v
\end{eqnarray}
と置くと, 
\begin{eqnarray}
p = \rho_g \langle R \rangle T
 \Deqlab{theta7}
\end{eqnarray}
となる. 
但し $R_d$, $R_v$ はそれぞれ非凝結性ガスの気体定数, 凝結性ガスの気体定数
である. 
$\langle R \rangle$ は気体定数の密度の重みつき平均であり, 本文書では平均
気体定数と呼ぶことにする. 



普遍気体定数を $R_{\ast}$ として, 
\begin{eqnarray}
\langle M \rangle 
 &\equiv& 
  \frac{R_{\ast}}{\langle R \rangle}
 = 
  \frac{q_d + \sum q_v}{q_d / M_d + \sum q_v / M_v} 
 \Deqlab{Tfaa} 
\end{eqnarray}
と置く. 
但し $M_d$, $M_v$ はそれぞれ非凝結性ガスの分子量, 凝結性ガスの分子量を表
す. 
また $\langle M \rangle$ を平均分子量と呼ぶことにする. 
本モデルでは $\langle R \rangle$, $\langle M \rangle$ を一定値とみなし, 
以降 $\langle R \rangle_0$, $\langle M \rangle_0$ と書くことにする. 
ここで
\begin{equation}
T_v \equiv \left( \frac{\langle M \rangle_0}{M_d} q_d
     + \sum \frac{\langle M \rangle_0}{ M_v } q_v
         \right) T \Deqlab{Tfb}
\end{equation}
を導入すると, 
\begin{eqnarray}
p 
&=& \rho_g \langle R \rangle T \nonumber \\
&\simeq& \rho_g \langle R \rangle_0 T \nonumber \\
&=& \rho \langle R \rangle_0 
    \frac{\rho_g \langle R \rangle_0}{\rho \langle R \rangle_0} T \nonumber
    \\
&=& \rho \langle R \rangle_0 
    \frac{ \rho_d R_d + \sum \rho_v R_v }
         {\rho \langle R \rangle_0 } T \nonumber \\
&=& \rho \langle R \rangle_0
    \left( \frac{\langle M \rangle_0}{M_d} q_d
     + \sum \frac{\langle M \rangle_0}{ M_v } q_v
         \right) T \nonumber \\
&=& \rho \langle R \rangle_0 T_v
  \Deqlab{TfA}
\end{eqnarray}
と近似される\footnote{
\Deqref{TfA} で行なった近似の精度については現在調査中である. 
}. 
更に温位
\begin{eqnarray}
\theta 
 \equiv 
  T \left( \frac{p_{0}}{p} \right)^{ \frac{ \langle R \rangle_0}
   {\langle c_{p} \rangle_0} } 
 \Deqlab{theta4}
\end{eqnarray}
及びエクスナー関数
\begin{eqnarray}
\Pi 
 \equiv \frac{T}{\theta}
 = \left( \frac{p}{p_{0}} \right)^{ \frac{ \langle R \rangle_0}
   {\langle c_{p} \rangle_0} }
 \Deqlab{theta5}
\end{eqnarray}
並びに
\begin{equation}
\theta_v 
 \equiv \frac{T_v}{\Pi} 
 =  \left( \frac{\langle M \rangle_0}{M_d} q_d
     + \sum \frac{\langle M \rangle_0}{ M_v } q_v
         \right) \theta
 \Deqlab{Tfd}
\end{equation}
を導入すると, 
\begin{eqnarray}
p &=& \rho \langle R \rangle_0 \Pi \theta_v \Deqlab{Tff}
\end{eqnarray}
となる. 
ここで
\begin{eqnarray}
\langle c_p \rangle &=& \frac{q_d c_{pd} + \sum q_v c_{pv}}{
 q_d + \sum q_v } 
 \Deqlab{cp-def}
\end{eqnarray}
であり, 本文書では $\langle c_{p} \rangle$ を平均定圧比熱と呼ぶことにす
る. 
$c_{pd}$, $c_{pv}$ はそれぞれ非凝結性ガスの定圧比熱, 凝結性ガスの定圧比
熱を表す. 
また本モデルでは $\langle c_{p} \rangle$ を一定値とみなし, $\langle
c_{p} \rangle_0$ と表すことにする.  
%
\Deqref{Tff}に\Deqref{theta5}を代入して $p$ を消去し, $\rho$ について整
理すると, 
\begin{eqnarray}
\rho 
 = \frac{
    p_{0} \Pi^{ \langle c_v \rangle_0 / \langle R \rangle_0}}
     { \langle R \rangle_0 \theta_v}
 = \frac{p_0 \Pi^{\langle c_v \rangle_0 / \langle R
   \rangle_0 } }{
    \langle R \rangle_0 \left( q_d  \langle M \rangle_0 / M_d
    + \sum q_v \langle M \rangle_0 / M_v \right) \theta }
 \Deqlab{Tfg}
\end{eqnarray}
となる. 


\subsection{
熱力学第一法則の書き換え -- 温位の式の導出 --
}


\Deqref{theta1} において, 凝結物が気体に比べて十分少なく, 
$q_d c_{pd} \gg q_c c_{pc}$, $q_d c_{pv} \gg q_c c_{pc}$ が成り立つとする
と, $\hat{c_p} \approx \langle c_p \rangle$ となるので, 
\begin{equation}
\langle c_p \rangle \DD{T}{t} - \Dinv{\rho} \DD{p}{t}
\approx \langle c_p \rangle \left( Q_{cond} + Q_{rad} + Q_{dis} \right)
 \Deqlab{theta2}
\end{equation}
と近似される. 
平均定圧比熱を定数で近似すると, 
\begin{equation}
\langle c_p \rangle_0 \DD{T}{t} - \Dinv{\rho} \DD{p}{t}
\approx \langle c_p \rangle_0 \left( Q_{cond} + Q_{rad} + Q_{dis} \right)
 \Deqlab{theta2-1}
\end{equation}
となる. 
さらに凝結物は全密度に比べて十分小さいとみなすと, 
\begin{equation}
\langle c_p \rangle_0 \DD{T}{t} - \Dinv{\rho_g} \DD{p}{t}
\approx \langle c_p \rangle_0 \left( Q_{cond} + Q_{rad} + Q_{dis} \right) 
 \Deqlab{theta3}
\end{equation}
と近似される. 
$\theta$ のラグランジュ微分をとると,  
\begin{eqnarray}
 \DD{\theta}{t}
  &=& 
   \left( \frac{p_{0}}{p} \right)^{ \frac{ \langle R \rangle_0}
   {\langle c_{p} \rangle} } \DD{T}{t}
   - \frac{ \langle R \rangle_0 }{ \langle c_p \rangle_0} 
     \frac{T}{p} 
    \left( \frac{p_{0}}{p} \right)^{ \frac{ \langle R \rangle_0}
    {\langle c_{p} \rangle} } \DD{p}{t}
    \nonumber \\
  &=& 
   \frac{\theta}{ \langle c_p \rangle_0 T }
    \left(
     \langle c_p \rangle_0 \DD{T}{t} - \Dinv{\rho_g} \DD{p}{t}
    \right)
     \nonumber \\
  &=&
   \frac{ Q_{cond} + Q_{rad} + Q_{dis} }{\Pi}
 \Deqlab{theta6}
\end{eqnarray}
となる. 
但し式変形の途中で $p = \rho_g \langle R \rangle_0 T$ を用いた. 
%非断熱加熱として凝結加熱 $Q_{cond}$, 放射加熱 $Q_{rad}$, 散逸加熱
%$Q_{dis}$ を考慮し, また
乱流拡散を考慮すると, 
\begin{eqnarray}
\DD{\theta}{t}
 = \Dinv{\Pi} \left( Q_{cond} + Q_{rad} + Q_{dis} \right) + D_{\theta}
 \Deqlab{theta8}
\end{eqnarray}
が得られる. 

予報変数として温位を用いる際には平均気体定数, 平均定圧比熱があまり大きく
変化することを許容しないことを前提とするので, 計算の適用範囲に制約が加わ
ることに常に注意しなければならない. 


\subsection{ 相当温位の導出 }

\Deqref{theta8}, \Deqref{sphm7} において, 放射加熱・散逸加熱・乱流拡散・
雲粒落下を無視すると, 
\begin{eqnarray}
d \theta 
 &=& \Dinv{\Pi} Q_{cond}
  = - \sum \frac{L}{\Pi \langle c_p \rangle_0} M_{src} (\rho_v), 
 \Deqlab{epot1} \\
d q_v
 &=& \Dinv{\rho} M_{src} (\rho_v)
 \Deqlab{epot2}
\end{eqnarray}
となる. 
ここで $d q_v / q_v \gg d \Pi / \Pi$ であると仮定すると
\footnote{
地球大気ではこの近似が成立するとされているが, 他の惑星大気で成立するかど
うかは常に確かめる必要がある. 
}, \Deqref{epot1}, \Deqref{epot2} より
\begin{eqnarray}
d \theta
 \simeq
  \sum \frac{L q_v}{\Pi \langle c_p \rangle}_0
\Deqlab{epot3}
\end{eqnarray}
となる. 
\Deqref{epot3}において
\begin{eqnarray}
\theta_e
 =
 \theta + \sum \frac{L q_v}{\Pi \langle c_p \rangle}_0
\Deqlab{epot4}
\end{eqnarray}
と置くと, $\theta_e$ は近似的に保存量となる. 


\subsection{圧力方程式の導出}


圧力方程式は密度の式と連続の式を組み合わせることで得られる. 
\Deqref{Tfg} のラグランジュ微分をとると, 
%
\begin{eqnarray}
\DD{\rho}{t}
&=&
\DD{}{t} \left[ \frac{p_0 \Pi^{\langle c_v \rangle_0 / \langle R
   \rangle}}{
    \langle R \rangle_0 \left( q_d  \langle M \rangle_0 / M_d
    + \sum q_v \langle M \rangle_0 / M_v \right) \theta }
\right] \nonumber \\
&=&
\frac{p_0 \Pi^{(\langle c_v \rangle_0 - 1) / \langle R
   \rangle}}{
    \langle R \rangle_0 \left( q_d  \langle M \rangle_0 / M_d
    + \sum q_v \langle M \rangle_0 / M_v \right) \theta }
    \frac{\langle c_v \rangle_0}{\langle R \rangle_0}
    \DD{\Pi}{t}
    \nonumber \\
& & - \frac{p_0 \Pi^{\langle c_v \rangle_0 / \langle R
   \rangle}}{
    \langle R \rangle_0 \left( q_d  \langle M \rangle_0 / M_d
    + \sum q_v \langle M \rangle_0 / M_v \right) \theta^2 } 
      \DD{\theta}{t}
   \nonumber \\
& & - \frac{p_0 \Pi^{\langle c_v \rangle_0 / \langle R
   \rangle}}{
    \langle R \rangle_0 \left( q_d  \langle M \rangle_0 / M_d
    + \sum q_v \langle M \rangle_0 / M_v \right)^2 \theta }
   \left( \frac{\langle M \rangle_0}{M_d} \DD{q_d}{t}
    + \sum \frac{\langle M \rangle_0}{M_v} \DD{q_v}{t}  \right) \nonumber \\
&=&
\langle c_p \rangle_0 \theta_v
\left( \frac{\langle c_v \rangle_0}{\langle c_p \rangle_0 \langle R \rangle_0
 \Pi \theta_v} \right)
\rho \DD{\Pi}{t}
- \frac{\rho}{\theta} \DD{\theta}{t} \nonumber \\
& & - \frac{\rho}{q_d  \langle M \rangle_0 / M_d
    + \sum q_v \langle M \rangle_0 / M_v}
\left( \frac{\langle M \rangle_0}{M_d} \DD{q_d}{t}
    + \sum \frac{\langle M \rangle_0}{M_v} \DD{q_v}{t}  \right) 
 \nonumber \\
&=&
\frac{\langle c_p \rangle_0 \rho \theta_v}{\langle c_s \rangle^2}
 \DD{\Pi}{t}
- \frac{\rho}{\theta} \DD{\theta}{t} \nonumber \\
& & - \frac{\rho}{q_d  \langle M \rangle_0 / M_d
    + \sum q_v \langle M \rangle_0 / M_v}
\left( \frac{\langle M \rangle_0}{M_d} \DD{q_d}{t}
    + \sum \frac{\langle M \rangle_0}{M_v} \DD{q_v}{t}  \right).
\Deqlab{pressure:theta-pi:drho}
\end{eqnarray}
%
となる. 
ここで $\langle c_s \rangle$ は音速であり, 
\begin{equation}
\langle c_s \rangle^2
= \frac{\langle c_p \rangle_0 \langle R \rangle_0
 \Pi \theta_v}{\langle c_v \rangle_0} \Deqlab{SoundVel}
\end{equation}
である. 
\Deqref{pressure:theta-pi:drho} に \Deqref{renzoku5} を適用すると, 
%
\begin{eqnarray}
\DD{\Pi}{t}
&=& \frac{\langle c_s \rangle^2}{\langle c_p \rangle_0 \theta_v}
\Bigg[ \Dinv{\rho} \DD{\rho}{t} + \Dinv{\theta} \DD{\theta}{t} \nonumber
\\
& &
+ \Dinv{q_d  \langle M \rangle_0 / M_d
    + \sum q_v \langle M \rangle_0 / M_v}
\left( \frac{\langle M \rangle_0}{M_d} \DD{ q_d}{t}
    + \sum \frac{\langle M \rangle_0}{M_v} \DD{ q_v}{t}  \right)
\Bigg] \nonumber \\
&=&
 \frac{\langle c_s \rangle^2}{\langle c_p \rangle_0 \theta_v}
\Bigg[ - \DP{u_j}{x_j} + \Dinv{\rho} \sum M_{fall} (\rho_r) +
\Dinv{\theta} \DD{\theta}{t} \nonumber \\
& &
+ \Dinv{q_d  \langle M \rangle_0 / M_d
    + \sum q_v \langle M \rangle_0 / M_v}
\left( \frac{\langle M \rangle_0}{M_d} \DD{ q_d}{t}
    + \sum \frac{\langle M \rangle_0}{M_v} \DD{ q_v}{t}  \right)
\Bigg]
%%%
\Deqlab{pressure:theta-p:p}
\end{eqnarray}
%
となり, 圧力方程式が得られる. 


\subsection{運動方程式の書き換え}


運動方程式の圧力勾配は, 温位とエクスナー関数を用いることで得られる. 
\Deqref{theta5}, \Deqref{Tff} より
%
\begin{eqnarray}
 \Dinv{\rho} \DP{p}{x_i} 
  &=& 
   \frac{ \langle R \rangle_0  \Pi \theta_v }{p} 
   \DP{ \left( 
     p_{0} \Pi^{ \langle c_{p} \rangle / \langle R \rangle_0}
    \right)}{x_i}
  \nonumber \\
  &=& 
   \frac{ \langle R \rangle_0 \Pi \theta_v }{p} 
   \left(
   \frac{p_{0} {c_{p}}_{d}}{R_{d}} \Pi^{ \langle c_{p} \rangle / \langle
   R \rangle_0 - 1} 
   \right) \DP{\Pi}{x_i}
  \nonumber \\
  &=& 
   \frac{ \langle R \rangle_0  \Pi \theta_v }{p} 
   \left(
   \frac{ \langle c_{p} \rangle}{ \langle R \rangle_0} p \Pi^{-1} 
   \right) \DP{\Pi}{x_i}
   \nonumber \\
  &=& 
   \langle c_{p} \rangle \theta_v  \DP{\Pi}{x_i}
   \nonumber \\
\end{eqnarray}
%
となる. 
以上より
\begin{eqnarray}
&& \DD{u}{t} 
    = - \langle c_{p} \rangle \theta_{v} \DP{\Pi}{x}
      + D_u 
  \Deqlab{equations:u2} \\
&& \DD{v}{t} 
    = - \langle c_{p} \rangle \theta_{v} \DP{\Pi}{y}
      + D_v 
  \Deqlab{equations:v2} \\
&& \DD{w}{t} 
    = - \langle c_{p} \rangle \theta_{v} \DP{\Pi}{z}
      - g + D_w
  \Deqlab{equations:w2} 
\end{eqnarray}
が得られる. 

\subsection{
温位 $\theta$, エクスナー関数 $\Pi$, 風速 $u, v, w$ を予報変数とする場合
の方程式系
}


以上より, 3 次元大気の状態を
温位 $\theta$, エクスナー関数 $\Pi$, 風速 $u, v, w$, 密度 $\rho$ 
で表現する場合, 基礎方程式系は以下のようになる. 

%
\begin{description}
%
 \item[運動方程式]~
%
\begin{eqnarray}
&&\DD{u}{t} = - \langle c_p \rangle_0 \theta_v \DP{\Pi}{x} + D_u
  \Deqlab{equations:u-w-pi-theta:u} \\
&&\DD{v}{t} = - \langle c_p \rangle_0 \theta_v \DP{\Pi}{y} + D_v
  \Deqlab{equations:u-w-pi-theta:v} \\
&&\DD{w}{t} = - \langle c_p \rangle_0 \theta_v \DP{\Pi}{z} - g + D_w
  \Deqlab{equations:u-w-pi-theta:w} 
\end{eqnarray}
%
 \item[圧力方程式]~
%
\begin{eqnarray}
\DP{\Pi}{t}
&=&
 \frac{\langle c_s \rangle^2}{\langle c_p \rangle_0 \theta_v}
\Bigg[ - \DP{u_j}{x_j} + \Dinv{\rho} \sum M_{fall} (\rho_r) +
\Dinv{\theta} \DD{\theta}{t} \nonumber \\
& &
+ \Dinv{q_d  \langle M \rangle_0 / M_d
    + \sum q_v \langle M \rangle_0 / M_v}
\left( \frac{\langle M \rangle_0}{M_d} \DD{ q_d}{t}
    + \sum \frac{\langle M \rangle_0}{M_v} \DD{ q_v}{t}  \right)
\Bigg] \nonumber \\
  \Deqlab{equations:u-w-pi-theta:continue}
\end{eqnarray}
%
 \item[状態方程式]~
%
\begin{eqnarray}
\rho
= \frac{p}{\langle R \rangle_0 \Pi \theta_v}
= \frac{p_0 \Pi^{\langle c_{v} \rangle / \langle R \rangle_0}}{\langle R
\rangle \theta_v}
  \Deqlab{equations:u-w-pi-theta:rho}
\end{eqnarray}
%
 \item[熱の式]~
%
\begin{eqnarray}
 \DD{\theta}{t} 
   = \Dinv{\Pi} \left( Q_{cond} + Q_{rad} + Q_{dis} \right)
     + D_{\theta}
  \Deqlab{equations:u-w-pi-theta:theta}
\end{eqnarray}
%
 \item[凝結性ガスおよび凝結物の比湿の式]~
%
\begin{eqnarray}
%&& \DD{q_{d}}{t}
%   = - \frac{q_{d}}{\rho} \sum M_{fall}(\rho_r)
%     + D_{q_d}
%  \Deqlab{equations:u-w-pi-theta:qd} \\
&& 
\DD{q_{v}}{t}
   = - \frac{q_{v}}{\rho} \sum M_{fall}(\rho_r)
     + \Dinv{\rho} M_{src}(\rho_v)
     + D_{q_v}
  \Deqlab{equations:u-w-pi-theta:qv} \\
&& 
\DD{q_{c}}{t}
   = - \frac{q_{c}}{\rho} \sum M_{fall}(\rho_r)
     + \Dinv{\rho} M_{src}(\rho_c)
     + D_{q_c}
  \Deqlab{equations:u-w-pi-theta:qc} \\
&& 
\DD{q_r}{t}
   = - \frac{q_r}{\rho} \sum M_{fall}(\rho_r)
     + \Dinv{\rho} M_{fall}(\rho_r) + \Dinv{\rho} M_{src}(\rho_r)
     + D_{q_r}
  \Deqlab{equations:u-w-pi-theta:qr}
\end{eqnarray}
%
\end{description}
%
ただし, エクスナー関数 $\Pi$ は, 
%
\begin{eqnarray}
\Pi \equiv \left( \frac{p}{p_0} \right)^{\langle R \rangle_0 / \langle
 c_p \rangle_0}
\end{eqnarray}
%
であり, 音速 $\langle c_{s} \rangle$ は
%
\begin{eqnarray}
\langle c_s \rangle^2
 &=& \frac{\langle c_p \rangle_0 \langle R \rangle_0
 \Pi \theta_v}{\langle c_v \rangle_0}
\end{eqnarray}
%
である. 



