はじめに
水文学，農業農村工学，気象学，かんがい排水，etc... 大学でこんな授業を受けたことのある大学生にとって，「ペンマン法」，「ペンマン式」，「蒸発散位」という言葉は一度は耳にしたことがあるのものではないでしょうか？
しかし， ・ 式の形が複雑 ・ 入力データをどう用意したらいいかわからない ・ 風速関数ってなんだよ... ・ 飽和水蒸気圧曲線の勾配？？ (゜Д゜)ハァ？？？
と意味不明なことが盛りだくさんで，理解するのを諦めた人も少なくないでしょう． 私もその一人でした．
だって，「ペンマン法」で検索しても堅苦しい論文しか出てこないんですもん(^ω^#)
というわけで，"ペンマン法をやってみたブログ第1号"として， ・ 学生目線で ・ 細かい話は置いといて ・ 計算の雰囲気が何となく伝わる感じの ペンマン式の計算方法の紹介記事にしたいと思います．
ペンマン法について
それでは，本題のペンマン式の登場です．
ここで，
ちなみに，↑の式は丸山先生・三野先生編の地域環境水文学(朝倉書店)という本から引っ張ってきたものです． この式の導出方法や各項の意味についてはこちらの記事↓で紹介しています． よかったらご覧ください．
【導出】ペンマン式の導出・証明に挑戦しよう！ https://nnhiroba.xsrv.jp/nnTechBlog/Article/Article?ID=20
$$ $$ さて，この記事の題目である「ペンマン法ってどうやって計算するの？」についてですが， 結論：この報文↓を読んだら出来ます！
三浦 健志，奥野 林太郎，1993．ペンマン式による蒸発散位計算方法の詳細，農業土木学会誌．(164)，157-163．
$$ $$
この報文の何がすごいってペンマン法の計算過程を順を追って丁寧に説明していることです． しかも教科書とか参考書よりもわかりやすいです． え？じゃあ，この記事書く意味ないんじゃないかって？？
まあ，意味があるかどうかはわかりませんが，私個人としては多少は意義があるのかなと思っています．具体的には， ・ 農業農村工学分野では，こういった報文をフランクに紹介している記事すら，ほぼ存在していない． ・ この報文に載っているプログラムは有益だけど，今はほとんど使われない言語(Basic)で書かれている．
というわけで，今の学生に馴染みがあるエクセルと農業農村工学と相性がいいVBAでペンマン法の計算をしてみまた！
入力データ
まずは入力データの用意ですね． 報文に載ってるデータを丸パクリしました．以上．
実はペンマン法に必要なデータはすべて気象庁のホームページからダウンロードできるので，皆さんもよかったらチャレンジしてみてください．
計算結果
はい，結果です． え？計算過程の説明は...？？
省略します！ というのも，報文の説明がめっちゃわかりやすいので，それを読むのが手っ取り早いです笑 本記事ではVBAで書いたプログラムで計算していますが，エクセルに標準搭載されている関数のみで計算可能なのでご安心ください． ※ペンマン法の計算用プログラム(VBA)は記事の最後に載せておきます．
蒸発散位は2〜8 mm/dくらいになりました． 右下のグラフ(各入力データと蒸発散位の相関係数)を見ると，気温，湿度，日照時間なんかが蒸発散位に効いてるようですね．
他の地域への適用
さて，報文のデータだけで計算するのもつまんないので，他の地域でも蒸発散位を求めてみることにします． 場所は彦根の気象台にしました． データは2019年のものを使用しました．
彦根を選んだ理由は風速計の設置高度が公表されているから．
※もちろん他にも風速計の設置高度が公表されている気象台はあります．
気象台の所在や風速計の設置高度は気象庁の資料「地域気象観測所一覧」にまとめられています． 彦根の気象台については35ページに書かれています．
気象庁：地域気象観測所一覧 https://www.jma.go.jp/jma/kishou/know/amedas/ame_master.pdf
彦根の蒸発散位(2019)
蒸発散位は1,120 mm/yearくらいでした． ちなみに降水量は1,400 mm/yearくらいでした． 彦根は気温とか日照時間が効いてるみたいですね．
計算用ファイルおよびソースコード
計算用ファイル：0_ペンマン法.xlsm
計算に使用したプログラム↓です． ご自由にお使いください．
Sub PenmmanPotentialEvapotranspiration()
'ペンマン法により蒸発散位を計算するプログラム(2019/12/19作成)
'圧力の単位mbarはhPaに等しいため，本プログラムではhPaと表記することにした
'参考文献1：三浦健志, 奥野林太郎, 1993. ペンマン式による蒸発散位の計算方法の詳細. 農業土木学会論文集. 164, 157-163.
'参考文献2：丸山利輔, 三野徹 編, 地域環境水文学, 朝倉書店, p.18, pp.36-37(2014)
'参考文献3：文字信貴, 平野高司, 高見晋一, 堀江武, 桜谷哲夫 編, 農学・生態学のための気象環境学, 丸善株式会社, p.8(2004)
'有効数字については特別考慮していないので実際に利用する場合は注意すること
Application.ScreenUpdating = False
'変数の定義
Dim i As Integer,j as Integer 'カウント用
Dim lat As Double '緯度(°)
Dim latR As Double '緯度(rad)
Dim H As Double '風速計地上高度(m)
Dim a As Double 'アルベド
Dim TD As Integer '積算日数(Total Days) (日)
Dim dd As Double '地球太陽間の距離(天文単位)
Dim Dcl As Double '赤緯(Declination) (°)
Dim DclR As Double '赤緯(Declination)をラジアンで表現
Dim w As Double '日没時の時角(rad)
Dim N as Double '可照時間(h)
Dim Esat As Double '飽和水蒸気圧(hPa)
Dim Ea As Double '水蒸気圧(hPa)
Dim dlt As Double '飽和水蒸気圧曲線の勾配(hPa/℃)
Dim l As Double '蒸発潜熱(MJ/kg)
Dim u2 As Double '地上2mにおける風速(m/s)
Dim fu As Double '風速関数
Dim ET0 As Double '蒸発散位(mm/d)
Dim Gamma As Double '乾湿計定数(hPa/℃)
Dim Sigma As Double 'ステファン-ボルツマン定数(W/m2/K4)
'データ出力欄の作成
Range("H1") = "出力データ"
Range("H1:M4").Merge
Range("H1:M4").HorizontalAlignment = xlCenter
Cells(5, 8) = "赤緯(°)"
Cells(5, 9) = "大気圏外日射量(MJ/m2/d)"
Cells(5, 10) = "純放射量(MJ/m2/d)"
Cells(5, 11) = "蒸発散位ET0第1項 (mm/d)"
Cells(5, 12) = "蒸発散位ET0第2項 (mm/d)"
Cells(5, 13) = "蒸発散位ET0 (mm/d)"
'地点パラメータの取得
lat = Cells(3, 3) '緯度(°)
latR = WorksheetFunction.Radians(Cells(3, 3)) '緯度(rad)
H = Cells(3, 4) '風速計地上高度(m)
a = Cells(3, 5) 'アルベド
'物理化学定数(参考文献1の値を使用)
Gamma = 0.66 '乾湿計定数(hPa/℃)
Sigma = 4.9 * 10 ^ (-9) 'ステファン-ボルツマン定数(W/m2/K4)
'ペンマン法の本体部分(参考文献1のpp.160-162)
i = 6
Do
'起算日(1/1)からの積算日数(日)
TD = Cells(i, 2) - DateSerial(Year(Cells(6, 2)), 1, 1) + 1
'太陽からの距離(天文単位)
dd = 1 + 0.01676 * Cos(WorksheetFunction.Radians(0.977 * (TD - 186)))
'赤緯(°)
Cells(i, 8) = 23.45 * Cos((TD - 173) * 0.966 * 3.1416 / 180)
'赤緯(rad)
DclR = WorksheetFunction.Radians(Cells(i, 8))
'日没時の時角(rad)
w = WorksheetFunction.Acos((-Tan(latR) * Tan(DclR)))
'可照時間(h)
N = 2 * WorksheetFunction.Degrees(w) / 15
'大気圏外日射量(MJ/m2)
Cells(i, 9) = 1.37 * 10 ^ (-3) / dd ^ 2 * 86400 / (4 * Atn(1)) _
* (w * Sin(latR) * Sin(DclR) + Sin(w) * Cos(latR) * Cos(DclR))
'飽和水蒸気圧(hPa)
Esat = 6.1078 * Exp(17.2694 * Cells(i, 3) / (Cells(i, 3) + 237.3))
'水蒸気圧の計算(hPa)
Ea = Esat * Cells(i, 4) / 100
'純放射量(MJ/m2)(参考文献1のステファン-ボルツマン定数の表現がちょっとおかしい)
Cells(i, 10) = (1 - a) * Cells(i, 9) * (0.18 + 0.55 * Cells(i, 6) / N) - Sigma * (Cells(i, 3) + 273.2) ^ 4 _
* (0.56 - 0.092 * 0.866 * Sqr(Ea)) * (0.1 + 0.9 * Cells(i, 6) / N)
'飽和水蒸気圧曲線の勾配(hPa/℃)
dlt = 0.4495 + 0.2721 * 10 ^ (-1) * Cells(i, 3) + 0.9873 * 10 ^ (-3) * Cells(i, 3) ^ 2 + 0.2907 * 10 ^ (-5) * Cells(i, 3) ^ 3 + 0.2538 * 10 ^ (-6) * Cells(i, 3) ^ 4
'蒸発潜熱(MJ/kg)
l = 2.5 - 0.0024 * Cells(i, 3)
'地上2m風速への換算と風速関数fu
u2 = Cells(i, 5) * Log(200) / Log(100 * H)
fu = 0.26 * (1 + 0.54 * u2)
'蒸発散位(mm/d)
Cells(i, 11) = dlt / (dlt + Gamma) * Cells(i, 10) / l
Cells(i, 12) = Gamma / (dlt + Gamma) * fu * (Esat - Ea)
Cells(i, 13) = Cells(i, 11) + Cells(i, 12)
i = i + 1
Loop Until Cells(i, 2) = ""
'合計と平均
Cells(i, 1) = "合計"
Cells(i + 1, 1) = "平均"
For j = 3 To 13
Cells(i, j) = WorksheetFunction.Sum(Range(Cells(6, j), Cells(i - 1, j)))
Cells(i + 1, j) = WorksheetFunction.Average(Range(Cells(6, j), Cells(i - 1, j)))
Next j
Application.ScreenUpdating = True
End Sub
おわりに
拙い記事ではありましたが，最後までお付き合い頂き，ありがとうございました． (実はこれが人生初ブログです．) 今後も農業農村工学(水文学，かんがい排水，土壌物理，水理学)を中心に記事を執筆していきたいと思います． リクエスト等も受け付けておりますので，ご遠慮なく連絡ください． Twitterアカウント：エビぐんかん@6LxAi9GCOmRigUI メール：nnCreatorCircle@gmail.com
なお，執筆にあたりデータのミス等には十分注意を払っておりますが，お気づきの点がございましたら，優しくご指摘いただきたく存じます．
本記事が農業農村工学に携わる学生の勉強を支える一助となれば幸いです．
