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

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


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

惑星表面に 1 層の板があるモデルにおいて, 
熱容量が有限の場合の熱収支の式 ... を整理すると, 
%
\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}
%
となる. 
これらを, 今後の式の整理の都合から, $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} 
                  + \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}
%
である
\footnote{
なお, バケツモデルを使う時には (表現が不正確だが), 
$b_{s,k,k-1} = 0$ とし, $b_{s,k,k}$ の 
$L \epsilon (TC)_{q,\frac{1}{2}} \frac{\partial q_s^*}{\partial T}$ を削除すれば良い.
}.

一方, 層の熱容量が無限大 (もしくは惑星表面温度が固定) の場合には, 同様にして
%
\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}} }
\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$ として), 
\footnote
{
こういう書き方しないよな.
}
%
\begin{eqnarray}
  F_{g,h,k-\frac{1}{2}} &=& F_{SR} + F_{LR} + F_{h,\frac{1}{2}} + L F_{q,\frac{1}{2}} 
\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}) を上部境界条件と考えることにし, 
... を $k = 0$ における式と考えることにする. 


土壌の熱拡散方程式を変形して整理すると, 
$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$ のときは, ... において ... の式を用いて変形すると, 
$k = 2$ とした場合に
%
\begin{eqnarray}
  T_{g,k-1} &=& T_s
\end{eqnarray}
%
とした式と同じ式となる. 


$k = 0$ のとき, 
この式は, 式の形としては, 地表面に熱容量ゼロの仮想的な層が存在すると
仮定することと等価である. 
%
\footnote
{
これはここに書くもんじゃないな.
}
%
そこで, ここではこの考えを拡張し, 一様な温度 $T_s$ を持ち, 単位面積当たりの熱容量が
$C_{s}$ である層が地表面直下にあると考えることにする. 
この層の熱収支の式は ... の式を拡張し, 以下のように書くことができる. 
%
\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}
%
$C_s = 0$ の場合には, 上に述べた ... の式を適応した場合に対応する
%
\footnote
{
このように定式化しておくと, 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{T_s^{t+\Delta t} - T_s^{t-\Delta t}}{ 2 \Delta 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}
  & & \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} 
\end{eqnarray}
%
となる. 
ただし, ここでは三重対角行列にするために, 潜熱フラックスは $t-\Delta t$ の時刻の
ものを用いる. 

これらをまとめると, 
%
\begin{eqnarray}
    \Dvect{B}_g \Dvect{x}_g = \Dvect{G}_g
    \label{eq:disc:soilheatdiff}
\end{eqnarray}
%
と書くことができる. 
ここで, 
%
\begin{eqnarray}
  \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), 
\\
  \Dvect{G}_g &=& \left( g_{g,0}, g_{g,1}, g_{g,2}, ..., g_{g,k_{s,max}} \right), 
\\
  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} 
\\
  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} 
\\
  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{海氷面上の熱収支}

海氷面上の熱収支式 ... を離散化する.

海氷面上の, 海氷に伝わる熱フラックスを
%
\begin{eqnarray}
  F_{g,\frac{1}{2}} &=& - \frac{\kappa_i}{h_i} \left( T_s - T_0 \right)
\end{eqnarray}
%
と書くことにする. 
ここで, $h_i$ は海氷の厚さであり, $T_0$ は海氷下の海水温である. 
このとき, 海氷面上の熱収支式 ... は下のように離散化される. 
%
\begin{eqnarray}
  & & \hspace{-1cm}
  \frac{C_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_{g,\frac{1}{2}}^{t+\Delta t}
\end{eqnarray}
%
ただし, 潜熱フラックスは $t-\Delta t$ の時刻の
ものを用いる
%
\footnote
{
これは, 整理した結果得られる行列を三重対角行列にするためである. 
これは, 水蒸気の式において惑星表面の水蒸気フラックスの値として $t - \Delta t$ の
時刻の値を使うことにしたことに起因しており, その場合には, ここでも $t - \Delta t$ の
時刻のフラックスを使わなければ水の質量が保存されない. 
もちろん, 他のやり方はあり得るだろう. 
}. 
%
これを整理すると, 
%
\begin{eqnarray}
  & & \hspace{-1cm}
    \left( 
        \frac{C_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_{g,\frac{1}{2}}^{t-\Delta t}
\end{eqnarray}
%
となる. 

以上を整理し, 今後の式の整理を念頭において, $k = 0$ として下の形に書くことにする. 
%
\begin{eqnarray}
  b_{i,k,k-1} 
    \left( T_s^{t+\Delta t} - T_s^{t-\Delta t} \right)
  + b_{i,k,k} 
    \left( T_1^{t+\Delta t} - T_1^{t-\Delta t} \right)
  = g_{i,k}
  \label{eq:disc:seaiceheatbudget}
\end{eqnarray}
%
ここで, 
%
\begin{eqnarray}
  b_{i,k,k  } &=& \frac{C_i}{2 \Delta t} 
                  + C_p (TC)_{h,\frac{1}{2}} 
                  + \frac{\partial F_{LR}}{\partial T_s}
                  + \frac{\kappa_i}{h_i} 
\\
  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} 
\\
  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_{g,\frac{1}{2}}^{t-\Delta t}
\end{eqnarray}
%
である. 


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

  ここはまだ整理していない.


融雪時の地表面の熱収支式 ... を離散化すると, 
%
\begin{eqnarray}
  F_g^{n+1} 
    &=& -\kappa \frac{ T_{s}^{n+1} - T_{g,1}^{n+1} } { z_\frac{1}{2} - z_1 } \\
    &=& F_s^{n+1} + F_L^{n+1} + F_c^{n+1} + L F_q^{n+1} + F_{SM}^{n+1}
\\
  C_{g,1} \frac{ T_{g,1}^{n+1} - T_{g,1}^{n-1} } { 2 \Delta t }
    &=& - \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}} 
          \left\{ 
            F_s^{n+1} + F_L^{n+1} + F_c^{n+1} + L F_q^{n+1} + F_{SM}^{n+1} - F_g^{n+1}
          \right\}, 
\end{eqnarray}
%
となる. 
ここで, 
$F_g^{n+1}, F_s^{n+1}, F_L^{n+1}, F_c^{n+1}, L F_q^{n+1}, F_{SM}^{n+1}$
はそれぞれ, 地下への熱伝導フラックス, 短波放射フラックス, 長波放射フラックス, 
顕熱フラックス, 潜熱フラックス, 融雪による熱フラックスである. 

$F_s^{n+1}, F_c^{n+1}, L F_q^{n+1}, F_g^{n+1}$ は, 前節までの方法で連立一次
方程式を解くことで得られた値を用いることとし, 地表面温度 $T_s^{n+1}, T_{g,1}^{n+1}$ 
を介して $F_L^{n+1}$ を調節することで融雪量, すなわち $F_{SM}^{n+1}$ を求める. 
ここで, 
%
\begin{eqnarray}
  F_L^{n+1} &=& F_L^{n-1} 
    + \frac{\partial F_L}{\partial T_1} \left( T_1^{n+1} - T_1^{n-1} \right)
    + \frac{\partial F_L}{\partial T_s} \left( T_s^{n+1} - T_s^{n-1} \right)
\end{eqnarray}
%
であるから, (...), (...) より, 
%
\begin{eqnarray}
  -\kappa \frac{ T_{s}^{n+1} - T_{g,1}^{n+1} }{ z_\frac{1}{2} - z_1 }
    &=& F_s^{n+1}
      + F_L^{n-1} 
        + \frac{\partial F_L}{\partial T_1} \left( T_1^{n+1} - T_1^{n-1} \right)
        + \frac{\partial F_L}{\partial T_s} \left( T_s^{n+1} - T_s^{n-1} \right) \nonumber
\\
    & & + F_c^{n+1} + L F_q^{n+1} + F_{SM}^{n+1}, 
\\
  C_{g,1} \frac{ T_{g,1}^{n+1} - T_{g,1}^{n-1} } { 2 \Delta t }
    &=& - \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}} 
          \biggl\{
            F_s^{n+1} 
          + F_L^{n-1} 
              + \frac{\partial F_L}{\partial T_1} \left( T_1^{n+1} - T_1^{n-1} \right)
              + \frac{\partial F_L}{\partial T_s} \left( T_s^{n+1} - T_s^{n-1} \right) \nonumber 
\\
    & &   + F_c^{n+1} + L F_q^{n+1} + F_{SM}^{n+1} - F_g^{n+1}
          \biggr\}, 
\end{eqnarray}
%
となる. 

これら 2 つの方程式に対して, 未知数は $T_s^{n+1}, T_{g,1}^{n+1}, F_{SM}^{n+1}$ 
の 3 つである. 
そこでまず, 積雪すべてが解けた場合を想定し, $T_s^{n+1}, T_{g,1}^{n+1}$ を求める
ことにする. それによって得られた $T_s^{n+1}, T_{g,1}^{n+1}$ が不適切であった場合
($T_{g,1}^{n+1} < T_{cond}$) には, $T_{g,1}^{n+1} = T_{cond}$ として, 
$T_s^{n+1}, F_{SM}^{n+1}$ を求める. 

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

このとき, $F_{SM}^{n+1} = \frac{M_{\rm snow}}{2 \Delta t}$ であり, 
$T_s^{n+1}$, $T_{g,1}^{n+1}$ が未知数である. 
ここで, $M_{\rm snow}$ は積雪量である. 

(...) より, 
%
\begin{eqnarray}
  - \frac{\kappa}{ z_\frac{1}{2} - z_1 } T_{s}^{n+1}
  + \frac{\kappa}{ z_\frac{1}{2} - z_1 } T_{g,1}^{n+1}
    &=& F_s^{n+1}
      + F_L^{n-1} 
        + \frac{\partial F_L}{\partial T_1} \left( T_1^{n+1} - T_1^{n-1} \right)
        + \frac{\partial F_L}{\partial T_s} T_s^{n+1}
        - \frac{\partial F_L}{\partial T_s} T_s^{n-1} \\
    & & + F_c^{n+1} + L F_q^{n+1} + F_{SM}^{n+1}, 
\\
  \frac{\kappa}{ z_\frac{1}{2} - z_1 } T_{g,1}^{n+1}
    &=& F_s^{n+1}
      + F_L^{n-1} 
        + \frac{\partial F_L}{\partial T_1} \left( T_1^{n+1} - T_1^{n-1} \right)
        + \frac{\partial F_L}{\partial T_s} T_s^{n+1}
        - \frac{\partial F_L}{\partial T_s} T_s^{n-1} \\
    & & + F_c^{n+1} + L F_q^{n+1} + F_{SM}^{n+1} \\
    & & + \frac{\kappa}{ z_\frac{1}{2} - z_1 } T_{s}^{n+1}, 
\\
  T_{g,1}^{n+1}
    &=& \left( 1 
           + \frac{ z_\frac{1}{2} - z_1 }{\kappa} \frac{\partial F_L}{\partial T_s}
        \right) T_s^{n+1} \\
    & & + \frac{ z_\frac{1}{2} - z_1 }{\kappa}
          \biggl\{
            F_s^{n+1}
          + F_L^{n-1} 
            + \frac{\partial F_L}{\partial T_1} \left( T_1^{n+1} - T_1^{n-1} \right)
            - \frac{\partial F_L}{\partial T_s} T_s^{n-1} \\
    & &   + F_c^{n+1} + L F_q^{n+1} + F_{SM}^{n+1}
          \biggr\}, 
\\
  T_{g,1}^{n+1}
    &=& \left( 1 
           + \frac{ z_\frac{1}{2} - z_1 }{\kappa} \frac{\partial F_L}{\partial T_s}
        \right) T_s^{n+1} \\
    & & + \frac{ z_\frac{1}{2} - z_1 }{\kappa} \alpha
\end{eqnarray}
%
ここで, 
%
\begin{eqnarray}
  \alpha 
    &=&   F_s^{n+1}
        + F_L^{n-1} 
          + \frac{\partial F_L}{\partial T_1} \left( T_1^{n+1} - T_1^{n-1} \right)
          - \frac{\partial F_L}{\partial T_s} T_s^{n-1} \\
    & & + F_c^{n+1} + L F_q^{n+1} + F_{SM}^{n+1}
\end{eqnarray}
%
である. 
(...) を (...) に代入すると, 
%
\begin{eqnarray}
  \frac{ C_{g,1} }{ 2 \Delta t } & & \hspace{-9mm}
    \left\{
      \left( 1 
         + \frac{ z_\frac{1}{2} - z_1 }{\kappa} \frac{\partial F_L}{\partial T_s}
      \right) T_s^{n+1}
      + \frac{ z_\frac{1}{2} - z_1 }{\kappa} \alpha
    \right\}
  - \frac{ C_{g,1} }{ 2 \Delta t } T_{g,1}^{n-1} \nonumber \\
    &=& - \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}} 
          \biggl\{
            F_s^{n+1}
          + F_L^{n-1} 
              + \frac{\partial F_L}{\partial T_1} \left( T_1^{n+1} - T_1^{n-1} \right)
              + \frac{\partial F_L}{\partial T_s} \left( T_s^{n+1} - T_s^{n-1} \right) \nonumber \\
    & &   + F_c^{n+1} + L F_q^{n+1} + F_{SM}^{n+1} - F_g^{n+1}
          \biggr\}, 
\\
    &=& - \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}} 
          \biggl\{
            F_s^{n+1}
          + F_L^{n-1} 
              + \frac{\partial F_L}{\partial T_1} \left( T_1^{n+1} - T_1^{n-1} \right)
              + \frac{\partial F_L}{\partial T_s} T_s^{n+1}
              - \frac{\partial F_L}{\partial T_s} T_s^{n-1} \nonumber \\
    & &   + F_c^{n+1} + L F_q^{n+1} + F_{SM}^{n+1}
          \biggr\}, \\
    & & + \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}} F_g^{n+1}, 
\\
    &=& - \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}} 
          \frac{\partial F_L}{\partial T_s} T_s^{n+1} \\
    & & - \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}} 
            \biggl\{
              F_s^{n+1}
            + F_L^{n-1} 
                + \frac{\partial F_L}{\partial T_1} \left( T_1^{n+1} - T_1^{n-1} \right)
                - \frac{\partial F_L}{\partial T_s} T_s^{n-1} \nonumber \\
    & &     + F_c^{n+1} + L F_q^{n+1} + F_{SM}^{n+1}
            \biggr\}, \\
    & & + \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}} F_g^{n+1}, 
\\
    &=& - \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}} 
          \frac{\partial F_L}{\partial T_s} T_s^{n+1} \\
    & & - \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}} \alpha \\
    & & + \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}} F_g^{n+1}, 
\end{eqnarray}


%
\begin{eqnarray}
    \frac{ C_{g,1} }{ 2 \Delta t } & & \hspace{-9mm}
      \left( 1 
         + \frac{ z_\frac{1}{2} - z_1 }{\kappa} \frac{\partial F_L}{\partial T_s}
      \right) T_s^{n+1}
  + \frac{ C_{g,1} }{ 2 \Delta t } 
      \frac{ z_\frac{1}{2} - z_1 }{\kappa} \alpha
  - \frac{ C_{g,1} }{ 2 \Delta t } T_{g,1}^{n-1} \nonumber \\
    &=& - \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}} 
          \frac{\partial F_L}{\partial T_s} T_s^{n+1}
        - \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}} \alpha
        + \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}} F_g^{n+1}, 
\\
  \biggl\{
    \frac{ C_{g,1} }{ 2 \Delta t } & & \hspace{-9mm}
      \left( 1 
         + \frac{ z_\frac{1}{2} - z_1 }{\kappa} \frac{\partial F_L}{\partial T_s}
      \right) 
  + \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}} 
           \frac{\partial F_L}{\partial T_s}
  \biggr\} T_s^{n+1} \nonumber \\
    &=&
        - \left( 
            \frac{ C_{g,1} }{ 2 \Delta t } \frac{ z_\frac{1}{2} - z_1 }{\kappa}
          + \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}}
          \right) \alpha
        + \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}} F_g^{n+1} 
        + \frac{ C_{g,1} }{ 2 \Delta t } T_{g,1}^{n-1}
\\
  T_s^{n+1} 
    &=&
      \biggl\{
        \frac{ C_{g,1} }{ 2 \Delta t }
          \left( 1 
             + \frac{ z_\frac{1}{2} - z_1 }{\kappa} \frac{\partial F_L}{\partial T_s}
          \right) 
      + \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}} 
               \frac{\partial F_L}{\partial T_s}
      \biggr\}^{-1} \nonumber \\
    & & \left\{
          - \left( 
              \frac{ C_{g,1} }{ 2 \Delta t } \frac{ z_\frac{1}{2} - z_1 }{\kappa}
            + \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}}
            \right) \alpha
          + \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}} F_g^{n+1} 
          + \frac{ C_{g,1} }{ 2 \Delta t } T_{g,1}^{n-1}
        \right\}
\end{eqnarray}


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

このとき, $T_{g,1}^{n+1} = T_{cond}$ であり, $F_{SM}$, $T_s^{n+1}$ が未知数である. 

(...) より, 
%
\begin{eqnarray}
  F_{SM}^{n+1}
    &=& - \kappa \frac{ T_{s}^{n+1} - T_{g,1}^{n+1} }{ z_\frac{1}{2} - z_1 }
        - F_s^{n+1}
        - F_L^{n-1} 
          - \frac{\partial F_L}{\partial T_1} \left( T_1^{n+1} - T_1^{n-1} \right)
          - \frac{\partial F_L}{\partial T_s} \left( T_s^{n+1} - T_s^{n-1} \right) \nonumber \\
    & & - F_c^{n+1} - L F_q^{n+1} 
\end{eqnarray}

(...) に代入すると, 
%
\begin{eqnarray}
  C_{g,1} \frac{ T_{g,1}^{n+1} - T_{g,1}^{n-1} } { 2 \Delta t }
    &=& - \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}} 
          \biggl\{
            F_s^{n+1} 
          + F_L^{n-1} 
              + \frac{\partial F_L}{\partial T_1} \left( T_1^{n+1} - T_1^{n-1} \right)
              + \frac{\partial F_L}{\partial T_s} \left( T_s^{n+1} - T_s^{n-1} \right) \nonumber \\
    & &   + F_c^{n+1} + L F_q^{n+1} \\
    & &   - \kappa \frac{ T_{s}^{n+1} - T_{g,1}^{n+1} }{ z_\frac{1}{2} - z_1 }
          - F_s^{n+1}
          - F_L^{n-1} 
            - \frac{\partial F_L}{\partial T_1} \left( T_1^{n+1} - T_1^{n-1} \right)
            - \frac{\partial F_L}{\partial T_s} \left( T_s^{n+1} - T_s^{n-1} \right) \nonumber \\
    & &   - F_c^{n+1} - L F_q^{n+1} \\
    & &   - F_g^{n+1}
          \biggr\}, 
\\
    &=& - \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}} 
          \biggl\{ 
            - \kappa \frac{ T_{s}^{n+1} - T_{g,1}^{n+1} }{ z_\frac{1}{2} - z_1 }
            - F_g^{n+1}
          \biggr\}
\\
    &=&   \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}} 
            \kappa \frac{ T_{s}^{n+1} - T_{g,1}^{n+1} }{ z_\frac{1}{2} - z_1 }
        + \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}} F_g^{n+1}
\\
    &=&   \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}} 
            \frac{\kappa}{ z_\frac{1}{2} - z_1 } T_{s}^{n+1}
        - \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}} 
            \frac{\kappa}{ z_\frac{1}{2} - z_1 } T_{g,1}^{n+1}
        + \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}} F_g^{n+1}
\end{eqnarray}




%
\begin{eqnarray}
  \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}} 
            \frac{\kappa}{ z_\frac{1}{2} - z_1 } T_{s}^{n+1}
    &=& C_{g,1} \frac{ T_{g,1}^{n+1} - T_{g,1}^{n-1} } { 2 \Delta t }
        + \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}} 
            \frac{\kappa}{ z_\frac{1}{2} - z_1 } T_{g,1}^{n+1}
        - \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}} F_g^{n+1}
\\
  T_{s}^{n+1}
    &=& 
      \left( z_\frac{1}{2} - z_\frac{3}{2} \right) \frac{ z_\frac{1}{2} - z_1 }{\kappa}
        \left( 
          C_{g,1} \frac{ T_{g,1}^{n+1} - T_{g,1}^{n-1} } { 2 \Delta t }
          + \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}} 
              \frac{\kappa}{ z_\frac{1}{2} - z_1 } T_{g,1}^{n+1}
          - \frac{1}{z_\frac{1}{2} - z_\frac{3}{2}} F_g^{n+1}
        \right)
\end{eqnarray}

