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

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

QRS=ncread('/Users/yangyang/Aerosol_Variaiton/IRUN/DTVERT/ncfile/DTVERT.nc','QRS_DJF');
QRS_weak_temp=zeros(96,30,sum(WIND_year));
n=0;
for yearID=1:length(year)
    if WIND_year(yearID)==1
        n=n+1;
        QRS_weak_temp(:,:,n)=QRS(:,:,yearID);
    end
end
QRS_weak=mean(QRS_weak_temp,3);
QRS_aveg=mean(QRS,3);
QRS_diff=QRS_weak-QRS_aveg;

QRL=ncread('/Users/yangyang/Aerosol_Variaiton/IRUN/DTVERT/ncfile/DTVERT.nc','QRL_DJF');
QRL_weak_temp=zeros(96,30,sum(WIND_year));
n=0;
for yearID=1:length(year)
    if WIND_year(yearID)==1
        n=n+1;
        QRL_weak_temp(:,:,n)=QRL(:,:,yearID);
       
    end
end
QRL_weak=mean(QRL_weak_temp,3);
QRL_aveg=mean(QRL,3);
QRL_diff=QRL_weak-QRL_aveg;

DTV=ncread('/Users/yangyang/Aerosol_Variaiton/IRUN/DTVERT/ncfile/DTVERT.nc','DTV_DJF');
DTV_weak_temp=zeros(96,30,sum(WIND_year));
n=0;
for yearID=1:length(year)
    if WIND_year(yearID)==1
        n=n+1;
        DTV_weak_temp(:,:,n)=DTV(:,:,yearID);
    end
end
DTV_weak=mean(DTV_weak_temp,3);
DTV_aveg=mean(DTV,3);
DTV_diff=DTV_weak-DTV_aveg;

DTCOND=ncread('/Users/yangyang/Aerosol_Variaiton/IRUN/DTVERT/ncfile/DTVERT.nc','DTCOND_DJF');
DTCOND_weak_temp=zeros(96,30,sum(WIND_year));
n=0;
for yearID=1:length(year)
    if WIND_year(yearID)==1
        n=n+1;
        DTCOND_weak_temp(:,:,n)=DTCOND(:,:,yearID);
    end
end
DTCOND_weak=mean(DTCOND_weak_temp,3);
DTCOND_aveg=mean(DTCOND,3);
DTCOND_diff=DTCOND_weak-DTCOND_aveg;

DTVERT_diff=QRS_diff+QRL_diff+DTV_diff+DTCOND_diff;

save('DTVERT_10th.mat','QRS_diff','QRL_diff','DTV_diff','DTCOND_diff','DTVERT_diff')

nccreate('DTVERT_10th.nc','lat','Dimensions',{'lat',96},'Format','classic')
ncwrite('DTVERT_10th.nc','lat', lat)
ncwriteatt('DTVERT_10th.nc','lat','units','degrees_north')
ncwriteatt('DTVERT_10th.nc','lat','long_name','Latitude')
nccreate('DTVERT_10th.nc','lev','Dimensions',{'lev',30},'Format','classic')
ncwrite('DTVERT_10th.nc','lev', lev)
nccreate('DTVERT_10th.nc','DTVERT_diff','Dimensions',{'lat',96,'lev',30},'Format','classic')
ncwrite('DTVERT_10th.nc','DTVERT_diff', DTVERT_diff)
nccreate('DTVERT_10th.nc','QRS_diff','Dimensions',{'lat',96,'lev',30},'Format','classic')
ncwrite('DTVERT_10th.nc','QRS_diff', QRS_diff)
nccreate('DTVERT_10th.nc','QRL_diff','Dimensions',{'lat',96,'lev',30},'Format','classic')
ncwrite('DTVERT_10th.nc','QRL_diff', QRL_diff)
nccreate('DTVERT_10th.nc','DTV_diff','Dimensions',{'lat',96,'lev',30},'Format','classic')
ncwrite('DTVERT_10th.nc','DTV_diff', DTV_diff)
nccreate('DTVERT_10th.nc','DTCOND_diff','Dimensions',{'lat',96,'lev',30},'Format','classic')
ncwrite('DTVERT_10th.nc','DTCOND_diff', DTCOND_diff)
ncdisp('DTVERT_10th.nc')