请教大家用EC数据画假相当位温的问题

lm8005507771

我现在有的EC的grib2格式数据是这个样子的，

                               
登录/注册后可看大图

，
里面的不同气象要素是在不同文件中的，算假相当位温要用的温度、相对湿度、气压的文件我都能找到，但怎么把它们放到一起计算，遇到问题了如下。
我想了两个办法。第一个办法是把三个数据文件解码以后，在脚本中对它们的ctl文件分别open，然后define不同要素不同的变量，带入公式计算，遇到的问题就是，好像只有第一个被open的参数才能被define，第二第三个被打开的就不行，请教各位大神，这个情况怎么弄，或者怎么百度，我连这种情况怎么搜索都不清楚。脚本用的是http://bbs.06climate.com/forum.p ... &extra=page%3D2 中给的，我把前后set的部分都是删了，想的是先出个结果，然后再调试具体内容。

'reinit'
'open f:/test1/echrxa85.ctl'
'open f:/test1/echtxa85.ctl'
'open f:/test1/echpxa89.ctl'

'define tk=tmp850mb' ;*某高度层的开氏温度K
'define tc=tk+273.16' ;*某高度层的摄氏温度C
'define rh=rh850mb' ;*某高度层的相对湿度Relative humidity%
'define prs=presmsl' ;*获得某层高度的气压

******************************************************************************************
*求饱和水汽压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'
*水汽压
'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))'
******************************************************************************************

'd se'

执行脚本后报错内容如下：

脚本出现的问题

第二个办法是根据兰溪大神的pdf中写的把三个文件复制成一个文件中再进行操作，利用dos下copy的命令来做的“copy /b ech*.032 test001.032” (我已经提前把三个参数的文件放到一个单独的文件夹里了)，但是这样操作之后，再进行!g2ctl和!gribmap后还是不行。
下面是合成以后文件的ctl内容。我的水平看起来这个好像都对着呢啊。
dset c:/short1/test001.032
index c:/short1/test001.032.idx
undef 9.999E+20
title c:/short1/test001.032
* produced by g2ctl v0.1.1
* command line options: c:/short1/test001.032
* griddef=1:0:(720 x 361):grid_template=0:winds(N/S): lat-lon grid:(720 x 361) units 1e-06 input WE:NS output WE:SN res 48 lat 90.000000 to -90.000000 by 0.500000 lon 0.000000 to 359.500000 by 0.500000 #points=259920:winds(N/S)

dtype grib2
ydef 361 linear -90.000000 0.5
xdef 720 linear 0.000000 0.500000
tdef 1 linear 12Z03jan2019 1mo
zdef 1 linear 1 1
vars 3
PRESmsl   0,101,0   0,3,0 ** mean sea level Pressure [Pa]
RH850mb   0,100,85000   0,1,1 ** 850 mb Relative Humidity [%]
TMP850mb   0,100,85000   0,0,0 ** 850 mb Temperature [K]
ENDVARS
EDEF 1
E2 1  12Z03jan2019 2
ENDEDEF


把合成的文件放到脚本里出不来图，我单独打开它的ctl，对三个参数画图检查一下有问题没，结果相对湿度有问题。

三合一文件出现的问题

请教大家，应该怎么处理？求大家指点，刚接触grabs，有些操作不当的地方还请大家不吝赐教！

lm8005507771

求各位大神指教！

伽蓝鸟

GrADS里，如果打开多个文件，比如先打开air.ctl, 第二个打开uwnd.ctl
需要使用其中变量的时候，d air.1   d uwnd.2
需要在变量后面加一个.第几个文件

lm8005507771

伽蓝鸟 发表于 2019-1-28 08:24
GrADS里，如果打开多个文件，比如先打开air.ctl, 第二个打开uwnd.ctl
需要使用其中变量的时候，d air.1    ...

意思是 第一个文件里的变量是t，第二个文件里的变量是rh，那么我在调用的时候应该写：t.1 和 rh.2 ？
感谢感谢！

伽蓝鸟

lm8005507771 发表于 2019-1-28 17:46
意思是 第一个文件里的变量是t，第二个文件里的变量是rh，那么我在调用的时候应该写：t.1 和 rh.2 ？
感 ...

是的。（为啥回复的总字数还得大于10呢。。）

lm8005507771

伽蓝鸟 发表于 2019-1-28 19:07
是的。（为啥回复的总字数还得大于10呢。。）

这个不敢吐槽（悄悄@版主们问问看）

lm8005507771

自己回复自己

问题基本都解决了，下面贴上我现在的脚本，方便以后我这样的小白交流。

'reinit'
*'open c:/test2/echhxa50.ctl '
*'open c:/test2/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' ;*某高度层的开氏温度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 cint 0.2'
'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、打开多个文件时，引用变量后面要加".1"、".2".... 按照打开文件的顺序来，打开的第一个文件的变量后面加“.1”，第二个文件的变量后面加".2"
2、注意自己编辑的公式，网上贴出来的一些公式会有小问题，所以要仔细核对公式，不然算出来的也不知道是啥，而且我一开始的原始数据还出现了问题，最后换了数据才搞好，这个也要注意