year=2003:2015;
year2=2003:2016;
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/disk/MODIS/MODIS_Aqua/',num2str(year(yearID)),'/MYD08_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



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

AOD_200212_temp=zeros(180,360);
AOD_200212_temp(:,:)=hdfread('/Volumes/disk/MODIS/MODIS_Aqua/2002/MYD08_M3.A200212.hdf', '/mod08/Data Fields/Aerosol_Optical_Depth_Land_Ocean_Mean_Mean', 'Index', {[1  1],[1  1],[180  360]});
AOD_200212=zeros(360,180);
for i=1:360
    for j=1:180
        AOD_200212(i,j)=AOD_200212_temp(181-j,i)*0.001;
    end
end

AOD_201601_temp=zeros(180,360);
AOD_201601_temp(:,:)=hdfread('/Volumes/disk/MODIS/MODIS_Aqua/2016/MYD08_M3.A201601.hdf', '/mod08/Data Fields/Aerosol_Optical_Depth_Land_Ocean_Mean_Mean', 'Index', {[1  1],[1  1],[180  360]});
AOD_201601=zeros(360,180);
for i=1:360
    for j=1:180
        AOD_201601(i,j)=AOD_201601_temp(181-j,i)*0.001;
    end
end

AOD_201602_temp=zeros(180,360);
AOD_201602_temp(:,:)=hdfread('/Volumes/disk/MODIS/MODIS_Aqua/2016/MYD08_M3.A201602.hdf', '/mod08/Data Fields/Aerosol_Optical_Depth_Land_Ocean_Mean_Mean', 'Index', {[1  1],[1  1],[180  360]});
AOD_201602=zeros(360,180);
for i=1:360
    for j=1:180
        AOD_201602(i,j)=AOD_201602_temp(181-j,i)*0.001;
    end
end


for i=1:360
    for j=1:180
        
        yearID=1;
        n=0;
        if AOD_200212(i,j)>0
            n=n+1;
            AOD_DJF(i,j,yearID)=AOD_DJF(i,j,yearID)+AOD_200212(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
        
        
        yearID=yearID+1;
        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_201601(i,j)>0
            n=n+1;
            AOD_DJF(i,j,yearID)=AOD_DJF(i,j,yearID)+AOD_201601(i,j);
        end
        if AOD_201602(i,j)>0
            n=n+1;
            AOD_DJF(i,j,yearID)=AOD_DJF(i,j,yearID)+AOD_201602(i,j);
        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




AOD_fine_MONTH_temp=zeros(180,360,length(month),length(year));
AOD_fine_MONTH=zeros(360,180,length(month),length(year));

AOD_fine_ANN_temp=zeros(180,360,length(year));
AOD_fine_ANN=zeros(360,180,length(year));

for yearID=1:length(year)
    for monthID=1:length(month)
        temp=hdfread(['/Volumes/disk/MODIS/MODIS_Aqua/',num2str(year(yearID)),'/MYD08_M3.A',num2str(year(yearID)),cell2mat(month(monthID)),'.hdf'], '/mod08/Data Fields/Aerosol_Optical_Depth_Small_Ocean_Mean_Mean');
        AOD_fine_MONTH_temp(:,:,monthID,yearID)=temp(2,:,:);
    end
end


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



AOD_fine_DJF=zeros(360,180,length(year)+1);

temp=hdfread('/Volumes/disk/MODIS/MODIS_Aqua/2002/MYD08_M3.A200212.hdf', '/mod08/Data Fields/Aerosol_Optical_Depth_Small_Ocean_Mean_Mean');
AOD_fine_200212_temp=zeros(180,360);
AOD_fine_200212_temp(:,:)=temp(2,:,:);
AOD_fine_200212=zeros(360,180);
for i=1:360
    for j=1:180
        AOD_fine_200212(i,j)=AOD_fine_200212_temp(181-j,i)*0.001;
    end
end

temp=hdfread('/Volumes/disk/MODIS/MODIS_Aqua/2016/MYD08_M3.A201601.hdf', '/mod08/Data Fields/Aerosol_Optical_Depth_Small_Ocean_Mean_Mean');
AOD_fine_201601_temp=zeros(180,360);
AOD_fine_201601_temp(:,:)=temp(2,:,:);
AOD_fine_201601=zeros(360,180);
for i=1:360
    for j=1:180
        AOD_fine_201601(i,j)=AOD_fine_201601_temp(181-j,i)*0.001;
    end
end

temp=hdfread('/Volumes/disk/MODIS/MODIS_Aqua/2016/MYD08_M3.A201602.hdf', '/mod08/Data Fields/Aerosol_Optical_Depth_Small_Ocean_Mean_Mean');
AOD_fine_201602_temp=zeros(180,360);
AOD_fine_201602_temp(:,:)=temp(2,:,:);
AOD_fine_201602=zeros(360,180);
for i=1:360
    for j=1:180
        AOD_fine_201602(i,j)=AOD_fine_201602_temp(181-j,i)*0.001;
    end
end


for i=1:360
    for j=1:180
        
        yearID=1;
        n=0;
        if AOD_fine_200212(i,j)>0
            n=n+1;
            AOD_fine_DJF(i,j,yearID)=AOD_fine_DJF(i,j,yearID)+AOD_fine_200212(i,j);
        end
        if AOD_fine_MONTH(i,j,1,yearID)>0
            n=n+1;
            AOD_fine_DJF(i,j,yearID)=AOD_fine_DJF(i,j,yearID)+AOD_fine_MONTH(i,j,1,yearID);
        end
        if AOD_fine_MONTH(i,j,2,yearID)>0
            n=n+1;
            AOD_fine_DJF(i,j,yearID)=AOD_fine_DJF(i,j,yearID)+AOD_fine_MONTH(i,j,2,yearID);
        end
            
        if n==0
            AOD_fine_DJF(i,j,yearID)=-999;
        else
            AOD_fine_DJF(i,j,yearID)=AOD_fine_DJF(i,j,yearID)/n;
        end
        
        
        for yearID=2:length(year)
            
            n=0;
            if AOD_fine_MONTH(i,j,12,yearID-1)>0
                n=n+1;
                AOD_fine_DJF(i,j,yearID)=AOD_fine_DJF(i,j,yearID)+AOD_fine_MONTH(i,j,12,yearID-1);
            end
            if AOD_fine_MONTH(i,j,1,yearID)>0
                n=n+1;
                AOD_fine_DJF(i,j,yearID)=AOD_fine_DJF(i,j,yearID)+AOD_fine_MONTH(i,j,1,yearID);
            end
            if AOD_fine_MONTH(i,j,2,yearID)>0
                n=n+1;
                AOD_fine_DJF(i,j,yearID)=AOD_fine_DJF(i,j,yearID)+AOD_fine_MONTH(i,j,2,yearID);
            end
            
            if n==0
                AOD_fine_DJF(i,j,yearID)=-999;
            else
                AOD_fine_DJF(i,j,yearID)=AOD_fine_DJF(i,j,yearID)/n;
            end
            
        end
        
        
        yearID=yearID+1;
        n=0;
        if AOD_fine_MONTH(i,j,12,yearID-1)>0
                n=n+1;
                AOD_fine_DJF(i,j,yearID)=AOD_fine_DJF(i,j,yearID)+AOD_fine_MONTH(i,j,12,yearID-1);
        end
        if AOD_fine_201601(i,j)>0
            n=n+1;
            AOD_fine_DJF(i,j,yearID)=AOD_fine_DJF(i,j,yearID)+AOD_fine_201601(i,j);
        end
        if AOD_fine_201602(i,j)>0
            n=n+1;
            AOD_fine_DJF(i,j,yearID)=AOD_fine_DJF(i,j,yearID)+AOD_fine_201602(i,j);
        end
            
        if n==0
            AOD_fine_DJF(i,j,yearID)=-999;
        else
            AOD_fine_DJF(i,j,yearID)=AOD_fine_DJF(i,j,yearID)/n;
        end
        
        
        
    end
end


AOD_coarse_DJF=AOD_DJF-AOD_fine_DJF;

for i=1:360
    for j=1:180
        for yearID=1:length(year2)
            if AOD_DJF(i,j,yearID)<0 || AOD_fine_DJF(i,j,yearID)<0
                AOD_coarse_DJF(i,j,yearID)=-999;
            end
        end
    end
end


save('AOD_MODIS_Aqua_2003_2016.mat','AOD_DJF','AOD_fine_DJF','AOD_coarse_DJF','lat','lon','year2')

nccreate('AOD_MODIS_Aqua_2003_2016.nc','Year','Dimensions',{'time',Inf},'Format','classic')
ncwrite('AOD_MODIS_Aqua_2003_2016.nc','Year', year2)
nccreate('AOD_MODIS_Aqua_2003_2016.nc','lat','Dimensions',{'lat',180},'Format','classic')
ncwrite('AOD_MODIS_Aqua_2003_2016.nc','lat', lat)
ncwriteatt('AOD_MODIS_Aqua_2003_2016.nc','lat','units','degrees_north')
ncwriteatt('AOD_MODIS_Aqua_2003_2016.nc','lat','long_name','Latitude')
nccreate('AOD_MODIS_Aqua_2003_2016.nc','lon','Dimensions',{'lon',360},'Format','classic')
ncwrite('AOD_MODIS_Aqua_2003_2016.nc','lon', lon)
ncwriteatt('AOD_MODIS_Aqua_2003_2016.nc','lon','units','degrees_east')
ncwriteatt('AOD_MODIS_Aqua_2003_2016.nc','lon','long_name','Longitude')
nccreate('AOD_MODIS_Aqua_2003_2016.nc','AOD_DJF','Dimensions',{'lon',360,'lat',180,'year',length(year2)},'Format','classic')
ncwrite('AOD_MODIS_Aqua_2003_2016.nc','AOD_DJF', AOD_DJF)
nccreate('AOD_MODIS_Aqua_2003_2016.nc','AOD_fine_DJF','Dimensions',{'lon',360,'lat',180,'year',length(year2)},'Format','classic')
ncwrite('AOD_MODIS_Aqua_2003_2016.nc','AOD_fine_DJF', AOD_fine_DJF)
nccreate('AOD_MODIS_Aqua_2003_2016.nc','AOD_coarse_DJF','Dimensions',{'lon',360,'lat',180,'year',length(year2)},'Format','classic')
ncwrite('AOD_MODIS_Aqua_2003_2016.nc','AOD_coarse_DJF', AOD_coarse_DJF)
ncdisp('AOD_MODIS_Aqua_2003_2016.nc')


