剩余速度计算问题

鸾依雨觞

  最近在算150-10hPa的剩余速度，但是画出来的图就是不对。所用的资料是NCEP/NCAR2的1980-2015年月平均v、w资料。本人已研究好多天，也检查了程序（程序是一个师兄给的，感谢师兄，自己做了修改），但还是不对，所以想发个帖请教各位大神们，希望能寻求到帮助。。。。本人将不甚感激，谢谢啦，挺急的。。。图及程序附后：
计算公式：

计算公式：

       程序：1、位温

1. %%计算春季位温（150-10层次的）
   
2. clc;
   
3. clear;
   
4. ncdisp('f:\data\air.mon.mean.nc');
   
5. lon=ncread('f:\data\air.mon.mean.nc','lon');
   
6. lat=ncread('f:\data\air.mon.mean.nc','lat');
   
7. time=ncread('f:\data\air.mon.mean.nc','time');
   
8. lev=ncread('f:\data\air.mon.mean.nc','level');
   
9. air=ncread('f:\data\air.mon.mean.nc','air');  %温度单位:K
   
10. air=air(:,:,11:17,:);                           %提取150-10层（共7层）
    
11. air=fliplr(air);                                %设置数据纬度为自南向北
    
12. air1=air;  
    
13. [nx,ny,nz,nt]=size(air);         
    
14. j=1;%air_monave=zeros(nx,ny,nz,24);
    
15. for i=15:12:444                           
    
16.     air_monave(:,:,:,j)=nanmean(air(:,:,:,i:i+2),4); %提取1980-2015春季数据
    
17.     j=j+1;
    
18. end
    
19. [nx,ny,nz,nt]=size(air_monave);
    
20. 
    
21. air1=nanmean(air_monave,4);
    
22. air2=nanmean(air1,1);
    
23. ave_T=reshape(air2,ny,nz);  %纬向平均温度
    
24. 
    
25. te=[];te=air_monave;
    
26. t=[];t=reshape(te,144,73,7,36);
    
27. p=[150 100 70 50 30 20 10];
    
28. m=0.2856;
    
29. th=[];
    
30. for n=1:36
    
31.    for z=1:7
    
32.       for j=1:73
    
33.           for i=1:144
    
34.               th(i,j,z,n)=t(i,j,z,n)*((1000.0/p(z))^m);
    
35.           end
    
36.       end
    
37.    end
    
38. end
    
39. %th:位温；ave_T：纬向平均温度
    
40. save F:\5\quan_bd\th.mat th ave_T;
    
41. fid=fopen('F:\5\quan_bd\th.dat','w');
    
42. fwrite(fid,th,'float32');     
    
43. fclose(fid);
    
44. 
    

复制代码

2、v

1. %计算春季剩余环流v*（150-10层次的）
   
2. clc;
   
3. clear;
   
4. p(1)=15000;
   
5. p(2)=10000;
   
6. p(3)=7000;
   
7. p(4)=5000;
   
8. p(5)=3000;
   
9. p(6)=2000;
   
10. p(7)=1000;
    
11. 
    
12. ncdisp('f:\data\vwnd.mon.mean.nc');
    
13. lon=ncread('f:\data\vwnd.mon.mean.nc','lon');
    
14. lat=ncread('f:\data\vwnd.mon.mean.nc','lat');
    
15. time=ncread('f:\data\vwnd.mon.mean.nc','time');
    
16. lev=ncread('f:\data\vwnd.mon.mean.nc','level');
    
17. b=ncread('f:\data\vwnd.mon.mean.nc','vwnd');
    
18. b=b(:,:,11:17,:);
    
19. b=fliplr(b);                                %设置数据纬度为自南向北
    
20. b_jian=b;  
    
21. [nx,ny,nz,nt]=size(b);
    
22. j=1;b_xuan=zeros(nx,ny,nz,36);
    
23. for i=15:12:444                              
    
24.     b_xuan(:,:,:,j)=nanmean(b(:,:,:,i:i+2),4); %提取1980-2015春季数据
    
25.     j=j+1;
    
26. end
    
27. %k=size(b_xuan,4);
    
28. b_mon=[];b_mon=b_xuan;
    
29. %k=size(b_mon,4);
    
30. b1=[];b1=reshape(b_mon,144,73,7,36);
    
31. b2=[];b2=nanmean(b1,4);
    
32. b3=[];b3=nanmean(b2,1);
    
33. b4=reshape(b3,ny,nz);                                  %平均经向风速
    
34. bb=[];
    
35. for j=1:ny
    
36.     for i=1:nz
    
37.         bb(:,j,i)=b2(:,j,i)-b4(j,i);                   %经向风速扰动量
    
38.     end
    
39. end
    
40. 
    
41. load('f:\5\quan_bd\th.mat');
    
42. T=th;                                                  %th:位温
    
43. [nx,ny,nz,nt]=size(T);
    
44. a1=[];a1=nanmean(T,4);
    
45. a4=[];a4=reshape(nanmean(a1,1),ny,nz);                      %经向平均位温
    
46. aa=[];
    
47. for j=1:ny
    
48.     for i=1:nz
    
49.         aa(:,j,i)=a1(:,j,i)-a4(j,i);                   %位温扰动量
    
50.     end
    
51. end
    
52. 
    
53. yy=[];yy=nanmean(bb.*aa,1);
    
54. y=zeros;y=reshape(yy,ny,nz);  %位温扰动量*经向风速扰动量的平均
    
55. for i=1:ny
    
56.     for j=1:nz
    
57.         if abs(y(i,j))<0.000001
    
58.            y(i,j)=0;
    
59.         end
    
60.     end
    
61. end
    
62. 
    
63. dx=[];fun=[];
    
64. dx(1)=-(p(2)-p(1)); dx(7)=-(p(7)-p(6));
    
65. for j=1:73
    
66.         fnn(j,1)=(a4(j,2)-a4(j,1))/dx(1);
    
67.         fnn(j,7)=(a4(j,7)-a4(j,6))/dx(7);
    
68.     for k=2:6
    
69.         dx(k)=-(p(k+1)-p(k-1));
    
70.         fnn(j,k)=(a4(j,k+1)-a4(j,k-1))/(2*dx(k));           %求平均位温偏导数
    
71.     end
    
72. end
    
73. 
    
74. fnn1=[];
    
75. for j=1:73
    
76.     fnn1(j,1)=(y(j,2)/fnn(j,2)-y(j,1)/fnn(j,1))/dx(1);
    
77.     fnn1(j,7)=(y(j,7)/fnn(j,7)-y(j,6)/fnn(j,6))/dx(7);
    
78.    for k=2:6;
    
79.       fnn1(j,k)=(y(j,k+1)/fnn(j,k+1)-y(j,k-1)/fnn(j,k-1))/(2*dx(k));
    
80.    end
    
81. end
    
82. 
    
83. svdata=[];
    
84. for j=1:73
    
85.     for k=1:7;
    
86.         svdata(j,k)=b4(j,k)-fnn1(j,k);       %剩余环流v*
    
87.         if abs(svdata(j,k))<=0.00001
    
88.             svdata(j,k)=0;                     
    
89.         end;
    
90.      end
    
91. end
    
92. %svdata:剩余环流v*;a4:经向平均位温；y:位温扰动量*经向风速扰动量的平均;fnn:平均位温偏导数
    
93. %b4:平均经向风速
    
94. save F:\5\quan_bd\v.mat svdata a4 y fnn b4;
    
95. fid=fopen('F:\5\quan_bd\v.dat','w');
    
96. fwrite(fid,svdata,'float32');     
    
97. fclose(fid);
    
98. 
    
99. 
    

复制代码

3、w

1. %计算春季剩余环流w*（150-10层次的）
   
2. clc;
   
3. clear;
   
4. p(1)=15000;
   
5. p(2)=10000;
   
6. p(3)=7000;
   
7. p(4)=5000;
   
8. p(5)=3000;
   
9. p(6)=2000;
   
10. p(7)=1000;
    
11. a=6371000;
    
12. g=9.80665;
    
13. R=287.;           %气体常数
    
14. load('F:\5\quan_bd\th.mat');
    
15. ave_T=ave_T;      %纬向平均温度
    
16. ncdisp('f:\data\omega.mon.mean.nc');
    
17. lon=ncread('f:\data\omega.mon.mean.nc','lon');
    
18. lat=ncread('f:\data\omega.mon.mean.nc','lat');
    
19. time=ncread('f:\data\omega.mon.mean.nc','time');
    
20. lev=ncread('f:\data\omega.mon.mean.nc','level');
    
21. omega=ncread('f:\data\omega.mon.mean.nc','omega');
    
22. omega1=omega(:,:,11:17,:);
    
23. omega2=fliplr(omega1);
    
24. [nx,ny,nz,nt]=size(omega2);
    
25. j=1;omega2_xuan=zeros(nx,ny,nz,36);
    
26. omega3=[];
    
27. for i=15:12:444                           
    
28.     omega3(:,:,:,j)=nanmean(omega2(:,:,:,i:i+2),4); %提取1980-2015春季数据
    
29.     j=j+1;
    
30. end
    
31. %k=size(omega3,4);
    
32. omega3_mon=[];omega3_mon=omega3;
    
33. omega4=reshape(omega3_mon,nx,ny,nz,36);
    
34. omega5=nanmean(omega4,4);
    
35. o1=nanmean(omega5,1);
    
36. o2=reshape(o1,ny,nz);                    %w*第一项
    
37. 
    
38. %svdata:剩余环流v*；y:位温扰动量*经向风速扰动量的平均;fnn:平均位温偏导数
    
39. load('F:\5\quan_bd\v.mat');
    
40. epf=y;t3=fnn;
    
41. [ny,nz]=size(epf);
    
42. 
    
43. lat1=flipud(lat);
    
44. for i=1:73
    
45.     coskk(i)=cos(lat1(i)/180*pi);
    
46.     if abs(lat1(i))==90
    
47.         coskk(i)=0;
    
48.     end
    
49. end
    
50. 
    
51. divx=[];
    
52. for  j=1:nz                               %w*第二项
    
53.     dx=2.5/180*pi;
    
54.     divx(1,j)=(epf(2,j)*coskk(2)/t3(2,j)-epf(1,j)...
    
55.         *coskk(1)/t3(1,j))/dx/(a*coskk(1));
    
56.     divx(ny,j)=(epf(ny,j)*coskk(ny)/t3(ny,j)-epf(ny-1,j)...
    
57.       *coskk(ny-1)/t3(ny-1,j))/dx/(a*coskk(ny));
    
58.     for i=2:ny-1
    
59.        divx(i,j)=(epf(i+1,j)*coskk(i+1)/t3(i+1,j)-epf(i-1,j)*...
    
60.      coskk(i-1)/t3(i-1,j))/(2*dx)/(a*coskk(i));
    
61.     end
    
62. end
    
63. 
    
64. womi=[];
    
65. for j=1:ny
    
66.    for k=1:nz;
    
67.        womi(j,k)=o2(j,k)+divx(j,k);
    
68.       if abs(womi(j,k))>=1000000
    
69.        womi(j,k)=0;
    
70.       end;
    
71.     end
    
72. end
    
73. 
    
74. w=[];
    
75. for k=1:nz
    
76.    for j=1:ny
    
77.        w(j,k)=-womi(j,k)*R*ave_T(j,k)/(g*p(k));       %转换为Z坐标,单位：m/s
    
78.    end
    
79. end
    
80. %womi:剩余环流w*；w：z坐标下w*
    
81. save F:\5\quan_bd\w.mat womi w;
    
82. 
    
83. fid=fopen('F:\5\quan_bd\w.dat','w');
    
84. fwrite(fid,w,'float32');     
    
85. fclose(fid);
    
86. 
    
87. 
    

复制代码

ctl文件：       gs文件：
画出来的图：

鸾依雨觞

怎么都没人理我呢？

伽蓝鸟

虽然是陈老师的学生，可我不做行星波那块。。

鸾依雨觞

那请问你有陈老师关于这个的程序吗？想对比看看是不是我的程序哪步弄错了，画不出来老师那个图啊

indigo_wang

请问楼主解决了嘛