%\subsection{離散表現: 氷相を含む Relaxed Arakawa-Schubert 積雲パラメタリゼーション}
%\subsection{離散表現 (風上差分版)}
\subsection{離散表現}

Relaxed Arakawa-Schubert 積雲パラメタリゼーションの離散化は, 
Moorthi and Suarez (1992) の方法を拡張することで導出する. 
また, この節では, 凝結物のエントレインメントも含む定式化について述べる.

\ref{sec:math:ras_with_ice:p-axis} 節の方程式を下のように離散化する.
%
\begin{align}
  \left( \DP{s_k}{t} \right)_{c} 
    &= \sum_{k'} \left( \DP{s_k}{t} \right)_{c,k'}
     = \sum_{k'} M_{B,k'} \Gamma_{s,k',k}
     = \sum_{k'} m_{B,k'} \Delta \lambda_{k'} \Gamma_{s,k',k}
\\
  \left( \DP{q_{v,k}}{t} \right)_{c}
    &= \sum_{k'} \left( \DP{q_{v,k}}{t} \right)_{c,k'}
     = \sum_{k'} M_{B,k'} \Gamma_{q_v,k',k}
     = \sum_{k'} m_{B,k'} \Delta \lambda_{k'} \Gamma_{q_v,k',k}
\\
  \left( \DP{q_{l,k}}{t} \right)_{c}
    &= \sum_{k'} \left( \DP{q_{l,k}}{t} \right)_{c,k'}
     = \sum_{k'} M_{B,k'} \Gamma_{q_l,k',k}
     = \sum_{k'} m_{B,k'} \Delta \lambda_{k'} \Gamma_{q_l,k',k}
\\
  \left( \DP{q_{i,k}}{t} \right)_{c}
    &= \sum_{k'} \left( \DP{q_{i,k}}{t} \right)_{c,k'}
     = \sum_{k'} M_{B,k'} \Gamma_{q_i,k',k}
     = \sum_{k'} m_{B,k'} \Delta \lambda_{k'} \Gamma_{q_i,k',k}
\end{align}
%
ここで, $k'$ は雲頂のインデックスである.
%
また, 
%
\begin{align}
  \Gamma_{s,k',k} & = 
    \begin{cases}
      \displaystyle
        - \frac{g}{\Delta p_k}
          \eta_{k',k+\frac{1}{2}} \left( s_k - s_{k+1} \right)
        \\ \hspace{5mm}
      \displaystyle
        - g L   D_{k',k} q_{l,k',k'}^c \left( 1 - r_{l,k'} \right)
        - g L_i D_{k',k} q_{i,k',k'}^c \left( 1 - r_{i,k'} \right)
    &
      ( k \le k' ) 
    \\
      0
    &
      ( k > k' ) 
    \end{cases}
\\
  \Gamma_{q_v,k',k} &=
    \begin{cases}
      \displaystyle
      - \frac{g}{\Delta p_k}
        \eta_{k',k+\frac{1}{2}} \left( q_{v,k} - q_{v,k+1} \right)
        \\ \hspace{5mm}
      \displaystyle
        + g D_{k',k} \left( q_{v,k'}^c - q_{v,k'} \right)
        + g D_{k',k} q_{l,k'}^c \left( 1 - r_{l,k'} \right)
    &
      ( k \le k' ) 
    \\
      0
    &
      ( k > k' ) 
    \end{cases}
\\
  \Gamma_{q_l,k',k} &=
    \begin{cases}
      \displaystyle
      - \frac{g}{\Delta p_k}
        \eta_{k',k+\frac{1}{2}} \left( q_{l,k} - q_{l,k+1} \right)
        \\ \hspace{5mm}
      \displaystyle
        + g D_{k',k} \left( q_{l,k'}^c - q_{l,k'} \right)
        - g D_{k',k} q_{l,k'}^c \left( 1 - r_{l,k'} \right)
        + g D_{k',k} q_{i,k'}^c \left( 1 - r_{i,k'} \right)
    &
      ( k \le k' ) 
    \\
      0
    &
      ( k > k' ) 
    \end{cases}
\\
  \Gamma_{q_i,k',k} &=
    \begin{cases}
      \displaystyle
      - \frac{g}{\Delta p_k}
        \eta_{k',k+\frac{1}{2}} \left( q_{i,k} - q_{i,k+1} \right)
        \\ \hspace{5mm}
      \displaystyle
        + g D_{k',k} \left( q_{i,k'}^c - q_{i,k'} \right)
        - g D_{k',k} q_{i,k'}^c \left( 1 - r_{i,k'} \right)
    &
      ( k \le k' ) 
    \\
      0
    &
      ( k > k' ) 
    \end{cases}
\\
  D_{k',k} &= \frac{1}{\Delta p_{k}} \eta_{k',k} \delta_{k,k'}
\end{align}
%
である
\footnote{
  Moorthi and Suarez (1992) では, この部分の離散化に中心差分を用いているが, 
  ここでは風上差分を用いている. これにより, 物理量が負になることを避ける
  ことができる. 
}.

\begin{align}
  \eta_{k',k-\frac{1}{2}} - \eta_{k',k+\frac{1}{2}}
    &= - \beta_{k} \theta_k \lambda_{k'}
\\
  \eta_{k',k'-\frac{1}{2}} - \eta_{k',k'}
    &= - \beta_{k'}' \theta_{k'} \lambda_{k'}
\\
  \beta_k &= \frac{C_p}{g} \left( P_{k-\frac{1}{2}} - P_{k+\frac{1}{2}} \right)
\\
  \beta_k' &= \frac{C_p}{g} \left( P_{k-\frac{1}{2}} - P_{k} \right)
\end{align}
\begin{align}
    \eta_{k',k-\frac{1}{2}} \left( h_{k',k-\frac{1}{2}}^c - L_i q_{i,k',k-\frac{1}{2}}^c \right)
 &- \eta_{k',k+\frac{1}{2}} \left( h_{k',k+\frac{1}{2}}^c - L_i q_{i,k',k+\frac{1}{2}}^c \right) \nonumber \\
    &=   \left( \eta_{k',k-\frac{1}{2}} - \eta_{k',k+\frac{1}{2}} \right)
         \left( h_{k} - L_i q_{i,k} \right)
\\
    \eta_{k',k'-\frac{1}{2}} \left( h_{k',k'-\frac{1}{2}}^c - L_i q_{i,k',k'-\frac{1}{2}}^c \right)
 &- \eta_{k',k'} \left( h_{k',k'}^c - L_i q_{i,k',k'}^c \right) \nonumber \\
    &=   \left( \eta_{k',k'-\frac{1}{2}} - \eta_{k',k'} \right)
         \left( h_{k'} - L_i q_{i,k'} \right)
\end{align}
\begin{align}
    \eta_{k',k-\frac{1}{2}} \left( q_{v,k',k-\frac{1}{2}}^c + q_{l,k',k-\frac{1}{2}}^c + q_{i,k',k-\frac{1}{2}}^c \right)
 &- \eta_{k',k+\frac{1}{2}} \left( q_{v,k',k+\frac{1}{2}}^c + q_{l,k',k+\frac{1}{2}}^c + q_{i,k',k+\frac{1}{2}}^c \right) \nonumber \\
    &=   \left( \eta_{k',k-\frac{1}{2}} - \eta_{k',k+\frac{1}{2}} \right)
         \left( q_{v,k} + q_{l,k} + q_{i,k} \right)
\\
    \eta_{k',k'-\frac{1}{2}} \left( q_{v,k',k'-\frac{1}{2}}^c + q_{l,k',k'-\frac{1}{2}}^c + q_{i,k',k'-\frac{1}{2}}^c \right)
 &- \eta_{k',k'} \left( q_{v,k',k'}^c + q_{l,k',k'}^c + q_{i,k',k'}^c \right) \nonumber \\
    &=   \left( \eta_{k',k'-\frac{1}{2}} - \eta_{k',k'} \right)
         \left( q_{v,k'} + q_{l,k'} + q_{i,k'} \right)
\end{align}

雲内の物理量, $T_{\lambda}^c(z)$, $q_{v,\lambda}^c(z)$, 
$q_{l,\lambda}^c(z)$, $q_{i,\lambda}^c(z)$, は, 
\ref{eq:math:ras_with_ice:cloud_eq_qt}, 
\ref{eq:math:ras_with_ice:cloud_eq_condstatene}, 
\ref{eq:math:ras_with_ice:cloud_prop_assumption1}, 
\ref{eq:math:ras_with_ice:cloud_prop_assumption2} 
を連立させて解くことで求めるが, これらは非線形方程式であり, 解析的には
解けないため, 繰り返し法により求める. 
また, 安定に解くために下のような仮定を置く.
%
\begin{align}
  \alpha_{k} &= \alpha\left( T_k^c \right)
\end{align}
%
... 続きはまた後で.


$k'$ 層目に雲頂を持つ雲のエントレインメントパラメータ $\lambda_{k'}$ は下のように
離散化する
\footnote{
  実際には, この $\lambda_{k'}$ は, 二次方程式の解の片方である.
  二つの解のうちのこちらが, 氷を含まない場合の定式化から得られる解 (この
  ときは, $\lambda_{k'}$ は一次方程式の解) と等しい. 
  なお, 二次方程式のもう一つの解は $\lambda_{k'} = - \frac{1}{B_{k'}} < 0$ である.
}.
%
\begin{align}
  \lambda_{k'}
    &= \frac{
         \left( h_{k_{mlt}} - L_i q_{i,k_{mlt}} \right) - h_{k',k'}^c + L_i \left( 1 - \alpha \right) \left\{ \left( q_{v,k_{mlt}} + q_{l,k_{mlt}} + q_{i,k_{mlt}} \right) - q_{v,k'}^c \right\} }
            { E_{k'} - L_i \left( 1 - \alpha \right) \left( C_{k'} - B_{k'} q_{v,k'}^c \right) }
\\
  B_{k'} &= \beta_k' \theta_{k'} + \sum_{l=k_{mlt}+1}^{k'-1} \beta_l \theta_l
\\
  C_{k'} &= \beta_k' \theta_{k'} \left( q_{v,k'} + q_{l,k'} + q_{i,k'} \right)
       + \sum_{l=k_{mlt}+1}^{k'-1} \beta_l \theta_l \left( q_{v,l} + q_{l,l} + q_{i,l} \right)
\\
  E_{k'} &= \beta_k' \theta_{k'} \left\{ h_{k'}^c - \left( h_{k'} - L_i q_{i,k'} \right) \right\}
  + \sum_{l=k_{mlt}+1}^{k'-1} \beta_l \theta_l \left\{ h_{k'}^c - \left( h_l - L_i q_{i,l} \right) \right\}
\\
  h_{k'}^c &= h_{k'}^*
\\
  q_{v,k'}^c &= q_{v,k'}^*
\end{align}
%
ここで, $k_{mlt}$ は混合層最上層のインデックスである. 

$k'$ 層目に雲頂を持つ雲の雲仕事関数 $A_{k'}$ は下のように離散化する.
%
\begin{align}
  A_{k'}
    &=   \sum_{l=l_{\rm mlt}+1}^{k'-1}
           \left\{
               \mu_l' \eta_{k',l+\frac{1}{2}}
                 \left( s_{k',l+\frac{1}{2}}^c - s_l \right)
             + \epsilon_l' \eta_{k',l-\frac{1}{2}}
                 \left( s_{k',l-\frac{1}{2}}^c - s_l \right)
           \right\}
       + \epsilon_{k'}' \eta_{k',k'-\frac{1}{2}}
           \left( s_{k',k'-\frac{1}{2}}^c - s_{k'} \right)
\\
  \mu_l' &= \frac{1}{P_l} \left( P_l - P_{l+\frac{1}{2}} \right)
\\
  \epsilon_l' &= \frac{1}{P_l} \left( P_{l-\frac{1}{2}} - P_l \right)
\end{align}
%
あるいは, 
%
\begin{align}
  A_{k'}
    &=   \sum_{l=l_{\rm mlt}+1}^{k'-1}
           \left\{
               \mu_l \eta_{k',l+\frac{1}{2}}
                 \left( h_{k',l+\frac{1}{2}}^c - h_l^* \right)
             + \epsilon_l \eta_{k',l-\frac{1}{2}}
                 \left( h_{k',l-\frac{1}{2}}^c - h_l^* \right)
           \right\}
       + \epsilon_{k'} \eta_{k',k'-\frac{1}{2}}
           \left( h_{k',k'-\frac{1}{2}}^c - h_{k'}^* \right)
\\
  \mu_l &= \frac{1}{1+\gamma_l}\frac{1}{P_l} \left( P_l - P_{l+\frac{1}{2}} \right)
\\
  \epsilon_l &= \frac{1}{1+\gamma_l}\frac{1}{P_l} \left( P_{l-\frac{1}{2}} - P_l \right)
\end{align}


雲仕事関数の積雲の効果による時間変化率は, 有限差分によって求める
\footnote{
  氷相を含まない場合には, $A_{k'}$ を時間微分することで求めていた. 
  しかし, 氷相を含む場合には雲の物理量, $h^c$, を環境場の物理量で表現する
  ことができない... と思われる... ため, ここでは有限差分によって力づくで
  求めている.
}. 
%
\begin{align}
  \left(\DD{A_{k'}}{t}\right)_c &= \frac{1}{m_{B,k'}^*} \frac{A_{k'}^* - A_{k'}}{\Delta t}
\end{align}
%
ここで, $m_{B,k'}^*$ は $\left(\DD{A_{k'}}{t}\right)_c$ を求めるために
与えた雲底での質量フラックスであり, 小さな値\footnote{具体的な値は試行錯誤により決める.}を与える. 
$A_{k'}^*$ は $m_{B,k'}^*$ を与えて積雲の効果を更新した上で更新した
温度, 比湿, 雲水と雲氷の混合比を使って計算した雲仕事関数である. 



なお, Relaxed Arakawa-Schubert パラメタリゼーションでは, 
$\displaystyle \left(\DP{A_\lambda}{t}\right)_{LS}$ は下のように与える.
%
\begin{align}
  \left(\DP{A_\lambda}{t}\right)_{LS} 
    &\sim \frac{1}{\tau_{RAS}} \left\{ A_{k'} - A_{eq} \right\}
\end{align}
%
ここで, $\tau_{RAS}$ は緩和時定数, 
$A_{eq}$ は雲仕事関数の「平衡値」である. 
また, $A_{k'}$ は格子スケール循環, 放射, 乱流混合などによる
寄与を経た後の雲仕事関数である.


メモ: 雲の物理量の計算について書いてない.