(Library) 断熱図作成のための熱力学関数集
作成者：辻野智紀
概要
エマグラムなどで必要とされる熱力学変数間の変換関数を集めたライブラリです. たとえば, ラジオゾンデで得られた気圧, 温度, 湿度から相当温位を計算するといったことができます.
(Fortran 90 で作成していた自作ライブラリを ruby に転用したものです.)
なお, <URL:http://ruby.gfd-dennou.org/products/met/> にはすでに堀之内さん作成の同様の処理をするライブラリがあります.
ソースファイル
- メインライブラリ 1 : thermo_function.rb. (機能を追加しました : 2014/05/26)
- メインライブラリ 2 : thermo_advanced_function.rb. (新規に追加しました : 2014/05/26)
- 定数ライブラリ 1 : phys_const.rb.
- 定数ライブラリ 2 : thermo_const.rb.
- 補助ライブラリ 1 : algebra.rb.
- 補助ライブラリ 2 : statistics.rb.
以上を tar でまとめたもの + αThermo_Function.tar.gz
使い方
- 上記 3 ファイルを同一ディレクトリもしくは, パスの通っているディレクトリに格納します.
require する場合は,
require "thermo_function" include Thermo_Function
でモジュールを読み込めます. ただし, 関数を使う際は必ず, narray を require し, NMath を使える状態にしておいてください.
require したのち, 関数を使用する場合は,
Thermo_Function.???
として使用して下さい.
依存モジュール
- NMath
関数一覧 (Thermo_Function)
- 一応, ソースファイルにもありますが, rdoc でうまく処理できない書き方をしている可能性があるので, 以下にまとめました.
- 引数が複数段にわたってある場合は, 上から第 1 引数, 第 2 引数という順になります.
- この順番を間違えると, とんでもない結果を生み出す可能性があります.
- また, 引数の単位は相対湿度 (単位 : %) を除いてすべて MKS 単位系で処理されますので, 気圧を hPa で代入するとか, 温度を degree Celsius で代入すると悲惨なことになります. (ありがちな間違いとして, 混合比を g/kg で代入する)
- 関数名の中で, "_2_" という形のものは, "hoge_2_fuga" としたとき, 変数 hoge を引数として, 結果 fuga を返すという決まりで命名しています.
- 引数が複数ある場合, その引数の与える順番も hoge と対応させてあります.
関数, 引数, 処理内容
tetens, 温度(K), tetens の実験式を用いて飽和水蒸気圧を計算する.(氷飽和を考慮)
es_Bolton, 温度(K), Bolton(1980) の手法を用いて飽和水蒸気圧を計算する.
es_TD, 水蒸気圧(Pa), es_Bolton を用いて飽和水蒸気圧の逆算から露点温度を計算する.
lh, 温度(K), キルヒホッフの式から得られる潜熱の温度依存性を考慮した潜熱評価式.
eP_2_qv, 水蒸気圧(Pa), 水蒸気圧と気圧から水蒸気混合比を計算する.
気圧(Pa)
qvP_2_e, 混合比(kg/kg), 混合比と気圧から水蒸気圧を計算する.
気圧(Pa)
theta_dry, 温度(K), 乾燥大気の温位を計算する.
気圧(Pa)
theta_moist, 温度(K), 湿潤大気における温位を計算する.
気圧(Pa)
混合比(kg/kg)
thetaP_2_T, 温位(K), 温位と気圧から温度を計算する(乾燥大気を仮定).
気圧(Pa)
thetaT_2_P, 温位(K), 温位と温度から温度を計算する(乾燥大気を仮定).
温度(K)
tqvP_2_TLCL, 温度(K), 温度, 混合比, 気圧から持ち上げ凝結(LCL)高度での温度を計算する.
混合比(kg/kg)
気圧(Pa)
thetae_Bolton, 温度(K), Bolton(1980) の手法を用いて相当温位を計算する.
混合比(kg/kg)
気圧(Pa)
thetaes_Bolton, 温度(K), Bolton(1980) の手法を用いて飽和相当温位を計算する.
気圧(Pa)
tqvP_2_thetae, 温度(K), 相当温位の定義式を用いて計算する.
混合比(kg/kg)
気圧(Pa)
tqvP_2_thetaes, 温度(K), 飽和相当温位の定義式を用いて計算する.
気圧(Pa)
rhT_2_e, 相対湿度(%), 相対湿度と温度から水蒸気圧を計算する.
温度(K)
eT_2_RH, 水蒸気圧(Pa), 水蒸気圧と温度から相対湿度を計算する.
温度(K)
rhTP_2_qv, 相対湿度(%), 相対湿度と温度, 気圧から混合比を計算する.
温度(K)
気圧(Pa)
qvTP_2_RH, 混合比(kg/kg), 混合比と温度, 気圧から相対湿度を計算する.
温度(K)
気圧(Pa)
exner_func_dry, 気圧(Pa), 乾燥大気のエクスナー関数を計算する.
hypsometric_form, 参照気圧(Pa), 気圧の高度補正. 海面更正式に基づく.
参照高度(m) 「参照」が基準となる高度での値, 「更正」が計算したい高度での値.
参照温度(K)
更正高度(m)
moist_laps_temp, 参照気圧(Pa), 参照気圧高度の値を用いて計算気圧高度での
参照温度(K) 湿潤断熱減率での温度を計算する.
計算高度(Pa)
関数一覧 (Thermo_Advanced_Function)
- 引数に "A" という記載があるものは鉛直 1 次元の配列として, 「地表面に近い方」を若い配列要素に格納して与える.
- 引数が [] でくくられているものはオプション引数である.
cape, 基準高度(m), 基準高度からパーセルを持ち上げたときの
高度(m;A) 対流有効位置エネルギーを計算する.
温度(K;A) 中立浮力高度の計算は
圧力(K;A) LFC から上向きに伸ばす方法 (1) と,
水蒸気混合比(kg/kg;A) 上空から下向きに伸ばす方法 (2) が選べる.
[計算方法 = 1]
[未定義値 = 999.0]
cin, 基準高度(m), 基準高度からの対流抑制エネルギーを計算する.
高度(m;A)
温度(K;A)
圧力(K;A)
水蒸気混合比(kg/kg;A)
[未定義値 = 999.0]
precip_water, 圧力(Pa;A), 可降水量を計算する.
水蒸気混合比(kg/kg;A)
[未定義値 = 999.0]
t_LNB, 基準高度(m), 基準高度からパーセルを持ち上げたときの
高度(m;A) 中立浮力高度における環境場の温度を計算する.
温度(K;A) 中立浮力高度の計算は
圧力(K;A) LFC から上向きに伸ばす方法 (1) と,
水蒸気混合比(kg/kg;A) 上空から下向きに伸ばす方法 (2) が選べる.
[計算方法 = 1]
[未定義値 = 999.0]
t_LFC, 基準高度(m), 基準高度からパーセルを持ち上げたときの
高度(m;A) 自由対流高度における環境場の温度を計算する.
温度(K;A)
圧力(K;A)
水蒸気混合比(kg/kg;A)
[未定義値 = 999.0]
z_LCL, 基準高度(m), 基準高度からパーセルを持ち上げたときの
高度(m;A) 持ち上げ凝結高度を計算する.
温度(K;A)
圧力(K;A)
水蒸気混合比(kg/kg;A)
[未定義値 = 999.0]
z_LFC, 基準高度(m), 基準高度からパーセルを持ち上げたときの
高度(m;A) 自由対流高度を計算する.
温度(K;A)
圧力(K;A)
水蒸気混合比(kg/kg;A)
[未定義値 = 999.0]
z_LNB, 基準高度(m), 基準高度からパーセルを持ち上げたときの
高度(m;A) 持ち上げ凝結高度を計算する.
温度(K;A) 中立浮力高度の計算は
圧力(K;A) LFC から上向きに伸ばす方法 (1) と,
水蒸気混合比(kg/kg;A) 上空から下向きに伸ばす方法 (2) が選べる.
[計算方法 = 1]
[未定義値 = 999.0]
使用例
- 実際のラジオゾンデ観測データから作成した Skew-T ダイアグラム.
- 青実線が温度プロファイル, 青破線が露点温度のプロファイル.
- 矢印は各高度における水平風 (紙面上向きが北).
(作図には DCL 使用)
サンプルスクリプト
- 上記の tar ファイルを解凍すると, 中に NetCDF 化されたゾンデデータが入っています.
- 同じく同梱の sample.rb を実行すると, ゾンデデータの鉛直プロファイルを用いて, 温位 (黒), 相当温位 (赤), 飽和相当温位 (緑) を計算し, グラフにすることができます.
- [注意]: 実行の際は, narray, netcdf, dcl が動作する環境で実行して下さい.
以下, sample.rb の中身:
require "numru/netcdf"
require "numru/dcl"
include NumRu
include NMath
require "thermo_function"
include Thermo_Function
require "thermo_advanced_function"
include Thermo_Advanced_Function
z_ref = 1000.0 # Reference height (using CAPE)
# ゾンデデータファイルオープン
file = NetCDF.open('F2011052814S7100069.nc')
temp = file.var('temp').get
pres = file.var('pres').get
rh = file.var('RH').get
z = file.var('z').get
# 描画用配列
qv = NArray.float(temp.total)
pt = NArray.float(temp.total)
ept = NArray.float(temp.total)
sept = NArray.float(temp.total)
nx = temp.total-2000
snd_linex = NArray.float(3)
snd_liney = NArray.float(3)
# 変数計算
for i in 0..nx
if temp[i] == 999.0 then # 欠損値処理
pt[i] = 999.0
ept[i] = 999.0
sept[i] = 999.0
else
# ここで関数を呼び出す. MKS 単位になっていることに注意.
qv[i] = Thermo_Function.rhTP_2_qv( rh[i], temp[i], pres[i] )
pt[i] = Thermo_Function.theta_dry( temp[i], pres[i] )
ept[i] = Thermo_Function.thetae_Bolton( temp[i], qv[i], pres[i] )
sept[i] = Thermo_Function.thetaes_Bolton( temp[i], pres[i] )
end
end
zlcl=Thermo_Advanced_Function.z_LCL( z_ref, z, temp, pres, qv )
zlfc=Thermo_Advanced_Function.z_LFC( z_ref, z, temp, pres, qv )
zlnb=Thermo_Advanced_Function.z_LNB( z_ref, z, temp, pres, qv, opt=2 )
pwv=Thermo_Advanced_Function.precip_water( pres, qv )
capev=Thermo_Advanced_Function.cape( pres, z, qv, temp, z_ref, opt=2, error=999.0 )
cinv=Thermo_Advanced_Function.cin( pres, z, qv, temp, z_ref, error=999.0 )
print "---------------------------------\n"
printf("%8s", "CAPE = ")
printf("%10f", capev)
printf("%8s\n", " [J/kg].")
printf("%8s", "CIN = ")
printf("%10f", cinv)
printf("%8s\n", " [J/kg].")
printf("%8s", "PW = ")
printf("%10f", pwv)
printf("%8s\n", " [mm].")
printf("%8s", "z_ref = ")
printf("%10f", z_ref)
printf("%8s\n", " [m].")
printf("%8s", "LCL = ")
printf("%10f", zlcl)
printf("%8s\n", " [m].")
printf("%8s", "LFC = ")
printf("%10f", zlfc)
printf("%8s\n", " [m].")
printf("%8s", "LNB = ")
printf("%10f", zlnb)
printf("%8s\n", " [m].")
print "---------------------------------\n"
iz_snd = Statistics.interpo_search_1d( z, z_ref )
snd_linex[0] = Statistics.interpolation_1d( z[iz_snd], z[iz_snd+1], ept[iz_snd], ept[iz_snd+1], z_ref )
iz_snd = Statistics.interpo_search_1d( z, zlfc )
snd_linex[1] = Statistics.interpolation_1d( z[iz_snd], z[iz_snd+1], sept[iz_snd], sept[iz_snd+1], zlfc )
iz_snd = Statistics.interpo_search_1d( z, zlnb )
snd_linex[2] = Statistics.interpolation_1d( z[iz_snd], z[iz_snd+1], sept[iz_snd], sept[iz_snd+1], zlnb )
snd_liney[0..2] = [z_ref, zlfc, zlnb]
# 以下 DCL 側の処理
DCL.gllset( 'lmiss', 'true' )
DCL.sgiset('IFONT', 2 )
puts 'device number [1,2,3,4 ?]'
num = gets
DCL.gropn(num)
DCL.grfrm
DCL.ussttl( 'temperature', 'K', 'altitude', 'm' )
DCL.grsvpt( 0.3, 0.75, 0.2, 0.8)
DCL.usspnt( pt[0..nx], z[0..nx] )
DCL.usspnt( ept[0..nx], z[0..nx] )
DCL.usspnt( sept[0..nx], z[0..nx] )
DCL.usspnt( snd_linex[0..1], snd_liney[0..1] )
DCL.usspnt( snd_linex[1..2], snd_liney[1..2] )
DCL.uspfit
DCL.grstrf
DCL.usdaxs
DCL.uulinz( pt[0..nx], z[0..nx], 1, 12 )
DCL.uulinz( ept[0..nx], z[0..nx], 1, 22 )
DCL.uulinz( sept[0..nx], z[0..nx], 1, 32 )
DCL.uumrkz( snd_linex[0..0], snd_liney[0..0], 3, 6, 0.02 )
DCL.uumrkz( snd_linex[1..1], snd_liney[1..1], 10, 6, 0.02 )
DCL.uumrkz( snd_linex[2..2], snd_liney[2..2], 16, 6, 0.02 )
DCL.grcls
結果
標準出力
--------------------------------- CAPE = 517.673395 [J/kg]. CIN = -9.807296 [J/kg]. PW = 67.571137 [mm]. z_ref = 1000.000000 [m]. LCL = 1103.846784 [m]. LFC = 1509.845477 [m]. LNB = 9746.783566 [m]. ---------------------------------
描画出力
- 黒・赤・緑実線はそれぞれ温位・相当温位・飽和相当温位を表す.
- 六花・黒丸・五芒星はそれぞれ CAPE 計算に用いた基準高度・自由対流高度・中立浮力高度を表す.
キーワード:[断熱図]
参照: