求取水汽通量散度的一些gs文件

mofangbao

首先感谢3群的一位朋友给我提供了很多他用过的gs文件，这些文件就是从那里面得到了，如果这些文件对比有帮助，不妨也可以贴出你自己的gs文件，分享给大家，谢谢~

1、700hpa水汽通量散度（标准）

1. 
   
2. 'reinit'
   
3. 'open d:\12090.ctl'
   
4. 
   
5. 'set lev 700'
   
6. #'set lev 850'
   
7. 'set lat 20 35'
   
8. 'set lon 100 120'
   
9. 'set t 5'
   
10. 'set mpdset cnworld'
    
11. 'set xlopts 1 4 0.15'
    
12. 'set ylopts 1 4 0.15'
    
13. 'set grads off'
    
14. 'set timelab off'
    
15. 'set grid off'
    
16. 
    
17. #'tc=(tmpprs-273.16)'
    
18. 'td=tc-((14.55+0.114*tc)*(1-0.01*rh) + pow((2.5+0.007*tc)*(1-0.01*rh),3) + (15.9+0.37*tc)*pow((1-0.01*rh),14))'
    
19. 'vapr=6.112*exp((17.67*td)/(td+243.5))'
    
20. 'e=vapr*1.001+(lev-100)/900*0.0034'
    
21. 'define mixr=0.62137*(e/(lev-e))*1000'
    
22. 'define qx=u*mixr'
    
23. 'define qy=v*mixr'
    
24. 'define mconv=hdivg(qx,qy)*1e6/9.8'
    
25. #水汽通量的计算方法就是hdivg（q*v)/g,负值区就是水汽的辐合区。单位是g/cm2.hpa.s. 量级为10-6到10-7.
    
26. #'define mconv=(-1)*hdivg(qx,qy)*1e6/9.8',这个在前面加-1就相反，正值区为辅合区
    
27. 'enable print d:\tlsandu1.gmf'
    
28. 'd mconv'
    
29. 'print'
    
30. 'disable print'
    

复制代码

2、水汽通量与等值线

1. 
   
2. #该gs绘出水汽通量是阴影图
   
3. 'reinit'
   
4. 'open f:\tt.ctl'
   
5. 'enable print d:\1243567788sqtl.gmf'
   
6. #经纬度设置
   
7. 'set lat 20 35 '
   
8. 'set lon 100 120'
   
9. 'set lev 700'
   
10. #先设定时间循环
    
11. t=1
    
12. while(t<=9)
    
13. 'c'
    
14. 'set t 't''
    
15. #设置图片角标底图
    
16. 'set grads off'
    
17. 'set grid off'
    
18. 'set mpdset cnworld'
    
19. #'set map  15 5 1'
    
20. #'draw map'
    
21. #设置划线格式
    
22. #'set cthick 9'
    
23. #'set xlopts 1 4 0.12'
    
24. #'set ylopts 1 4 0.12'
    
25. #'set cint 0.5'
    
26. #★重点：编程水气通量的算法
    
27. 'define prs=lev'
    
28. #饱和水汽压的计算，单位：毫巴
    
29. 'define es=6.112*exp(17.67*(Tk-273.15)/(Tk-29.65))'
    
30. #饱和比湿的计算  单位：克/克
    
31. 'define qs=0.622*es/(lev-0.278*es)'
    
32. #rh是比湿，利用比湿求水汽
    
33. 'define q=rh*qs/100'
    
34. #a为|V|
    
35. 'define a=sqrt(u*u+v*v)'
    
36. #地球引力系数
    
37. 'define g=9.8'
    
38. #以阴影图显示，命令一直到run cbarn.gs
    
39. 'set gxout shaded'
    
40. 'set cmin 10'
    
41. #10以上才画阴影图
    
42. 'd a*q*1000/9.8'
    
43. #本来是q*|V|/g，但这里q的单位为：克/克，而水汽通量中的q一般取克、千克，所以乘以1000
    
44. #显示一下尺度标,画阴影图用
    
45. 'run cbarn.gs'
    
46. #set gxout contour是重新设定图形输出格式为二维等值线图，叠加在阴影上，所以一般先画阴影，再画等值线，不然阴影覆盖了等值线
    
47. 'set gxout contour'
    
48. 'set ccolor 0'
    
49. 'set cmin 6'
    
50. #只画值为6以上的等值线
    
51. 'd a*q*1000/9.8'
    
52. 
    
53. 'print'
    
54. pull dummy
    
55. t=t+1
    
56. endwhile
    
57. 'disable print'
    

复制代码

3、水汽通量-标准

1. 
   
2. * This is a script for displaying moisture convergence
   
3. * Written by Michael Maxwell
   
4. *   
   
5. * rh    = Relative Humidity in %
   
6. * t     = Temp at *set level in degrees Kelvin
   
7. * tc    = Temp in degrees C
   
8. * td    = Dewpoint at *set level in degrees C
   
9. * e     = Vapor pressure
   
10. * mixr  = Mixing ratio
    
11. * u     = U-wind in m/s
    
12. * v     = V-wind in m/s
    
13. * mconv = moisture convergence/divergence. convergence is positive and divergence is negative.
    
14. 'reinit'
    
15. 'open d:\12090.ctl'
    
16. 
    
17. 'set lev 700'
    
18. 'set lat 20 35'
    
19. 'set lon 100 120'
    
20. 'set t 5'
    
21. 'set mpdset cnworld'
    
22. 'set xlopts 1 4 0.15'
    
23. 'set ylopts 1 4 0.15'
    
24. 'set grads off'
    
25. 'set timelab off'
    
26. 'set grid off'
    
27. 
    
28. #'tc=(tmpprs-273.16)'
    
29. 'td=tc-((14.55+0.114*tc)*(1-0.01*rh) + pow((2.5+0.007*tc)*(1-0.01*rh),3) + (15.9+0.37*tc)*pow((1-0.01*rh),14))'
    
30. 'vapr=6.112*exp((17.67*td)/(td+243.5))'
    
31. 'e=vapr*1.001+(lev-100)/900*0.0034'
    
32. 'define mixr=0.62137*(e/(lev-e))*1000/9.8'
    
33. #通过此gs文件算出的水汽通量的单位是g/(cm.s),一般为 十几到几十
    
34. 'define qx=u*mixr'
    
35. 'define qy=v*mixr'
    
36. 'define mconv=hdivg(qx,qy)*1e6'
    
37. 'enable print d:\2lb-shuiqi900.gmf'
    
38. 'set arrscl 1 500'
    
39. #'set lon 124'
    
40. #'set cmin 30'
    
41. 'set gxout shaded'
    
42. 'd mag(qx,qy)'
    
43. #'set cmin 30'
    
44. 'set gxout contour'
    
45. 'd mag(qx,qy)'
    
46. #'d theta'
    
47. #'d skip(u,8);skip(v,8)'
    
48. 'print'
    
49. 'disable print'
    

复制代码

mofangbao

liantai4567 发表于 2011-10-7 20:57

                               
登录/注册后可看大图

版权啊！这些个gs文件不都是我的吗，连数据的描述文件“12090.ctl“和”tt.ctl”都是我的哎，直接没改就贴上 ...

对大家当然会有帮助啦，感谢你提供，当时找到这个脚本的时候直接就贴上来了

嘟嘟鱼

在计算这些物理量的时候，NCL有比GrADS更大的优势，希望敢于尝试新事物的同志多多专研，一起多多交流！！！

传说中的谁

建议每种代码都说明一下用的是什么资料:lol

zhaolanxiai

好东东，谢谢分享了！

飞雪残剑

楼主，你的.ctl文件所描述的是什么格式的文件，有没有NECP FNL格式文件的gs呀？

mofangbao

飞雪残剑 发表于 2011-7-2 00:44

                               
登录/注册后可看大图

楼主，你的.ctl文件所描述的是什么格式的文件，有没有NECP FNL格式文件的gs呀？

fnl我没用过，如果是一般的nc是不需要ctl的，fnl的好像要导出ctl的，不知道是不是这样的

269984288

分享了~~~:)谢谢楼主

传说中的谁

用的时候要注意，有的变量没有定义，例如rh，可以自己定义一下，'define rh=rhprs'

aaayyy110119

有用的东东~~谢谢

MeteoInfo

不知能否提供一个具体的例子，包括数据。我想用MeteoInfo写script试试能否做这项工作。