私にはＣ＋＋は難しい、関数を作ってもファイルからデータをstringで読み込んで、charに変換して、さらにdoubleに変換してやっと計算ができるみたいです。
ハードルが高すぎるのでサンプルプログラムを作るのはあきらめました。
SSIは５００hPaの飽和相当温位から８５０hPaの相当温位の差を求めればよいので作成しませんでした。
SSIは温度差ですが、飽和相当温位と相当温位の差はＣpをかければエネルギーの差になります。
単位系は普通の「物理」にある単位系にしました。
状態方程式PV＝ｎRTでRが窒素とヘリウムで値が変わらないようにしました。
気象では１ｋｇあたりの状態方程式を使いますのでRの値が違ってしまう不具合を避けました。
いろいろ書きたいのですが・・
まず、関数をみましょう。
不具合があっても責任は取りません。個人の責任で使用するのはかまいません。
template
<typename TYPE>
TYPE
Abs(TYPE a);
//絶対値を返す関数
double
e(double td);
//露点温度（ｔｄ）から水蒸気圧（e）を求める
double
Td(double rh,double t);
//相対湿度（）と気温（）から露点温度を求める
double
g(double hight);
//高さ（hight）毎の重力加速度（g）を求める
double
Tw(double temp,double td,double press);
//気温（temp）露点（td）気圧（press）から湿急温度（tw）を求める
double
pt(double hight,double temp);
//高さ（hight）と気温（temp）から温位（pt）を求める cf気温と気圧から求めるられる温位より正確です
double
ept(double hight,double temp,double td,double press);
//高さ（hight）,気温（temp）,露点温度（）と気圧から相当温位（ept）を求める cf普通求める相当温位より正確です
int
main(){
/*????
ゴメン
????*/
}
template
<typename TYPE>
TYPE
Abs(TYPE a)
{
return
a<0 ? -a : a;
}
double
e(double td){
static
double const e0=6.11;//０℃の水の水蒸気圧
static
double const LR0=5374.48884;// L/R0でLは昇華熱R0は水蒸気のガス定数 水蒸気は理想気体ではない
static
double const ek=2.791233047;
static
double const t0=273.15;//0℃の絶対温度
double
LR=LR0-ek*(td-20);//水は理想気体ではないのでL/Rに温度依存性がある
double
result=e0*exp(LR*(1/t0-1/(td+t0)));
return
result;
}
//水蒸気圧表からekを決定している
double
Td(double rh,double t)
{
double
t1,t2,tr;
double*
t3;
double
d,a;
t3=&tr;
cout<<"
rh : "<<rh<<" t : "<<t<<endl;
t1=t-100.0;//十分大きな差をとること
t2=t;
*t3=0.5*(t1+t2);
//
cout<<"*t : "<<*t3<<endl;
a=rh*e(t)/100;
d=e(*t3)-a;
//cout<<"d
: "<<d<<endl;
while(Abs(d)>0.0001){
if(
d<0 ){
t1=*t3;
}else{
t2=*t3;
}
*t3=0.5*(t1+t2);
d=e(*t3)-a;
//
cout<<"*t3 : "<<*t3<<" d : "<<d<<"
a : "<<a<<endl;
}
return
tr;
}
//水蒸気圧からTdを追い込んでいる この追い込みは別に関数化できると思うが・・；
double
g(double hight)
{
static
double const k=0.0000031;//一般気象学の表からkを求めた高度１０ｋｍ以上だと影響がでそう
double
result =9.807-k*hight;
return
result;
}
//ｇを固定すると２０〜３０ｋｍで不自然に温位が高くなる
double
Tw(double temp,double td,double press){
static
double const L=51012;
static
double const cp=29.093;
double
t1=temp;
double
t2=td;
double
*t3;
double
tr;
t3=&tr;
*t3=0.5*(t1+t2);
cout<<*t3<<endl;
double
d;
double
a=cp*(temp+273.15)+L*e(td)/press;
d=a-cp*(*t3+273.15)-L*e(*t3)/press;
cout<<"a
: "<<a<<endl;
cout<<"d
: "<<d<<endl;
while(Abs(d)>0.01){
//
for(int i=0;i<10;++i){
cout<<"d
: "<<d<<endl;
if(
d<0 ){
t1=*t3;
//
cout<<"t2 : "<<t2<<endl;
}else{
t2=*t3;
//
cout<<"t1 : "<<t1<<endl;
}
*t3=0.5*(t1+t2);
//cout<<"a
: "<<a<<endl;
d=a-cp*(*t3+273.15)-L*e(*t3)/press;
cout<<"*t3
: "<<*t3<<" L*e(*t3)/press : "<<
L*e(*t3)/press<<endl;
}
return
tr;
}
//
ＣpT+Le(Td)/P=ＣpTw+Le(Tw)/PからTwを追い込んでいる 追い込みはTdと同じ方法
double
pt(double hight,double temp){
static
double const cp=29.093;//N2:28.01 78.088% O2:32.00 20.949% Ar:39.94
0.93% 単原子のArまで考慮した
static
double const m=0.028957; //N2:28.01 78.088% O2:32.00 20.949%
Ar:39.94 0.93% 単原子のArまで考慮した
MKSに統一したほうがよいのだけど・・yoku
matigaeru
static
double const t0=273.15;//平成２３年理科年表P390
//cout<<"m:
"<<m<<" g: "<<g(hight)<<"
T: "<<temp+t0<< endl;
double
result=m*g(hight)*hight/cp + temp+t0;
return
result;
}
//エネルギー ＣpT+mgh をＣpで割っている普通の温位より物理的意味は明確
double
ept(double hight,double temp,double td,double press){
static
double const L=51012.0;//気象関係の値を採用
static
double const cp=29.093;
double
result=pt(hight,temp)+L*e(td)/press/cp;
return
result;
}
//エネルギー ＣpT+mgh+Le(td)/P をＣpで割っている普通の相当温位より物理的意味は明確
//普通の温位は1000Phaまで圧縮するエネルギーを無視するミスを犯している結果SSIがプラスでも不安定なんて変な観測結果になる。
１．重力加速度
一般気象学 P２３
表２ 諸物理量の各高度における値(米国標準大気モデル,
1976)
(tama は対象を３００００ｍ程度までにした）
高度（Km) 重力加速度（ｍｓ−２）
0 9.807
5
9.791
10
9.776
15
9.761
20
9.745
25
9.730
30
9.715
0mと２００００ｍを基準とする
ｇ（ｈ）＝９．８０７−ｋｈ
ｇ（０）＝９．８０７
ｇ（２００００）＝９．８０７−２００００ｋ＝９．７４５
結果 ｋ＝０．０００００３１
２．水蒸気圧
Eｓ＝E0Exp{（L/R)・(１/T0
ー１/T)}
大気科学講座２ P１４
（注 Rは理想気体のガス状数とは異なる）
Rは温度に依存する可能性があり
L/Rが温度に依存すると考える
E0＝6.11
hPa
T=ー20℃=257.15Kの場合 水蒸気圧1．25hPa
1.25=6.11Exp{L/R・(１/273.15-1/253.15)}
ln(1.25)-ln(6.11)=L/R・(１/273.15-1/253.15
)
ー1.586783222=ーL/R(-0.000289235)
L/R=5486.138163
T=20℃=297.15Kの場合 水蒸気圧23.39hPa
L/R=5374.488841
グラフにすると (温度は本来絶対温度で示すべきだった)
飽和水蒸気圧の値は 大気科学講座２ P１３から
LはTに依存するL=L(T)
T=20℃=297.15Kを基準として（L／R)の計算式（近似式）を作ると
L/R=5374.48884-2.791233047*(T-297.15)
Eｓ＝E0Exp{（L/R)・(１/T0
ー１/T)}
高層観測の蒸気圧の値は相対湿度から計算することになる
氷の飽和蒸気圧か水の飽和蒸気圧を使うで値がことなってしまう
軽井沢の地上データを見ると水の飽和蒸気圧を使っている。
天気相談所で確かめておいたほうがよさそうだが・・；
マッイイカ（^^
３．定圧比熱
一般気象学P52
定積比熱Cv＝717JK−１ｋｇ−１
定積比熱Cp＝1004JK−１ｋｇ−１
大気科学講座P17
定積比熱Cv＝716JK−１ｋｇ−１
定積比熱Cp＝1005JK−１ｋｇ−１
現代熱力学 イリヤ・プリゴジンP３３
T=２９８．１５K P=1atm
N2
Cp=29.12 Cv=20.74 (Jmol-1)
O2
Cp=29.36 Cv=20.95 (Jmol-1)
Ar
Cp=20.79 Cv=12.47 (Jmol-1)
CO2
Cp=37.11 Cv=28.46 (Jmol-1)
一般気象学P13
地表付近の大気組成（８０ｋｍまでほとんど変わらない）
N2 分子量 ２８．０１ 存在比 ７８．０８８
O2 ３２．００ ２０．９４９
Ar
３９．９４ ０．９３
CO2
４４．０１ ０．０３
CO2等を無視すると
N2,O2,Arで９９．９５７%
乾燥空気をN2,O2,Arと考えるとCp＝２９．０９３ Cv＝２０．７０７を採用
平均分子量は２８.９５７１２６６５ｇ/mole
これで計算すると
1kgのCpは１００４．６９
1kgのCvは７１５．０９
多分これが正しい値
４．融解熱
一般気象学P58
氷の融解 ３３４ ｊ／ｇ→ ６０１２ ｊ／mole
液体の蒸発
２５００ j／g→ ４５０００ j／mole
大気講座２ P12
液体の蒸発
２５０１ j／g→ ４５０１８ j／mole
氷を水蒸気にする
２８３４ j／g→ ５１０１２ j／mole
現代熱力学イリヤ・プリゴジン P１３１
P＝１０００００Pa＝１０００ｈPa
氷の融解 T＝２７３．１５
６００８ j/mole
水の蒸発 T＝３７３．１５ ４０６５６ j/mole
気象関係と水の蒸発が大きく違う
理科年表 平成２３年
氷の融解
６０１０ j／mole
飽和蒸気圧下
水の蒸発
４４０００ j／mole 25度 P=２３．７５Torr＝３１６５．８７５Pa＝３１．６５８７７５ｈPa
Pa３１．６６は２５℃の飽和蒸気圧３１．６９に近い（多分誤差範囲）
水の蒸発 ４０７００ j／mole
100度 飽和蒸気圧下
氷の融解熱は共通している
水の蒸発圧力や温度に左右されるようだ
根拠ははっきりしないが気象関係の教科書に従ったほうが無難と考える
欲しいのは上空の冷たい環境で水蒸気が氷になる時のエネルギー（潜熱）
潜熱は水蒸気が上空で氷になる潜熱L＝５１０１２ ｊ／moleを採用する
クラジュウス・クラペイロンの式（熱力学・統計力学 W．グライナー P６５）
平衡である条件
T液体 ＝ T気体
P液体 ＝ P気体
μ液体 ＝ μ気体
ギブスーデュエムの関係式（P59）から独立ではない
もし状態方程式が知られいればTとPからμを求めることができる
μ液体(P,T) ＝ μ気体(P,T)
(3.12)
(3.12)が成り立つように温度をdTだけ変えると、平衡を保つためには蒸気圧もある値dPだけ変える必要がある。
そうした変化では
dμ液体(P,T) ＝ dμ気体(P,T)
が成り立っていないといけないギブスーデュエムの関係式
SdT ー VdP ＋ Ndμ ＝ ０
を液体と気体について考えると
dμ液体(P,T) ＝ ー（S液体/N液体）dT ＋ （V液体/N液体）dP
dμ気体(P,T) ＝ ー（S気体/N気体）dT ＋ （V気体/N気体）dP
ｓ液体＝（S液体/N液体）
ｖ液体＝（V液体/N液体）
ｓ気体＝（S気体/N気体）
ｖ気体＝（V気体/N気体）
とすれば
dP(v液体 ー v気体）＝dT(s液体 ー s気体）
s液体 ー s気体
は液体同じ温度で全部気体した時のエントロピーの変化とみなせるから
s液体 ー s気体＝⊿Q’（液体→気体）／T
⊿Q’は液体を気体にするために必要な熱
dP/dT=⊿Q’（液体→気体）／T(v液体 ー v気体）〜ー⊿Q’（液体→気体）／T(v気体）
⊿Q’’を気体を液体にすると出てくる熱とするとー⊿Q’＝⊿Q’’だから
dP/dT〜⊿Q’’（液体→気体）／T(v気体）
1moleの水蒸気を考えると⊿Q’’は水にしたの熱量L
１moleの水蒸気はPV＝RvTに従ったとして
dP/dT=（L／T）（P／RvT）
dP/P=(L/Rv)dT/T2
lnP(T)=(L/Rv)(ー１／T)＋const
lnP(273.15=0℃）＝ｌｎ（6.11）＝(L/Rv)（−１／２７３．１５）＋const
lnP(T)=(L/Rv)(ー１／T) ＋ ｌｎ（6.11） ー (L/Rv)（−１／２７３．１５）
ln(P(T)/6.11） ＝ ー (L/Rv)(1／Tー ／２７３．１５）
P（T)＝6.11exp{ー (L/Rv)(1／Tー ／２７３．１５）}
水蒸気の飽和蒸気圧はeやesであらわされるようなので
e(T)=6.11exp{ー (L/Rv)(1／Tー ／２７３．１５）}
Rvは水蒸気の気体定数で大気科学講座２P１４でRv＝０．４６１５ Jｇ−１K−１
水の分子量は１８だからRv=8.307
Jmole−１K−１と気体定数(H23理科年表P３６７
)R＝8.314472に近い
Rvを使ってと大気科学講座２P１３の氷の蒸気圧を使ってLを求めて見よう
Td＝ー10 L＝51017 j/mole
Td＝ー20 L＝51133 j/mole
Td＝ー30 L＝51080 j/mole
Rvが一定である根拠もないしLも51012にかなり近い
大気講座２ P12
液体の蒸発
２５０１ j／g→ ４５０１８ j／mole
氷を水蒸気にする
２８３４ j／g→ ５１０１２ j／mole
上空の飽和蒸気圧に対応する露点温度はマイナス
下層の水蒸気が上空に行って氷になると仮定してL＝５１０１２を採用して良さそう