【[注(2006-12-16)] 2006年12月21日の授業では、 実習課題としてではなく、 データの応用例として紹介する。】

大気・流域水収支

エネルギー循環・水循環にかかわる量を求めるデータ処理の例として、 大気・流域水収支([講義])の 計算をしてみたい。 直接観測されていないが収支計算の結果として求まる量は、 蒸発量Eと、陸水貯留量の時間変化dS/dtである。 ただし、ここでは、水蒸気収束量C (鉛直積算した水蒸気の水平輸送の収束量)を求めるのに、 ラジオゾンデ観測データをそのまま使うことにする。

【近ごろの研究では、Cを求めるのに気象データ同化 ([講義])で得られた 格子点データを使うことが多い。 GAME再解析 ([別ページ参照])も使われている。 温帯の空間スケール数百km以上の領域では、 複数の同化データから求めたCがほぼ一致し、 おそらく現実の値をよく近似していると思われる。 (空間スケールの制約は、 観測データの代表的時間間隔の間に空気が動く距離の内側の 情報は混ざってしまうと考えればよいようだ。 水蒸気と風の鉛直分布に関する主要な観測データである ラジオゾンデの観測が多くの国で12時間間隔である。) しかし、熱帯では、同化データによるCは、 観測値よりも予報モデルの降水・対流過程の表現の影響を強く受けるようで、 複数の同化データから得られたCが大きくくいちがうことが多い。 同化データに頼らなければならないことも多いのだが、 その質を判断するために、 同化に頼らず観測データだけで求めてみることも重要である。 観測データだけによる場合、観測地点の間は単純内挿によるしかなく、 それで得られたCが同化データによるものよりも現実に近いと、 必ずしも言い切れないという問題は残るが。】

同化データを使う計算を、 [実習10]として追加する。

対象地域・時期

水蒸気収束量Cを求めるのに、 GAME-Tプロジェクトのラジオゾンデ強化観測 ([実習3]) から求めることにする。 対象地域と時期はこれに制約されて決まってくる。 対象地域は、Bangkok, Ubon (正式にはUbon Ratchathani), Nong Khaiの 3地点でかこまれた三角形領域とする ([地図]参照; 3地点の緯度経度は、[表のページ]参照)。 この領域の大部分は、Mekong川の支流であるMun川(Nam Mun)の流域である (Bangkokはよく知られているようにChao Phraya川(Mae Nam Chao Phraya)の流域、 Nong KhaiはMekong川本流沿いだが) ので、三角形領域の大気水収支で求められたP-Eがこの流域をも代表するとして、 流域水収支も考えてみることにする。 対象時期はラジオゾンデの現業観測地点でないNong Khaiでの観測が制約条件になる。 ここでは、1998年5月後半から6月前半の強化観測期間のうちで、 3地点のゾンデ観測がそろった最初の日時から最後の日時までとする。 (前後に観測のある場合の欠損は時間内挿で補えるが、時間外挿はしないという方針)。

水蒸気収束量C、鉛直積算した水蒸気量W

対象時期の全部の時点について、ゾンデデータから、 鉛直積算した水蒸気量と水蒸気輸送ベクトル(東西成分・南北成分)を計算する プログラムはすでに用意してある ([実習3の作業3])。 時間平均の計算と、ベクトルの収束(発散の符号を逆にしたもの)の計算は いずれも線形演算なので、原理的にはどちらを先にしても結果は同じになる (演算誤差や、データ欠損の埋めかたからくる違いはあるが)。 ここでは、 3地点の観測がそろっていることが確認された時刻のデータだけを使うことにし、 各地点の水蒸気輸送ベクトルの東西成分・南北成分の期間平均をまず計算しておき、 それから期間平均の収束量を求めることにする。

観測値の品質検査。

データがないところのほかに、 ゾンデデータから鉛直積算した空気の(面積あたり)質量が (たとえば)7000 kg/m²に達しないところは欠損とみなす。 できれば、観測値の鉛直分布をみて、 かなり長い高さ区間で観測が欠けているものや、 機器あるいは通信の異常と思われる値が出ているものは除いたほうがよい。 また、積算した結果の時系列のグラフをみて、 異常と思われる値を除いたほうがよいかもしれない。

期間平均。

ここでは、次のいずれかの方法をとることにする。

• 有効な(欠損でない)観測値がそろった時刻のデータを平均したものを 期間平均とみなす。
• データ欠損の時刻の値を、前後のデータのある時刻の値から直線 (1次関数)内挿して埋め、それから時間平均する。

熱帯では降水の日変化が大きい (雨が1日のうちのある時間帯に集中しやすい)ことが知られており、 この対象地域でも確かめられている (気象研究ノート「東南アジアのモンスーン気候学」 [読書ノート]の大沢の報告、 里村の報告など参照)。 したがって水蒸気収束量の日変化も大きいと思われ、 データ欠損の時刻の値を推定するには、このことを考慮に入れるべきかもしれない。 しかしその方法を考えるのは次の段階としたい。

座標系。

グローバルや高緯度の地域の計算では、球面幾何で計算すべきだが (さらに精密に楕円体を考慮することもあるが) 、 熱帯の地域では平面で近似しても比較的誤差が小さいので、 ここでは平面幾何で計算することにする。 地球の周の長さを40000 km、緯度1度の長さはこれの1/360、 代表緯度を北緯15度とし、経度1度の長さは緯度1度の長さにcos(15度)をかけた ものとする。 北緯15度、東経100度を原点とし、東向きにx、北向きにyをとって、 適当な単位(m, kmまたは緯度1度の長さ)の座標系で計算するとよい。 ただし収束量の定義で分母にくる長さの単位はmなので、必要に応じて換算すること。 [地点表(緯度、経度の情報を含む)のファイルを読んで、 ここで採用した平面x,y座標値(長さの単位は緯度1度の長さ)で 書き出すプログラム例を /pub/sfc/earthobs/example/sonde/wb/vertex_xy.awk として用意した。 入力ファイル例は同じディレクトリのvertex_lola.txt。 使い方は
awk -f vertex_xy.awk vertex_lola.txt > vertex_xy.txt 。 なお、プログラム中のpiは円周率で、角度の「度」とラジアンの 換算のために必要となる。 数値の書きまちがいを避けるため、 awkの組みこみ関数atan2を使って 求めるようにした。

収束量の計算。

• 三角形の辺上で水蒸気輸送ベクトルが(その東西成分・南北成分がそれぞれ) 直線的に(1次関数で)変化していると仮定する。 すると、辺上の水蒸気輸送ベクトルは両端の値の単純平均で代表することができる。 (水蒸気量と風とがそれぞれ直線的に変化していると仮定した場合とは 結果が一致せず、どちらがよい仮定かは検討が必要な問題だが、 ここでは計算手順の便宜上これでよいとしておく。)
• 水蒸気輸送ベクトルの、辺に垂直な方向の成分を求める。
  • 辺の方向の単位ベクトルは、1頂点から他頂点に向かうベクトルをその長さで割れば 求められる。
  • 辺に垂直な方向の単位ベクトルは、辺の方向の単位ベクトルを90度回転すればよい。 平面幾何では、x成分= -旧y成分、y成分= 旧x成分。
  • ベクトルQの、辺に垂直な方向の成分は、ベクトルQと辺に垂直な方向の単位ベクトルとの内積をとればよい。 (ここでは2次元直交直線座標で考えているので、 二つのベクトルの内積を求めるには、 x成分どうし、y成分どうしの積を求めて合計すればよい。) 符号の約束を、三角形の中向きか外向きのどちらかを常に正にするように決めておくこと。
• これに辺の長さをかければ、(単位時間あたりに)辺をとおって三角形領域に出入りする水蒸気の質量が求められる。
• これを、三角形の中向きを正にとって集計し、三角形の面積で割れば、 水蒸気収束量が求められる。
• プログラム例を /pub/sfc/earthobs/example/sonde/wb/inflow.awk として用意した。 ただし、最後の三角形の面積で割るところは行なっていない。 また、この中で使われている辺の長さの単位は緯度1度の長さなので、 結果をSI単位で示すには換算(緯度1度の長さをmで表わしたものをかける)が必要である。 計算手順は上の記述に対応するが、 辺の長さで割ってまたかける演算を省略したので中間変数は少し違う。 水蒸気輸送ベクトルのデータファイルたとえばqvec.txtとする)と、
  awk -f inflow.awk qvec.txt > inflow.txt
  のように起動する。ただし、vertex_lola.txtも 同じディレクトリにあることを想定している。 入力ファイル(例qvec.txt)は、 各行に空白くぎりで水蒸気輸送ベクトルの東西成分と南北成分の2つの数値、 1行が1地点に対応し、地点の順序はvertex_lola.txtと 同じでなければならない。 出力ファイル(例inflow.txt)は1行で、 空白くぎりで、まず正味中向きに入る水蒸気質量の合計、 それに続いて3辺それぞれの寄与を順に示す。

領域を代表する水蒸気量Wの計算。

領域の中にも観測点があればそれも使って内挿すべきだが、 ここでは頂点にしかないので、3頂点の単純平均でよいとしておく。 なお、GPS観測によるデータ ( [講義] [実習6] )を使うこともできるが、 ここでは観測値の一貫性からゾンデ観測を主の材料としておく。

水蒸気量の時間変化dW/dt。

ここでは対象期間の平均的時間変化をみる。次のどちらかの方法をとる。

• 有効なデータのある最後と最初の時刻の値の差を、期間の長さで割る。
• 対象全時刻のデータに直線を最小2乗法であてはめて傾きを採用する。

時間で割る際の時間の単位は秒であることに注意。

降水量P

空間内挿してから集計するのがふつうだが、 ここでは、 対象地域の三角形の中に含まれている雨量計観測をもとに、 どの地点も同じ程度に代表性があるとみなして単純平均で求めることにする。 まずは、[実習2]で扱った タイ気象庁(TMD)の日降水量データを使ってほしい。 メコン川委員会(MRC)のデータもあるが、 そのうちにはまだ、TMDとの重複やまちがいと思われるものがあるので、 増田の判断でそれを省いたものを用意する予定だが、 レポートの期限に間に合わない可能性が高い。 なお、日降水量の定義は、TMD, MRC(少なくともタイの地点で)のいずれも、 世界時0時(地方標準時7時)までの24時間累積値を前日の降水量としている。

雨量計の観測地点一覧(緯度、経度の情報を含む)のファイルを読んで、 三角形の内側にある地点だけ、この課題で採用した平面x,y座標値で 書き出すプログラム例を /pub/sfc/earthobs/example/precip/triangle/inside_stn_xy.awk として用意した。 このプログラムと、同じディレクトリにあるvertex_lola.txt (三角形の頂点の緯度経度の表、上記sonde/wbにあるものと同じ)を作業ディレクトリに置いて
awk -f inside_stn_xy.awk /pub/sfc/earthobs/data/gamet/precip/st_list/station_TMD.txt > inside_stn_xy.txt
のように実行すると、 得られたinside_stn_xy.txtは三角形の内側の地点の表であり、 空白区切りで左からx座標、y座標、地点ID、地点名である。 なお、vertex_lola.txtの点の順序は、上から見て反時計まわりに しておく必要がある。 もし逆まわりだと、三角形の外側の地点を選択することになってしまう。

流出量R

MRCのLower Mekong Hydrologic Yearbookに収録された 河川流量データから、 Mun川のUbonでの1998年の毎日の流量(単位、m³/s)のデータを、 /pub/sfc/earthobs/data/river/mun/mun_ubon_1998.txt におく。 データ形式は空白くぎりで左の列から年、月、日、通しの日、流量とした。 水位(単位 m、基準はローカルに決められている)のデータも levelで始まるファイル名で置く。 収支式のRは面積あたりの流出量(流出高ともいう)であり、 (流量観測点から下流に出ていく以外の水の出入りが無視できることを前提として) 流量を流域面積で割ることによって求められる。 この地点から上流の流域面積は104000 km²である。 (単位がm²でないことに注意)。

計算結果の例 (別ページへ)