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

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

QRS=ncread('/Users/yangyang/Aerosol_Variaiton/DRUN/DTSF/ncfile/DTSF.nc','QRS_DJF');
QRS_ECHINA=zeros(length(year),1);
QRS_ECHINA(:)=mean(mean(QRS(45:50,60:72,:),2),1);
QRS_weak_temp=zeros(144,96,sum(WIND_year));
QRS_ECHINA_weak=0;
n=0;
for yearID=1:length(year)
    if WIND_year(yearID)==1
        n=n+1;
        QRS_weak_temp(:,:,n)=QRS(:,:,yearID);
        QRS_ECHINA_weak=QRS_ECHINA_weak+QRS_ECHINA(yearID)/sum(WIND_year);
    end
end
QRS_weak=mean(QRS_weak_temp,3);
QRS_aveg=mean(QRS,3);
QRS_diff=QRS_weak-QRS_aveg;
QRS_ECHINA_aveg=mean(QRS_ECHINA);
QRS_ECHINA_diff=QRS_ECHINA_weak-QRS_ECHINA_aveg;

QRL=ncread('/Users/yangyang/Aerosol_Variaiton/DRUN/DTSF/ncfile/DTSF.nc','QRL_DJF');
QRL_ECHINA=zeros(length(year),1);
QRL_ECHINA(:)=mean(mean(QRL(45:50,60:72,:),2),1);
QRL_weak_temp=zeros(144,96,sum(WIND_year));
QRL_ECHINA_weak=0;
n=0;
for yearID=1:length(year)
    if WIND_year(yearID)==1
        n=n+1;
        QRL_weak_temp(:,:,n)=QRL(:,:,yearID);
        QRL_ECHINA_weak=QRL_ECHINA_weak+QRL_ECHINA(yearID)/sum(WIND_year);
    end
end
QRL_weak=mean(QRL_weak_temp,3);
QRL_aveg=mean(QRL,3);
QRL_diff=QRL_weak-QRL_aveg;
QRL_ECHINA_aveg=mean(QRL_ECHINA);
QRL_ECHINA_diff=QRL_ECHINA_weak-QRL_ECHINA_aveg;

DTV=ncread('/Users/yangyang/Aerosol_Variaiton/DRUN/DTSF/ncfile/DTSF.nc','DTV_DJF');
DTV_ECHINA=zeros(length(year),1);
DTV_ECHINA(:)=mean(mean(DTV(45:50,60:72,:),2),1);
DTV_weak_temp=zeros(144,96,sum(WIND_year));
DTV_ECHINA_weak=0;
n=0;
for yearID=1:length(year)
    if WIND_year(yearID)==1
        n=n+1;
        DTV_weak_temp(:,:,n)=DTV(:,:,yearID);
        DTV_ECHINA_weak=DTV_ECHINA_weak+DTV_ECHINA(yearID)/sum(WIND_year);
    end
end
DTV_weak=mean(DTV_weak_temp,3);
DTV_aveg=mean(DTV,3);
DTV_diff=DTV_weak-DTV_aveg;
DTV_ECHINA_aveg=mean(DTV_ECHINA);
DTV_ECHINA_diff=DTV_ECHINA_weak-DTV_ECHINA_aveg;

DTCOND=ncread('/Users/yangyang/Aerosol_Variaiton/DRUN/DTSF/ncfile/DTSF.nc','DTCOND_DJF');
DTCOND_ECHINA=zeros(length(year),1);
DTCOND_ECHINA(:)=mean(mean(DTCOND(45:50,60:72,:),2),1);
DTCOND_weak_temp=zeros(144,96,sum(WIND_year));
DTCOND_ECHINA_weak=0;
n=0;
for yearID=1:length(year)
    if WIND_year(yearID)==1
        n=n+1;
        DTCOND_weak_temp(:,:,n)=DTCOND(:,:,yearID);
        DTCOND_ECHINA_weak=DTCOND_ECHINA_weak+DTCOND_ECHINA(yearID)/sum(WIND_year);
    end
end
DTCOND_weak=mean(DTCOND_weak_temp,3);
DTCOND_aveg=mean(DTCOND,3);
DTCOND_diff=DTCOND_weak-DTCOND_aveg;
DTCOND_ECHINA_aveg=mean(DTCOND_ECHINA);
DTCOND_ECHINA_diff=DTCOND_ECHINA_weak-DTCOND_ECHINA_aveg;

DTSF_diff=QRS_diff+QRL_diff+DTV_diff+DTCOND_diff;

save('DTSF_10th.mat','QRS_diff','QRL_diff','DTV_diff','DTCOND_diff','DTSF_diff')

nccreate('DTSF_10th.nc','lat','Dimensions',{'lat',96},'Format','classic')
ncwrite('DTSF_10th.nc','lat', lat)
ncwriteatt('DTSF_10th.nc','lat','units','degrees_north')
ncwriteatt('DTSF_10th.nc','lat','long_name','Latitude')
nccreate('DTSF_10th.nc','lon','Dimensions',{'lon',144},'Format','classic')
ncwrite('DTSF_10th.nc','lon', lon)
ncwriteatt('DTSF_10th.nc','lon','units','degrees_east')
ncwriteatt('DTSF_10th.nc','lon','long_name','Longitude')
nccreate('DTSF_10th.nc','DTSF_diff','Dimensions',{'lon',144,'lat',96},'Format','classic')
ncwrite('DTSF_10th.nc','DTSF_diff', DTSF_diff)
nccreate('DTSF_10th.nc','QRS_diff','Dimensions',{'lon',144,'lat',96},'Format','classic')
ncwrite('DTSF_10th.nc','QRS_diff', QRS_diff)
nccreate('DTSF_10th.nc','QRL_diff','Dimensions',{'lon',144,'lat',96},'Format','classic')
ncwrite('DTSF_10th.nc','QRL_diff', QRL_diff)
nccreate('DTSF_10th.nc','DTV_diff','Dimensions',{'lon',144,'lat',96},'Format','classic')
ncwrite('DTSF_10th.nc','DTV_diff', DTV_diff)
nccreate('DTSF_10th.nc','DTCOND_diff','Dimensions',{'lon',144,'lat',96},'Format','classic')
ncwrite('DTSF_10th.nc','DTCOND_diff', DTCOND_diff)
ncdisp('DTSF_10th.nc')