load ~/NC/grid
year=1:150;

WIND_year=ncread('/Users/yangyang/Aerosol_Variaiton/IRUN/WIND/ncfile/WIND_10th.nc','WIND_year');
TS=ncread('/Users/yangyang/Aerosol_Variaiton/IRUN/TS/ncfile/TS.nc','TS_DJF');

TS_ECHINA=zeros(length(year),1);
TS_ECHINA(:)=mean(mean(TS(45:50,60:72,:),2),1);


TS_weak_temp=zeros(144,96,sum(WIND_year));
TS_ECHINA_weak=0;
n=0;
for yearID=1:length(year)
    if WIND_year(yearID)==1
        n=n+1;
        TS_weak_temp(:,:,n)=TS(:,:,yearID);
        TS_ECHINA_weak=TS_ECHINA_weak+TS_ECHINA(yearID)/sum(WIND_year);
    end
end
TS_weak=mean(TS_weak_temp,3);


TS_aveg=mean(TS,3);
TS_diff=TS_weak-TS_aveg;

TS_ECHINA_aveg=mean(TS_ECHINA);
TS_ECHINA_diff=TS_ECHINA_weak-TS_ECHINA_aveg;

save('TS_10th.mat','TS_weak','TS_aveg','TS_diff')

nccreate('TS_10th.nc','lat','Dimensions',{'lat',96},'Format','classic')
ncwrite('TS_10th.nc','lat', lat)
ncwriteatt('TS_10th.nc','lat','units','degrees_north')
ncwriteatt('TS_10th.nc','lat','long_name','Latitude')
nccreate('TS_10th.nc','lon','Dimensions',{'lon',144},'Format','classic')
ncwrite('TS_10th.nc','lon', lon)
ncwriteatt('TS_10th.nc','lon','units','degrees_east')
ncwriteatt('TS_10th.nc','lon','long_name','Longitude')
nccreate('TS_10th.nc','TS_weak','Dimensions',{'lon',144,'lat',96},'Format','classic')
ncwrite('TS_10th.nc','TS_weak', TS_weak)
nccreate('TS_10th.nc','TS_aveg','Dimensions',{'lon',144,'lat',96},'Format','classic')
ncwrite('TS_10th.nc','TS_aveg', TS_aveg)
nccreate('TS_10th.nc','TS_diff','Dimensions',{'lon',144,'lat',96},'Format','classic')
ncwrite('TS_10th.nc','TS_diff', TS_diff)
ncdisp('TS_10th.nc')