エマグラムを描いてみようと思った。まあ実際には skew-T log p diagram か温位エマグラムにするのだろうが、いずれにしても目盛り代わりの乾燥断熱線と湿潤断熱線が必要になる。
乾燥断熱線は温位の等値線、湿潤断熱線は相当温位の等値線であるから、要は相当温位の（温度の関数としての）逆関数が必要になる。
まあしょせん単調増加関数なので二分法で充分高速に計算できた。
ポカミスをしていないかチェックのため、ワイオミング大学のサイトに http://weather.uwyo.edu/upperair/sounding.html にならって
等温位線を -60 ... 180℃、等相当温位線を -50 ... 130℃ について計算してみると次のようになる。
図と見合わせても問題ない。
なお相当温位の算法は配信資料に関する技術情報第371号 https://www.data.jma.go.jp/suishin/jyouhou/pdf/371.pdf による。
※2024-01-31修正：コメント指摘をうけて動作を再確認したところ、ept_bolton() の正常動作範囲外の引数を与えて不適切な結果を得ているところがあったので修正しています。
#include <stdio.h>
#include <math.h>
// 気温t[K], 気圧p[hPa]に対して温位 [K] を与える
double
potemp(double t, double p)
{
double th = t * pow(1.0e3 / p, 0.2854);
return th;
}
// 温位t[K], 気圧p[hPa]に対して気温 [K] を与える
double
inv_potemp(double th, double p)
{
double t = th * pow(p / 1.0e3, 0.2854);
return t;
}
// 気温t[K], 相対湿度rh[%], 気圧p[hPa]に対して温位 [K] を与える
// Bolton(1980) による。
// https://doi.org/10.1175/1520-0493(1980)108%3C1046:TCOEPT%3E2.0.CO;2
// 動作範囲の注意:
// Bolton が精度を検証しているのは -35℃ .. 35℃ まで。かなり高温まで使えるが、
// 96℃ 100%1atmで 2.6e50 を返し、98℃までに爆発する（infやnanを返す）
double
ept_bolton(double t, double rh, double p)
{
double es = 6.112 * exp(17.67 * (t - 273.15) / (t - 29.65));
double e = es * rh * 0.01;
double ke = log(e/6.112) / 17.67;
double td = 0.05 * (-5463.0 + 593.0 * ke) / (-1.0 + ke);
double tlcl = 1.0 / (1.0 / (td - 56.0) + log(t / td) / 800.0) + 56.0;
double x = 0.622 * e / (p - e);
double thdl = t * pow(1.0e3 / (p - e), 0.2854) * pow(t / tlcl, 0.28 * x);
double ept = thdl * exp((3036.0 / tlcl - 1.78) * x * (1.0 + 0.448 * x));
return ept;
}
// 相当温位 ept[K] の（気温の関数としての）逆関数。
// 気圧は p[hPa] で与える。相対湿度は 100% とする。
// 第一推定値 t1 からステップ tstep で振って二分法で求める。
double
inv_ept_core(double ept, double p, double t1, double tstep)
{
const double rh = 100.0;
const double TSTEP_MIN = 0.1;
double ept1, ept2;
// 上側の関数値。
// 上側になっていなければ tstep だけ上がったところでやりなおし
SWEEP_UP:
ept1 = ept_bolton(t1, rh, p);
if (ept1 < ept) {
t1 += tstep;
goto SWEEP_UP;
} else if (isnan(ept1)) {
// 高温すぎるばあいは NaN を返す
return ept1;
}
// 下側の関数値。
// 適当なステップ tstep だけ下がったところで関数値を試算して、
// 下側になっていなければ t1-tstep を上側としてやりなおし
SWEEP_DOWN:
ept2 = ept_bolton(t1-tstep, rh, p);
if (ept2 > ept) {
t1 -= tstep;
ept1 = ept2;
goto SWEEP_DOWN;
}
if (tstep < TSTEP_MIN) {
// 収束していれば ept1, ept2 で線形補間して返す
return t1-((ept1-ept)/(ept1-ept2))*tstep;
} else {
// 収束していなければステップを半分にしてやりなおし
return inv_ept_core(ept, p, t1, 0.5*tstep);
}
}
// 相当温位 ept の（気温の関数としての）逆関数。
// 気圧は p [hPa] で与える。相対湿度は 100% とする。
double
inv_ept(double ept, double p)
{
// 第一推定値として、過小な値が望ましく、273.15 K を与える。
// ステップは 32K とする。float64 的に切りの良い数字であり、
// 32*6 = 96℃までの範囲がアルゴリズム的に破綻せず動くことを狙ったもの。
return inv_ept_core(ept, p, 273.15, 32.0);
}
double plevs[] = { 1000.0, 925.0, 850.0, 700.0, 500.0,
300.0, 250.0, 200.0, 150.0, 100.0 };
unsigned nplevs = sizeof(plevs)/sizeof(double);
int
main(int argc, char **argv)
{
// 相当温位近似式の適用限界のチェック
for (double tc = 70.0; tc < 200.0; tc+=1.0) {
double t = tc + 273.15;
double ept = ept_bolton(t, 100.0, 1013.25);
printf("tc=%g t=%g ept=%g\n", tc, t, ept);
if (isnan(ept)) break;
}
// 温位の逆関数のテスト
for (double th_c = -60.0; th_c <= 180.0; th_c += 20.0) {
printf("= TH %5.1f %5.1f\n", th_c, th_c+273.15);
for (unsigned iz = 0; iz < nplevs; iz++) {
double t = inv_potemp(th_c+273.15, plevs[iz]);
printf("%6.1f\t%5.1f\t%5.1f\n", plevs[iz], t, t-273.15);
}
}
// 相当温位の逆関数のテスト
for (double ept_c = -50.0; ept_c <= 130.0; ept_c += 20.0) {
printf("= EPT %5.1f %5.1f\n", ept_c, ept_c+273.15);
for (unsigned iz = 0; iz < nplevs; iz++) {
double t = inv_ept(ept_c+273.15, plevs[iz]);
printf("%6.1f\t%5.1f\t%5.1f\n", plevs[iz], t, t-273.15);
}
}
return 0;
}
実行結果
tc=70 t=343.15 ept=3586.2
tc=71 t=344.15 ept=4234.44
tc=72 t=345.15 ept=5076.69
tc=73 t=346.15 ept=6192.74
tc=74 t=347.15 ept=7704.84
tc=75 t=348.15 ept=9805.92
tc=76 t=349.15 ept=12810.9
tc=77 t=350.15 ept=17253.7
tc=78 t=351.15 ept=24079.9
tc=79 t=352.15 ept=35049.5
tc=80 t=353.15 ept=53634.4
tc=81 t=354.15 ept=87162.8
tc=82 t=355.15 ept=152393
tc=83 t=356.15 ept=291501
tc=84 t=357.15 ept=623702
tc=85 t=358.15 ept=1.53781e+06
tc=86 t=359.15 ept=4.55143e+06
tc=87 t=360.15 ept=1.71242e+07
tc=88 t=361.15 ept=8.89702e+07
tc=89 t=362.15 ept=7.22202e+08
tc=90 t=363.15 ept=1.10896e+10
tc=91 t=364.15 ept=4.39366e+11
tc=92 t=365.15 ept=7.65658e+13
tc=93 t=366.15 ept=1.57206e+17
tc=94 t=367.15 ept=2.81278e+22
tc=95 t=368.15 ept=4.47968e+31
tc=96 t=369.15 ept=2.63753e+50
tc=97 t=370.15 ept=1.60054e+100
tc=98 t=371.15 ept=inf
tc=99 t=372.15 ept=inf
tc=100 t=373.15 ept=-nan
= TH -60.0 213.1
1000.0 213.1 -60.0
925.0 208.5 -64.7
850.0 203.5 -69.7
700.0 192.5 -80.6
500.0 174.9 -98.3
300.0 151.2 -122.0
250.0 143.5 -129.6
200.0 134.6 -138.5
150.0 124.0 -149.1
100.0 110.5 -162.7
= TH -40.0 233.1
1000.0 233.1 -40.0
925.0 228.0 -45.1
850.0 222.6 -50.6
700.0 210.6 -62.6
500.0 191.3 -81.8
300.0 165.4 -107.8
250.0 157.0 -116.2
200.0 147.3 -125.9
150.0 135.7 -137.5
100.0 120.8 -152.3
= TH -20.0 253.1
1000.0 253.1 -20.0
925.0 247.6 -25.6
850.0 241.7 -31.5
700.0 228.6 -44.5
500.0 207.7 -65.4
300.0 179.5 -93.6
250.0 170.4 -102.7
200.0 159.9 -113.2
150.0 147.3 -125.8
100.0 131.2 -141.9
= TH 0.0 273.1
1000.0 273.1 0.0
925.0 267.1 -6.0
850.0 260.8 -12.4
700.0 246.7 -26.4
500.0 224.1 -49.0
300.0 193.7 -79.4
250.0 183.9 -89.3
200.0 172.6 -100.6
150.0 158.9 -114.2
100.0 141.6 -131.6
= TH 20.0 293.1
1000.0 293.1 20.0
925.0 286.7 13.5
850.0 279.9 6.7
700.0 264.8 -8.4
500.0 240.5 -32.6
300.0 207.9 -65.2
250.0 197.4 -75.8
200.0 185.2 -88.0
150.0 170.6 -102.6
100.0 151.9 -121.2
= TH 40.0 313.1
1000.0 313.1 40.0
925.0 306.3 33.1
850.0 299.0 25.8
700.0 282.8 9.7
500.0 256.9 -16.2
300.0 222.1 -51.1
250.0 210.8 -62.3
200.0 197.8 -75.3
150.0 182.2 -90.9
100.0 162.3 -110.8
= TH 60.0 333.1
1000.0 333.1 60.0
925.0 325.8 52.7
850.0 318.1 44.9
700.0 300.9 27.8
500.0 273.4 0.2
300.0 236.3 -36.9
250.0 224.3 -48.9
200.0 210.5 -62.7
150.0 193.9 -79.3
100.0 172.7 -100.5
= TH 80.0 353.1
1000.0 353.1 80.0
925.0 345.4 72.2
850.0 337.1 64.0
700.0 319.0 45.8
500.0 289.8 16.6
300.0 250.5 -22.7
250.0 237.8 -35.4
200.0 223.1 -50.1
150.0 205.5 -67.6
100.0 183.0 -90.1
= TH 100.0 373.1
1000.0 373.1 100.0
925.0 364.9 91.8
850.0 356.2 83.1
700.0 337.0 63.9
500.0 306.2 33.0
300.0 264.6 -8.5
250.0 251.2 -21.9
200.0 235.7 -37.4
150.0 217.1 -56.0
100.0 193.4 -79.7
= TH 120.0 393.1
1000.0 393.1 120.0
925.0 384.5 111.3
850.0 375.3 102.2
700.0 355.1 81.9
500.0 322.6 49.4
300.0 278.8 5.7
250.0 264.7 -8.5
200.0 248.4 -24.8
150.0 228.8 -44.4
100.0 203.8 -69.4
= TH 140.0 413.1
1000.0 413.1 140.0
925.0 404.1 130.9
850.0 394.4 121.3
700.0 373.2 100.0
500.0 339.0 65.8
300.0 293.0 19.9
250.0 278.2 5.0
200.0 261.0 -12.2
150.0 240.4 -32.7
100.0 214.1 -59.0
= TH 160.0 433.1
1000.0 433.1 160.0
925.0 423.6 150.5
850.0 413.5 140.4
700.0 391.2 118.1
500.0 355.4 82.3
300.0 307.2 34.0
250.0 291.6 18.5
200.0 273.6 0.5
150.0 252.1 -21.1
100.0 224.5 -48.6
= TH 180.0 453.1
1000.0 453.1 180.0
925.0 443.2 170.0
850.0 432.6 159.5
700.0 409.3 136.1
500.0 371.8 98.7
300.0 321.4 48.2
250.0 305.1 31.9
200.0 286.3 13.1
150.0 263.7 -9.5
100.0 234.9 -38.3
= EPT -50.0 223.1
1000.0 223.0 -50.1
925.0 218.2 -55.0
850.0 213.0 -60.2
700.0 201.5 -71.6
500.0 183.1 -90.1
300.0 158.3 -114.9
250.0 150.2 -122.9
200.0 141.0 -132.2
150.0 129.9 -143.3
100.0 115.7 -157.5
= EPT -30.0 243.1
1000.0 242.3 -30.8
925.0 237.3 -35.9
850.0 231.8 -41.3
700.0 219.5 -53.6
500.0 199.5 -73.7
300.0 172.4 -100.7
250.0 163.7 -109.5
200.0 153.6 -119.6
150.0 141.5 -131.7
100.0 126.0 -147.1
= EPT -10.0 263.1
1000.0 259.5 -13.6
925.0 254.7 -18.4
850.0 249.4 -23.7
700.0 237.0 -36.1
500.0 215.8 -57.3
300.0 186.6 -86.5
250.0 177.2 -96.0
200.0 166.2 -106.9
150.0 153.1 -120.0
100.0 136.4 -136.8
= EPT 10.0 283.1
1000.0 272.9 -0.3
925.0 268.8 -4.4
850.0 264.1 -9.0
700.0 252.8 -20.4
500.0 231.8 -41.4
300.0 200.8 -72.4
250.0 190.6 -82.5
200.0 178.9 -94.3
150.0 164.8 -108.4
100.0 146.8 -126.4
= EPT 30.0 303.1
1000.0 282.6 9.4
925.0 279.1 5.9
850.0 275.1 2.0
700.0 265.4 -7.7
500.0 246.4 -26.8
300.0 214.9 -58.3
250.0 204.1 -69.1
200.0 191.5 -81.7
150.0 176.4 -96.7
100.0 157.1 -116.0
= EPT 50.0 323.1
1000.0 289.7 16.5
925.0 286.6 13.4
850.0 283.1 10.0
700.0 274.8 1.7
500.0 258.4 -14.7
300.0 228.5 -44.6
250.0 217.3 -55.8
200.0 204.1 -69.1
150.0 188.0 -85.1
100.0 167.5 -105.7
= EPT 70.0 343.1
1000.0 295.0 21.9
925.0 292.2 19.1
850.0 289.1 16.0
700.0 281.7 8.6
500.0 267.6 -5.6
300.0 241.0 -32.1
250.0 230.1 -43.1
200.0 216.5 -56.6
150.0 199.6 -73.5
100.0 177.9 -95.3
= EPT 90.0 363.1
1000.0 299.2 26.1
925.0 296.6 23.5
850.0 293.7 20.6
700.0 287.0 13.8
500.0 274.4 1.2
300.0 251.4 -21.8
250.0 241.5 -31.6
200.0 228.4 -44.7
150.0 211.2 -62.0
100.0 188.2 -84.9
= EPT 110.0 383.1
1000.0 302.7 29.5
925.0 300.2 27.0
850.0 297.5 24.3
700.0 291.1 18.0
500.0 279.6 6.5
300.0 259.4 -13.7
250.0 250.9 -22.3
200.0 239.1 -34.1
150.0 222.3 -50.8
100.0 198.6 -74.6
= EPT 130.0 403.1
1000.0 305.5 32.4
925.0 303.1 30.0
850.0 300.5 27.4
700.0 294.5 21.3
500.0 283.7 10.6
300.0 265.6 -7.6
250.0 258.1 -15.0
200.0 247.9 -25.3
150.0 232.6 -40.6
100.0 208.8 -64.4