在这里学到了挺多知识，一直没有做过贡献，今天我也发个画整层水汽通量和散度的经验贴

568679112

1，提取资料
从necp上下载下来shum，uwnd，vwnd和pres的月平均资料，然后提取资料如下：
提取shum：

1. 'reinit'
   
2. mon.1=Jan;mon.2=Feb;mon.3=Mar;mon.4=Apr;mon.5=May;mon.6=Jun;
   
3. mon.7=Jul;mon.8=Aug;mon.9=Sep;mon.10=Oct;mon.11=Nov;mon.12=Dec
   
4. 'sdfopen e:/lunwen/ziliao/shum.mon.mean.nc'
   
5. 'set gxout fwrite'
   
6. imon=1
   
7. while(imon<=12)
   
8. 'set fwrite e:/sst/ex1/sqtl/'mon.imon'_shum.grd'
   
9. y=1948
   
10. y1=1971
    
11. y2=2010
    
12. while(y1<=y2)
    
13. ta=(y1-y)*12+imon
    
14. 'set t 'ta''
    
15. nz=1
    
16. while(nz<=8)
    
17. 'set z 'nz''
    
18. 'set x 1 144'
    
19. 'set y 1 73'
    
20. 'd shum'
    
21. nz=nz+1
    
22. endwhile
    
23. y1=y1+1
    
24. endwhile
    
25. 'disable fwrite'
    
26. imon=imon+1
    
27. endwhile
    
28. ;
    

'reinit'<br /> mon.1=Jan;mon.2=Feb;mon.3=Mar;mon.4=Apr;mon.5=May;mon.6=Jun;<br /> mon.7=Jul;mon.8=Aug;mon.9=Sep;mon.10=Oct;mon.11=Nov;mon.12=Dec<br /> 'sdfopen e:/lunwen/ziliao/shum.mon.mean.nc'<br /> 'set gxout fwrite'<br /> imon=1<br /> while(imon<=12)<br /> 'set fwrite e:/sst/ex1/sqtl/'mon.imon'_shum.grd'<br /> y=1948<br /> y1=1971<br /> y2=2010<br /> while(y1<=y2)<br /> ta=(y1-y)*12+imon<br /> 'set t 'ta''<br /> nz=1<br /> while(nz<=8)<br /> 'set z 'nz''<br /> 'set x 1 144'<br /> 'set y 1 73'<br /> 'd shum'<br /> nz=nz+1<br /> endwhile<br /> y1=y1+1<br /> endwhile<br /> 'disable fwrite'<br /> imon=imon+1<br /> endwhile<br /> ;<br />

提取u，v和pres类似，只是pres只有一层，所以不用z循环了。

2.写这些提取出来资料的ctl，由于我提取的为每种资料各个月份的，所以ctl为单个月份40年的：
dset e:/sst/ex1/sqtl/jan_shum.grd
undef  -9.99e+8
xdef           144 linear  0 2.5
ydef           73 linear  -90 2.5
zdef 8 levels 1000  925 850 700 600 500 400 300
tdef           40 linear DEC1971 1yr
vars 1
shum  8  99  eigenvactors matrix of EOF
endvars
（要注意里面的层数，是8层，vars那里的8要对应层数，我一开始没注意，在这里吃力好多苦头，整了2个小时没弄出来，后来想放弃了，在论坛里不经意看到才弄出来）

u和v的ctl类似的，pres的ctl为：
dset e:/sst/ex1/sqtl/jan_pres.grd
undef  -9.99e+8
xdef           144 linear  0 2.5
ydef           73 linear  -90 2.5
zdef 1 linear 1 1
tdef           40 linear DEC1971 1yr
vars 1
p  0  99  eigenvactors matrix of EOF
endvars
（这里的vars下面的层数要写0代表1层，我一开始写了1，也是怎么弄也弄不出来！！！）

3.逐层水汽通量，是根据论坛里大一个大侠的帖子写的：

1. *flying_nuist
   
2. 'reinit'
   
3. 'open e:\sst\ex1\sqtl\jan_uwnd.ctl'
   
4. 'open e:\sst\ex1\sqtl\jan_vwnd.ctl'
   
5. 'open e:\sst\ex1\sqtl\jan_shum.ctl'
   
6. 'set fwrite e:\sst\ex1\sqtl\quv1.grd'
   
7. 'set gxout fwrite'
   
8. y1=1
   
9. while(y1<=40)
   
10. nz=1
    
11. while(nz<=8)
    
12. 'set z 'nz''
    
13. 'set t 'y1''
    
14. 'set y 1 73'
    
15. 'set x 1 144'
    
16. 'define qu=u*shum.3/9.8'
    
17. 'd qu'
    
18. nz=nz+1
    
19. endwhile
    
20. 
    
21. nz=1
    
22. while(nz<=8)
    
23. 'set z 'nz''
    
24. 'set t 'y1''
    
25. 'set y 1 73'
    
26. 'set x 1 144'
    
27. 'define qv=v.2*shum.3/9.8'
    
28. 'd qv'
    
29. nz=nz+1
    
30. endwhile
    
31. 
    
32. nz=1
    
33. while(nz<=8)
    
34. 'set z 'nz''
    
35. 'set t 'y1''
    
36. 'set y 1 73'
    
37. 'set x 1 144'
    
38. 'define qu=u*shum.3/9.8'
    
39. 'define qv=v.2*shum.3/9.8'
    
40. 'd hdivg(qu,qv)'
    
41. nz=nz+1
    
42. endwhile
    
43. y1=y1+1
    
44. endwhile
    
45. 'disable fwrite'
    
46. 'reinit'
    
47. ;
    
48. 
    
49. 
    

*flying_nuist<br /> 'reinit'<br /> 'open e:\sst\ex1\sqtl\jan_uwnd.ctl'<br /> 'open e:\sst\ex1\sqtl\jan_vwnd.ctl'<br /> 'open e:\sst\ex1\sqtl\jan_shum.ctl'<br /> 'set fwrite e:\sst\ex1\sqtl\quv1.grd'<br /> 'set gxout fwrite'<br /> y1=1<br /> while(y1<=40)<br /> nz=1<br /> while(nz<=8)<br /> 'set z 'nz''<br /> 'set t 'y1''<br /> 'set y 1 73'<br /> 'set x 1 144'<br /> 'define qu=u*shum.3/9.8'<br /> 'd qu'<br /> nz=nz+1<br /> endwhile<br /> <br /> nz=1<br /> while(nz<=8)<br /> 'set z 'nz''<br /> 'set t 'y1''<br /> 'set y 1 73'<br /> 'set x 1 144'<br /> 'define qv=v.2*shum.3/9.8'<br /> 'd qv'<br /> nz=nz+1<br /> endwhile <br /> <br /> nz=1<br /> while(nz<=8)<br /> 'set z 'nz''<br /> 'set t 'y1''<br /> 'set y 1 73'<br /> 'set x 1 144'<br /> 'define qu=u*shum.3/9.8'<br /> 'define qv=v.2*shum.3/9.8'<br /> 'd hdivg(qu,qv)'<br /> nz=nz+1<br /> endwhile<br /> y1=y1+1<br /> endwhile<br /> 'disable fwrite'<br /> 'reinit'<br /> ;<br /> <br /> <br />

对应的ctl为：
dset e:/sst/ex1/sqtl/quv1.grd
undef  -9.99e+8
xdef           144 linear  0 2.5
ydef           73 linear  -90 2.5
zdef 8 levels 1000  925 850 700 600 500 400 300
tdef           40 linear DEC1971 1yr
vars 3
qu   8  99  eigenvactors matrix of EOF
qv   8  99  eigenvactors matrix of EOF
quv  8  99  eigenvactors matrix of EOF
endvars
（注意层数）

4.整层水汽通量和散度为（也是根据上面提到的那位大侠写得）：

1. *flying_nuist
   
2. 'reinit'
   
3. 'open e:\sst\ex1\sqtl\quv1.ctl'
   
4. 'open e:\sst\ex1\sqtl\jan_pres.ctl'
   
5. 'set fwrite e:\sst\ex1\sqtl\quv1_all.grd'
   
6. 'set gxout fwrite'
   
7. ta=1
   
8. while(ta<=40)
   
9. 'set t 'ta''
   
10. 'set z 1'
    
11. 'set y 1 73'
    
12. 'set x 1 144'
    
13. 'define qua=vint(p.2,qu,300)*9.8'
    
14. 'define qva=vint(p.2,qv,300)*9.8'
    
15. 'd qua'
    
16. 
    
17. 'd qva'
    
18. 
    
19. 'd hdivg(qua,qva)'
    
20. ta=ta+1
    
21. endwhile
    
22. 'disable fwrite'
    
23. 'reinit'
    
24. ;
    

*flying_nuist<br /> 'reinit' <br /> 'open e:\sst\ex1\sqtl\quv1.ctl' <br /> 'open e:\sst\ex1\sqtl\jan_pres.ctl' <br /> 'set fwrite e:\sst\ex1\sqtl\quv1_all.grd'<br /> 'set gxout fwrite'<br /> ta=1<br /> while(ta<=40)<br /> 'set t 'ta''<br /> 'set z 1'<br /> 'set y 1 73'<br /> 'set x 1 144'<br /> 'define qua=vint(p.2,qu,300)*9.8' <br /> 'define qva=vint(p.2,qv,300)*9.8'<br /> 'd qua'<br /> <br /> 'd qva'<br /> <br /> 'd hdivg(qua,qva)'<br /> ta=ta+1<br /> endwhile<br /> 'disable fwrite'<br /> 'reinit'<br /> ;<br />

对应ctl：
dset e:/sst/ex1/sqtl/quv1_all.grd
undef  -9.99e+8
xdef           144 linear  0 2.5
ydef           73 linear  -90 2.5
zdef 1 linear 0 1
tdef           40 linear DEC1971 1yr
vars 3
qua   0  99  eigenvactors matrix of EOF
qva   0  99  eigenvactors matrix of EOF
quva  0  99  eigenvactors matrix of EOF
endvars
（也是注意层数，今天主要被ctl里的层数正懵了，所以老提醒层数！！！！！！！！！！）

5.出图gs：
'reinit'
'open e:/sst/ex1/sqtl/quv1_all.ctl'
*'set gxout stream'

'set lat 0 60'
'set lon 20 150'
'set gxout shaded'
'define a=ave(qua,t=1,t=40)'
'define b=ave(qva,t=1,t=40)'
'd sqrt(a*a+b*b)'
'set gxout vector'
'd a;b'
;
贴个图看看（还没修饰）：

言深深

新春快乐，万事如意！

568679112

言深深 发表于 2014-1-28 14:42

                               
登录/注册后可看大图

新春快乐，万事如意！

哈哈哈，3q，也祝超级版主新春快乐呀！！！！！

Aires

经验帖要顶！

568679112

Aires 发表于 2014-1-28 15:04

                               
登录/注册后可看大图

经验帖要顶！

ayzqs

顶一个，多多共享！

一方通行

ding

snowie

好东西啊！！

kongfeng0824

新年快乐

果子

顶～～