湿位涡的GS程序

蓝月亮

不知道有没有人发过，从动力论坛的一个回复里面翻出来的，比较难找，转载后希望有更多的人看到，原帖子地址：
http://bbs.lasg.ac.cn/bbs/thread-37279-5-4.html

1. 
   
2. title.1=1000;title.2=925;title.3=850;title.4=700;title.5=500
   
3. title.6=400;title.7=300;title.8=250;title.9=200;title.10=150;title.11=100;
   
4. t.1=09073008;t.2=09073020;t.3=09073108;t.4=09073120
   
5. 'open e:\swwsz\111.ctl'
   
6. 'set grads off'
   
7. 'set lon 110 135'
   
8. 'set lat 30 50'
   
9. 'set lev 850'
   
10. 'set t 1'
    
11. 'smth9 (mpv2)'
    
12. *'set mpdset cnworld'
    
13. i=1
    
14. while(i<=11)
    
15. 'set z 'i
    
16. j=1
    
17. while(j<=4)
    
18. 'set t 'j
    
19. *'set z 'z''
    
20. 'define f=2*7.292*sin(lat*3.14159/180)'
    
21. 'define g=9.8'
    
22. 'define dthetse=thetse(z+1)-thetse(z-1)'
    
23. 'define du=u(z+1)-u(z-1)'
    
24. 'define dv=v(z+1)-v(z-1)'
    
25. 'define dp=100*(lev(z+1)-lev(z-1))'
    
26. 'define dx=cdiff(lon,x)*3.1416/180'
    
27. 'define dy=cdiff(lat,y)*3.1416/180'
    
28. 'define mpv1=-g*(vor*1e2+f)*dthetse/dp'
    
29. 'define mpv2=g*((dv/dp)*(dthetse/dx)-(du/dp)*(dthetse/dy))'
    
30. *'define mpv2=g*((dv/dp)*cdiff(thetse,x)-(du/dp)*cdiff(thetse,y))'
    
31. 'define mpv=mpv1+mpv2'
    
32. *'d mpv1'
    
33. 'd mpv2'
    
34. 'draw title 't.j' 'title.i'hPa mpv2'
    
35. 'printim E:\swwsz\2\'t.j''title.i'hPa_mpv2.gif white'
    
36. 'c'
    
37. j=j+1
    
38. endwhile
    
39. i=i+1
    
40. endwhile
    
41. *'d mpv'
    

复制代码

没有测试过...

传说中的谁

x方向差分的计算有点问题，木有考虑纬圈的半径是随纬度变化而变化的（实际上计算结果差别好像不大）。个人感觉应该写成：'define dx=6.37e6*cos(lat*3.1415926/180)*cdiff(lon,x)*3.1416/180'

传说中的谁

嗯，需要先算出theta_se，要不在做theta_se的垂直差分时会出错的。感谢蓝月亮的分享。

zsmnantong

感谢小篮分享，有时间调一下看看

善人/jw

感谢蓝月亮~~~

godenflame135

有意思！多谢楼主分享。

meehooqq

顶一个,收藏了!

kazami2003

感謝，測試看看

iloveplmmok

以前算过这个，呵呵

wutongyeyu

拿走  去试试