matlab求南海夏季风指数差别很大

陈伯英

求助，想问问这个程序哪里出了问题，最后算出的南海夏季风指数与李崇银老师的差别过大了程序参考了这个帖子↓
李崇银南海夏季风指数matlab程序
http://bbs.06climate.com/forum.p ... &fromuid=128407
(出处: 气象家园)

1. clc
   
2. clear
   
3. close all;
   
4. 
   
5. div1=zeros(9,5);
   
6. div2=zeros(9,5);
   
7. D=zeros(9,5,9);
   
8. Id=zeros(9,5,9);
   
9. u=ncread('I:\NCEPdata\pressure\uwnd.1984.nc','uwnd');
   
10. v=ncread('I:\NCEPdata\pressure\vwnd.1984.nc','vwnd');
    
11. ncdisp I:\NCEPdata\pressure\uwnd.1984.nc
    
12. ncdisp I:\NCEPdata\pressure\vwnd.1984.nc
    
13. u_lon=ncread('I:\NCEPdata\pressure\uwnd.1984.nc','lon');
    
14. u_lat=ncread('I:\NCEPdata\pressure\uwnd.1984.nc','lat');
    
15. 
    
16. U=squeeze(u(43:51,40:44,2,:));
    
17. V=squeeze(v(43:51,40:44,2,:));
    
18. U1=squeeze(u(43:51,40:44,8,:));
    
19. V1=squeeze(v(43:51,40:44,8,:));
    
20. [X,Y]=meshgrid(u_lon(43:51),u_lat(40:44));
    
21. for time=1:365
    
22.     div1(:,:,time)=divh_atmos(Y',X',U(:,:,time),V(:,:,time));%850hap
    
23. end
    
24. for time=1:365
    
25.     div2(:,:,time)=divh_atmos(Y',X',U1(:,:,time),V1(:,:,time));%200hap
    
26. end
    
27. for i=1:365
    
28.    D(i)=nanmean(nanmean((div1(:,:,i)-div2(:,:,i)),1),2);
    
29. end
    
30. for i=2:365
    
31.     Id(i)=D(i)/sqrt(sum(D(1:i).^2)/i);
    
32. end
    
33. Id(1)=0;
    
34. Id(:);
    
35. xlswrite('ceshi.xlsx',Id(:),'sheet1','C2')

复制代码

1. function divh=divh_atmos(longitude, latitude, u, v)

复制代码

1. function dy=dy_atmos(latitude)
   
2.     R=6.3781e6; % earth's radius
   
3.     [~, dy]=gradient(latitude);
   
4.     dy=dy.*(pi./180).*R;
   
5. end

复制代码

1. function dx=dx_atmos(longitude, latitude)
   
2.     R=6.3781e6; % earth's radius
   
3.     [dx, ~]=gradient(longitude);
   
4.     dx=dx.*(pi./180).*R.*cos(latitude*pi./180);
   
5. end

复制代码

                               
登录/注册后可看大图

算出的指数全部都小于0，出图的形状也不相似

膘膘

你的散度如果是月资料求的，就不对。

陈伯英

是算空间平均的地方出了问题吗

陈伯英

膘膘 发表于 2023-3-15 11:22
你的散度如果是月资料求的，就不对。

用的是ncep再分析资料逐天的

下面这个是算散度的程序

1. function divh=divh_atmos(longitude, latitude, u, v)
   
2.     R=6.3781e6; % earth's radius
   
3.     dx=dx_atmos(longitude, latitude);
   
4.     dy=dy_atmos(latitude);
   
5.     [du, ~]=gradient(u);
   
6.     [~, dv]=gradient(v);
   
7.     divh=du./dx+dv./dy-v.*tan(latitude.*pi./180)./R;
   
8.     divh(abs(latitude)==90)=NaN;
   
9. end

复制代码