广义位温和广义湿位涡

平流层的萝卜

最近要用到位涡，分享一个自己编写的广义位温和广义湿位涡的gs，所用资料ncep1*1fnl, 请交流指正。

吴国雄的位涡文章（用位温算的是位涡）：


把位温换成广义位温后，位涡就变成广义湿位涡：

gs语句：
'open D:\selection\mpv_test\ss.ctl'
'open D:\selection\quality_rain\rain.ctl'
    i=1
    while(i<=4)
    'set grads off'
    'set mpdset cnworld'
    'set t 'i''
    'set lon 70 150'
    'set lat 15 55'
    'set z 5 7'
*计算广义位温
   'define p=lev'
   'define theta=tmpprs*pow(1000/p,0.286)'
   'define Cp=0.24'
****cp=0.24卡/克*度=1.004焦耳/克*度
   'define L=597.3-0.566*(tmpprs-273.16)'
****L为凝结函数
   'define a=17.27'
   'define b=35.86'
   'define es=6.1078*exp(a*(tmpprs-273.16)/(tmpprs-b))'
*****此处计算了es,比较了es和p=850之间的大小，q/qs约等于e/es。
    'define qs=622*es/(p-0.378*es)'
    'define theta1=theta*exp(L/Cp*qs/1000/tmpprs*pow(rhprs/100,9))'
****theta1即为广义位温

*    'd theta1'
*    'cbarn.gs'

*****计算850hPa的广义湿位涡：
    'define Rd=287'
    'define P=85000'
    'define a=Rd/P*theta*pow(P/100000,0.286)'
***算好比容
    'define vor=hcurl(ugrdprs,vgrdprs)'
    'define rr=6.37e6'
***地球半径
    'define dx=cdiff(lon,x)*3.14/180*cos(lat*3.14/180)*rr'
    'define dy=cdiff(lat,y)*3.14/180*rr'
***差分时，x和y方向的实际距离
     'define g=9.8'
     'define zizhuan=7.292e-5'
***地球自转角速度
     'define ee=2*zizhuan*cos(lat*3.14/180)'
     'define f=2*zizhuan*sin(lat*3.14/180)'

****差分时，z方向的距离
      'set z 6'
***算850hPa广义湿位涡
     'define wz=-a*vvelprs/g'
***z坐标下的垂直速度**********此处进行了修改，dz需要放在固定层的后面
     'define dz=(hgtprs(z+1)-hgtprs(z-1))'
     'define pvx=a*(cdiff(wz,y)/dy-(vgrdprs(z+1)-vgrdprs(z-1))/dz)*cdiff(theta1,x)/dx'
     'define pvy=a*(ee+(ugrdprs(z+1)-ugrdprs(z-1))/dz-cdiff(wz,x)/dx)*cdiff(theta1,y)/dy'
     'define pvz=a*(f+cdiff(vgrdprs,x)/dx-cdiff(ugrdprs,y)/dy)*(theta1(z+1)-theta1(z-1))/dz'
****上边是广义湿位涡的x、y、z上的分量，算好相加即为广义湿位涡
    'set lon 102 130'
    'set lat 20 41'
    'define_colors.gs'
    'set gxout shaded'   
    'set clevs -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3'
     'set ccols 47 46 45 44 43 42 41 21 22 23 24 25 26 27'         
*    'd (pvx+pvy+pvz)*1e6'        
    'cnbasemap (pvx+pvy+pvz)*1e6'
*****draw广义湿位涡
    'cbarn.gs'

***画风场看看
     'set gxout vector'
     'set arrscl 0.5 20'
     'd ugrdprs;vgrdprs'

*****画对应时次的降水

    'set z 1'
     'set gxout stnmark'
     'set clevs 50 100'
     'set ccols 59 69 39'
     'set cmark 3'
    'set digsiz 0.15'

      'd maskout(rain.2(t='i+7152'),rain.2(t='i+7152')-10)'

'printim D:\selection\mpv_test\'i'.gif white'

    'c'
i=i+1
    'disable print'
endwhile

效果图：（该日20时24小时降水与该日的6h一次的风场/广义湿位涡的叠加图），左侧色标是降水值，右侧色标是广义湿位涡（单位：PVU, ）


可见广义湿位涡异常区与降水大值区对应关系较好。

平流层的萝卜

星雨 发表于 2015-4-15 21:20
楼主你好，我想问问位涡分析用在中小尺度天气过程中合适吗？位涡守恒的条件是静力平衡，但中心尺度天气过程 ...

绝热无摩擦的大气的个别位涡变化是守恒的，静力平衡怎么成为条件了。建议你搜集位涡诊断这方面的文献看一下吧，有很多的

river

很好很强大，支持！

ooo

很好，谢谢分享!

才仔

两张GIF图不显示{:eb328:}

Aires

支持一下~虽然不懂= =

topmad

很详实的记录  是我们的榜样

葱儿

多上论坛，总会有收获

quentin011099

lz，请问一下SS.ctl里面是什么变量啊，可以共享一下这个ctl文件吗？

又是那隻貓

谢谢分享！！

平流层的萝卜

quentin011099 发表于 2013-10-17 04:04

                               
登录/注册后可看大图

lz，请问一下SS.ctl里面是什么变量啊，可以共享一下这个ctl文件吗？

ss是我随便起的，就是用命令行g2ctl等命令生成的一个描述4个fnl1*1资料的ctl文件