year=2001:2014;
lat=88.75:-2.5:-88.75;
lon=1.25:2.5:358.75;
month=1:168;

AOD_temp=ncread('/Users/yangyang/Aerosol_Variaiton/IRUN/AOD/ncfile/AOD_MODIS_25x25.nc','AOD_MONTH');

AOD_MONTH=zeros(length(lon),length(lat),168);

AOD_temp(AOD_temp<0)=nan;

% Deseasonal
for i=1:length(lon)
    for j=1:length(lat)
        for monthID=1:12
            temp=nanmean(AOD_temp(i,j,monthID,:),4);
            for yearID=1:14
                AOD_temp(i,j,monthID,yearID)=AOD_temp(i,j,monthID,yearID)-temp;
            end
        end
    end
end

% Change size
for i=1:length(lon)
    for j=1:length(lat)
        n=0;
        for yearID=1:length(year)
            for monthID=1:12
                n=n+1;
                AOD_MONTH(i,j,n)=AOD_temp(i,j,monthID,yearID);
            end
        end
    end
end




%for i=1:length(lon)
%    for j=1:length(lat)
%        for n=1:168
%            if AOD_MONTH(i,j,n)<0
%                AOD_MONTH(i,j,:)=-999;
%            end
%        end
%    end
%end
                



save('AOD_MODIS_monthly_25x25.mat','AOD_MONTH')

nccreate('AOD_MODIS_monthly_25x25.nc','Month','Dimensions',{'time',Inf},'Format','classic')
ncwrite('AOD_MODIS_monthly_25x25.nc','Month', month)
ncwriteatt('AOD_MODIS_monthly_25x25.nc','Month','long_name','Month')
nccreate('AOD_MODIS_monthly_25x25.nc','lat','Dimensions',{'lat',length(lat)},'Format','classic')
ncwrite('AOD_MODIS_monthly_25x25.nc','lat', lat)
ncwriteatt('AOD_MODIS_monthly_25x25.nc','lat','units','degrees_north')
ncwriteatt('AOD_MODIS_monthly_25x25.nc','lat','long_name','Latitude')
nccreate('AOD_MODIS_monthly_25x25.nc','lon','Dimensions',{'lon',length(lon)},'Format','classic')
ncwrite('AOD_MODIS_monthly_25x25.nc','lon', lon)
ncwriteatt('AOD_MODIS_monthly_25x25.nc','lon','units','degrees_east')
ncwriteatt('AOD_MODIS_monthly_25x25.nc','lon','long_name','Longitude')
nccreate('AOD_MODIS_monthly_25x25.nc','AOD_MONTH','Dimensions',{'lon',length(lon),'lat',length(lat),'month',length(month)},'Format','classic')
ncwrite('AOD_MODIS_monthly_25x25.nc','AOD_MONTH', AOD_MONTH)
ncdisp('AOD_MODIS_monthly_25x25.nc')