1．はじめに
【流出解析】シリーズの第6回ではタンクモデルを紹介しました． その記事の最後の方でも簡単に触れましたが，タンクモデルを使うときは積雪量や融雪量の計算を求められることが多いです．
というのも，冬に降った雪は山の上に蓄えられ，春先に気温が上がるまで山の上にずーっと留まったままです． そして，春になると融け出して，雪解け水が川に流れ込みます． つまり，降雪が数か月遅れて流出に寄与するわけです．
普通の4段タンクモデルでは，この遅れをうまく表せません． 積雪量のデータが手に入れば計算に組み込めるのですが，残念ながら日本における積雪量データの空間的な充実度はあまり高くはありません．
そこで，積雪量や融雪量を推定するために『積雪・融雪モデル』というものが作られました． いくつかある積雪・融雪モデルのうち，今回は菅原正巳さんが開発した計算方法で降雪量と融雪量を推定してみたいと思います．
菅原正巳さんといえば，タンクモデルを開発した人です． さすが開発者が同じだけあってタンクモデルとは相性が良く，そして何より使いやすいです．
なお，積雪・融雪モデルには，気温日数法や熱収支法といったモデルもあります． 熱収支法であれば，物理的根拠がはっきりしているという利点がありますが，計算に必要なデータが揃わない場合もあります． 一方で，菅原さんの方法は気温と降水量のデータのみで積雪量・融雪量の計算が可能です．
計算用のワークシートは下記URLからダウンロードできます． ご自由にお使いください． 積雪・融雪モデル.xlsm
【流出解析】シリーズ 一覧
【流出解析】シリーズ 第1回 流域界の決定と入力データの用意 【流出解析】シリーズ 第2回 貯留関数法【前編】 【流出解析】シリーズ 第3回 貯留関数法【後編】 【流出解析】シリーズ 第4回 流域実蒸発散量の推定―蒸発散比法― 【流出解析】シリーズ 第5回 実蒸発散量の推定―補完法― 【流出解析】シリーズ 第6回 タンクモデル 【流出解析】シリーズ 第7回 積雪・融雪モデル ―降雪量・融雪量の推定― 【流出解析】シリーズ 第8回 タンクモデルの最適化 【流出解析】シリーズ 第9回 単位図法 【流出解析】シリーズ 第10回 流出量・流出高・比流量の違い 【流出解析】シリーズ 第11回 表面流モデル【定常降雨：前編】 【流出解析】シリーズ 第12回 表面流モデル【定常降雨：後編】
2．積雪・融雪モデル
積雪・融雪モデルの作業は大きく5つのステップに分けられます．
ステップ1．流域のいくつかの地帯に分割 ステップ2．地帯別に日平均気温を推定 ステップ3．地帯別に降雨と降雪を判別し，降雨量，積雪量を計算 ステップ4．地帯別に融雪量を計算 ステップ5．流域全体での降雨量，降雪量，積雪量，融雪量を計算
ただし，ステップ1は最初に1回やるだけでいいので，実質的にはステップ2〜5の計算を日単位で繰り返します．
なお，改良菅原雪モデルと呼ばれるモデルもあり，改良版では降雨や融雪水の一部が積雪層に貯留されたり，再凍結する過程を加味してたモデルとなっています． 今回は学生や若手技術者向けの記事ですので，改良版ではなくオリジナルの積雪・融雪モデルについて解説しています．
タンクモデルや積雪・融雪モデルについては以下の書籍でも紹介されているのでよかったらご覧ください．
(7) 丸山利輔，石川重雄，大槻恭一，高瀬恵次，永井明博，田中丸治哉，駒村正治，赤江剛夫，堀野治彦，三野徹，武田育郎，金木亮一，渡辺紹裕：地域環境水文学，朝倉書店，2014，pp.67-70． (10) 田中丸治哉，大槻恭一，近森秀高，諸泉利嗣：地域環境水文学，朝倉書店，p.131，2019．
それでは，各ステップについて説明していきます．
2.1．流域のいくつかの地帯に分割
はじめに流域を標高別にN個の地帯に分割し，各地帯の面積\( A_{i} \)（i=1〜N）を求めます． Nは地帯分割数と呼ばれ，4程度とすることが多いようです． 菅原さんも4程度で十分みたいなことを言ってました．
例えば，標高0〜1,000mの範囲にある流域を4つの地帯に分けるとすると以下のようになります．
| 地帯番号 | 標高(m) | 地帯代表標高(m) ※ |
|---|---|---|
| 1 | 0〜250 | 125 |
| 2 | 250〜500 | 375 |
| 3 | 500〜750 | 625 |
| 4 | 750〜1000 | 875 |
※各地帯を代表する標高をここでは地帯代表標高と呼び，各地帯の最高標高と最低標高の平均値とします．
詳細は省略しますが，GIS(地理情報システム)と国土地理院のDEMを使えば各地帯の面積を求められます． (要望があれば記事にしたいと思います．)
2.2．地帯別の日平均気温を推定
次に地帯別の日平均気温(以下，地帯代表気温)を求めます． 地帯\( i \)の地帯代表気温\( T_{i} \)（℃）は，気象観測所の標高\( H_{0} \)（m），気象観測所の日平均気温\( T_{0} \)（℃），地帯代表標高\( H_{i} \)（m），および気温逓減率\( \alpha \)（℃/100m）から求められます．
なお，気温逓減率は-0.6 (℃/100m)※を使うことが多いです．
※標高が100m上がれば気温が0.6℃下がる．
例えば，上図のオレンジ色のプロット(標高200m，日平均気温20℃)が気象観測所だとします． 青色のプロットが地帯\( i \)の代表標高(500m)とすると，地帯iの地帯代表気温\( T_{i} \)は
2.3．降雨と降雪の判別および積雪量の計算
世間では，降水と降雨という言葉は混同されがちですが，正確には異なります． 『降水』は空から降ってくる水のことであり，それが液体なら『降雨』，固体なら『降雪』です． なので，降水量とは降雨量と降雪量の合計になります．
菅原さんの積雪・融雪モデルでは，降雨と降雪の判別を地帯ごとに行います． 判別には地帯代表気温を用い，気温が0℃以上なら降雨であると判断し，0℃未満なら降雪と判断します．
なお，標高が高い地帯に気象観測所がない場合は観測所の雨量に割増係数をかけて降水量を推定するそうです． (どれくらいの割増係数がいいかは記事の最後の参考文献などをご覧ください．)
降雨の場合はそのままタンクモデルへの入力になりますが，降雪の場合は各地帯ごとに積雪として貯めておきます． 計算上は積雪量をカウントするタンクを設け，降雪の場合はタンクに降雪量を追加していきます．
2.4．融雪量の計算
地帯iの日融雪量\( SM_{i} \)(mm/d)は次式で表されます．
ただし，
第1項が気温による雪の融解，第2項が雨による雪の融解を表しています． 式中の80という数値は，氷の融解熱80 (cal/g)に由来しています． 気温融雪率は6(mm/d/℃)が採用されることが多いようです．
積雪があれば積雪量から融雪量を引きます． ただし，積雪量は0 mm未満にはなりません．
2.5．流域全体での降雨量，降雪量，積雪量，融雪量を計算
流域全体での降雨量\( R \)(mm/d)，降雪量\( SF \)(mm/d)，積雪量\( SC \)(mm)，融雪量\( SM \)(mm/d)を求めます．
ただし，
3．使用データ
今回使用したデータは下表に示す10地点の気象データです． 期間は2010/1/1〜2019/12/31の10年間です．
| 気象観測所名 | 上田 | 佐久 | 立科 | 鹿教湯 | 東御 | 軽井沢 | 北相木 | 野辺山 | 白樺湖 | 乙女湖 |
|---|---|---|---|---|---|---|---|---|---|---|
| 標高(m) | 502 | 683 | 715 | 721 | 958 | 999 | 1185 | 1350 | 1440 | 1465 |
| 気温データ | 〇 | 〇 | 〇 | × | 〇 | 〇 | × | 〇 | × | × |
| 降水量データ | 〇 | 〇 | 〇 | 〇 | 〇 | 〇 | 〇 | 〇 | 〇 | 〇 |
気象データは気象庁のホームページからダウンロードしました． 気象観測所の標高は，「地域気象観測所一覧」から読み取りました．
なお，降雪量，積雪量の計算結果が妥当であるかを確認するため，同じく2010/1/1〜2019/12/31の生田観測所の河川流量(流出量)を水文水質データベースからダウンロードしました．
気象庁：地域気象観測所一覧 【Excel VBA】気象庁のデータをスクレイピング 【水文水質データベース】【使い方】流量データを自動で簡単にダウンロードしてみよう！【スクレイピング】
4．具体的な計算手順
生田流域に積雪・融雪モデルを適用してみます．
4.1．流域のいくつかの地帯に分割
生田流域は，標高約400〜3,000mの範囲にあります．
下に示す地図のように流域を325mずつ8地帯に分割することにしました．
地帯面積はおよそ以下の通りです．
| 地帯番号 | 地帯面積 (km2) | 流域面積全体に占める割合(%) | 地帯代表標高(m) |
|---|---|---|---|
| 1 | 255 | 12.56 | 562.5 |
| 2 | 634 | 31.20 | 887.5 |
| 3 | 561 | 27.62 | 1212.5 |
| 4 | 367 | 18.04 | 1537.5 |
| 5 | 144 | 7.08 | 1862.5 |
| 6 | 60 | 2.96 | 2187.5 |
| 7 | 10 | 0.51 | 2512.5 |
| 8 | 1 | 0.03 | 2837.5 |
4.2．地帯別の日平均気温を推定
はじめに，気温を測定している6つの観測所(上田，佐久，立科，東御，軽井沢，野辺山)におけるそれぞれの日平均気温Tkのデータから，次式を使って地帯代表気温の一次近似値\( T_{i}(k) \)を推定しました． なお，k=1〜6は上田〜野辺山のそれぞれの観測所を表すこととします．
ただし，
地帯代表気温の一次近似値\( T_{i}(k) \)は観測所の数(6)だけ求められるので，これを平均して地帯代表気温\( T_{i} \)としました．
なお，気温逓減率\( \alpha \)は-0.6 (℃/100m)としました．
4.3．降雨と降雪の判別および積雪量の計算
地帯1〜8それぞれで日降水量を求めます．
地帯1〜4には気象観測所があり降水量を観測しています．
地帯1〜4では，各地帯の日降水量のデータを算術平均して，その地帯の日降水量としました．
上図に示すように地帯1〜4だけで流域面積の約9割を占めています． 地帯5〜8には降水量データがありませんが，地帯5〜8の降水量が全体の計算に与える影響は小さいと判断し，地帯4の降水量で代用することにしました．
また，同様の理由で割増係数も今回は使用していません．
4.3で求めた地帯代表気温を使い，降雨と降雪の閾値を0℃として降雨と降雪を判別しました．
降雪の場合は積雪量をカウントするタンクに降水量を追加していきました．
4.4．融雪量の計算
地帯iの融雪量は2.4で示した日融雪量の推定式で求めました． なお，気温融雪率は6 (mm/d/℃)としました．
4.5．流域全体での降雨量，降雪量，積雪量，融雪量を計算
2.5に示した式で計算しました． 地帯分割数は8です．地帯面積は4.1で示した値になります．
5．計算結果
春先に融雪している様子がよくわかりますね！
また，生田流域では1月1日よりも前から雪が降り始めるため，蒸発散比法(※)で1年の境目を1月1日とするのは問題があるということが確認できましたね．
8月とか10月とか雪の影響のなさそうな時期を選んだら良さそうですかね．
※【流出解析】シリーズ 第4回 参照
6．検証
今回の計算が妥当な結果になっているかを確認するため，6つの気温観測所の気温データ，降雪量を観測してる軽井沢の降雪量，ハイドログラフという3つの観点からデータを眺めてみたいと思います．
6.1．地帯代表気温の推定結果の妥当性
6つある気温観測所の1つを選び，他の5つの気温観測所のデータで気温を推定しました． 日平均気温の実測値(1つ選んだ観測所のデータ)と推定値(他の5つの観測所のデータから推定した気温の平均値)を比較したところ，おおむね一致しているため，気温の推定はだいたい問題ないと判断しました．
6.2．降雪量の推定結果の妥当性
気温に比べれば，予想通り推定精度はいまいちですが，概ね良いのではないでしょうか？
ただ，実測では降雪があるのに，モデルでは降雪量が0となっている点(グラフのx軸上のプロット)があることが気にはなりますが．．．
6.3．冬季〜春季におけるハイドログラフとハイエトグラフの比較
ハイエトグラフ(降雨量+融雪量の時系列グラフ)とハイドログラフ(流出量の時系列グラフ)を比較してみます．
グラフの波形は概ね一致しており，まあまあいい推定結果だと思います． もちろん，降雨と降雪の閾値や気温融雪率などを最適化したり，改良菅原雪モデルを使うなどすればもっと精度は良くなると思います． その他にも，降水量の割増係数についての検討も必要になるかもしれません．
（余談）最近は高解像度の衛星画像が無料で手に入りますから，衛星画像と比較して積雪をチェックしてみるのも面白いと思います．
7．計算用ワークシート
計算用のプログラムを記事の最後に載せています．
プログラムを正常に動かすために下図のような「生データ」と「結果」のワークシートを作成してください．
観測所の数や地帯分割数は自由に変えてもらって構いません．
(プログラム実行時に表示されるポップアップウィンドウで観測所の数や地帯分割数を聞かれます．)
ちなみに計算結果はこんな感じで表示されます．
8．おわりに
最後までお付き合いいただきありがとうございました．
積雪・融雪モデルの計算の具体的なイメージは掴めましたか？ 説明に関してわかりくい箇所などありましたらご指摘いただけると幸いです．
それでは，また次回の記事でお会いしましょう． 次回はタンクモデルの最適化についてです． お楽しみに！
今後も農業農村工学(水文学，かんがい排水，土壌物理，水理学)を中心に記事を執筆していきたいと思います． リクエスト等も受け付けておりますので，ご遠慮なく連絡ください． Twitterアカウント：エビぐんかん@6LxAi9GCOmRigUI メール：nnCreatorCircle@gmail.com
【流出解析】シリーズ 一覧
【流出解析】シリーズ 第1回 流域界の決定と入力データの用意 【流出解析】シリーズ 第2回 貯留関数法【前編】 【流出解析】シリーズ 第3回 貯留関数法【後編】 【流出解析】シリーズ 第4回 流域実蒸発散量の推定―蒸発散比法― 【流出解析】シリーズ 第5回 実蒸発散量の推定―補完法― 【流出解析】シリーズ 第6回 タンクモデル 【流出解析】シリーズ 第7回 積雪・融雪モデル ―降雪量・融雪量の推定― 【流出解析】シリーズ 第8回 タンクモデルの最適化 【流出解析】シリーズ 第9回 単位図法 【流出解析】シリーズ 第10回 流出量・流出高・比流量の違い 【流出解析】シリーズ 第11回 表面流モデル【定常降雨：前編】 【流出解析】シリーズ 第12回 表面流モデル【定常降雨：後編】
(余談) 参考文献(5)の論文では，3つの気象観測所がある流域において，各観測所が代表する地区を高度別に4地帯に分けていました． つまり，流域を3×4=12分割しています．ただし，12に分割された地帯全てが異なる標高の範囲にあるわけではなく，ある地帯の標高の範囲と別の地帯の標高の範囲は被っていることもあるようです．
実は菅原さんの本(参考文献(1))では，4地帯に分割するだけで十分と書いてありました． 僕も今回の計算は8地帯に分割しましたが，やはり4地帯で十分だと感じました． おもしろいものですね．
検索用Key Words： 降雪，積雪，融雪，菅原正巳，タンクモデル，高地，標高，高度，気温，気温融雪率，水文学，流出解析，流出モデル，長短期両用モデル，短期流出解析，長期流出解析，水循環，気温日数法，積雪融雪解析
引用・参考文献
(1) 菅原正巳：水文学講座7 流出解析法，共立出版株式会社，p.99-105，1972． (2) 菅原正巳：水文学講座 別巻 続・流出解析法，共立出版株式会社，1979． (3) 安藤義久，知久岳史，安池宏之，1968．山地流域の積雪・融雪モデルを含む日単位の長期流出解析．水文・水資源学会．1-1，69-74． (4) 佐藤晃三，倉島栄一，松山幹男，1989．4種類の積雪・融雪変換モデルの推算精度の比較．農土論集．141，45-52． (5) 菅原正巳，1965．長良川流域の水収支について．水利科学．44，26-40． (6) 水津重雄，2002．広域に適用可能な融雪・積雪水量モデル．日本雪氷学会誌 雪氷．64-6，617-630． (7) 丸山利輔，石川重雄，大槻恭一，高瀬恵次，永井明博，田中丸治哉，駒村正治，赤江剛夫，堀野治彦，三野徹，武田育郎，金木亮一，渡辺紹裕：地域環境水文学，朝倉書店，2014，pp.67-70． (8) 菅原正巳，1985．タンク・モデル ―河川の流量を雨量から算出する一つのモデルについて―．地学雑誌．94-4，209-221． (9) 安富裕二，1983．算術平均法による面積平均雨量の推定誤差．天気．30-7，319-325． (10) 田中丸治哉，大槻恭一，近森秀高，諸泉利嗣：地域環境水文学，朝倉書店，p.131，2019．
9．計算プログラム
ご自由にお使いください． VBAの標準モジュールに下記のプログラムをコピペしてF5を押せばプログラムが起動します． 表示されるポップアップウィンドウに従って，流域面積，降雨と降雪を判別する気温の閾値，気象観測所の数，流域の最高標高，流域の最低標高，地帯分割数を入力すれば，計算が実行されます．
VBAの使い方はこちら↓
【VBA】Excelでプログラミングできる環境を構築し簡単なプログラムを組んでみる【プログラミング】 VBA の使い方、開発環境を整える
Option Explicit
Function DataInputBox(ByRef InputVariable As Variant, ByVal MessageVariable As String)
'インプットボックスを使って入力値を取得する
Dim ans As String 'InputBoxの入力値を受け取る
ans = Application.InputBox(MessageVariable & "を入力してください．")
If ans = False Then
InputVariable = -999 'キャンセルされたことを示すフラグ
Exit Function
End If
Do Until Val(ans) >= 1 And IsNumeric(ans) = True
ans = Application.InputBox("1以上の数値を入力してください．" & vbCrLf & MessageVariable & "を入力してください．")
If ans = False Then
InputVariable = -999 'キャンセルされたことを示すフラグ
Exit Function
End If
Loop
InputVariable = ans
End Function
Sub SnowFall_SnowMelt_Estimation()
'プログラムの内容：菅原正巳先生の積雪融雪モデルを使って気温と降水量データから，降雨量，降雪量，積雪量，融雪量を求める
'作成年月日：20201203
'最終更新日時：20201208
'入力データ：日平均気温(℃)，日降水量(mm/d)，観測所の標高(m)
'画面更新の停止
Application.ScreenUpdating = False
'計算開始メッセージ
Debug.Print "計算開始 "; Now()
'エラー防止
On Error GoTo ErrorFlag
'************変数の定義************
Dim i As Long, j As Integer, k As Integer, Cnt As Integer 'カウント用
Dim SR As Long, ER As Long '開始行，終了行
Dim SheetRaw As String '生データのワークシート
Dim SheetResult As String '計算結果を出力するワークシート
Dim StaNum As Integer '気象観測所の数
Dim StaTemperature() As Double '気温(℃)
Dim StaPrecipitation() As Double '降水量(mm/d)
Dim StaElevation() As Double '気象観測所の標高(m)
Dim StaZoneNo() As Double '気象観測所が属する地帯番号
Dim TempFlag() As Double '気温データフラグ(0：気温データなし，1：気温データあり)
Dim TempStaNum '気温を測定している観測所の数
Dim ZoneElevation() As Double '地帯代表標高(m)
Dim ZoneArea() As Double '地帯面積(km2)
Dim ZoneSnowTank() As Double '地帯別積雪タンク(mm)
Dim ZonePrecipitation() As Double '地帯平均降水量(mm)
Dim ZoneTemperature() As Double '地帯代表気温(℃)
Dim ZoneSnowMelt() As Double '地帯融雪量(mm)
Dim ArealRain As Double '流域全体での日降雨量(m3/d)
Dim TotalSnowFall As Double '流域全体での日降雪量(m3/d)
Dim TotalSnowCover As Double '流域全体での積雪量(m3)
Dim TotalSnowMelt As Double '流域全体での日融雪量(m3)
Dim m As Double '気温融雪率(mm/℃/d)
Dim A As Double '流域面積(km2)
Dim Threshold As Double '降雨と降雪を判別する気温の閾値(℃)
Dim DivisionNum As Integer '地帯数
Dim MaximumElevation As Double '流域の最高標高(m)
Dim MinimumElevation As Double '流域の最低標高(m)
Dim HeighestZoneSta As Integer '一番標高が高い気象観測所がある地番号
Dim DecRate As Double '標高の変化に対する気温低下率(℃/100m)
'************変数の初期化************
'ワークシート名
SheetRaw = "生データ"
SheetResult = "結果"
'物理定数など
TempStaNum = 0 '気温を測定している観測所の数
A = 2031 '流域面積(km2)
m = 6 '気温融雪率(mm/℃/d)
DecRate = 0.6 '標高の変化に対する気温低下率(℃/100m)
Threshold = 0 '降雨と降雪を判別する気温の閾値(℃)
'流域面積を入力
Call DataInputBox(A, "流域面積(km2)")
If A < 0 Then
MsgBox ("キャンセルされました．計算を終了します．")
Exit Sub
End If
'降雨と降雪を判別する気温の閾値を入力
Threshold = Application.InputBox("降雨と降雪を判別する気温の閾値(℃)を入力してください．")
'気象観測所の数を入力
Call DataInputBox(StaNum, "気象観測所の地点数")
If StaNum < 0 Then
MsgBox ("キャンセルされました．計算を終了します．")
Exit Sub
End If
'流域の最高標高を入力
Call DataInputBox(MaximumElevation, "流域の最高標高(m)")
If MaximumElevation < 0 Then
MsgBox ("キャンセルされました．計算を終了します．")
Exit Sub
End If
'流域の最低標高を入力
Call DataInputBox(MinimumElevation, "流域の最低標高(m)")
If MinimumElevation < 0 Then
MsgBox ("キャンセルされました．計算を終了します．")
Exit Sub
End If
'地帯分割数を入力
Call DataInputBox(DivisionNum, "地帯分割数(流域の分割数)")
If DivisionNum < 0 Then
MsgBox ("キャンセルされました．計算を終了します．")
Exit Sub
End If
'配列の再定義
ReDim StaTemperature(1 To StaNum), StaPrecipitation(1 To StaNum), StaElevation(1 To StaNum)
ReDim StaZoneNo(1 To StaNum), TempFlag(1 To StaNum)
ReDim ZoneElevation(1 To DivisionNum), ZoneArea(1 To DivisionNum), ZoneSnowTank(1 To DivisionNum)
ReDim ZonePrecipitation(1 To DivisionNum), ZoneTemperature(1 To DivisionNum), ZoneSnowMelt(1 To DivisionNum)
'観測所のパラメータを取得
For j = 1 To StaNum
StaElevation(j) = Worksheets(SheetRaw).Cells(3, j + 1) '標高の取得
StaZoneNo(j) = Worksheets(SheetRaw).Cells(4, j + 1 + StaNum) '気象観測所が属する地帯番号
TempFlag(j) = Worksheets(SheetRaw).Cells(4, j + 1) '気温フラグの取得
'気温を測定している観測所の数
If TempFlag(j) = 1 Then
TempStaNum = TempStaNum + 1
End If
Next j
'一番標高が高い気象観測所がある地番号
HeighestZoneSta = WorksheetFunction.Max(Worksheets(SheetRaw).Range( _
Worksheets(SheetRaw).Cells(4, StaNum + 2), Worksheets(SheetRaw).Cells(4, StaNum * 2 + 1)))
Debug.Print "流域面積(km2)="; A
Debug.Print "降雨と降雪を判別する気温の閾値(℃)="; Threshold
Debug.Print "気象観測所の数="; StaNum
Debug.Print "流域の最高標高(m)="; MaximumElevation
Debug.Print "流域の最低標高(m)="; MinimumElevation
Debug.Print "地帯数(流域の分割数)="; DivisionNum
Debug.Print "気温を測定している気象観測所の数="; TempStaNum
Debug.Print "一番標高が高い気象観測所がある地番号="; HeighestZoneSta
For j = 1 To DivisionNum
ZoneArea(j) = Worksheets(SheetRaw).Cells(3, StaNum * 2 + 2 + j) '地帯面積を取得
ZoneSnowTank(j) = 0 '積雪タンクをはじめは空と仮定
Next j
'************ラベル作成************
Worksheets(SheetResult).Cells(1, 1) = "地帯番号"
Worksheets(SheetResult).Cells(2, 1) = "地帯代表標高 (m)"
Worksheets(SheetResult).Cells(3, 1) = "地帯面積 (km2)"
For j = 1 To DivisionNum
Worksheets(SheetResult).Cells(4, j + 1) = "地帯代表気温(℃)"
If j = 1 Then
ZoneElevation(j) = (MaximumElevation - MinimumElevation) / 2 / DivisionNum + MinimumElevation
Else
ZoneElevation(j) = (MaximumElevation - MinimumElevation) / DivisionNum + ZoneElevation(j - 1)
End If
Worksheets(SheetResult).Cells(1, j + 1) = j
Worksheets(SheetResult).Cells(2, j + 1) = ZoneElevation(j)
Worksheets(SheetResult).Cells(3, j + 1) = ZoneArea(j)
Next j
Worksheets(SheetResult).Cells(4, 10) = "流域全体での日降雨量(m3/d)"
Worksheets(SheetResult).Cells(4, 11) = "流域全体での日降雪量(m3/d)"
Worksheets(SheetResult).Cells(4, 12) = "流域全体での積雪量(m3)"
Worksheets(SheetResult).Cells(4, 13) = "流域全体での日融雪量(m3/d)"
Worksheets(SheetResult).Cells(4, 14) = "流域全体での日降雨量+日融雪量(m3/d)"
Worksheets(SheetResult).Cells(4, 15) = "流域全体での日降雨量(mm/d)"
Worksheets(SheetResult).Cells(4, 16) = "流域全体での日降雪量(mm/d)"
Worksheets(SheetResult).Cells(4, 17) = "流域全体での積雪量(mm)"
Worksheets(SheetResult).Cells(4, 18) = "流域全体での日融雪量(mm/d)"
Worksheets(SheetResult).Cells(4, 19) = "流域全体での日降雨量+日融雪量(mm/d)"
'************日ごとの降雨量，降雪量，積雪量，融雪量を求める************
i = 5
Do
'日付を出力
Worksheets(SheetResult).Cells(i, 1).NumberFormatLocal = "yyyym/d"
Worksheets(SheetResult).Cells(i, 1) = Worksheets(SheetRaw).Cells(i, 1)
'i日目の入力データ(気温と降水量)を取得
For j = 1 To StaNum
StaTemperature(j) = Worksheets(SheetRaw).Cells(i, j + 1)
StaPrecipitation(j) = Worksheets(SheetRaw).Cells(i, j + 1 + StaNum)
Next j
'************地帯代表気温を求める************
For j = 1 To DivisionNum
ZoneTemperature(j) = 0
For k = 1 To StaNum
If TempFlag(k) = 1 Then
ZoneTemperature(j) = ZoneTemperature(j) - (ZoneElevation(j) - StaElevation(k)) / 100 * DecRate + StaTemperature(k)
End If
Next k
ZoneTemperature(j) = ZoneTemperature(j) / TempStaNum
Worksheets(SheetResult).Cells(i, j + 1) = ZoneTemperature(j)
Next j
'初期化
ArealRain = 0 'i日目の流域全体での日降雨量を0として初期化
TotalSnowCover = 0 'i日目の流域全体での積雪量を0として初期化
TotalSnowFall = 0 'i日目の流域全体での日降雪量を0として初期化
TotalSnowMelt = 0 'i日目の流域全体での日融雪量を0として初期化
For j = 1 To 8
ZonePrecipitation(j) = 0 'i日目の地帯降水量を0として初期化
ZoneSnowMelt(j) = 0 'i日目の地帯融雪量を0として初期化
Next j
'************地帯降水量を求める************
For j = 1 To DivisionNum
Cnt = 0
'雨量計がある地帯では，地帯降水量(mm/d)を各観測点の算術平均とする
If j <= HeighestZoneSta Then
For k = 1 To StaNum
If j = StaZoneNo(k) Then
ZonePrecipitation(j) = ZonePrecipitation(j) + StaPrecipitation(k)
Cnt = Cnt + 1
End If
Next k
ZonePrecipitation(j) = ZonePrecipitation(j) / Cnt
'雨量計がない地帯は地帯4の雨量で代用
ElseIf j > HeighestZoneSta Then
For k = 1 To StaNum
If StaZoneNo(k) = HeighestZoneSta Then
ZonePrecipitation(j) = ZonePrecipitation(j) + StaPrecipitation(k)
Cnt = Cnt + 1
End If
Next k
ZonePrecipitation(j) = ZonePrecipitation(j) / Cnt
End If
Next j
'************流域全体での日降雨量，積雪量，流域全体での日降雪量を求める************
For j = 1 To DivisionNum
'地帯代表気温が0℃以上のとき降水を降雨とし，流域全体での降水量(地帯降雨量(mm)×地帯面積(km2)の合計)を求める
If ZoneTemperature(j) >= Threshold Then
ArealRain = ArealRain + ZonePrecipitation(j) / 1000 * (ZoneArea(j) * 1000000)
'地帯代表気温が0℃未満のとき降水を降雪とし，各地帯の降雪タンクSnowTankに貯める
ElseIf ZoneTemperature(j) < Threshold Then
ZoneSnowTank(j) = ZoneSnowTank(j) + ZonePrecipitation(j)
TotalSnowFall = TotalSnowFall + ZonePrecipitation(j) / 1000 * (ZoneArea(j) * 1000000)
End If
Next j
'************流域全体での日融雪量，融雪量を差し引いた地帯積雪量・流域積雪量を求める************
For j = 1 To DivisionNum
'気温0℃より高いなら融雪すると仮定
If ZoneTemperature(j) > 0 Then
'融雪量を計算
ZoneSnowMelt(j) = m * ZoneTemperature(j) + ZonePrecipitation(j) * ZoneTemperature(j) / 80
'融雪量の上限は積雪量とする
If ZoneSnowMelt(j) > ZoneSnowTank(j) Then
ZoneSnowMelt(j) = ZoneSnowTank(j)
ZoneSnowTank(j) = 0
Else
ZoneSnowTank(j) = ZoneSnowTank(j) - ZoneSnowMelt(j)
End If
'流域全体での融雪量(地帯融雪量(mm)×地帯面積(km2)の合計)
TotalSnowMelt = TotalSnowMelt + ZoneSnowMelt(j) / 1000 * (ZoneArea(j) * 1000000)
End If
'i日目の流域全体での積雪量
TotalSnowCover = TotalSnowCover + ZoneSnowTank(j) / 1000 * (ZoneArea(j) * 1000000)
Next j
Worksheets(SheetResult).Cells(i, 10) = ArealRain
Worksheets(SheetResult).Cells(i, 11) = TotalSnowFall
Worksheets(SheetResult).Cells(i, 12) = TotalSnowCover
Worksheets(SheetResult).Cells(i, 13) = TotalSnowMelt
Worksheets(SheetResult).Cells(i, 14) = ArealRain + TotalSnowMelt
Worksheets(SheetResult).Cells(i, 15) = ArealRain / (A * 1000000) * 1000
Worksheets(SheetResult).Cells(i, 16) = TotalSnowFall / (A * 1000000) * 1000
Worksheets(SheetResult).Cells(i, 17) = TotalSnowCover / (A * 1000000) * 1000
Worksheets(SheetResult).Cells(i, 18) = TotalSnowMelt / (A * 1000000) * 1000
Worksheets(SheetResult).Cells(i, 19) = (ArealRain + TotalSnowMelt) / (A * 1000000) * 1000
i = i + 1
Loop Until Worksheets(SheetRaw).Cells(i, 1) = ""
'画面更新の再開
Application.ScreenUpdating = True
'計算終了メッセージ
Debug.Print "計算終了 "; Now()
MsgBox ("計算が終了しました．")
Exit Sub
'エラー防止
ErrorFlag:
MsgBox ("エラーが発生しました")
End Sub
