year=2001:2015;
year2=2001: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'};

AODBD_MONTH_temp=zeros(180,360,length(month),length(year));
AODBD_MONTH=zeros(360,180,length(month),length(year));

AODBD_ANN_temp=zeros(180,360,length(year));
AODBD_ANN=zeros(360,180,length(year));

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


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



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

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

AODBD_201601_temp=zeros(180,360);
AODBD_201601_temp(:,:)=hdfread('/Volumes/disk/MODIS/MODIS_Terra/2016/MOD08_M3.A201601.hdf', '/mod08/Data Fields/Deep_Blue_Aerosol_Optical_Depth_550_Land_Mean_Mean', 'Index', {[1  1],[1  1],[180  360]});
AODBD_201601=zeros(360,180);
for i=1:360
    for j=1:180
        AODBD_201601(i,j)=AODBD_201601_temp(181-j,i)*0.001;
    end
end

AODBD_201602_temp=zeros(180,360);
AODBD_201602_temp(:,:)=hdfread('/Volumes/disk/MODIS/MODIS_Terra/2016/MOD08_M3.A201602.hdf', '/mod08/Data Fields/Deep_Blue_Aerosol_Optical_Depth_550_Land_Mean_Mean', 'Index', {[1  1],[1  1],[180  360]});
AODBD_201602=zeros(360,180);
for i=1:360
    for j=1:180
        AODBD_201602(i,j)=AODBD_201602_temp(181-j,i)*0.001;
    end
end


for i=1:360
    for j=1:180
        
        yearID=1;
        n=0;
        if AODBD_200012(i,j)>0
            n=n+1;
            AODBD_DJF(i,j,yearID)=AODBD_DJF(i,j,yearID)+AODBD_200012(i,j);
        end
        if AODBD_MONTH(i,j,1,yearID)>0
            n=n+1;
            AODBD_DJF(i,j,yearID)=AODBD_DJF(i,j,yearID)+AODBD_MONTH(i,j,1,yearID);
        end
        if AODBD_MONTH(i,j,2,yearID)>0
            n=n+1;
            AODBD_DJF(i,j,yearID)=AODBD_DJF(i,j,yearID)+AODBD_MONTH(i,j,2,yearID);
        end
            
        if n==0
            AODBD_DJF(i,j,yearID)=-999;
        else
            AODBD_DJF(i,j,yearID)=AODBD_DJF(i,j,yearID)/n;
        end
        
        
        for yearID=2:length(year)
            
            n=0;
            if AODBD_MONTH(i,j,12,yearID-1)>0
                n=n+1;
                AODBD_DJF(i,j,yearID)=AODBD_DJF(i,j,yearID)+AODBD_MONTH(i,j,12,yearID-1);
            end
            if AODBD_MONTH(i,j,1,yearID)>0
                n=n+1;
                AODBD_DJF(i,j,yearID)=AODBD_DJF(i,j,yearID)+AODBD_MONTH(i,j,1,yearID);
            end
            if AODBD_MONTH(i,j,2,yearID)>0
                n=n+1;
                AODBD_DJF(i,j,yearID)=AODBD_DJF(i,j,yearID)+AODBD_MONTH(i,j,2,yearID);
            end
            
            if n==0
                AODBD_DJF(i,j,yearID)=-999;
            else
                AODBD_DJF(i,j,yearID)=AODBD_DJF(i,j,yearID)/n;
            end
            
        end
        
        
        yearID=yearID+1;
        n=0;
        if AODBD_MONTH(i,j,12,yearID-1)>0
                n=n+1;
                AODBD_DJF(i,j,yearID)=AODBD_DJF(i,j,yearID)+AODBD_MONTH(i,j,12,yearID-1);
        end
        if AODBD_201601(i,j)>0
            n=n+1;
            AODBD_DJF(i,j,yearID)=AODBD_DJF(i,j,yearID)+AODBD_201601(i,j);
        end
        if AODBD_201602(i,j)>0
            n=n+1;
            AODBD_DJF(i,j,yearID)=AODBD_DJF(i,j,yearID)+AODBD_201602(i,j);
        end
            
        if n==0
            AODBD_DJF(i,j,yearID)=-999;
        else
            AODBD_DJF(i,j,yearID)=AODBD_DJF(i,j,yearID)/n;
        end
        
        
        
    end
end





save('AODBD_MODIS_2001_2016.mat','AODBD_DJF','lat','lon','year2')

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


