MATLAB/Octaveで海面フラックス計算
2022年6月 川合義美
MATLABまたはGNU Octaveで海面フラックスを計算します。COARE3の公開プログラムを使用します。GNU Octaveでも動きます。
論文等を執筆するのに使用した場合は参考文献やプログラムのダウンロード元を明記し謝意を表しましょう。参考文献やバージョン間の違いはご自分で調べてください。
以下のFTPサイトから、cor3_0 またはcor3_5 またはcor3_6をディレクトリごとダウンロードします。
ftp://ftp1.esrl.noaa.gov/BLO/Air-Sea/bulkalg/
1)coare3.0の場合
入力データの行列(例:dat)を作成し、関数cor3_0ah.mに渡し、出力データの行列(例:flx)を得ます。
>> flx = cor3_0ah(dat);
ファイルcor3_0ah.mの中で、風速計、気温計、湿度計の海面下からの高度(zu, zt, zq)、及びバルク水温の深度(ts_depth)を指定する必要があります(84-87行目)。
Skin layerの影響を考慮したい場合は jcool = 1(デフォルト), 考慮したくない場合は0にしてください。
Warm layerの影響を考慮したい場合は jwarm = 1(デフォルト), 考慮したくない場合は0にしてください。
入力データの行列は以下の通りです。全ての列は同じ長さのベクトルです。
1列目/ 日時:YYYYMMDDhhmmss の形で与える。例)2022年6月19日12時30分00秒(国際標準時)の場合は、20220619123000 とする(整数)。
2列目/ 風速 (m/s)
3列目/ サブスキン水温 ※使用しないので、0など適当なダミーデータで可
4列目/ 気温 (°C)
5列目/ 比湿 (g/kg) ※単位に注意
6列目/ 下向き短波放射 (W/m2)
7列目/ 下向き長波放射 (W/m2)
8列目/ 降水量 (mm/hr)
9列目/ 緯度 (°N)
10列目/ 経度 (°E)
11列目/ バルク水温 (°C)
出力データの3, 4, 6列目がそれぞれhsb (顕熱フラックス), hlb (潜熱フラックス), taub (風応力)です。それ以外の説明は割愛します(プログラムファイルを読んでください)。
2) coare3.5の場合
Warm layerを考慮しない場合には、
>> flx = coare35vn(u, zu, t, zt, rh, zq, P, ts, Rs, Rl, lat, zi, rain, cp, sigH);
ここで、u(風速:m/s), t(気温:°C), rh(相対湿度:%), ts(水温:°C), P(気圧:hPa), ts(水温:°C), Rs(下向き短波放射:W/m2), Rl(下向き長波放射:W/m2), rain(降水量:mm/hr), cp(卓越波高の位相速度:m/s), sigH(有義波高:m)はベクトルで与えます。Warm layerを考慮しないので、tsにバルク水温を与えた場合でもサブスキン水温と見なされます。
zu(風速計の高度), zt(気温計の高度), zq(湿度計の高度), lat(緯度), zi(境界層高度)はスカラーで与えます。
P, Rs, Rl, lat, ziのデータがない場合には、NaNを入力します。その場合、デフォルト値が使われます(P=1015 hPa, Rs=150 W/m2, Rl=370 W/m2, lat=45°N, zi=600 m)。
cp, sigHのデータがない場合には、NaNを入力します。この場合、ベクトルでもスカラーでもかまいません。
Skin layerの影響を考慮したい場合は jcool = 1(デフォルト), 考慮しない場合は0にしてください。
出力値の行列flxの中身は、前からusr (摩擦速度), tau (風応力), hsb (顕熱フラックス), hlb (潜熱フラックス), hbb (浮力フラックス), … です。coare3.0と少し異なるので注意しましょう。
Warm layerを考慮したい場合には、
>> flx = coare35vnWarm(yday, ur, zu, t, zt, rh, zq, P, ts, Rs, Rl, lat, lon, zi, rain , ts_depth);
yday: 1月1日からの日数(国際標準時で1月1日0時を基準とする)
ur: 相対風速(m/s) ※風速-海面流速、流速がわからない場合は風速で代用
ts: 深度ts_depthにおける水温(°C)
lon: 経度(°E)
ts_depth: 水温を測定した深度(m)
ここではlatとlonはベクトルで与えます。ts_depthはスカラーです。
!!! 注意 !!! そのままだと、grvという変数がないというエラーが出ます。ディレクトリcor3_6にあるgrv.mをcor3_5にコピーして置くと解決します。
3) coare3.6の場合
Warm layerを考慮しない場合には、
>> flx = coare36vn_zrf(ur, zu, t, zt, rh, zq, P, ts, Rs, Rl, lat, zi, rain, Ss, cp, sigH, zrf_u, zrf_t, zrf_q);
ここで、urは相対風速(風速-海面流速、m/s)、Ssは海面の実用塩分(psu)です。
zrf_u, zrf_t, zrf_qは、それぞれ風速、気温、湿度の基準高度(m)です。この高度における風速(Urf)、気温(Trf)、湿度(Qrf)の推定値が出力されます。必要ない場合は適当な値(10とか)を入れておけばいいでしょう。
zu, zt, zq, lat, zi, zrf_u, zrf_t, zrf_qはスカラー、それ以外はベクトルで与えます。cp, sigHはNaNでもベクトルで与える必要があります。Warm layerを考慮しないので、tsにバルク水温を与えた場合でもサブスキン水温と見なされます。
出力値flxの前から5列は、coare3.5と同じです。
Warm layerを考慮したい場合には、
>> flx = coare36vnWarm(yday, ur, zu, t, zt, rh, zq, P, ts, Rs, Rl, lat, lon, zi, rain , ts_depth, Ss, cp, sigH, zrf_u, zrf_t, zrf_q);
ここでは、zu, zt, zq, lat, lon, cp, sigH, zrf_u, zrf_t, zrf_qもベクトルで与える必要があるようなので、注意してください。ziとts_depthはスカラーでよいようです。