GrADS计算湿位涡全套程序

左手疆环

清风大师，我参考了您分享的gs,然后修改了自己的gs文件，可是到最后算出来的mpv值不管什么时次、层次全是0。检查好久，实在是不知道哪一块出错了，想向您请教，麻烦您抽空帮我看一下，万分谢谢！~'reinit'
'set gxout fwrite'
'set fwrite D:\keyan\fnl\vars\05\mpv.dat'
'open D:\keyan\fnl\vars\05\eqt.ctl'
'open D:\keyan\fnl\vars\05\uwind.ctl'
'open D:\keyan\fnl\vars\05\vwind.ctl'
it=1
while(it<=124)
'set x 1 181'
'set y 1 181'
'set t 'it''
'define thetse=eqt.1'
'define u=ugrdprs.2'
'define v=vgrdprs.3'
****************************1000等压面
'set lev 1000'
lev1=1000
lev2=975
'define f=2*7.292*sin(lat*3.14159/180)*1e-05'
'define g=9.8'
'define dp=2500.0'
'define dthetse=thetse(lev='lev1')-thetse(lev='lev2')'
'define du=u(lev='lev1')-u(lev='lev2')'
'define dv=v(lev='lev1')-v(lev='lev2')'
'define dx=2*6.37e6*cos(lat*3.1415926/180.0)*3.1416/180.0'
'define dy=2*6.37e6*3.1416/180.0'
'define dtx=cdiff(thetse,x)'
'define dty=cdiff(thetse,y)'
'define mpv1=-g*(hcurl(u,v)+f)*dthetse/dp'
'define mpv2=g*((dv/dp)*(dtx/dx)-(du/dp)*(dty/dy))'
'define mpv=mpv1+mpv2'
'd mpv'
******************************此段为间隔为25等压面层次
iz=1
while(iz<=3)
level=1000-iz*25
lev1=1000-(iz-1)*25
lev2=1000-(iz+1)*25
'set lev 'level''
'define f=2*7.292*sin(lat*3.14159/180)*1e-05'
'define g=9.8'
'define dp=5000.0'
'define dthetse=thetse(lev='lev1')-thetse(lev='lev2')'
'define du=u(lev='lev1')-u(lev='lev2')'
'define dv=v(lev='lev1')-v(lev='lev2')'
'define dx=2*6.37e6*cos(lat*3.1415926/180.0)*3.1416/180.0'
'define dy=2*6.37e6*3.1416/180.0'
'define dtx=cdiff(thetse,x)'
'define dty=cdiff(thetse,y)'
'define mpv1=-g*(hcurl(u,v)+f)*dthetse/dp'
'define mpv2=g*((dv/dp)*(dtx/dx)-(du/dp)*(dty/dy))'
'define mpv=mpv1+mpv2'
'd mpv'
iz=iz+1
endwhile
******************此段为间隔为50等压面层次
iz=2
while(iz<=17)
level=1000-iz*50
lev1=1000-(iz-1)*50
lev2=1000-(iz+1)*50
'set lev 'level''
'define f=2*7.292*sin(lat*3.14159/180)*1e-05'
'define g=9.8'
'define dp=10000.0'
'define dthetse=thetse(lev='lev1')-thetse(lev='lev2')'
'define du=u(lev='lev1')-u(lev='lev2')'
'define dv=v(lev='lev1')-v(lev='lev2')'
'define dx=2*6.37e6*cos(lat*3.1415926/180.0)*3.1416/180.0'
'define dy=2*6.37e6*3.1416/180.0'
'define dtx=cdiff(thetse,x)'
'define dty=cdiff(thetse,y)'
'define mpv1=-g*(hcurl(u,v)+f)*dthetse/dp'
'define mpv2=g*((dv/dp)*(dtx/dx)-(du/dp)*(dty/dy))'
'define mpv=mpv1+mpv2'
'd mpv'
iz=iz+1
endwhile
******************************100等压面
'set lev 100'
lev1=150
lev2=100
'define f=2*7.292*sin(lat*3.14159/180)*1e-05'
'define g=9.8'
'define dp=5000.0'
'define dthetse=thetse(lev='lev1')-thetse(lev='lev2')'
'define du=u(lev='lev1')-u(lev='lev2')'
'define dv=v(lev='lev1')-v(lev='lev2')'
'define dx=2*6.37e6*cos(lat*3.1415926/180.0)*3.1416/180.0'
'define dy=2*6.37e6*3.1416/180.0'
'define dtx=cdiff(thetse,x)'
'define dty=cdiff(thetse,y)'
'define mpv1=-g*(hcurl(u,v)+f)*dthetse/dp'
'define mpv2=g*((dv/dp)*(dtx/dx)-(du/dp)*(dty/dy))'
'define mpv=mpv1+mpv2'
'd mpv'
it=it+1
endwhile
'disable fwrite'
;

mofangbao

左手疆环 发表于 2013-11-26 15:24

                               
登录/注册后可看大图

清风大师，我参考了您分享的gs,然后修改了自己的gs文件，可是到最后算出来的mpv值不管什么时次、层次全是0。 ...

这个程序也不是我写的  我也是分享一下  抱歉不能给你更多的帮助

yuzhanyou09

很好啊，可以学习一下

左手疆环

mofangbao 发表于 2013-11-26 19:30

                               
登录/注册后可看大图

这个程序也不是我写的  我也是分享一下  抱歉不能给你更多的帮助

没事，谢谢您啊~~~

1092995269

厄，我醉了

jinlili_001

apple973 发表于 2013-11-28 14:47

                               
登录/注册后可看大图

谢谢分享，打算带走

谢谢分享，可惜我没权限哪

大果子

请问PV1和PV2分别代表什么呀？湿位涡是将它俩相加吗？谢谢！

wuyu

好东西，收了

hust420

谢谢

shan

谢谢分享，下载学习啊