假相当位温gs！个人总结修改增值版

ljh110011 · 发表于 2016-5-17 20:58:29

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

本帖最后由 ljh110011 于 2016-5-18 17:28 编辑

grads计算并制作假相当位温的脚本
说明：
其实气象家园里面有关假相当位温的grads脚本已经很多很多了。
但是，当我这几天正好要作这种图，找了几个准备用的时候，发现好多个脚本是没有说明的，甚至有些脚本有错。
看了一轮脚本以后莫名其妙啊，只好重新看书，学习怎么计算假相当位温，顺便写了个gs。并把计算假相当位温，控制出图的一些命令用途都注释清楚了，现在共享一下。不能保证我的gs完全没有错，如果有错请大家指出，共同学习。

计算假相当位温的方法过程：

1.计算饱和水汽压es（利用Tetens经验公式）
1）水面公式
2）冰面公式
2.计算饱和比湿qs
3.计算比湿q（利用相对湿度rh和饱和比湿qs）
4.计算水汽压e
5.计算假相当位温（利用Bolton公式）

登录/注册后可看大图

好的，先放效果图：

登录/注册后可看大图

file:///C:\Users\LJH110~1\AppData\Local\Temp\msohtmlclip1\01\clip_image002.jpg

多情况的面积平均垂直假相当位温图

登录/注册后可看大图

file:///C:\Users\LJH110~1\AppData\Local\Temp\msohtmlclip1\01\clip_image004.jpg

纬向假相当位温剖面图

登录/注册后可看大图

file:///C:\Users\LJH110~1\AppData\Local\Temp\msohtmlclip1\01\clip_image006.jpg

某一点的垂直假相当位温图

登录/注册后可看大图

经向向假相当位温剖面图

脚本gs程序：（*****************之间的部分不需要修改****************）

'reinit'
'sdfopen h:/Ryan/ERA-Int_pl_19950917_19950923.nc'
'set t 5' ;*注意时间
'set lat 10 29'
'set lon 108 128'
'set lev 1000 100'
'set zlog on' ;*Z轴取对数坐标（气压不等距）
'set ylint 100' ;*纵坐标间隔
'set xlint 2' ;*横坐标间隔
'set clopts -1 -1 0.14' ;*等值线标值属性颜色粗细大小
'set xlopts 1 4 0.14' ;*坐标刻度和标值的属性颜色粗细大小
'set ylopts 1 4 0.14'
'set csmooth on' ;*光滑开关
'define tc=t-273.16' ;*某高度层的摄氏温度C
'define tk=t' ;*某高度层的开氏温度K
'define rh=r' ;*某高度层的相对湿度Relative humidity%
'define prs=lev' ;*获得某层高度的气压
******************************************************************************************
*求饱和水汽压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-35.86))'
endif
*饱和比湿
'define qs=0.622*es/(prs-0.378*es)'
*用相对湿度等求比湿
'q=rh*qs/100'
*水汽压
'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/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 gxout line'
'set lat 14.7' ;*台风中心经纬度，如果下面画的是区域平均，这里不起作用，但是还是要设置
'set lon 114.1'
'set lev 1000 150'
'set cint 5'
*'set clopts -1 -1 0.12'
'set xlint 5' ;*设置x坐标间隔
'set axlim 345 365' ;*设置x坐标值的范围
*设置线条颜色，标记等
'set ccolor 1'
'set cmark 2'
'set cstyle 1'
'set grads off'
'set grid off'
'd ave(ave(se,lon=113.17,lon=115.03),lat=13.80,lat=15.60)' ;*台风中心半径20km范围的经纬度范围
*'d se'