% 表題   DCPAM5 第1部 物理過程 放射
%
% 履歴 
%\Drireki{93/03/18 沼口敦・保坂征宏}
%\Drireki{2010/04/14 高橋芳幸}
%\Drireki{2010/11/04 石渡正樹}
%

\section{大気上端での恒星の放射フラックス}

惑星上の 1 点における大気上端での恒星の放射フラックス $F_{0}(\chi)$ は 
%
\begin{eqnarray}
  F_{0}(\chi) &=& 
 \left\{
  \begin{array}{ll}
   F_{00} \left( \frac{1}{r_{S}} \right)^2 \cos \chi & (\cos \chi >0) \\
   0 & (\cos \chi \leq 0)
  \end{array}
\right.
  \Deqlab{大気上端入射放射フラックス}
\end{eqnarray}
%
と書くことができる. 
$F_{00}$ は軌道長半径における恒星の放射フラックスであり, 太陽定
数に相当する\footnote{太陽定数のより正確な定義は
恒星-惑星間の平均距離における恒星の放射フラックスである.
ケプラー運動する惑星の恒星からの平均距離は
\begin{eqnarray}
  a \left( 1 + \frac{1}{2} e^2 \right)  \nonumber
\end{eqnarray}
となる. ($a$ は軌道長半径, $e$ は離心率). 
$e \ll 1$ であれば平均距離は $a$ とほぼ等しい.
}.
$r_{S}$ は惑星の軌道長半径で規格化した恒星-惑星間距離である.
$\chi$ は考えている地点における恒星の天頂角である.


\vspace{5mm}

$\cos \chi$ は, 
%
\begin{eqnarray}
  \cos \chi &=& \cos \phi \cos \delta_S \cos H + \sin \phi \sin \delta_S
  \Deqlab{天頂角余弦}
\end{eqnarray}
%
と表わされる. 
ここで, $\phi$ は緯度, $\delta_S$ は恒星の傾斜角
(惑星の赤道面から測った太陽の角度である. 
Liou, 2002 によれば declination of the sun.  
%Hartmann によれば declination angle of the sun, 
%Oort and Peixiot によれば solar declination
天球上における恒星の赤緯とも等しい
)
である.
$H$ は時角 (hour angle)であり, 
考える点の現在の位置と正午になる時の位置との経度の差
(もしくは恒星直下点の経度を基準にした経度)である. 
\Deqref{天頂角余弦}は球面三角関数の公式を使えば導くことができる
(Liou, 2002 の2.2 節および Appendix C 参照).
各種の角度の関係を \Dfigref{惑星恒星系における各種角度の関係} と
\Dfigref{天頂角・緯度・恒星の傾斜角・時角の関係} に示しておく.
なお, \Dmodel の放射計算においては, $\cos \chi$よりも $\sec \chi$ の形の方が
便利であるので, 変数としては $\sec \chi$ の値を格納したものを
用意している.

\begin{figure}
      \begin{center}
       \Depsf[120mm]{radiation/radiation-math-images/Liou_2002_Fig2.5.ps} 
      \end{center}
      \caption[惑星恒星系における各種角度の関係]{
            {\bf 惑星恒星系における各種角度の関係}.
            原図は Liou (2002) の Figure 2.5. いずれ自分で書きかえないと
            いけないだろう.
            図中の $\delta$ が $\delta_S$ に, 
            $\nu$ が$\Phi$ に対応する.
            }
      \Dfiglab{惑星恒星系における各種角度の関係}
\end{figure}

\begin{figure}
      \begin{center}
       \Depsf[120mm]{radiation/radiation-math-images/Liou_2002_Fig2.6.ps}
      \end{center}
      \caption[天頂角・緯度・恒星の傾斜角・時角の関係]{
            {\bf 天頂角・緯度・恒星の傾斜角・時角の関係}.
            原図は Liou (2002) の Figure 2.6. いずれ自分で書きかえないと
            いけないだろう.
            図中の $\theta_0$ が $\chi$ に, 
            $\delta$ が$\delta_S$ に対応する.
            }
      \Dfiglab{天頂角・緯度・恒星の傾斜角・時角の関係}
\end{figure}

$r_S$ は惑星の軌道要素から次のように計算される 
(ランダウ・リフシッツ「力学」\S 15 参照).
\begin{eqnarray}
  r_S &=& ( 1 - e \cos \xi ) 
  \Deqlab{惑星恒星間距離}
\end{eqnarray}
ここで, $e$ は離心率である.
$\xi$  は離心近点角 (または離心近点離角, eccentric anomaly) であり, 
楕円上の位置を与えるパラメータであり
\footnote{楕円の中心を原点に, 近日点の方向を $x$ 軸にとった時に
          楕円上の点 $(x,y)$ と離心近点角 $\xi$ との関係は以下のように
          なる.
          \begin{eqnarray}
	    x &=& a (\cos \xi - e), \nonumber \\
	    y &=& a \sqrt{1-e^2}\sin \xi, \nonumber
	  \end{eqnarray}
          ただし, $a$ は軌道長半径, $e$ は離心率である.
 }, 
ケプラーの方程式
\begin{eqnarray}
 \xi - e \sin \xi = l
 \Deqlab{ケプラーの方程式}
\end{eqnarray}
を解くことによって求められる.
ここで $l$  は平均近点角 (mean anomaly)であり, 
惑星軌道を円で近似して, 軌道中心を原点として近日点の方向から測った
惑星の位置までの角度である.
時刻を $t$, 公転周期を $T_{orb}$ とすると
\begin{eqnarray}
 l =   \frac{2 \pi (t- t_0)}{T_{orb}}
     + ( \Phi_{Epoch} - \Phi_0 ) \frac{\pi}{180}
\end{eqnarray}
となる.
$t_0$ は元期 (げんき, 天体位置計算の基準となる時刻. 
天体観測では軌道要素が観測された日を元期とする) を示す時刻である.
$\Phi_{Epoch}$ は
元期における惑星の平均近点角(単位は degree, 元期における太陽の黄経に
$\pi$ を足したもの), 
$\Phi_0$ は近日点黄経
である.
\Dmodel では,  各時刻における $l$ を求めた後に 
Newton 法によって \Deqref{ケプラーの方程式} を
$\xi$ について解いている.

\vspace{5mm}

$\delta_S$ は以下の式で計算される. 
\begin{eqnarray}
  \sin \delta_S = - \sin \theta_p \sin (\Phi_0 + \Phi)
  \Deqlab{恒星の傾斜角}
\end{eqnarray}
$\theta_p$ は赤道傾斜角 (天体の軌道面と赤道面のなす角.
Liou, 2002 によれば oblique angle of the earth's axis.
自転軸と公転軸のなす角にも等しい). 
$\Phi$  は真近点角 (true anomaly) であり, 
恒星を原点として, 近日点から測る惑星の軌道上の位置を表す角度である.
$\Phi$ は以下の式から決定される
\begin{eqnarray}
  \tan \frac{\Phi}{2}
    &=& \sqrt{ \frac{1 + e}{1 - e}}
        \tan \frac{\xi}{2}
\end{eqnarray}
$\Phi_0$ は近日点黄経であり, 春分点の方向と近日点のなす角である.
$\Phi_0 + \Phi$ は
恒星を原点として, 春分点の方向から惑星の位置まで測った角度
となっている.

\vspace{5mm}

時角 $H$ は以下の式で決定する.
\begin{eqnarray}
   H = 2 \pi t_{ByDay} - \pi + \lambda
\end{eqnarray}
$t_{ByDay}$ は時刻を日単位で表現したものである
($t_{ByDay}=0$ が深夜 0 時に, $t_{ByDay}=0.5$ が正午に対応する).
日の出と日の入りの時の時角 $H_0$ は
\begin{equation}
   \cos H_0 = - \tan \phi \tan \delta_S 
\end{equation}
となる.

\vspace{5mm}

以下では, 現在 \Dmodel に実装されている日射分布の計算方法に
関する記述を行う.
\begin{enumerate}
\item 年変化(季節変化)を日変化する場合.

      この場合には, 上記の \Deqref{大気上端入射放射フラックス}, 
      \Deqref{天頂角余弦}, 
      \Deqref{惑星恒星間距離}, \Deqref{恒星の傾斜角}
      を用いて緯度・経度, 時間を与えた場合の大気上端における
      恒星からの放射フラックス分布を計算する.

      \Dmodel のデフォルトのパラメータ設定を使って
      計算した大気上端における日平均日射量の時間--緯度分布を
      \Dfigref{日平均太陽放射(デフォルト)} に, 
      現実の惑星にあわせたパラメータ設定を使って計算した
      日平均日射量の時間--緯度分布を
      \Dfigref{日平均太陽放射(地球・火星)} に示す.
      現実の地球の場合の結果 (\Dfigref{日平均太陽放射(地球・火星)}b)
      は, Liou (2002) の Figure2.8 で示された結果と同じパターンになっている.

      \begin{figure}
      \begin{center}
       \Depsf[60mm]{radiation/radiation-math-images/DailyMeanInsolation-Default.eps}
      \end{center}
      \caption[日平均太陽放射時間-緯度分布]{
            {\bf \Dmodel の入射太陽放射ルーチンのデフォルト設定を用いて得られる
            大気上端における日平均太陽放射の時間-緯度分布}.
            横軸は1 年の開始日からの日数, 縦軸は緯度.
            大気上端における太陽放射分布を 1 時間ごとに計算し, 
            日平均をとった値を示している.
            \Dmodel の大気上端での恒星の放射フラックスを与える
            サブルーチンを用いて計算した. 水平解像度は T21.
            $F_{00} = 1380$ W/m$^{-2}$, $\theta_p = 23.5^{\circ}$, 
            $\Phi_0= 0.0$, 
            $\epsilon = 0.0$, 
            元期における惑星の黄経は $280.0$ とした場合.
            1 年の長さは 365 日.
            }
      \Dfiglab{日平均太陽放射(デフォルト)}
      \end{figure}

      \begin{figure}
      \begin{center}
      \begin{tabular}{cc}
       \raisebox{80mm}{\large (a)} 
       \Depsf[60mm]{radiation/radiation-math-images/DailyMeanInsolation-Earth.eps} &
       \raisebox{80mm}{\large (b)} 
       \Depsf[60mm]{radiation/radiation-math-images/DailyMeanInsolation-Mars.eps}
      \end{tabular}
      \end{center}
      \caption[日平均太陽放射時間-緯度分布(地球・火星)]{
            {\bf 現実の惑星の設定を用いて得られる
            大気上端における日平均太陽放射の時間-緯度分布}.
            横軸は1 年の開始日からの日数, 縦軸は緯度.
            大気上端における太陽放射分布を 1 時間ごとに計算し, 
            日平均をとった値を示している.
            \Dmodel の大気上端での恒星の放射フラックスを与える
            サブルーチンを用いて計算した. 水平解像度は T21.
            (a) 現実の地球を模したパラメータ設定を用いた場合.
            $F_{00} = 1369$ W/m$^{-2}$, $\theta_p = 23.44^{\circ}$, 
            $\Phi_0= 102.768413 + 180.0$, 
            $\epsilon = 0.016713$, 
            元期における惑星の黄経は $99.403308 + 180.0$. 
            1 年の長さは 365 日.
            (b) 現実の火星を模したパラメータ設定を用いた場合.
            $F_{00} = 588.98$ W/m$^{-2}$, $\theta_p = 25.19^{\circ}$, 
            $\Phi_0= 258.98$, 
            $\epsilon = 0.0934$, 
            元期における惑星の黄経は $-10.342$,
            1 年の長さは 669 日.
            }
      \Dfiglab{日平均太陽放射(地球・火星)}
      \end{figure}


\item 年平均・日平均日射分布を用いる場合

      %年平均, 日平均日射分布に関する記述はいずれまた
      %\footnote
      %{
      %  地球のパラメータの場合の式・数値はあるが, あまり書く気にならないな.
      %}.

      \Dmodel においては, 年平均入射量および年平均入射角は, 
      以下の近似式を用いて計算している.
      \begin{equation}
	\overline{F_S^I} (\varphi)
	\simeq - S_0 ( A_{ins} + B_{ins} \cos^2 \varphi ) ,
	     \Deqlab{年平均入射放射分布}
      \end{equation}
      \begin{equation}
	\overline{\cos \chi} \simeq A_{\chi} + B_{\chi} \cos^2 \varphi .
      \end{equation}
      なお, 
      短波放射の放射伝達方程式で必要となる $\sec \zeta$ は 
      \begin{eqnarray}
	\overline{\sec \chi} = \frac{1}{A_{\chi} + B_{\chi} \cos^2 \varphi} 
      \end{eqnarray}
      として計算する.

      $A_{ins}$, $B_{ins}$, $A_{\chi}$, $B_{\chi}$ の値
      を \Dtabref{Ains} に示す. 
      これらの値は AGCM5 で使用されていたものであり, 
      どのように決定されたのかについては正確なところは確認されていない.
      しかし, \Deqref{大気上端入射放射フラックス} によって
      入射放射量を計算し日平均・年平均した結果を用いて, 
      最小二乗法で \Deqref{年平均入射放射分布} へのフィッティングを行うと
      \Dtabref{Ains}に示した $A_{ins}$ および $B_{ins}$ の値とほぼ等しい数値が得られる
      (T42 で計算した場合 $A_{ins}=0.12756$, $B_{ins}=0.18340$ となる).
      $A_{\chi}$, $B_{\chi}$ については, \Deqref{大気上端入射放射フラックス}
      の時間平均を取ったものが \Deqref{年平均入射放射分布} であると考えれば
      $A_{ins}$, $B_{ins}$ をそれぞれ定数倍したものが
      $A_{\chi}$, $B_{\chi}$ になるはずである.
      赤道での$\overline{\cos \chi}$ の値が 1 になるように定数を決めると
      \Dtabref{Ains} の $A_{\chi}$ および $B_{\chi}$ と等しい値が得られる
      (T42 の計算で得られた $A_{ins}=0.12756$, $B_{ins}=0.18340$ を用いると
      $A_{\chi}=0.41021$, $B_{\chi}=0.58979$ となる).

      \begin{table}
	\begin{center}
	  \begin{tabular}{|c|c|c|c|}
	  $A_{ins}$ & $B_{ins}$ & $A_{\chi}$ & $B_{\chi}$ \\
	  0.127 & 0.183 & 0.410 & 0.590   
	  \end{tabular}
	\end{center}
	\caption{現実の地球を想定した場合の
		$A_{ins}$, $B_{ins}$, $A_{\chi}$, $B_{\chi}$ の値}
	\Dtablab{Ains}
      \end{table}

      %日平均日射の場合
      %\begin{eqnarray}
      % h= 0
      %\end{eqnarray}
      %と設定している. これ必要なのか???

      ちなみに, 日平均日射分布・年平均日射分布に関する正確な表式は以下の通りで
      ある. これらの式に基づき年平均・日平均日射分布の表式を構成できるはずであ
      るがやっていない.
      \begin{itemize}
      \item 日平均日射分布の正確な表式

	    Liou (2002) によれば, 
	    日平均放射量 $F_{d}$ は次のように計算される. 
	    \begin{eqnarray}
		  \overline{F_{day}}(\phi)
	      &=& F_{00} \left( \frac{1}{r_S} \right)^2 
		  \frac{S(r)}{\pi} \nonumber \\
              &&  \quad \times ( \cos \phi \sin h_{0} \cos \delta_S
				  + h_{0} \sin \phi \sin \delta_S ).
	    \end{eqnarray}
	    ここで, $\delta$ と $S(r)$ の 1 日の間での
	    変化量は小さいとする近似を用いている.

      \item 日平均・年平均日射分布の正確な表式

	    Liou (2002) によれば, 日平均・年平均日射分布は
	    \begin{eqnarray}
		   \overline{F_y} (\phi) 
	       &=& \frac{F_{00} T_{orb} \tilde{S} (\phi, \epsilon)}{\pi (1-e^2)^{1/2}}, \\
      %
		    \tilde{S} (\phi, \epsilon) 
	       &\equiv&  \frac{\sin \phi \sin \epsilon}{2 \pi}
		    \int^{2 \pi}_{0} (h_0 - \tan h_0) \sin \lambda d \lambda
	    \end{eqnarray}
	    で与えられる.
      \end{itemize}

      参考として, North (1975) で用いられている式も挙げておく.
      North (1975) では, 大気上端における恒星の放射フラックスの年平均・
      日平均分布を
      \begin{eqnarray}
         F_0(x) &=& \frac{F_{00}}{4} \left\{ 1 + S_2  P_2 (x) \right\}, \\
         S_2 &=&  - 0.482
      \end{eqnarray}
      として, エネルギーバランスモデルによる計算を行っている.
      ここで $x = \sin \phi$ である. $S_2$ の値は, 
      Ch\'{y}lek and Coakley (1975)
      の地球における長波放射の吸収量の観測値に基づき決定したものである
      (当時はまだ短波放射に関する衛星観測が無かったのだと思われる).
      これから, $A_{ins}$, $B_{ins}$ に対応する量を計算すると
      \begin{eqnarray}
          A_{ins} &=& 0.1295, \\
          B_{ins} &=& 0.1808
      \end{eqnarray}
      となる.

\item 特定の日の日射分布を用いる (perpetual run).

      この場合には, $\sin \delta_S$ と $r_{S}$ に定数値を与えて 
      \Deqref{大気上端入射放射フラックス}, 
      \Deqref{天頂角余弦}を用いて, 大気上端における
      恒星からの放射フラックス分布を計算する.


\item 昼半球・夜半球固定の日射分布 (同期回転惑星設定)

      この場合には, 太陽直下点の経度 $\lambda_{subsolar}$ と経度
      $\phi_{subsolar}$ を与える. 
      $\phi_{subsolar}=0$ の場合を考えて, 天頂角を
      \begin{eqnarray}
	\cos \chi  = \cos \phi \cos (\lambda - \lambda_{subsolar})
      \end{eqnarray}
      とする.
      これにより, 恒星からの放射フラックス分布を
      \begin{eqnarray}
         F_0 (\phi) = F_{00} \cos \chi
      \end{eqnarray}
      で決定する. 

\end{enumerate}


\section{放射計算で用いるパラメータ}

放射計算で指定するべきパラメータの主なもの
(いずれ網羅する予定)を以下に挙げる.

\begin{itemize}
\item 短波に対する大気アルベド:
      \Dmodel のソースコードに記述されている値は 0.2 .

      現実の地球の場合では, 短波に対する大気アルベドは 0.225 である.
      Kiehl and Trenberth (1997) による地球全体の熱収支の
      見積りでは, 全球平均放射量 342 W/m$^2$ のうち, 
      77 W/m$^2$ が大気により反射される.

\item 軌道の離心率:
      \Dmodel のソースコードに記述されている値は 0.0.

      現実の地球の場合, 軌道離心率は 0.0167 (理科年表による).

      現実の火星の場合, 軌道離心率は 0.0934 (Allison, 1997).

\item 近日点黄経:
      \Dmodel のソースコードに記述されている値は 0.0 .

      現実の地球の場合, 近日点黄経は 102.924$^{\circ}$ (理科年表による)
      \footnote{Duffett-Smith (1988) も確認するべし.}
      \Dmodel で使用する場合には 102.768413 + 180.0 を与える.

      現実の火星の場合, 近日点黄経は 250.98$^{\circ}$ (Allison, 1997).

\item 太陽定数:
      \Dmodel のデフォルト値は 1380 W/m$^2$. 
      この値は Ishiwatari et al. (2002) による. 以下に述べる
      ように現実の地球の太陽定数よりもやや大きい\footnote{デフォルト値も
      1370 W/m$^2$ にした方が良いかも.}

      現実の地球の場合, 太陽定数は 1367 W/m$^2$ である (Hartmann, 1994).
      太陽が単位時間に発する放射エネルギー $L_0$ として
      $3.85 \times 10^{26}$ W (理科年表), 太陽地球間の平均距離 
      $\overline{r}_S$ として 1 天文単位 ($1.496 \times 10^{11}$ m; 理科年表, 1995) 
      を用いると, 
      太陽定数は
      \begin{eqnarray}
        \frac{L_0}{4 \pi \overline{r}^2_S }
       = \frac{3.85 \times 10^{26}}{4 \times 3.142 \times (1.496 \times 10^{11})^2}
       = 1368.8 {\rm W}/{\rm m}^2
      \end{eqnarray}
      と計算される.

     現実の火星の場合, 太陽定数は 588.98 W/m$^2$ (Kieffer et al., 1992)

\end{itemize}
