year=1986:2006;
lon=-180:2.5:177.5;
lat=-90:2:90;

WIND_day=ncread('/Users/yangyang/Aerosol_Variaiton/GEOS_Chem/ncfile/WIND_10th.nc','WIND_day');
PM25_temp=ncread('/Users/yangyang/Extreme_Aerosol/MET/ncfile/Aerosol_daily.nc','pm25');

PM25=zeros(144,91,1890);
for yearID=1:length(year)
    PM25(:,:,yearID*90-89:yearID*90)=PM25_temp(:,:,:,yearID);
end

PM25_ECHINA=zeros(1890,1);
PM25_ECHINA(:)=mean(mean(PM25(117:122,57:68,:),2),1);


PM25_weak_temp=zeros(144,91,sum(WIND_day));
PM25_ECHINA_weak=0;
n=0;
for dayID=1:1890
    if WIND_day(dayID)==1
        n=n+1;
        PM25_weak_temp(:,:,n)=PM25(:,:,dayID);
        PM25_ECHINA_weak=PM25_ECHINA_weak+PM25_ECHINA(dayID)/sum(WIND_day);
    end
end
PM25_weak=mean(PM25_weak_temp,3);


PM25_aveg=mean(PM25,3);
PM25_diff=PM25_weak-PM25_aveg;

PM25_ECHINA_aveg=mean(PM25_ECHINA);
PM25_ECHINA_diff=PM25_ECHINA_weak-PM25_ECHINA_aveg;

save('PM25_10th.mat','PM25_weak','PM25_aveg','PM25_diff')

nccreate('PM25_10th.nc','lat','Dimensions',{'lat',91},'Format','classic')
ncwrite('PM25_10th.nc','lat', lat)
ncwriteatt('PM25_10th.nc','lat','units','degrees_north')
ncwriteatt('PM25_10th.nc','lat','long_name','Latitude')
nccreate('PM25_10th.nc','lon','Dimensions',{'lon',144},'Format','classic')
ncwrite('PM25_10th.nc','lon', lon)
ncwriteatt('PM25_10th.nc','lon','units','degrees_east')
ncwriteatt('PM25_10th.nc','lon','long_name','Longitude')
nccreate('PM25_10th.nc','PM25_weak','Dimensions',{'lon',144,'lat',91},'Format','classic')
ncwrite('PM25_10th.nc','PM25_weak', PM25_weak)
nccreate('PM25_10th.nc','PM25_aveg','Dimensions',{'lon',144,'lat',91},'Format','classic')
ncwrite('PM25_10th.nc','PM25_aveg', PM25_aveg)
nccreate('PM25_10th.nc','PM25_diff','Dimensions',{'lon',144,'lat',91},'Format','classic')
ncwrite('PM25_10th.nc','PM25_diff', PM25_diff)
ncdisp('PM25_10th.nc')