grads自定义gs脚本函数的应用心得：常用物理量的脚本函数

平流层的萝卜

      由于总是要用到那么几个变量（水汽通量及其散度、假相当位温、位涡等），自然想到如果把它们的计算公式写成一个脚本gs，下次再计算它们的时候直接调用该gs，而不必用到一次重新写一次，岂不方便很多！好比想算散度，直接define sandu=hdivg(u,v),sandu就这么算出来了。研究了很长时间function的用法，终于有点小进展，事实证明，用自定义的gs脚本函数，也可以像define sandu=hdivg(u,v)一样方便地计算物理量。
     举个例子，如果想算假相当位温，可以直接“thetase 温度 湿度 气压”，就可以算出来theatase。方法就是用到了自定义gs脚本函数。具体做法如下：先上图：


测试的gs脚本如下：
'open H:\paper\fnl.ctl'
'enable print H:\paper\cross\thetase_test.gmf'

    'set grads off'
    'set mpdset cnworld'
    'set t 5'
    'set lon 106 126'
    'set lat 24 38'
    'set lev 850'
     'thetase tmpprs rhprs lev'
     'd thetase'

    'drawtime'
    'print'
    'disable print'
    ;
    thetase即假相当位温，是温度、湿度和气压的函数，建立一种函数关系：thetase  输入tmpprs温度，rhprs湿度和lev气压，就可以算出来thetase. 当然，不同的情况下，温度，湿度和气压的变量名不同，需要简单修改一下。
thetase是自己写的一个gs脚本函数文件的文件名，theatase.gs的内容如下：
function thetase(args)
t=subwrd(args,1)
rh=subwrd(args,2)
p=subwrd(args,3)


'define es=(6.112*exp(17.67*('t'-273.15)/('t'-29.65)))'
'define q='rh'*(0.62197*es/('p'-0.378*es))/100' （ 0.378*）应该是要加上去的，计算比湿所用，第一次给的gs此处未加，现已修改。
'define e='p'*q/(0.62197+q)+1e-10'
'define tlcl=55.0+2840.0/(3.5*log('t')-log(e)-4.805)'
'define theta='t'*pow((1000/'p'),(0.2854*(1.0-0.28*q)))'
'define thetase=theta*exp(((3376./tlcl)-2.54)*q*(1.0+0.81*q))'
return
   把上边的这些保存在thetase.gs里，放在grads的script文件夹下边，下次调用的时候就可以直接用了！
  函数用法不用多说了，东西不多，很简单。
  输入几个参量，就能得到想要的结果，这正是函数得以存在的原因。
  所有的物理量计算应该都可以这样弄，这样可以搞一个物理量gs脚本函数的大集合！方便所有气象人！！
  待续，欢迎交流指正。 下边两个附件，是我重新用其他数据做的，可以正常出图的。但是脚本到底是否通用，还请各位用自己的数据一试，做的不对的地方请指正！！！多交流！！
test1.rar (44.48 KB, 下载次数: 54)

2013-11-22 14:06 上传

点击文件名下载附件
用来测试的grd,ctl和gs 测试用的文件
thetase.gs (460 Bytes, 下载次数: 83)

2013-11-23 19:43 上传

点击文件名下载附件
gs脚本函数文件 gs脚本函数文件（放在script文件夹下）
续2：
根据上述的想法，继续写一个水汽通量散度的gs脚本函数：
脚本函数文件内容如下：
function qvdiv(args)
u=subwrd(args,1)
v=subwrd(args,2)
t=subwrd(args,3)
r=subwrd(args,4)
p=subwrd(args,5)
    'define es=6.112*exp(17.67*('t'-273.15)/('t'-29.65))'    计算饱和水汽压，Tetan公式
    'define qs=0.62197*es/('p'-0.378*es)*1000'           饱和比湿
    'define q=qs*'r'/100'     比湿
    'define qx=q*'u'/9.8'     水汽通量的x方向分量
    'define qy=q*'v'/9.8'     水汽通量的y方向分量
    'define qvdiv=hdivg(qx,qy)'    水汽通量散度qvdiv
    'undefine es'
    'undefine qs'
    'undefine q'
    'undefine qx'
    'undefine qy'
return
用法：
qvdiv u风速 v风速 温度（单位：K） 湿度（单位：0.01） 气压（单位：百帕）
函数算得qvdiv即是水汽通量散度。
用fnl资料测试之：
'open H:\paper\fnl.ctl'
'enable print H:\paper\water_vapor_stream\water.gmf'
i=1
while(i<=13)
    'set t 'i''
    'set lon 80 140'
    'set lat 0 45'
    'set lev 850'
    'set grads off'
    'set mpdset mres cnriver '
    'qvdiv ugrdprs vgrdprs tmpprs rhprs 850'
     'set gxout shaded'
     'set lon 81 139'
    'set lat 1 44'
     'define_colors'
     'set clevs -100 -80 -60 -40 -20 0 20 40 60 80'
     'set ccols 26 25 24 23 22 21 41 42 43 44 45'
    'd qvdiv*1e6'    
    'cbarn'
    ’drawtime'
    'print'
    'c'
i=i+1
endwhile
'disable print'
;

部分结果如下：


水汽通量散度的gs脚本函数： qvdiv.gs (433 Bytes, 下载次数: 78)

2013-11-24 19:53 上传

点击文件名下载附件 ，放在script文件夹下边，以后算水汽通量散度的时候只需要输入几个参数就好了。

kongfeng0824

谢谢分享  

斥鷃

学grads这么久了，都没好好用过函数的说~楼主要多分享经验哈~

善人/jw

很不错啊，感谢分享

river

很不错，谢谢分享

lqouc

很不错，收藏了

mofangbao

楼主能否提供一个小示例数据供测试，按照官方的描述，默认情况下主文件是取不到脚本函数里面的变量的，我看你直接'd thetase'了，这个应该是脚本函数的变量吧，而你在脚本函数也没有将thetase这个变量return，所以有些怀疑，不过还是非常鼓励你这种分享精神啦，可以继续讨论一下，也许是GrADS的规则更新了说不定。

dmjdmj11

赞！

ningning

楼主威武！

ayzqs

回去试一试，谢谢分享啦