% 表題   DCPAM5  惑星表面と地下の熱収支
%
% 履歴
%\Drireki{1993/03/18 保坂征宏}
%\Drireki{2010/04/14 高橋芳幸}
%
%  \Dchapterhead
\section{離散表現}
\label{sec:disc:heatbudget}

ここでは, 惑星表面・地下の熱収支の離散化について述べる. 


\subsection{惑星表面 1 層モデル}

惑星表面に 1 層の板があるモデルにおいて, 
熱容量が有限の場合の熱収支の式 \Deqref{Surf1LayModelHeatBudget} を, 
\Dchapref{乱流過程} に示した惑星表面におけるフラックスを用いて整理すると, 
%
\begin{eqnarray}
  & & \hspace{-1cm}
  - L \epsilon (TC)_{q,\frac{1}{2}} \left ( q_1^{t+\Delta t} - q_1^{t-\Delta t} \right) \nonumber
\\
  & &
  + \left ( 
        \frac{C_s}{2 \Delta t} 
      + C_p (TC)_{h,\frac{1}{2}} 
      + \frac{\partial F_{LR}}{\partial T_s} 
      + L \epsilon (TC)_{q,\frac{1}{2}} \frac{\partial q_s^*}{\partial T} 
    \right) 
    \left ( T_s^{t+\Delta t} - T_s^{t-\Delta t} \right) \nonumber
\\
  & &
  + \left ( 
      - C_p \frac{P_{\frac{1}{2}}}{P_{1}} (TC)_{h,\frac{1}{2}} 
      + \frac{\partial F_{LR}}{\partial T_1} 
    \right) 
    \left ( T_1^{t+\Delta t} - T_1^{t-\Delta t} \right) \nonumber
\\
    &=& - F_{SR}^{t+\Delta t}
        - F_{LR}\left (T_s^{t-\Delta t}, T_1^{t-\Delta t}\right) 
        - F_{h,\frac{1}{2}}^{t-\Delta t}
        - L F_{q,\frac{1}{2}}^{t-\Delta t}
        + F_g^{t+\Delta t}
\end{eqnarray}
%
となる. 
なお, この変形においては, 惑星表面におけるフラックスとして $t+\Delta t$ の時刻の
値を用いている. 
もし, 惑星表面における水蒸気フラックスとして $t-\Delta t$ の値を用いる場合には, 
左辺第一項がなくなり, また 
$\displaystyle L \epsilon (TC)_{q,\frac{1}{2}} \frac{\partial q_s^*}{\partial T}$ 
を削除すれば良い.
これらを, 今後の式の整理の都合から, $k = 0$ として下の様に書き直す. 
%
\begin{eqnarray}
  b_{s,k,k-1} \left ( q_1^{t+\Delta t} - q_1^{t-\Delta t} \right)
  + b_{s,k,k} 
    \left ( T_s^{t+\Delta t} - T_s^{t-\Delta t} \right) \nonumber
\\
  + b_{s,k,k+1} 
    \left ( T_1^{t+\Delta t} - T_1^{t-\Delta t} \right)
    = g_{s,k}
    \label{eq:disc:surface1layermodelheatbudget}
\end{eqnarray}
%
ここで, 
%
\begin{eqnarray}
  b_{s,k,k-1} &=& - L \epsilon (TC)_{q,\frac{1}{2}}
\\
  b_{s,k,k  } &=& \frac{C_s}{2 \Delta t} 
                  + C_p (TC)_{h,\frac{1}{2}} 
                  + L \epsilon (TC)_{q,\frac{1}{2}} \frac{\partial q_s^*}{\partial T_s} 
                  + \frac{\partial F_{LR}}{\partial T_s}
\\
  b_{s,k,k+1} &=& - C_p \frac{P_{\frac{1}{2}}}{P_{1}} (TC)_{h,\frac{1}{2}}
                + \frac{\partial F_{LR}}{\partial T_1} 
\\
  g_{s,k  } 
    &=& - F_{SR}^{t+\Delta t}
        - F_{LR}\left (T_s^{t-\Delta t}, T_1^{t-\Delta t}\right) 
        - F_{h,\frac{1}{2}}^{t-\Delta t}
        - L F_{q,\frac{1}{2}}^{t-\Delta t}
        + F_g^{t+\Delta t}
\end{eqnarray}
%
である. 

一方, 層の熱容量が無限大 (もしくは惑星表面温度が固定) の場合には, 同様にして
%
\begin{eqnarray}
  b_{s,k,k-1} &=& 0
\\
  b_{s,k,k  } &=& 1
\\
  b_{s,k,k+1} &=& 0
\\
  g_{s,k  } &=& 0.
\end{eqnarray}
%
となる. 


%
%for $k = 0$
%
%\begin{eqnarray}
%  b_{k,k} &=& \frac{ C_s }{ 2 \Delta t }
%\\
%  b_{k,k} &=& 1 ( for ocean )
%\end{eqnarray}
%



\subsection{地表面における熱収支と地下における熱伝導方程式}

土壌の熱伝導方程式は下のように離散化される. 
%
\begin{eqnarray}
  C_{g} \frac{ T_{g,k}^{t+\Delta t} - T_{g,k}^{t-\Delta t} }{ 2 \Delta t }
    &=& - \frac{ F_{g,h,k+\frac{1}{2}}^{t+\Delta t} - F_{g,h,k-\frac{1}{2}}^{t+\Delta t} }{ z_{k+\frac{1}{2}} - z_{k-\frac{1}{2}} }
  \Deqlab{energybudget:disc:heatdiffeq}
\end{eqnarray}
%
ここで, $1 \le k \le k_{s,max}-1$ のとき, 
%
\begin{eqnarray}
  F_{g,h,k+\frac{1}{2}} 
    &=& - (TC)_{g,k+\frac{1}{2}} \left ( T_{g,k+1} - T_{g,k} \right), 
\\
  (TC)_{g,k+\frac{1}{2}} &=& \kappa_{g,k+\frac{1}{2}} \frac{1}{ z_{k+1} - z_k }
\end{eqnarray}
%
であり, 上部境界条件は ($k = 1$ のとき), 
%
\begin{eqnarray}
  F_{g,h,k-\frac{1}{2}} &=& F_{SR} + F_{LR} + F_{h,\frac{1}{2}} + L F_{q,\frac{1}{2}} 
  \Deqlab{disc:surface energy budget with subsurface diffusion}
\end{eqnarray}
%
であり, 下部境界条件は ($k = k_{s,max}$ のとき),
%
\begin{eqnarray}
  F_{g,h,k+\frac{1}{2}} &=& 0
\end{eqnarray}
%
である. 


しかし, このままでは方程式の数よりも未知数 (大気温度 $T$ ($k_{max}$ 個), 
地表面温度 $T_s$ (1 個), 土壌温度 $T_{g}$ ($k_{s,max}$ 個)) の数の方が
多いために解けない. 
そこで以下の式を導入する. 
%
\begin{eqnarray}
  F_{g,h,\frac{1}{2}} &=& - (TC)_{g,\frac{1}{2}} \left ( T_{g,1} - T_{s} \right)
  \label{eq:diff_soil_heat_diff_ubc}
\\
  (TC)_{g,\frac{1}{2}} &=& \kappa_{g,\frac{1}{2}} \frac{1}{ z_1 - 0 }
\end{eqnarray}
%
今後, (\ref{eq:diff_soil_heat_diff_ubc}) を上部境界条件と考え, 
同時に, \Deqref{disc:surface energy budget with subsurface diffusion} を 
$k = 0$ における式と考えることで, 大気と土壌の熱収支を仲介させる
\footnote{
  ここでは, 
  $-k_{max} \le k \le -1$ が大気中の層のインデクスであり, 
  $k = 0$ が, 言わば, 地表面のインデクスであり, 
  $1 \le k \le k_{s,max}$ が土壌中のインデクスとなる.
}. 


土壌の熱拡散方程式を変形して整理すると, 
$2 \le k \le k_{s,max}-1$ のとき, 
%
\begin{eqnarray}
  & & \hspace{-1cm}
  - (TC)_{g,k-\frac{1}{2}} 
      \left ( T_{g,k-1}^{t+\Delta t} - T_{g,k-1}^{t-\Delta t} \right) \nonumber
\\
  & & + \left\{ 
          \frac{ 1 }{ 2 \Delta t } 
            C_{g,k}, \left ( z_{k+\frac{1}{2}} - z_{k-\frac{1}{2}} \right)
          + (TC)_{g,k-\frac{1}{2}}
          + (TC)_{g,k+\frac{1}{2}}
        \right\}
        \left ( T_{g,k  }^{t+\Delta t} - T_{g,k  }^{t-\Delta t} \right) \nonumber
\\
   & & - (TC)_{g,k+\frac{1}{2}} 
           \left ( T_{g,k+1}^{t+\Delta t} - T_{g,k+1}^{t-\Delta t} \right) \nonumber
\\
   &=& - \left ( F_{g,h,k+\frac{1}{2}}^{t-\Delta t} - F_{g,h,k-\frac{1}{2}}^{t-\Delta t} \right)
\end{eqnarray}
%
となり, 
%
$k = k_{s,max}$ のとき, 
%
\begin{eqnarray}
  & & \hspace{-1cm}
  - (TC)_{g,k-\frac{1}{2}} \left ( T_{g,k-1}^{t+\Delta t} - T_{g,k-1}^{t-\Delta t} \right) \nonumber
\\
  & & + \left\{ 
          \frac{ 1 }{ 2 \Delta t } 
            C_{g,k}, \left ( z_{k+\frac{1}{2}} - z_{k-\frac{1}{2}} \right)
          + (TC)_{g,k-\frac{1}{2}}
        \right\}
        \left ( T_{g,k  }^{t+\Delta t} - T_{g,k  }^{t-\Delta t} \right) \nonumber
\\
   &=& - \left ( F_{g,h,k+\frac{1}{2}}^{t-\Delta t} - F_{g,h,k-\frac{1}{2}}^{t-\Delta t} \right)
\end{eqnarray}
%
となる. 

$k = 1$ における式は, \Deqref{energybudget:disc:heatdiffeq} を
(\ref{eq:diff_soil_heat_diff_ubc}) の式を用いて変形することで得られる. 
得られる式は, $k = 2$ の式において, 
%
\begin{eqnarray}
  T_{g,k-1} &=& T_s
\end{eqnarray}
%
とした式と同じである. 


$k = 0$ のとき, 
この式は, 式の形としては, 地表面に熱容量ゼロの仮想的な層が存在すると
仮定することと等価である. 
そこで, ここではこの考えを拡張し, 一様な温度 $T_s$ を持ち, 単位面積当たりの熱容量が
$C_{s}$ である層が地表面直下にあると考えることにする. 
この層の熱収支の式は \Deqref{disc:surface energy budget with subsurface diffusion} を
拡張し, 以下のように書くことができる
%
\footnote{
  $C_s = 0$ の場合に, \Deqref{disc:surface energy budget with subsurface diffusion} と
  等しくなることは容易に確認できる.
  このように定式化しておくと, slub ocean の条件に適応できる. 
  例えば, $C_s \ne 0$, $F_{g,h,k-\frac{1}{2}} = 0, (TC)_{g,\frac{1}{2}} = 0$ の場合には, 
  slub ocean に対応する. 
  しかし, これは単に計算上 / モデル開発上の工夫であるが, 実際にどの程度役に立つかは
  未知数. 
}. 
%
\begin{eqnarray}
  C_s \frac{\partial T_s}{\partial t} 
    &=& - F_{SR} - F_{LR} - F_{h,\frac{1}{2}} - L F_{q,\frac{1}{2}} + F_{g,h,\frac{1}{2}} 
\end{eqnarray}

この式を時間に関して離散化すると, 
%
\begin{eqnarray}
  \Deqlab{surface energy budget:disc:surface energy equation}
  C_s \frac{T_s^{t+\Delta t} - T_s^{t-\Delta t}}{ 2 \Delta t } 
    &=& - F_{SR}^{t+\Delta t} - F_{LR}^{t+\Delta t} - F_{h,\frac{1}{2}}^{t+\Delta t} - L F_{q,\frac{1}{2}}^{t+\Delta t/*} + F_{g,h,\frac{1}{2}}^{t+\Delta t/*}
\end{eqnarray}
%
となる.
$L F_{q,\frac{1}{2}}^{t+\Delta t/*}$ は, 
$L F_{q,\frac{1}{2}}^{t+\Delta t}$ か,
$L F_{q,\frac{1}{2}}^{*}$ のどちらかの値を表し,
%
\begin{align}
  L F_{q,\frac{1}{2}}^{t+\Delta t} 
  &=   L F_{q,\frac{1}{2}}^{t-\Delta t} 
     - L \epsilon (TC)_{q,\frac{1}{2}} ( q_1^{t+\Delta t} - q_1^{t-\Delta t} )
     + L \epsilon (TC)_{q,\frac{1}{2}} \frac{\partial q_s^*}{\partial T}
       ( T_s^{t+\Delta t} - T_s^{t-\Delta t} )
\\
  L F_{q,\frac{1}{2}}^{*}
  &=   L F_{q,\frac{1}{2}}^{t-\Delta t} 
     - L \epsilon (TC)_{q,\frac{1}{2}} ( q_1^{\dagger} - q_1^{t-\Delta t} )
     + L \epsilon (TC)_{q,\frac{1}{2}} \frac{\partial q_s^*}{\partial T}
       ( T_s^{t+\Delta t} - T_s^{t-\Delta t} )
\\
  &=   L F_{q,\frac{1}{2}}^{\dagger} 
     + L \epsilon (TC)_{q,\frac{1}{2}} \frac{\partial q_s^*}{\partial T}
       ( T_s^{t+\Delta t} - T_s^{t-\Delta t} )
\end{align}
%
である. ここで, $q_1^\dagger$ は仮の値である. 
同様に, 
$F_{g,h,\frac{1}{2}}^{t+\Delta t/*}$ は
$F_{g,h,\frac{1}{2}}^{t+\Delta t}$ か
$F_{g,h,\frac{1}{2}}^{*}$ のどちらかを表し,
%
\begin{align}
  F_{g,h,\frac{1}{2}}^{t+\Delta t}
  &=   F_{g,h,\frac{1}{2}}^{t-\Delta t}
     - (TC)_{g,\frac{1}{2}} ( T_{g,1}^{t+\Delta t} - T_{g,1}^{t-\Delta t} )
     + (TC)_{g,\frac{1}{2}} ( T_s^{t+\Delta t} - T_s^{t-\Delta t} )
\\
  F_{g,h,\frac{1}{2}}^{*}
  &=   F_{g,h,\frac{1}{2}}^{t-\Delta t}
     - (TC)_{g,\frac{1}{2}} ( T_{g,1}^{\dagger} - T_{g,1}^{t-\Delta t} )
     + (TC)_{g,\frac{1}{2}} ( T_s^{t+\Delta t} - T_s^{t-\Delta t} )
\\
  &=   F_{g,h,\frac{1}{2}}^{\dagger}
     + (TC)_{g,\frac{1}{2}} ( T_s^{t+\Delta t} - T_s^{t-\Delta t} )
\end{align}
%
である. ここで, $T_{g,1}^\dagger$ は仮の値である. 

なお, 
$L F_{q,\frac{1}{2}}^{t+\Delta t}$, $F_{g,h,\frac{1}{2}}^{t+\Delta t}$ を
用いて離散化して整理すると, 
%この式を整理すると, 
%#####
\begin{eqnarray}
  \Deqlab{surface energy budget:disc:surface energy equation rearranged}
  & & \hspace{-1cm}
%
  - L \epsilon (TC)_{q,\frac{1}{2}} ( q_1^{t+\Delta t} - q_1^{t-\Delta t} )
\\
  & &
  + \left ( 
      - C_p \frac{P_{\frac{1}{2}}}{P_{1}} (TC)_{h,\frac{1}{2}} 
      + \frac{\partial F_{LR}}{\partial T_1} 
    \right) 
    \left ( T_1^{t+\Delta t} - T_1^{t-\Delta t} \right)  \nonumber
\\
  & &
  + \left ( 
        \frac{C_s}{ 2 \Delta t } 
      - (TC)_{g,\frac{1}{2}}
%        \frac{C_{g}}{2 \Delta t} 
      + C_p (TC)_{h,\frac{1}{2}} 
      + \frac{\partial F_{LR}}{\partial T_s} 
      + L \epsilon (TC)_{q,\frac{1}{2}} \frac{\partial q_s^*}{\partial T} 
    \right) 
    \left ( T_s^{t+\Delta t} - T_s^{t-\Delta t} \right)  \nonumber
\\
  & &
  + (TC)_{g,\frac{1}{2}} \left ( T_{g,1}^{t+\Delta t} - T_{g,1}^{t-\Delta t} \right)  \nonumber
\\
    &=& - F_{SR}^{t+\Delta t}
        - F_{LR}\left (T_s^{t-\Delta t}, T_1^{t-\Delta t}\right) 
        - F_{h,\frac{1}{2}}^{t-\Delta t}
        - L F_{q,\frac{1}{2}}^{t-\Delta t}
        + F_{g,h,\frac{1}{2}}^{t-\Delta t} 
\end{eqnarray}
%
となる. 
\footnote{
  メモ:
  %
\begin{eqnarray}
  F_{q,k-\frac{1}{2}} &=& - \epsilon (TC)_{q,k-\frac{1}{2}} ( q_k - q_s^* )
\\
  F_{q,k-\frac{1}{2}}^{n+1}
    &=& - \epsilon (TC)_{q,k-\frac{1}{2}} ( q_k^{n+1} - (q_s^*)^{n+1} )
\\
    &=& - \epsilon (TC)_{q,k-\frac{1}{2}} \left[ q_k^{n+1} - \left\{ (q_s^*)^{n-1} + \left( \frac{\partial q^*}{\partial T} \right) ( T_s^{n+1} - T_s^{n-1} ) \right\} \right]
\\
    &=& - \epsilon (TC)_{q,k-\frac{1}{2}} \left\{ q_k^{n+1} - (q_s^*)^{n-1} \right\}
        + \epsilon (TC)_{q,k-\frac{1}{2}} \left( \frac{\partial q^*}{\partial T} \right) ( T_s^{n+1} - T_s^{n-1} )
\\
    &=& - \epsilon (TC)_{q,k-\frac{1}{2}} \left\{ q_k^{n+1} - q_k^{n-1} + q_k^{n-1} - (q_s^*)^{n-1} \right\}
        + \epsilon (TC)_{q,k-\frac{1}{2}} \left( \frac{\partial q^*}{\partial T} \right) ( T_s^{n+1} - T_s^{n-1} )
\\
    &=& - \epsilon (TC)_{q,k-\frac{1}{2}} \left\{ q_k^{n+1} - q_k^{n-1} \right\}
        - \epsilon (TC)_{q,k-\frac{1}{2}} \left\{ q_k^{n-1} - (q_s^*)^{n-1} \right\}
        + \epsilon (TC)_{q,k-\frac{1}{2}} \left( \frac{\partial q^*}{\partial T} \right) ( T_s^{n+1} - T_s^{n-1} )
\\
    &=& - \epsilon (TC)_{q,k-\frac{1}{2}} \left\{ q_k^{n+1} - q_k^{n-1} \right\}
        + F_{q,k-\frac{1}{2}}^{n+1}
        + \epsilon (TC)_{q,k-\frac{1}{2}} \left( \frac{\partial q^*}{\partial T} \right) ( T_s^{n+1} - T_s^{n-1} )
\end{eqnarray}
}
%
しかし, この式では, 他の式と連立させて解くときに, 三重対角行列に収まらない. 
そこで, 仮の値, $q_1^\dagger$, $T_{g,1}^\dagger$ を用いて, 連立一次方程式
を分割して計算する. 


$L F_{q,\frac{1}{2}}^{t+\Delta t}$, $F_{g,h,\frac{1}{2}}^{*}$
を用いて離散化する時, 整理すると, 
%また, 土壌熱伝導フラックスに $t-\Delta t$ の時刻の値を用いる場合には, 
%#####
\begin{eqnarray}
%  \Deqlab{surface energy budget:disc:surface energy equation rearranged}
  & & \hspace{-1cm}
%
  - L \epsilon (TC)_{q,\frac{1}{2}} ( q_1^{t+\Delta t} - q_1^{t-\Delta t} )
\\
  & &
  + \left ( 
      - C_p \frac{P_{\frac{1}{2}}}{P_{1}} (TC)_{h,\frac{1}{2}} 
      + \frac{\partial F_{LR}}{\partial T_1} 
    \right) 
    \left ( T_1^{t+\Delta t} - T_1^{t-\Delta t} \right)  \nonumber
\\
  & &
  + \left ( 
        \frac{C_s}{ 2 \Delta t } 
%        \frac{C_{g}}{2 \Delta t} 
      + C_p (TC)_{h,\frac{1}{2}} 
      + \frac{\partial F_{LR}}{\partial T_s} 
      + L \epsilon (TC)_{q,\frac{1}{2}} \frac{\partial q_s^*}{\partial T} 
      - (TC)_{g,\frac{1}{2}} 
    \right) 
    \left ( T_s^{t+\Delta t} - T_s^{t-\Delta t} \right)  \nonumber
\\
    &=& - F_{SR}^{t+\Delta t}
        - F_{LR}\left (T_s^{t-\Delta t}, T_1^{t-\Delta t}\right) 
        - F_{h,\frac{1}{2}}^{t-\Delta t}
        - L F_{q,\frac{1}{2}}^{t-\Delta t}
        + F_{g,h,\frac{1}{2}}^{\dagger} 
\end{eqnarray}
%
となる. 


また,
$L F_{q,\frac{1}{2}}^{*}$, $F_{g,h,\frac{1}{2}}^{t+\Delta t}$
を用いて離散化する時, 整理すると
\footnote{
  実際には $L F_{q,\frac{1}{2}}^{t-\Delta t}$ を用いている.
  このため, $\left ( T_s^{t+\Delta t} - T_s^{t-\Delta t} \right)$ の係数の
  $+ L \epsilon (TC)_{q,\frac{1}{2}} \frac{\partial q_s^*}{\partial T}$
  は存在せず, 右辺の $- L F_{q,\frac{1}{2}}^{\dagger}$ は
  $- L F_{q,\frac{1}{2}}^{t-\Delta t}$ となっている.
  (yot, 2016/05/19)
}, 
%また, 潜熱フラックスに $t-\Delta t$ の時刻の値を用いる場合には, 
%
\begin{eqnarray}
%  \Deqlab{surface energy budget:disc:surface energy equation rearranged}
  & & \hspace{-1cm}
%
    \left ( 
      - C_p \frac{P_{\frac{1}{2}}}{P_{1}} (TC)_{h,\frac{1}{2}} 
      + \frac{\partial F_{LR}}{\partial T_1} 
    \right) 
    \left ( T_1^{t+\Delta t} - T_1^{t-\Delta t} \right)  \nonumber
\\
  & &
  + \left ( 
        \frac{C_s}{ 2 \Delta t } 
      - (TC)_{g,\frac{1}{2}}
%        \frac{C_{g}}{2 \Delta t} 
      + C_p (TC)_{h,\frac{1}{2}} 
      + \frac{\partial F_{LR}}{\partial T_s} 
      + L \epsilon (TC)_{q,\frac{1}{2}} \frac{\partial q_s^*}{\partial T} 
    \right) 
    \left ( T_s^{t+\Delta t} - T_s^{t-\Delta t} \right)  \nonumber
\\
  & &
  + (TC)_{g,\frac{1}{2}} \left ( T_{g,1}^{t+\Delta t} - T_{g,1}^{t-\Delta t} \right)  \nonumber
\\
    &=& - F_{SR}^{t+\Delta t}
        - F_{LR}\left (T_s^{t-\Delta t}, T_1^{t-\Delta t}\right) 
        - F_{h,\frac{1}{2}}^{t-\Delta t}
        - L F_{q,\frac{1}{2}}^{\dagger}
        + F_{g,h,\frac{1}{2}}^{t-\Delta t} 
\end{eqnarray}
%
となる. 

これらをまとめると, 
%
\begin{eqnarray}
\Deqlab{energybudget-disc:soilheatdiff}
    \Dvect{B}_g \Dvect{x}_g = \Dvect{G}_g
\end{eqnarray}
%
と書くことができる. 
ここで, 
%
\begin{eqnarray}
\Deqlab{surface energy budget:disc:simul. lin. eq. explanation of array and vector 1}
  \Dvect{x}_g &=& \left ( 
      T_1^{t+\Delta t} - T_1^{t-\Delta t}, 
      T_s^{t+\Delta t} - T_s^{t-\Delta t},  \right. \nonumber
\\
      & & \left.
      T_{g,1}^{t+\Delta t} - T_{g,1}^{t-\Delta t}, 
      T_{g,2}^{t+\Delta t} - T_{g,2}^{t-\Delta t}, 
      ..., 
      T_{g,k_{s,max}}^{t+\Delta t} - T_{g,k_{s,max}}^{t-\Delta t} \right), 
\\
\Deqlab{surface energy budget:disc:simul. lin. eq. explanation of array and vector 2}
  \Dvect{G}_g &=& \left ( g_{g,0}, g_{g,1}, g_{g,2}, ..., g_{g,k_{s,max}} \right), 
\\
\Deqlab{surface energy budget:disc:simul. lin. eq. explanation of array and vector 3}
  g_{g,0  } 
    &=& - F_{SR}^{t+\Delta t}
        - F_{LR}\left (T_s^{t-\Delta t}, T_1^{t-\Delta t}\right) 
        - F_{h,\frac{1}{2}}^{t-\Delta t} \nonumber
\\
    & & - L F_{q,\frac{1}{2}}^{t-\Delta t}
        + F_{g,h,\frac{1}{2}}^{t-\Delta t} 
\\
\Deqlab{surface energy budget:disc:simul. lin. eq. explanation of array and vector 4}
  g_{g,k \ge 1} &=& - \left ( F_{g,k+\frac{1}{2}}^{t-\Delta t} - F_{g,k-\frac{1}{2}}^{t-\Delta t} \right)
\end{eqnarray}
%
ここで, 
$1 \le k \le k_{s,max}-1$ のとき, 
%
\begin{eqnarray}
  b_{g,k, k-1} &=& - (TC)_{g,k-\frac{1}{2}} 
\\
  b_{g,k, k  } 
    &=& \frac{ 1 }{ 2 \Delta t } 
          C_{g,k}, \left ( z_{k+\frac{1}{2}} - z_{k-\frac{1}{2}} \right)
      + (TC)_{g,k-\frac{1}{2}}
      + (TC)_{g,k+\frac{1}{2}}
\\
  b_{g,k, k+1} &=& - (TC)_{g,k+\frac{1}{2}} 
\end{eqnarray}
%
であり, $k = 0$ のとき, 
%
\begin{eqnarray}
  b_{g,k,k-1} &=& - C_p \frac{P_{\frac{1}{2}}}{P_{1}} (TC)_{h,\frac{1}{2}}
                + \frac{\partial F_{LR}}{\partial T_1} 
\\
\Deqlab{surface energy budget:disc:explanation of b_{g,k,k} at k = 0 for simul. lin. eq.}
  b_{g,k,k  } &=&     \frac{C_s}{ 2 \Delta t } 
                  - (TC)_{g,\frac{1}{2}}
%                    \frac{C_{g}}{2 \Delta t} 
                  + C_p (TC)_{h,\frac{1}{2}} 
%                  + L \epsilon (TC)_{q,\frac{1}{2}} \frac{\partial q_s^*}{\partial T} 
                  + \frac{\partial F_{LR}}{\partial T_s}
\\
  b_{g,k,k+1} &=& (TC)_{g,\frac{1}{2}}
\end{eqnarray}
%
であり
\footnote
{
  ここでは, $T_g$ (土壌温度), $T_s$ (地表面温度), $T$ (大気温度) の順番に書いているが, 
  2010/02/20 時点のコードでは逆の順番になっている.
}
, $k = k_{s,max}$ のとき, 
%
\begin{eqnarray}
  b_{g,k, k-1} &=& - (TC)_{g,k-\frac{1}{2}} 
\\
  b_{g,k, k  } &=& \frac{ 1 }{ 2 \Delta t } 
                   C_{g,k}, \left ( z_{k+\frac{1}{2}} - z_{k-\frac{1}{2}} \right)
              + (TC)_{g,k-\frac{1}{2}}
\end{eqnarray}
%
である. 

ただし, $\Dvect{B}_g$ は $k_{s,max}+1$ 行 $k_{s,max}+2$ 列の行列であり, この式だけでは
未知数が方程式数よりも多いために閉じない. 
方程式を閉じるために, 熱の鉛直拡散の式と同時に解く. 


なお, 層の熱容量が無限大 (もしくは惑星表面温度が固定) の場合には, 
%
\begin{eqnarray}
  b_{g,k,k-1} &=& 0
\\
  b_{g,k,k  } &=& 1
\\
  b_{g,k,k+1} &=& 0
\\
  g_{g,k  } &=& 0
\end{eqnarray}
%
である. 


\subsection{氷の融解・融雪による熱収支の修正}


氷の融解および融雪時の地表面の熱収支式,
\Deqref{surface energy budget:math:surface energy equation when snow melts}
を離散化すると,
%
\begin{eqnarray}
  F_g^{t+\Delta t} 
    &=& -\kappa \frac{ T_{s}^{t+\Delta t} - T_{g,1}^{t+\Delta t} } { z_\frac{1}{2} - z_1 } \nonumber \\
\Deqlab{surface energy budget:disc:surface energy equation when snow melts}
    &=& F_{SR}^{t+\Delta t} + F_{LR}^{t+\Delta t} + F_{h,\frac{1}{2}}^{t+\Delta t} + L F_{q,\frac{1}{2}}^{t+\Delta t} (+ F_{IM}^{t+\Delta t}) + F_{SM}^{t+\Delta t}
\end{eqnarray}
%
となる. 
さらに, 
\Deqref{surface energy budget:disc:surface energy equation} と同様の方法に基づき, 
\Deqref{surface energy budget:disc:surface energy equation when snow melts} を, 
地表面に熱容量が $C_s$ である層があると考えて離散化し直すと, 
%
\begin{eqnarray}
\Deqlab{surface energy budget:disc:surface energy equation when snow melts with Cs}
  C_s \frac{T_s^{t+\Delta t} - T_s^{t-\Delta t}}{ 2 \Delta t } 
    &=& - F_{SR}^{t+\Delta t} - F_{LR}^{t+\Delta t} - F_{h,\frac{1}{2}}^{t+\Delta t} - L F_{q,\frac{1}{2}}^{t+\Delta t} + F_{g,h,\frac{1}{2}}^{t+\Delta t} (- F_{IM}^{t+\Delta t}) - F_{SM}^{t+\Delta t}
\end{eqnarray}
%
となる. 
ここで, 
$F_{g,\frac{3}{2}}^{t+\Delta t}, F_{SR}^{t+\Delta t}, F_{LR}^{t+\Delta t}, F_{h,\frac{1}{2}}^{t+\Delta t}, L F_{q,\frac{1}{2}}^{t+\Delta t}, F_{IM}^{t+\Delta t}, F_{SM}^{t+\Delta t}$
はそれぞれ, 地下への熱伝導フラックス, 短波放射フラックス, 長波放射フラックス, 
顕熱フラックス, 潜熱フラックス, 氷の融解による熱フラックス, 融雪による熱フラックスである.
$\kappa$ は土壌の熱拡散係数である.

ここでは, まず $T_s^{t+\Delta t} = T_c$ となると仮定して, 
氷の融解・融雪による熱フラックスを未知数として求め, ... やめた. 後で書く (yot, 2011/12/28).


\subsubsection{積雪がすべて解ける場合}

このとき, 積雪量 $M_{\rm snow}$ を用いると, 
$\displaystyle F_{SM}^{t+\Delta t} = \frac{L_{fusion} M_{\rm snow}}{2 \Delta t}$ 
であるから, 
この潜熱を与えて連立一次方程式を解き直せばよい. 
ここで, $L_{fusion}$ は融解潜熱である.

\Deqref{surface energy budget:disc:surface energy equation when snow melts with Cs} 
を整理すると, 
%
\begin{eqnarray}
  & & \hspace{-1cm}
%
  (TC)_{g,\frac{1}{2}} \left ( T_{g,1}^{t+\Delta t} - T_{g,1}^{t-\Delta t} \right)  \nonumber
\\
%
%  - L \epsilon (TC)_{q,\frac{1}{2}} \left ( q_1^{t+\Delta t} - q_1^{t-\Delta t} \right) \\
  & &
  + \left ( 
        \frac{C_s}{ 2 \Delta t } 
      - (TC)_{g,\frac{1}{2}}
%        \frac{C_{g}}{2 \Delta t} 
      + C_p (TC)_{h,\frac{1}{2}} 
      + \frac{\partial F_{LR}}{\partial T_s} 
%      + L \epsilon (TC)_{q,\frac{1}{2}} \frac{\partial q_s^*}{\partial T} 
    \right) 
    \left ( T_s^{t+\Delta t} - T_s^{t-\Delta t} \right)  \nonumber
\\
  & &
  + \left ( 
      - C_p \frac{P_{\frac{1}{2}}}{P_{1}} (TC)_{h,\frac{1}{2}} 
      + \frac{\partial F_{LR}}{\partial T_1} 
    \right) 
    \left ( T_1^{t+\Delta t} - T_1^{t-\Delta t} \right)  \nonumber
\\
    &=& - F_{SR}^{t+\Delta t}
        - F_{LR}\left (T_s^{t-\Delta t}, T_1^{t-\Delta t}\right) 
        - F_{h,\frac{1}{2}}^{t-\Delta t}
        - L F_{q,\frac{1}{2}}^{t-\Delta t}
        + F_{g,h,\frac{1}{2}}^{t-\Delta t} 
       (- F_{IM}^{t+\Delta t}) - F_{SM}^{t+\Delta t}
\end{eqnarray}
%
となる. 
ただし, ここでは三重対角行列にするために, 潜熱フラックスは $t-\Delta t$ の時刻の
ものを用いる. 

これらを \Deqref{energybudget-disc:soilheatdiff} の形式にまとめると, 
\Deqref{surface energy budget:disc:simul. lin. eq. explanation of array and vector 3} 
において
%
\begin{eqnarray}
  g_{g,0  } 
    &=& - F_{SR}^{t+\Delta t}
        - F_{LR}\left (T_s^{t-\Delta t}, T_1^{t-\Delta t}\right) 
        - F_{h,\frac{1}{2}}^{t-\Delta t} \nonumber
\\
    & & - L F_{q,\frac{1}{2}}^{t-\Delta t}
        + F_{g,h,\frac{1}{2}}^{t-\Delta t} 
       (- F_{IM}^{t+\Delta t}) - F_{SM}^{t+\Delta t}
\end{eqnarray}
%
とすれば良い. 



\subsubsection{積雪がすべて解けない場合}

このとき, $T_s^{t+\Delta t} = T_{cond}$ として連立一次方程式を解き直せばよい. 

\Deqref{energybudget-disc:soilheatdiff} の形式にまとめると, 
\Deqref{surface energy budget:disc:simul. lin. eq. explanation of array and vector 3}, 
\Deqref{surface energy budget:disc:explanation of b_{g,k,k} at k = 0 for simul. lin. eq.}
において
%
\begin{eqnarray}
  b_{g,k,k-1} &=& 0
\\
  b_{g,k,k  } &=& 1
\\
  b_{g,k,k+1} &=& 0
\\
  g_{g,0  } 
    &=& T_c - T_s^{t-\Delta t}
\end{eqnarray}
%
とすれば良い. 

$F_{SM}$ は, $T_s^{t+\Delta t} = T_{cond}$ として, 惑星表面の熱収支式, 
\Deqref{surface energy budget:disc:surface energy equation when snow melts with Cs},
から求める.



\subsection{海氷面上の熱収支}

海氷面上の, 海氷に伝わる熱フラックスを
%
\begin{eqnarray}
  F_b &=& - \kappa_I \frac{ T_s - T_0 }{ h_I }
\end{eqnarray}
%
と書くことにする. 
ここで, $\kappa_I$ は海氷の熱伝導率, $h_I$ は海氷の厚さであり, $T_0$ は海氷下の
海水温である. 
このとき, 海氷面上の熱収支式
\Deqref{surface energy budget:math:energy budget of 1 layer sea ice}
は下のように離散化される. 
%
\begin{eqnarray}
  & & \hspace{-1cm}
  \frac{C_I h_I}{2 \Delta t} \left ( T_s^{t+\Delta t} - T_s^{t-\Delta t} \right) \nonumber
\\
    &=& - F_{SR}^{t+\Delta t}
        - F_{LR}^{t+\Delta t}
        - F_{h,\frac{1}{2}}^{t+\Delta t}
        - L F_{q,\frac{1}{2}}^{t+\Delta t}
        + F_b^{t+\Delta t}
\end{eqnarray}
%
これを整理すると, 
%############
\begin{eqnarray}
  & & \hspace{-1cm}
  - L \epsilon (TC)_{q,\frac{1}{2}} ( q_1^{t+\Delta t} - q_1^{t-\Delta t} )
\\
  & &
  + \left ( 
        \frac{C_I h_I}{2 \Delta t} 
      + C_p (TC)_{h,\frac{1}{2}} 
      + \frac{\partial F_{LR}}{\partial T_s} 
      + \frac{\kappa_I}{h_I} 
      + L \epsilon (TC)_{q,\frac{1}{2}} \frac{\partial q_s^*}{\partial T} 
    \right) 
    \left ( T_s^{t+\Delta t} - T_s^{t-\Delta t} \right) \nonumber
\\
  & &
  + \left ( 
      - C_p \frac{P_{\frac{1}{2}}}{P_{1}} (TC)_{h,\frac{1}{2}} 
      + \frac{\partial F_{LR}}{\partial T_1} 
    \right) 
    \left ( T_1^{t+\Delta t} - T_1^{t-\Delta t} \right) \nonumber
\\
    &=& - F_{SR}^{t+\Delta t}
        - F_{LR}\left (T_s^{t-\Delta t}, T_1^{t-\Delta t}\right) 
        - F_{h,\frac{1}{2}}^{t-\Delta t}
        - L F_{q,\frac{1}{2}}^{t-\Delta t}
        + F_{b}^{t-\Delta t}
\end{eqnarray}
%
となる. 


また, 潜熱フラックスに $t-\Delta t$ の時刻の値を用いる場合には, 
%
\begin{eqnarray}
  & & \hspace{-1cm}
  \frac{C_I h_I}{2 \Delta t} \left ( T_s^{t+\Delta t} - T_s^{t-\Delta t} \right) \nonumber
\\
    &=& - F_{SR}^{t+\Delta t}
        - F_{LR}^{t+\Delta t}
        - F_{h,\frac{1}{2}}^{t+\Delta t}
        - L F_{q,\frac{1}{2}}^{t-\Delta t}
        + F_b^{t+\Delta t}
\end{eqnarray}
%
\footnote
{
これは, 整理した結果得られる行列を三重対角行列にするためである. 
これは, 水蒸気の式において惑星表面の水蒸気フラックスの値として $t - \Delta t$ の
時刻の値を使うことにしたことに起因しており, その場合には, ここでも $t - \Delta t$ の
時刻のフラックスを使わなければ水の質量が保存されない. 
もちろん, 他のやり方はあり得るだろう. 
}. 
%
これを整理すると, 
%
\begin{eqnarray}
  & & \hspace{-1cm}
    \left ( 
        \frac{C_I h_I}{2 \Delta t} 
      + C_p (TC)_{h,\frac{1}{2}} 
      + \frac{\partial F_{LR}}{\partial T_s} 
      + \frac{\kappa_I}{h_I} 
    \right) 
    \left ( T_s^{t+\Delta t} - T_s^{t-\Delta t} \right) \nonumber
\\
  & &
  + \left ( 
      - C_p \frac{P_{\frac{1}{2}}}{P_{1}} (TC)_{h,\frac{1}{2}} 
      + \frac{\partial F_{LR}}{\partial T_1} 
    \right) 
    \left ( T_1^{t+\Delta t} - T_1^{t-\Delta t} \right) \nonumber
\\
    &=& - F_{SR}^{t+\Delta t}
        - F_{LR}\left (T_s^{t-\Delta t}, T_1^{t-\Delta t}\right) 
        - F_{h,\frac{1}{2}}^{t-\Delta t}
        - L F_{q,\frac{1}{2}}^{t-\Delta t}
        + F_{b}^{t-\Delta t}
\end{eqnarray}
%
となる. 

以上を整理し, 今後の式の整理を念頭において, $k = 0$ として下の形に書くことにする. 
%
\begin{eqnarray}
  b_{i,k,k} 
    \left ( T_s^{t+\Delta t} - T_s^{t-\Delta t} \right)
  + b_{i,k,k+1} 
    \left ( T_1^{t+\Delta t} - T_1^{t-\Delta t} \right)
  = g_{i,k}
  \label{eq:disc:seaiceheatbudget}
\end{eqnarray}
%
ここで, 
%
\begin{eqnarray}
\Deqlab{surface energy budget:disc:sea ice energy equation coefficients 1}
  b_{i,k,k  } &=& \frac{C_I h_I}{2 \Delta t} 
                  + C_p (TC)_{h,\frac{1}{2}} 
                  + \frac{\partial F_{LR}}{\partial T_s}
                  + \frac{\kappa_I}{h_I} 
\\
\Deqlab{surface energy budget:disc:sea ice energy equation coefficients 2}
  b_{i,k,k+1} &=& - C_p \frac{P_{\frac{1}{2}}}{P_{1}} (TC)_{h,\frac{1}{2}}
                + \frac{\partial F_{LR}}{\partial T_1} 
\\
\Deqlab{surface energy budget:disc:sea ice energy equation coefficients 3}
  g_{i,k  } 
    &=& - F_{SR}^{t+\Delta t}
        - F_{LR}\left (T_s^{t-\Delta t}, T_1^{t-\Delta t}\right) 
        - F_{h,\frac{1}{2}}^{t-\Delta t}
        - L F_{q,\frac{1}{2}}^{t-\Delta t}  \nonumber
\\
    & & + F_{b}^{t-\Delta t}
\end{eqnarray}
%
である. 



\subsection{海氷の融解による熱収支の修正}


海氷の融解時の熱収支式
\Deqref{surface energy budget:math:energy budget of 1 layer sea ice when ice melts}
を離散化すると,
%
\begin{eqnarray}
\Deqlab{surface energy budget:disc:surface energy equation when sea ice melts}
  \frac{C_I h_I}{2 \Delta t} \left ( T_s^{t+\Delta t} - T_s^{t-\Delta t} \right)
  &=& - F_{SR}^{t+\Delta t}
        - F_{LR}^{t+\Delta t}
        - F_{h,\frac{1}{2}}^{t+\Delta t}
        - L F_{q,\frac{1}{2}}^{t-\Delta t}
        + F_b^{t+\Delta t}
        - F_{IM}^{t+\Delta t}
\\
  T_s^{t+\Delta t} &=& T_c
\end{eqnarray}
%
となる. 
ここで, $T_c$ は凝結温度である. 
$T_c$ は既知であるから, $T_s = T_c$ として連立方程式を解き直す. 
$F_{IM}$ は, $T_s = T_c$ として, 
\Deqref{surface energy budget:disc:surface energy equation when sea ice melts}
から求める. 

なお, 連立方程式を解く際には, 
\Deqref{surface energy budget:disc:sea ice energy equation coefficients 1},
\Deqref{surface energy budget:disc:sea ice energy equation coefficients 2},
\Deqref{surface energy budget:disc:sea ice energy equation coefficients 3}
を下のように置き直して用いる. 
%
\begin{eqnarray}
  b_{i,k,k  } &=& 1
\\
  b_{i,k,k+1} &=& 0
\\
  g_{i,k  } 
    &=& T_c - T_s^{t+\Delta t}
\end{eqnarray}


%$F_{IM}$ を未知数とするため, 大気温度の連立一次方程式には, 
%\Deqref{vdiff-disc:elements of g and B_a with a given T_s} において
%$F_{SM}$ を $F_{IM}$ と置き換えたもの
%を用いる.