湿位涡编程出图_假相当位温（两个问题）

和清谈

最近在研究湿位涡编程出图，涉及到假相当位温。后期根据情况会接着更新，期待大家多多指导！
因为比较懒，所以想直接用GRADS出图，在此也谢谢论坛里各位大牛的分享。我编程有个臭习惯，喜欢把每个步骤都尽量要求做到看懂，看不懂就很抓狂。今天又遇到了难题，如果有知晓的，烦请指导指导。
【1】之前论坛里好多用的是NCEP资料，里面包含相对湿度rhum物理量，编写的程序基本如下：
（详见http://bbs.06climate.com/forum.php?mod=viewthread&tid=9884&highlight=%CA%AA%CE%BB%CE%D0）

*************   计算假相当位温   **********（里面注释是我自己摸索着增加的，有问题的地方用红色标注）
'set lev 1000'
**定义：气压prs，单位:mb=Pa）
'define prs=lev*100'
*'define ms=17.67*(air-273.15)/(air-29.65)'**这一项基本没用
**定义：饱和水气压es，单位:Pa（air-温度，单位:K）
**应用修正的Tetens公式计算：
'define es=100*(6.112*exp(17.67*(air-273.15)/(air-29.65)))'
**定义：饱和比湿qs
*'define qs=0.6212*es/(prs-0.378*es)'             **公式好像有些不对，貌似是这个，少了个系数？？
'define qs=0.62197*es/(prs-es)'
**定义：比湿qv（rhum是相对湿度，单位:%）
'define qv=rhum.2*qs/100'
**定义：水汽压e（单位:hPa）                          **这里单位为啥要变成hPa，是为了代入后面tlcl公式？?
'define e=(prs/100.)*qv/(0.62197+qv)+1e-10'    **为什么要加1e-10？？
**定义：凝结高度的绝对温度tlcl（单位:K）
'define tlcl=55.0+2840.0/(3.5*log(air)-log(e)-4.775)'
**定义：位温theta（单位:K）
'define theta=air*pow((100000./prs),(0.2854*(1.0-0.28*qv)))'*最大的问题？和找到的公式有出入
**定义：假相当位温these（单位:K）
'define thse=theta*exp(((3376./tlcl)-2.54)*qv*(1.0+0.81*qv))'*最大的问题？和找到的公式有出入
'd thse'
公式如下图：里面r指的是混合比mixing ratio，也就是水汽质量与干空气质量之比，而程序中qv明显就是比湿，即水汽质量与湿空气总质量之比。我就不懂了，虽然英语比较差，但mixing ratio怎么翻译也应该是混合比的意思吧？论坛内基本上所有计算假相当位温的程序都是这样的，是大家一直没注意到这个问题，还是我自己疏漏了哪一项？请大家指正。

【2】不知什么时候开始，Ncep资料已经无法下载（自动回复：请不要使用迅雷等下载工具，点我查看下载帮助）（自动回复：请不要使用迅雷等下载工具，点我查看下载帮助）了，迫于无奈，我只能去下载EC数据了，其中数据和Ncep有所不同，里面直接有Specific humidity，也就是比湿数据。由此推出水汽压e和混合比mr，再代入到以上假相当位温公式中去。              
'reinit'
'sdfopen E:\ECWMF\201801.nc'
i=1
while(i<=977)
'set lon 50 160'
'set lat 10 70'
'set t 'i'''set lev 1000'
'define prs=lev'
'define e=(prs*q)/(0.622+0.378*q)+1e-10'   
**水汽压e(单位:hPa)
'define mr=q/(1-q)'
**混合比mr(单位:Kg/Kg)
'define tlcl=55.0+2840.0/(3.5*log(t)-log(e)-4.805)'
'define theta=t*pow((1000./prs),(0.2854*(1.0-0.28*mr)))'
'define thse=theta*exp(((3376./tlcl)-2.54)*mr*(1.0+0.81*mr))'
'd thse'
i=i+1
endwhile
'disable fwrite'
'reinit'
运行结果如下：
       log()里面的温度t和水汽压e都不应该是小于零的数啊，为啥呢？
      而且如果将这些数据存到dat文件，最大只有12931kb，应该也是由问题的，但具体是什么就不清楚了。
      有谁能告知一下吗？

river

给你两个参考吧，实在没时间在仔细研究了
'reinit'
'sdfopen g:\2012\uwnd.2012.nc'
'sdfopen g:\2012\vwnd.2012.nc'
'sdfopen g:\2012\rhum.2012.nc'
'sdfopen g:\2012\shum.2012.nc'
'sdfopen g:\2012\omega.2012.nc'
'sdfopen g:\2012\air.2012.nc'

'set t 709 721'
*-----水平散度------
'set lev 1000 100'
'define div=hdivg(uwnd.1,vwnd.2)*1e5'
*---------水汽通量--------
'set lev 1000 300'
'define qu=(uwnd.1*shum.4*1e3/9.8)'
'define qv=(vwnd.2*shum.4*1e3/9.8)'
*-----水汽通量散度------
'define qdiv=hdivg(qu,qv)*(1e5)'
*--------计算假相当位温----------
'define br=(0.1158*log(rhum.3)-1.5332787)*air.6-789.92816'                              
'define cr=789.92816*air.6'                                                         
'define tc2=(-br-sqrt(br*br-4*cr))/2.0'                                               
'define pc2=lev*pow(tc2/air.6,3.5005574)'                                                        
'define etc2=6.1078*pow(273.16/tc2,5.1347779)*exp(3148.0973*(tc2-273.16)/(126.06334*tc2))'   
'define wc2=0.622*etc2/(pc2-etc2)'                                                        
'define ed2=tc2*pow(1000.0/(pc2-etc2),0.28765924)'                                          
'define lc2=2500.79-2.3697*(tc2-273.16)'                                                  
'define eqt=ed2*exp(wc2*lc2/1.0048/tc2)'
*--------GRADS设置------------
'set grads off'
'set grid off'
'set csmooth on'
'set map 15 1 9'
'set parea 1 10.2 1.5 8'
'set xlopts 1 6 0.14'
'set ylopts 1 6 0.14'
'set font 0'

*---------输出图形1------------
*'enable print g:\2012\qdiv.gmf'
'set lon 110'
'set lat 41'
'set zlog on'
*'set xlint 4'
'set ylevs 1000 925 850 700 600 500 400 300'
'set xlabs 08|12|16|20|00|04|08|12|16|20|00|04|08|12|16|20|00|04|08'
'set gxout shaded'
'set cint 0.3'
'set cmax 0'
'd qdiv'
'run cbar.gs'
'set ylpos 0 r'
'set gxout contour'
'set cint 0.3'
*'set cthick 11'
*'set clskip 1'
'd qdiv'
'set gxout vector'
'set arrowhead 0.1'
'set arrscl 0.5 20'
'd qu;qv'

*'draw title  0000UTC 'i' JAN 2008'
*'writehz 0.4 4 气压（hPa） 1 2 1 1.3 0.6 90 11'
*'writehz 1 8.1 水汽通量散度 1 2.5 1 1.3 0.6 0 11'
*'writehz 0.7 1.0 6月26日 1 2 1 1.3 0.6 0 7'
*'writehz 2.8 1.0 6月27日 1 2 1 1.3 0.6 0 7'
*'writehz 5.9 1.0 6月28日 1 2 1 1.3 0.6 0 7'
*'writehz 9.1 1.0 6月29日 1 2 1 1.3 0.6 0 7'

*'print'
*'disable print'
'gxprint g:/2012/qdiv.png white'
'c'
*---------输出图形2------------
*'enable print g:\2012\jia.gmf'
'set grads off'
'set lon 110'
'set lat 41'
'set zlog on'
'set ylevs 1000 925 850 700 600 500 400 300'
'set xlabs 08|12|16|20|00|04|08|12|16|20|00|04|08|12|16|20|00|04|08'
'set gxout contour'
'set cint 3'
*'set cthick 6'
'set clskip 1'
'd eqt'
'set ylpos 0 r'
'set gxout contour'
'set cint 3'
*'set cthick 11'
*'set clskip 1'
'd eqt'

*'draw title  0000UTC 'i' JAN 2008'
*'writehz 0.4 4 气压（hPa） 1 2 1 1.3 0.6 90 7'
*'writehz 1 8.1 假相当位温 1 2.5 1 1.3 0.6 0 7'
*'writehz 0.7 1.0 6月26日 1 2 1 1.3 0.6 0 7'
*'writehz 2.8 1.0 6月27日 1 2 1 1.3 0.6 0 7'
*'writehz 5.9 1.0 6月28日 1 2 1 1.3 0.6 0 7'
'writehz 9.1 1.0 6月29日 1 2 1 1.3 0.6 0 7'

*'print'
*'disable print'
'gxprint g:/2012/jia.png white'
'c'
*---------输出图形3------------
*'enable print g:\2012\omega.gmf'
'set lev 1000 100'
'set grads off'
'set lon 110'
'set lat 41'
'set zlog on'
'set ylevs 1000 925 850 700 600 500 400 300 250 200 150 100'
'set xlabs 08|12|16|20|00|04|08|12|16|20|00|04|08|12|16|20|00|04|08'
'set gxout shaded'
'set cint 0.05'
'set cmax -0.15'
'd omega.5'
'run cbar.gs'
'set ylpos 0 r'
'set gxout contour'
'set cint 0.05'
*'set cthick 11'
*'set clskip 1'
'd omega.5'
*'draw title  0000UTC 'i' JAN 2008'
*'writehz 0.4 4 气压（hPa） 1 2 1 1.3 0.6 90 7'
*'writehz 1 8.1 垂直速度 1 2.5 1 1.3 0.6 0 7'
*'writehz 0.7 1.0 6月26日 1 2 1 1.3 0.6 0 7'
*'writehz 2.8 1.0 6月27日 1 2 1 1.3 0.6 0 7'
*'writehz 5.9 1.0 6月28日 1 2 1 1.3 0.6 0 7'
*'writehz 9.1 1.0 6月29日 1 2 1 1.3 0.6 0 7'
*'print'
*'disable print'
'gxprint g:/2012/omega.png white'
'c'
*---------输出图形4------------
*'enable print g:\2012\div.gmf'
'set grads off'
'set lon 110'
'set lat 41'
'set zlog on'
'set ylevs 1000 925 850 700 600 500 400 300 250 200 150 100'
*'set xlabs on'
*'set xaxis 0 72 4'
'set xlabs 08|12|16|20|00|04|08|12|16|20|00|04|08|12|16|20|00|04|08'
'set gxout shaded'
'set cint 0.5'
'set cmax 0'
'd div'
'run cbar.gs'
'set ylpos 0 r'
'set gxout contour'
'set cint 0.5'
*'set cthick 11'
*'set clskip 1'
'd div'
'set gxout vector'
*'set arrowhead 0.1'
*'set arrscl 0.5 20'
'd mag(uwnd.1,vwnd.2);omega.5*100'
*'draw title  0000UTC 'i' JAN 2008'
*'writehz 0.4 4 气压（hPa） 1 2 1 1.3 0.6 90 7'
*'writehz 1 8.1 散度 1 2.5 1 1.3 0.6 0 7'
*'writehz 0.7 1.0 6月26日 1 2 1 1.3 0.6 0 7'
*'writehz 2.8 1.0 6月27日 1 2 1 1.3 0.6 0 7'
*'writehz 5.9 1.0 6月28日 1 2 1 1.3 0.6 0 7'
*'writehz 9.1 1.0 6月29日 1 2 1 1.3 0.6 0 7'
*'print'
*'disable print'
'gxprint g:/2012/div.png white'
;

另一个去看这个帖子http://bbs.06climate.com/forum.php?mod=viewthread&tid=51560

和清谈

river 发表于 2019-1-23 22:58
给你两个参考吧，实在没时间在仔细研究了
'reinit'
'sdfopen g:\2012%uwnd.2012.nc'

您好，谢谢您的分享，但有一事不明，能告诉我这个程序中计算假相当位温的公式吗？

river

和清谈 发表于 2019-1-25 08:24
您好，谢谢您的分享，但有一事不明，能告诉我这个程序中计算假相当位温的公式吗？

啊，不记得了。好像是大气物理书上的