year=2001:2014;
lat=-89.5:1:89.5;
lon=-179.5:1:179.5;
month={'01','02','03','04','05','06','07','08','09','10','11','12'};

AOD_MONTH_temp=zeros(180,360,length(month),length(year));
AOD_MONTH=zeros(360,180,length(month),length(year));

AOD_ANN_temp=zeros(180,360,length(year));
AOD_ANN=zeros(360,180,length(year));

for yearID=1:length(year)
    for monthID=1:length(month)
        AOD_MONTH_temp(:,:,monthID,yearID)=hdfread(['/Volumes/WORK5T/MODIS/',num2str(year(yearID)),'/MOD08_M3.A',num2str(year(yearID)),cell2mat(month(monthID)),'.hdf'], '/mod08/Data Fields/Aerosol_Optical_Depth_Land_Ocean_Mean_Mean', 'Index', {[1  1],[1  1],[180  360]});
    end
end


for i=1:360
    for j=1:180
        AOD_MONTH(i,j,:,:)=AOD_MONTH_temp(181-j,i,:,:)*0.001;
    end
end


for i=1:360
    for j=1:180
        for yearID=1:length(year)
            
            n=0;
            for monthID=1:length(month)
                if AOD_MONTH(i,j,monthID,yearID)>0
                    n=n+1;
                    AOD_ANN(i,j,yearID)=AOD_ANN(i,j,yearID)+AOD_MONTH(i,j,monthID,yearID);
                else
                    AOD_MONTH(i,j,monthID,yearID)=-999;
                end
            end
            if n==0
                AOD_ANN(i,j,yearID)=-999;
            else
                AOD_ANN(i,j,yearID)=AOD_ANN(i,j,yearID)/n;
            end
        
        end
    end
end

AOD_DJF=zeros(360,180,length(year));

AOD_200012_temp=zeros(180,360);
AOD_200012_temp(:,:)=hdfread('/Volumes/WORK5T/MODIS/2000/MOD08_M3.A200012.hdf', '/mod08/Data Fields/Aerosol_Optical_Depth_Land_Ocean_Mean_Mean', 'Index', {[1  1],[1  1],[180  360]});
AOD_200012=zeros(360,180);
for i=1:360
    for j=1:180
        AOD_200012(i,j)=AOD_200012_temp(181-j,i)*0.001;
    end
end


for i=1:360
    for j=1:180
        
        yearID=1;
        n=0;
        if AOD_200012(i,j)>0
            n=n+1;
            AOD_DJF(i,j,yearID)=AOD_DJF(i,j,yearID)+AOD_200012(i,j);
        end
        if AOD_MONTH(i,j,1,yearID)>0
            n=n+1;
            AOD_DJF(i,j,yearID)=AOD_DJF(i,j,yearID)+AOD_MONTH(i,j,1,yearID);
        end
        if AOD_MONTH(i,j,2,yearID)>0
            n=n+1;
            AOD_DJF(i,j,yearID)=AOD_DJF(i,j,yearID)+AOD_MONTH(i,j,2,yearID);
        end
            
        if n==0
            AOD_DJF(i,j,yearID)=-999;
        else
            AOD_DJF(i,j,yearID)=AOD_DJF(i,j,yearID)/n;
        end
        
        
        for yearID=2:length(year)
            
            n=0;
            if AOD_MONTH(i,j,12,yearID-1)>0
                n=n+1;
                AOD_DJF(i,j,yearID)=AOD_DJF(i,j,yearID)+AOD_MONTH(i,j,12,yearID-1);
            end
            if AOD_MONTH(i,j,1,yearID)>0
                n=n+1;
                AOD_DJF(i,j,yearID)=AOD_DJF(i,j,yearID)+AOD_MONTH(i,j,1,yearID);
            end
            if AOD_MONTH(i,j,2,yearID)>0
                n=n+1;
                AOD_DJF(i,j,yearID)=AOD_DJF(i,j,yearID)+AOD_MONTH(i,j,2,yearID);
            end
            
            if n==0
                AOD_DJF(i,j,yearID)=-999;
            else
                AOD_DJF(i,j,yearID)=AOD_DJF(i,j,yearID)/n;
            end
            
        end
        
        
        
    end
end


save('AOD_MODIS.mat','AOD_MONTH','AOD_ANN','lat','lon','year')

nccreate('AOD_MODIS.nc','Year','Dimensions',{'time',Inf},'Format','classic')
ncwrite('AOD_MODIS.nc','Year', year)
nccreate('AOD_MODIS.nc','lat','Dimensions',{'lat',180},'Format','classic')
ncwrite('AOD_MODIS.nc','lat', lat)
ncwriteatt('AOD_MODIS.nc','lat','units','degrees_north')
ncwriteatt('AOD_MODIS.nc','lat','long_name','Latitude')
nccreate('AOD_MODIS.nc','lon','Dimensions',{'lon',360},'Format','classic')
ncwrite('AOD_MODIS.nc','lon', lon)
ncwriteatt('AOD_MODIS.nc','lon','units','degrees_east')
ncwriteatt('AOD_MODIS.nc','lon','long_name','Longitude')
nccreate('AOD_MODIS.nc','AOD_MONTH','Dimensions',{'lon',360,'lat',180,'month',12,'year',length(year)},'Format','classic')
ncwrite('AOD_MODIS.nc','AOD_MONTH', AOD_MONTH)
nccreate('AOD_MODIS.nc','AOD_ANN','Dimensions',{'lon',360,'lat',180,'year',length(year)},'Format','classic')
ncwrite('AOD_MODIS.nc','AOD_ANN', AOD_ANN)
ncdisp('AOD_MODIS.nc')


