% 表題   DCPAM5  力学過程  鉛直離散化
%
% 履歴
%\Drireki{2009/03/03 森川靖大}
%\Drireki{2008/06/14 森川靖大}
%\Drireki{1993/03/18 沼口敦・保坂征宏}
%

\section{離散表現: 鉛直離散化}

ここでは支配方程式を鉛直方向に離散化する. 
Arakawa and Suarez(1983) に従って,
\Deqref{連続の式}〜\Deqref{水蒸気の式}
を鉛直方向に差分によって離散化する. 
各方程式の離散化表現は次のようになる.


\subsection{連続の式, 鉛直速度}
\vspace*{-5ex}
\begin{align}
  \DP{\pi}{t}
 &= - \sum_{k=1}^{K} ( D_k + \Dvect{v}_k \cdot \nabla \pi ) 
       \Delta  \sigma_k, \\
  \dot{\sigma}_{k-1/2}
 &= - \sigma_{k-1/2} \DP{\pi}{t}
   - \sum_{l=k}^{K} ( D_l + \Dvect{v}_l \cdot \nabla \pi )          
       \Delta  \sigma_l \qquad (k = 2, \cdots, K),  \\
 \dot{\sigma}_{1/2} &= \dot{\sigma}_{K+1/2} = 0. 
\end{align}
%
ここで,
%
\begin{align}
 \Dvect{v}_k \cdot \nabla \pi =
  \frac{U_k}{a (1-\mu^2)} \DP{\pi}{\lambda}
  + \frac{V_k}{a (1-\mu^2)} (1-\mu^2) \DP{\pi}{\mu}. 
\end{align}
%

\subsection{静水圧の式}
\vspace*{-5ex}
\begin{align}
 \begin{split}
  \Phi_{1} & = \Phi_{s} + C_{p} ( \sigma_{1}^{-\kappa} - 1  ) T_{v,1} \\
           & = \Phi_{s} + C_{p} \alpha_{1} T_{v,1}.
 \end{split}
\end{align}
%
\begin{align}
 \begin{split}
  \Phi_k - \Phi_{k-1} 
    & = C_{p}
      \left[ \left( \frac{ \sigma_{k-1/2} }{ \sigma_k } \right)^{\kappa}
             - 1
      \right] T_{v,k} 
      + C_{p}
      \left[ 1- 
         \left( \frac{ \sigma_{k-1/2} }{ \sigma_{k-1} } \right)^{\kappa}
      \right] T_{v,k-1} \\
  & = C_{p} \alpha_k T_{v,k} + C_{p} \beta_{k-1} T_{v,k-1}.
 \end{split}
\end{align}
%
ここで,
%
\begin{align}
 \Deqlab{静水圧係数α}
 \alpha_k  & = \left( \frac{ \sigma_{k-1/2} }
                               { \sigma_k } \right)^{\kappa} -1 , \\
 \Deqlab{静水圧係数β}
 \beta_k   & = 1- \left( \frac{ \sigma_{k+1/2} }
                               { \sigma_k } \right)^{\kappa} , \\
 \Deqlab{地表ジオポテンシャル}
 \Phi_{s}  & = gz_{s}
\end{align}
%
であり, $z_{s}$は地表面高度である. 


\subsection{運動方程式}
\vspace*{-5ex}
\begin{align}
  \Deqlab{渦度方程式の鉛直差分}
  \DP{\zeta_k}{t}
    & =   \Dinv{a}
          \left(
              \Dinv{1 - \mu^2} \DP{{V_A}_{,k}}{\lambda}
            - \DP{{U_A}_{,k}}{\mu}
          \right)
        + {\cal D}(\zeta_k), 
   \\
  \Deqlab{発散方程式の鉛直差分}
  \DP{D_k}{t} 
    & =   \Dinv{a}
          \left(
              \Dinv{1 - \mu^2} \DP{{U_A}_{,k}}{\lambda}
            + \DP{{V_A}_{,k}}{\mu}
          \right)
        - \nabla^{2}_{\sigma}
         ( \Phi_k + C_{p} \hat{\kappa}_k \overline{T}_k \pi 
           + (\mbox{\sl KE})_k )
        + {\cal D}(D_k). 
\end{align}
%
ここで,
%
\begin{align}
 \begin{split}
  {U_A}_{,1}
    & =   ( \zeta_1 + f ) V_1 
        - \frac{1}{2 \Delta \sigma_1} 
                 \dot{\sigma}_{3/2} ( U_1     - U_2 ) 
        - \frac{C_{p} \hat{\kappa}_1 T_{v,1}'}{a}
            \DP{\pi}{\lambda} 
        + {\cal F}_{\lambda, 1}  \cos \varphi, \\
  {U_A}_{,k}
    & =   ( \zeta_k + f ) V_k 
        - \frac{1}{2 \Delta \sigma_k} 
             [   \dot{\sigma}_{k-1/2} ( U_{k-1} - U_k   )
               + \dot{\sigma}_{k+1/2} ( U_k     - U_{k+1} ) ] \\
    & \quad
        - \frac{C_{p} \hat{\kappa}_k T_{v,k}'}{a} 
            \DP{\pi}{\lambda} 
        + {\cal F}_{\lambda, k} \cos \varphi, \qquad (k=2, \cdots, K-1) \\
  {U_A}_{,K}
    & =   ( \zeta_K + f ) V_K 
        - \frac{1}{2 \Delta \sigma_K} 
                 \dot{\sigma}_{K-1/2} ( U_{K-1} - U_K   )
        - \frac{C_{p} \hat{\kappa}_K T_{v,K}'}{a} 
            \DP{\pi}{\lambda} 
        + {\cal F}_{\lambda, K} \cos \varphi, 
 \end{split}
\end{align}
%
\begin{align}
 \begin{split}
  {V_A}_{,1}
    & = - ( \zeta_1 + f ) U_1 
        - \frac{1}{2 \Delta \sigma_1} 
                 \dot{\sigma}_{3/2} ( V_1 - V_2 ) 
        - \frac{C_{p} \hat{\kappa}_1 T_{v,1}'}{a} 
             (1-\mu^2) \DP{\pi}{\mu} 
        + {\cal F}_{\varphi, 1} \cos \varphi, \\
  {V_A}_{,k}
    & = - ( \zeta_k + f ) U_k 
        - \frac{1}{2 \Delta \sigma_k} 
             [   \dot{\sigma}_{k-1/2} ( V_{k-1} - V_k   )
               + \dot{\sigma}_{k+1/2} ( V_k   - V_{k+1} ) ] \\
    & \quad
        - \frac{C_{p} \hat{\kappa}_k T_{v,k}'}{a} 
             (1-\mu^2) \DP{\pi}{\mu} 
        + {\cal F}_{\varphi, k} \cos \varphi,  \qquad (k=2, \cdots, K-1) \\
  {V_A}_{,K}
    & = - ( \zeta_K + f ) U_K 
        - \frac{1}{2 \Delta \sigma_K} 
                 \dot{\sigma}_{K-1/2} ( V_{K-1} - V_K   )  \\
    & \quad
        - \frac{C_{p} \hat{\kappa}_K T_{v,K}'}{a} 
             (1-\mu^2) \DP{\pi}{\mu} 
        + {\cal F}_{\varphi, K} \cos \varphi, 
 \end{split}
\end{align}
%
\begin{align}
 \begin{split}
  \Deqlab{鉛直差分のκの定義}
   \hat{\kappa}_k 
     & = \frac{  \sigma_{k-1/2}(   \sigma^{\kappa}_{k-1/2} 
                                 - \sigma^{\kappa}_k      ) 
               + \sigma_{k+1/2}(   \sigma^{\kappa}_k 
                                 - \sigma^{\kappa}_{k+1/2}  ) }
              { \sigma^{\kappa}_k ( \sigma_{k-1/2} - \sigma_{k+1/2} ) } \\
     & = \frac{ \sigma_{k-1/2} \alpha_k + \sigma_{k+1/2} \beta_k }
              { \Delta \sigma_k                                  },  
 \end{split} \\
 %
 T_{v,k}' & = T_{v,k} - \overline{T}_k,  \\
 %
 (\mbox{\sl KE})_k &= \frac{U^{2}_k + V^{2}_k}{2 (1-\mu^2)}.
\end{align}

\subsection{熱力学の式}
\vspace*{-5ex}
\begin{align}
 \begin{split}
  \DP{T_k}{t}
    & = - \Dinv{a \cos \varphi}
            \left(    \Dinv{1 - \mu^2} \DP{U_k T^{\prime}_k}{\lambda}
                    + \DP{V_k T^{\prime}_k}{\mu}
            \right)
        + H_k  \\
    & \quad
        + \frac{Q_k}{C_{p}}
        + {\cal D}(T_k) 
        + {\cal D}'(\Dvect{v}).
 \end{split}
\end{align}
%
ここで,
\begin{align}
 \begin{split}
   H_k 
    & \equiv
           T_k' D_k
         - \frac{1}{\Delta \sigma_k} 
             [   \dot{\sigma}_{k-1/2} ( \hat{T}_{k-1/2} - T_k   )
               + \dot{\sigma}_{k+1/2} ( T_k   - \hat{T}_{k+1/2} ) ]
               \\
    & \quad
         + \left\{ \alpha_k
                    \left[   \sigma_{k-1/2} \Dvect{v}_k \cdot \nabla \pi
                           - \sum_{l=k}^{K} 
                               ( D_l + \Dvect{v}_l \cdot \nabla \pi )
                                 \Delta  \sigma_l
                    \right]
           \right.
               \\
    & \quad
         + \quad
           \left. \beta_k
                   \left[   \sigma_{k+1/2} \Dvect{v}_k \cdot \nabla \pi
                          - \sum_{l=k+1}^{K} 
                              ( D_l + \Dvect{v}_l \cdot \nabla \pi )
                                \Delta  \sigma_l
                   \right]
           \right\} 
             \frac{1}{\Delta \sigma_k} T_{v,k} \\
%
    & =    T_k' D_k 
         - \frac{1}{\Delta \sigma_k} 
             [   \dot{\sigma}_{k-1/2} ( \hat{T}_{k-1/2} - T_k   )
               + \dot{\sigma}_{k+1/2} ( T_k   - \hat{T}_{k+1/2} ) ]
               \\
    & \quad
         + \hat{\kappa}_k \Dvect{v}_k \cdot \nabla \pi T_{v,k} 
               \\
    & \quad
         - \alpha_k \sum_{l=k}^{K} 
                      ( D_l + \Dvect{v}_l \cdot \nabla \pi )
                        \Delta  \sigma_l 
                          \frac{T_{v,k}}{\Delta \sigma_k} 
               \\
    & \quad
         - \beta_k \sum_{l=k+1}^{K} 
                     ( D_l + \Dvect{v}_l \cdot \nabla \pi )
                       \Delta  \sigma_l 
                         \frac{T_{v,k}}{\Delta \sigma_k}
      \qquad (k=1, \cdots, K-1), \\
%
   H_K
    & \equiv
           T_K' D_K 
         - \frac{1}{\Delta \sigma_K} 
             [   \dot{\sigma}_{K-1/2} ( \hat{T}_{K-1/2} - T_K   )
               + \dot{\sigma}_{K+1/2} ( T_K   - \hat{T}_{K+1/2} ) ]
               \\
    & \quad
         + \hat{\kappa}_K \Dvect{v}_K \cdot \nabla \pi T_{v,K} 
               \\
    & \quad
         - \alpha_K 
                      ( D_K + \Dvect{v}_K \cdot \nabla \pi )
                        \Delta  \sigma_K 
                          \frac{T_{v,K}}{\Delta \sigma_K} 
 \end{split}
\end{align}
%
であり, 
%
\begin{align} 
 \begin{split}
  \hat{T}_{k-1/2}
    & = \frac{
          \left[
            \left( \displaystyle
                \frac{ \sigma_{k-1/2} }{ \sigma_k } \right)^{\kappa}
              - 1
            \right] \sigma_{k-1}^{\kappa} T_k
        + \left[
              1
            - \left( \displaystyle
                  \frac{ \sigma_{k-1/2} }{ \sigma_{k-1} }
              \right)^{\kappa}
          \right] \sigma_k^{\kappa} T_{k-1}                   }
        { \sigma_{k-1}^{\kappa} - \sigma_k^{\kappa}           }
        \\
    & = a_k T_k + b_{k-1} T_{k-1}
        \qquad  \quad (k = 2, \cdots, K), \\
  \hat{T}_{1/2} & = 0, \\
  \hat{T}_{K + 1/2} & = 0, 
 \end{split}
\end{align}
%
\begin{align}
  \Deqlab{温度補間係数a}
  a_k & = \alpha_k 
            \left[ 1- \left( \frac{ \sigma_k }{ \sigma_{k-1} }
                        \right)^{\kappa} \right]^{-1},    \\
  \Deqlab{温度補間係数b}
  b_k & = \beta_k 
            \left[ \left( \frac{ \sigma_k }{ \sigma_{k+1} } 
                     \right)^{\kappa} - 1 \right]^{-1} .  
\end{align}

\subsection{水蒸気の式}
\vspace*{-5ex}
\begin{align}
 \Deqlab{水蒸気の式の鉛直離散化}
  \DP{q_k}{t}
    & =   - \Dinv{a}
              \left(    \Dinv{1 - \mu^2} \DP{U_k q_k}{\lambda}
                      + \DP{V_k q_k}{\mu}
              \right)
          + R_k 
          + S_{q,k}
          + {\cal D}(q_k). 
\end{align}
%
ここで, 
%
\begin{align}
 \Deqlab{水蒸気の変化Ｒ}
 \begin{split}
  R_1 &=    q_1 D_1 
          - \frac{1}{2 \Delta \sigma_1} 
                       \dot{\sigma}_{  3/2} ( q_1     - q_{2} )          , \\
  R_k &=    q_k D_k 
          - \frac{1}{2 \Delta \sigma_k} 
              \left[   \dot{\sigma}_{k-1/2} ( q_{k-1} - q_k   )
                     + \dot{\sigma}_{k+1/2} ( q_k     - q_{k+1} ) \right],
            \qquad (k=2,\cdots,K-1) \\
  R_K &=    q_K D_K 
          - \frac{1}{2 \Delta \sigma_K} 
                       \dot{\sigma}_{K-1/2} ( q_{K-1} - q_K   )          .
 \end{split}
\end{align}

\subsection{鉛直差分に関する補足}
\footnote{(2011-09-07 石渡)  この節はまだ編集中である.}

Arakawa and Suarez (1983) による鉛直差分の取り方に関して, 説明を
追加しておく.

Arakawa and Suarez (1983) の鉛直差分式は以下の特徴を持つ.
\begin{itemize}
 \item 圧力勾配項の鉛直積分を地形に沿った閉曲線に沿って線積分
       すると 0 になる
       (Arakawa and Suarez, 1983 は循環が 0 になるという書き方をしている)

       これは運動方程式の圧力勾配項に関する条件である.
       (地形無しの場合には角運動量が保存することを意味する).

 \item 熱エネルギーと運動エネルギーとの変換項の有限差分近似は
       運動エネルギーの式でも熱エネルギーの式でも
       同じ形をとる.

 \item $\theta$ の全球質量積分(?) が断熱過程で保存されることが保証されている
       (保存, が指す具体的な内容は? flux form に準拠した差分式に
        なっている, で正しい?)

 \item 静水圧の式から得られる最下層の層厚は local な形をしている
 \item 静水圧の式の差分形は鉛直方向にエントロピー一様な大気において
       正確な形をしている.
 \item $p_{top}=0$, 3 次元等エントロピー大気の場合
       ($T= a + b (p/p_0)^{\kappa}$ の場合?)
       には, 圧力勾配項は正確な形になっている.
 \item しかし, $\theta$ の関数の global mass integral は断熱過程
       で保存されない.
\end{itemize}

%彼らによって作られた差分式は以下の通りである.
%\begin{itemize}
%\item 連続の式
%\begin{equation}
% \DP{}{t} \tilde{p_s} + \Ddiv (\tilde{p_s} \Dvect{v}_k) 
% + \left[ \frac{(\delta \tilde{p_s} \dot{\sigma})_{k}}{(\delta \sigma)_{k}}\right]
% = 0
%\end{equation}
%
%\item 熱力学方程式
%\begin{eqnarray}
%   \DP{}{t} (\tilde{p_s} \theta_{k})
% + \nabla_{\sigma} (\tilde{p_s} \Dvect{v}_{k} \theta_{k})
% + \left[ \frac{\delta(\tilde{p_s} \dot{\sigma} \theta)}{\delta \sigma} 
%   \right]_{k}
%= 0
%\ \ \ \ \ \ \ AS1983 の (8.1) \Deqlab{(8.1)}
%\end{eqnarray}
%$\theta_{k-1/2}$ の計算には \Deqref{(8.5)} を使う.
%
%\item $\theta_{k-1/2}$ の補間式
%\begin{eqnarray}
%        \theta_{k-1/2} 
%& = &  
%     \frac{
%              ( \Pi_{k-1} - \Pi_{k-1/2} ) \theta_{k-1}
%            + ( \Pi_{k-1/2} - \Pi_{k}  ) \theta_{k}
%         }{
%            \Pi_{k-1} - \Pi_{k}
%          } \ \ \ \ \ AS1983 の (8.5) \Deqlab{(8.5)}
%\end{eqnarray}
%
%\item 圧力勾配項
%\begin{equation}
%    - (\nabla_{p} \Phi)_{k}
% = - \nabla_{\sigma} \Phi_{k}
%   - c_p \theta_{k} \nabla_{\sigma} \Pi_{k}
%   \ \ \ \ \ AS1983 の  (8.2) \Deqlab{(8.2)}
%\end{equation}
%
%\item 静水圧の式
%\begin{eqnarray}
% \Phi_{1} - \Phi_{1/2} & = & c_p \theta_{1} (\Pi_{1/2} - \Pi_{ 1 })
%  \ \ \ \ \  AS1983 の (8.3) \Deqlab{(8.3)}, \\
%   \Phi_{k} - \Phi_{k-1/2} 
%& = & c_p 
%   \left[     (\Pi_{k-1/2} - \Pi_{k}) \theta_{k} 
%            + (\Pi_{k-1} - \Pi_{k-1/2}) \theta_{k-1} 
%   \right] \nonumber \\
%& = & c_p \theta_{k-1/2} (\Pi_{k-1} - \Pi_{k}),     \nonumber \\
%& &  \ \ \  (2 \le k \le K),  \ \ \ \ AS1983 の (8.4) \Deqlab{(8.4)}
%\end{eqnarray}
%ただし, 最後の変形では \Deqref{(8.5)} を使った.
%
%
%\item 半整数グリッドの $\Pi$
%\begin{eqnarray}
% \Pi_{k-1/2} = \left( \frac{p_{k-1/2}}{p_0} \right)^{\kappa}
%\end{eqnarray}
%
%\item 整数グリッドの $\Pi$
%\Deqref{(6.28)}:
%\begin{eqnarray}
%         \Pi_{k} 
%& = & \frac{1}{\kappa+1} 
%      \left[   \frac{   p_{k-1/2} \Pi_{k-1/2}
%                      - p_{k+1/2} \Pi_{k+1/2}
%                   }{p_{k-1/2} - p_{k+1/2} }
%      \right],
% \ \ \ \  \ AS1983 の (8.6) \Deqlab{(8.6)}
%\end{eqnarray}
%
%\item $\DD{\Pi}{\tilde{p_s}}$
%
%\Deqref{(6.4)} より, あるいは \Deqref{(8.6)} を微分して\footnotemark
%\begin{eqnarray}
%    \DD{\Pi_{k}} {\tilde{p_s}}
%  = \frac{   \sigma_{k-1/2} (\Pi_{k-1/2} - \Pi_{k})
%           - \sigma_{k+1/2} (\Pi_{k} - \Pi_{k+1/2})
%                   }{p_{k-1/2} - p_{k+1/2} }
%  \ \ \ \ \ AS1983 の (8.7) \Deqlab{(8.7)}
%\end{eqnarray}
%\footnotetext{
%この式, 未確認.
%}
%\end{itemize}


\subsection{鉛直レベルの決め方}
\footnote{(2011-02-22 石渡)  この節はまだ編集中である.}

AGCM5 では, L16 の設定として以下のものが使われてきた.

\begin{table}[hbtp]
 \begin{center}
 \begin{tabular}{cc}                \hline \hline
   $k$   & $\sigma_{k+1/2}$    \\ \hline
    0    &   1       \\
    1    &   0.990   \\
    2    &   0.970   \\
    3    &   0.930   \\
    4    &   0.870   \\
    5    &   0.790   \\
    6    &   0.700   \\
    7    &   0.600   \\
    8    &   0.500   \\
    9    &   0.410   \\
    10   &   0.330   \\
    11   &   0.260   \\
    12   &   0.200   \\
    13   &   0.150   \\
    14   &   0.100   \\
    15   &   0.050   \\
    16   &   0.000   \\
    \hline
   \end{tabular}
   \end{center}
   \caption{AGCM5 の L16 計算で習慣的に使われてきた鉛直格子点の位置. 
            $\sigma_{k+1/2}$ は半整数グリッドにおける $\sigma$ 座標値.}
  \Dtablab{vlevel_AGCM5_L16}
\end{table}


