关于最近画假相当位温脚本的总结

lm8005507771

刚开始学用grads，总结一下关于画假相当位温的方法，最后一共写了三个脚本，跟大家分享一下！
我用的是欧洲中心全球中期数值天气预报产品，名字长这样echtxa85.032的奇怪数据，问了很多同学和师兄师姐，都没见过这种数据，反正是grib2的数据，就跟着网上grib2 的处理方法来处理了。

下面说一下三个方法的脚本

方法一：
用了师弟帮我从大气物理书里翻出来的公式，还贴心的帮我推倒了，感谢师弟！！感谢师弟！！感谢师弟！！！
脚本内容如下：


'reinit'
'open c:/test3/echtxa85.ctl '
'set lon 70 135'
'set lat 15 55'
'set mpdset cnworld'
'set grid on'
'define tk=tmp850mb' ;*某高度层的开氏温度K
'define tc=tk-273.16' ;*某高度层的摄氏温度C
'define prs=850.0'
******************************************************************************************
*水面es,tk开氏温度,tc摄氏温度
if(tk>273.16)
'define es=6.1078*exp(17.2693882*tc/(tk-35.86))'
endif
*冰面es,tk开氏温度,tc摄氏温度
if(tk<=273.16)
'define es=6.1078*exp(21.8745584*tc/(tk-7.66))'
endif
'define rs=(0.622*es)/(prs-es)';*rs为混合比
'define theta=tk*pow((1000.0/prs),(287.05/1004.07))';*theta为位温
'define thetase=theta*exp((rs*2501.0)/(1004.07*tk))';*thetase为假相当位温
******************************************************************************************
'set gxout contour'
'set cthick 5'
'd theta'
'printim c:\test2\theta1.png white'
'c'
'd thetase'
'printim c:\test2\thetase51.png white'
'c'
;



公式出处：
我也不知道这样推倒有问题没，我大气物理学的是渣... 我相信师弟！

1

（这个图为啥转不过来...大家凑合看吧）

2

3

4

5

6

最后的结果图如下：

5.0

方法二和方法三都是借用网上的公式来，具体参考了哪些我也记不太清了，也十分感谢大神们的分享！
主要参考内容来自以下：
http://bbs.06climate.com/forum.php?mod=viewthread&tid=45075
http://bbs.06climate.com/forum.p ... &extra=page%3D3
http://bbs.06climate.com/forum.p ... &extra=page%3D4

方法二：
脚本如下：

'reinit'
'open c:/test3/echtxa85.ctl '
'open c:/test3/echrxa85.ctl '
'define tk=tmp850mb.1' ;*某高度层的开氏温度K
'define tc=tk-273.16' ;*某高度层的摄氏温度C
'define rh=rh850mb.2' ;*某高度层的相对湿度Relative humidity%
'define prs=850.0' ;*获得某层高度的气压
'set lon 70 135'
'set lat 15 55'
'set mpdset cnworld'
'set grid on'
*'set csmooth on'  ; *光滑开关
******************************************************************************************
*求饱和水汽压Tetens经验公式
*水面es,tk开氏温度,tc摄氏温度
if(tk>273.16)
'define es=6.1078*exp(17.2693882*tc/(tk-35.86))'
endif
*冰面es,tk开氏温度,tc摄氏温度
if(tk<=273.16)
'define es=6.1078*exp(21.8745584*tc/(tk-7.66))'
endif
*饱和比湿
'define qs=0.622*es/(prs-0.378*es)'
*用相对湿度等求比湿
'q=rh*qs/100.0'
*水汽压
'e=prs*q/(0.622+q)'
*凝结高度的绝对温度,tk起始面上绝对温度K，
'define tlcl=55.0+2840.0/(3.5*log(tk)-log(e)-4.805)'
*求假相当位温Bolton公式,se为开氏温度K
'define theta=tk*pow((1000.0/prs),(0.2854*(1.0-0.28*q)))'
'define se=theta*exp(((3376./tlcl)-2.54)*q*(1.0+0.81*q))'
******************************************************************************************
'set gxout contour'
'set cthick 5'
'd se'
'printim c:\test2\thetase2.png white'
'c'
;


结果图如下：

2.0

方法三：
脚本如下：

'reinit'
'open c:/test3/echrxa85.ctl '
'open c:/test3/echtxa85.ctl '
'set lon 70 135'
'set lat 15 55'
'set mpdset cnworld'
'set grid on'
'define rh=rh850mb.1' ;*某高度层的相对湿度Relative humidity%
'define tk=tmp850mb.2' ;*某高度层的开氏温度K
'define tc=tk-273.16' ;*某高度层的摄氏温度C
******************************************************************************************
'define es=exp((17.26*tc)/(tk-35.86))';*饱和水汽压
'define C=((log10(es)-log10(6.11))/log10(10.0))';
'define td=(235.0*C)/(7.45-C)';*露点温度
'define qs=0.622*es/(850.0-0.378*es)'
'define q=rh*qs/100.0'
'define thetase=tk*exp(0.28586*log(1000.0/850.0)+(2500.0*q)/(338.52-0.24*tk+1.24*td))'
******************************************************************************************
'set gxout contour'
'set cthick 5'
'd thetase'
'printim c:\test2\thetase4.png white'
'c'
;


结果图如下：

4.0

对比结果发现，方法一和方法三的结果比较接近，和方法二的差别比较大。

最后也恳请各位大神帮我看看我现在的方法里还有什么问题没？还有什么不对的地方？
感谢大家！

伽蓝鸟

对比一下文献里的计算结果，比如丁一汇老师2018年大气科学上那篇综述？

lm8005507771

伽蓝鸟 发表于 2019-2-1 22:39
对比一下文献里的计算结果，比如丁一汇老师2018年大气科学上那篇综述？

好的，多谢，我学习一下

yuyumage

哪种更准确呢？？？？

lm8005507771

yuyumage 发表于 2019-4-19 21:13
哪种更准确呢？？？？

最后也没有验证哪个更准

大仓鼠

请问假相当位温平面图是用来分析什么的呢?我知道剖面图可以看假相当位温的垂直变化率从而看对流不稳定.

lm8005507771

大仓鼠 发表于 2019-5-9 13:49
请问假相当位温平面图是用来分析什么的呢?我知道剖面图可以看假相当位温的垂直变化率从而看对流不稳定.

我也不知道这个平面图能干啥，我只是用来验证我们另一个项目做出来的假相当位温是否准确的

布拉风666

谢谢楼主

逐风

方法一的计算有点问题，书上（6.3.9）公式里面的pd是干空气的分压，不是脚本里面的prs。按照方法一中的编写的方法来看假相当位温只和当前的温度、气压有关了，与当前湿度无关，我觉得这应该是有问题的。仅代表个人理解和观点，希望大家交流指正

Itsljw

谢谢楼主分享！