MAKEPROXY;
clear all; close all;
load /rootdir/s1/directory1/hadcru3v/infillannt_3d/instfilled_3d
%load /rootdir/s1/directory1/hadcru3v/enough_3d
%load /rootdir/s1/directory1/hadcru3v/toofew_3d

% standardize sst/sat based on interval 1850-1995 for the full series(1850-2006)
keepers=instfilled_3d(47:146,:);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
means=nanmean(keepers(1:100,:));
stdevs=nanstd(keepers(1:100,:));
stdkeepers=keepers-repmat(means,100,1);
stdkeepers=stdkeepers./repmat(stdevs,100,1);
[U,S,V] = svd(stdkeepers,0); % do the decomposition, economy style. U's are the PCs
instpcs=U;
insteofs=V; insteofs=insteofs';
instsingvals=diag(S);
instpcmeans=nanmean(instpcs);
instpcs=instpcs-repmat(instpcmeans,100,1);
instpcstds=nanstd(instpcs(:,:));
instpcs=instpcs./repmat(instpcstds,100,1);

save('/rootdir/s1/directory1/work1/temanncfr/input/instmeans','means','-ASCII');
save('/rootdir/s1/directory1/work1/temanncfr/input/inststds','stdevs','-ASCII');

save('/rootdir/s1/directory1/work1/temanncfr/input/instpcmeans','instpcmeans','-ASCII');
save('/rootdir/s1/directory1/work1/temanncfr/input/instpcstds','instpcstds','-ASCII');
save('/rootdir/s1/directory1/work1/temanncfr/input/u','U','-ASCII');
save('/rootdir/s1/directory1/work1/temanncfr/input/v','V','-ASCII');
save('/rootdir/s1/directory1/work1/temanncfr/input/s','S','-ASCII');

% get high-f and low-f instpcs
stdkeepers=instpcs;
[nrows,ncols]=size(stdkeepers)
    cutoff=20; %years
	lowlim=1./(cutoff./2)

[nrows,ncols]=size(stdkeepers) 
gotstart=zeros(1,ncols);
gotstop=zeros(1,ncols);
start=zeros(1,ncols);
stop=zeros(1,ncols);
stop(:)=nrows;
[lowa,lowb]=butter(5,lowlim,'low');


for i=1:nrows
	for j=1:ncols
		if ~isnan(stdkeepers(i,j)) & gotstart(j)~=1
			start(j)=i;
			gotstart(j)=1;
		elseif isnan(stdkeepers(i,j)) & gotstart(j)==1 & gotstop(j)~=1
			stop(j)=i-1;
			gotstop(j)=1;	
		end
	end
end
%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%
%%%%	Apply the zero-phase butterworth filter but only to the proxies
	highf=zeros(nrows,ncols);
	lowf=zeros(nrows,ncols); 
	lowf=lowf./lowf;
	highf=highf./highf;
	for k=1:ncols
		clear y y2
		data=stdkeepers(start(k):stop(k),k);
%		lowf(start(k):stop(k),k)=lowpass(data,0.05,2,2);
        frequency = 0.05;
        [smoothed0,icb,ice,mse0] = lowpassmin(data,frequency);
        lowf(start(k):stop(k),k)= smoothed0;
%		lowf(start(k):stop(k),k)=filtfilt(lowa,lowb,data);
		highf(start(k):stop(k),k)=data-lowf(start(k):stop(k),k);
	end
%%%%%
   save('/rootdir/s1/directory1/work1/temanncfr/input/instpchigh','highf','-ASCII');
   save('/rootdir/s1/directory1/work1/temanncfr/input/instpclow','lowf','-ASCII');

%--------------------------------------------------------------------------
%%%%%%%%%%%%%%
   % here, do the EOF on the instrumental and on the proxynework to decide how many
 % EOFs to keep.

clear all; close all;

z=load('/rootdir/s1/directory1/work1/temanncfr/data/rtable1209');

%-----------  r-table description --------------
%    row 1:  longitude
%    row 2:  latitude
%    row 3:  proxy ID (e.g. 8000 annually-resolved ice)
%    row 4:  r -value for raw data over 1850-1995(entire),      n=145 dof.. |r|=0.162
%    row 5:  r -value for low-pass data over 1850-1995,         n=14 dof... |r|=0.496
%    row 6:  r -value for raw data over 1896-1995(late100)       n=99 dof... |r|=0.195
%    row 7:  r -value for low-pass data over 1896-1995          n=9 dof.... |r|=0.602
%---------------------------------------------------
load /rootdir/s1/directory1/work1/temanncfr/data/itrdbmatrix
x=itrdbmatrix;
[n1,m1]=size(x)

%-- select all these chronologies that located within a certain area --
%                   (e.g. -120W~-60W, 20N~60N)


%%-- global proxy network---
ilon1 = -180; ilon2 = 180; 
ilat1 = -90; ilat2 = 90;   
%%--------------------------

%%-- North hemisphere proxy network---
%ilon1 = -180; ilon2 = 180; 
%ilat1 = 0; ilat2 = 90;   
%%--------------------------

%%-- South hemisphere proxy network---
%ilon1 = -180; ilon2 = 180; 
%ilat1 = -90; ilat2 = 0;   
%%--------------------------
%------------------------------------------------------------------

%------------------------------------------------------------------
%% below is for selecting full/screened/1856-1925 screened/1926-1995 screened proxy-network
%% replacing "abs(z(4,i))>=0.165"/"abs(z(5,i))>=0.513" in line 75/84 with the followings for your expected proxy-network
%%        abs(z(4,i))>=0   / abs(z(5,i))>=0          (full proxy-network)
%%        abs(z(4,i))>=0.162 / abs(z(5,i))>=0.496    (screening over 1850-1995)
%%        abs(z(6,i))>=0.195 / abs(z(7,i))>=0.602      (screening over 1896-1995)
%------------------------------------------------------------------


for istep=1:19
%istep=1;

byr=1800-(istep-1)*100;
eyr=1995;

kk=find(x(:,1)==byr);
kkk=find(x(:,1)==eyr); 

x(kk,1:3)
x(kkk,1:3)

nc=0;
for i=1:m1-1 % 
%    (z(3,i)==9000 |  z(3,i)==8000 |  z(3,i)==7500 | z(3,i)==7000 |  z(3,i)==6000 |  z(3,i)==5000 |  z(3,i)==4000 |  z(3,i)==3000 |  z(3,i)==2000)  &...
    if (z(3,i)==9000 |  z(3,i)==8000 |  z(3,i)==7500 | z(3,i)==7000 |  z(3,i)==6000 |  z(3,i)==5000 |  z(3,i)==4000 |  z(3,i)==3000 |  z(3,i)==2000)  &...
         x(kk,i+1)>-99999 & x(kkk,i+1)>-99999 &...
            z(1,i)>=ilon1 & z(1,i)<=ilon2 & z(2,i)>=ilat1 & z(2,i)<=ilat2 & abs(z(6,i))>=0
        nc=nc+1;
        yc(1:kkk-kk+1,nc)=x(kk:kkk,i+1);
        zzc(:,nc)=z(:,i);
    end
end
nd=0; yd=[]; zzd=[];
for i=1:m1-1 % 
%    if (z(3,i)==9001 |  z(3,i)==8001 |  z(3,i)==7501 | z(3,i)==7001 |  z(3,i)==6001 |  z(3,i)==5001 |  z(3,i)==4001 |  z(3,i)==3001)  &...%no Lut data
    if (z(3,i)==9001 |  z(3,i)==8001 |  z(3,i)==7501 | z(3,i)==7001 |  z(3,i)==6001 |  z(3,i)==5001 |  z(3,i)==4001 |  z(3,i)==3001)  &...
         x(kk,i+1)>-99999 & x(kkk,i+1)>-99999 &...
            z(1,i)>=ilon1 & z(1,i)<=ilon2 & z(2,i)>=ilat1 & z(2,i)<=ilat2 & abs(z(7,i))>=0
        nd=nd+1;
        yd(1:kkk-kk+1,nd)=x(kk:kkk,i+1);
        zzd(:,nd)=z(:,i);
    end
end


%num(istep,1)=nc+nd;
num(istep,1)=nc;
num(istep,2)=nd;
y=[yc yd];
zz=[zzc zzd];

% only use lag 0-------------------------------------
[mmm,nnnc]=size(yc); [mmm,nnnd]=size(yd);
yyy(1:mmm,1:ceil(nnnc/3)+ceil(nnnd/3))=NaN;
yyy(1:mmm,1:ceil(nnnc/3))=y(1:mmm,1:ceil(nnnc/3));
yyy(1:mmm,ceil(nnnc/3)+1:ceil(nnnc/3)+ceil(nnnd/3))=y(1:mmm,nnnc+1:nnnc+ceil(nnnd/3));
%---------------------------------------------------------------

save(strcat('/rootdir/s1/directory1/work1/temanncfr/input/proxy',num2str(istep)),'yyy','-ASCII'); % no normalize again this time
clear y zz yc yd zzc zzd yyy mmm nd nc
end

clear all;

for istep=1:19
    
    z=load(strcat('/rootdir/s1/directory1/work1/temanncfr/input/proxy',num2str(istep)));
    keepers=z;
    [nrows,ncols]=size(keepers)

    means=nanmean(keepers(nrows-99:nrows,:)); % normalization over 1850-1995/ or whatever intervals
    stdevs=nanstd(keepers(nrows-99:nrows,:));
    stdkeepers=keepers-repmat(means,nrows,1);
    stdkeepers=stdkeepers./repmat(stdevs,nrows,1);
    
    save(strcat('/rootdir/s1/directory1/work1/temanncfr/input/proxy',num2str(istep)),'stdkeepers','-ASCII');
end


for istep=1:19
    
    z=load(strcat('/rootdir/s1/directory1/work1/temanncfr/input/proxy',num2str(istep)));
    keepers=z;
    [nrows,ncols]=size(keepers)

    stepproxies=keepers(nrows-99:nrows,:); 
    [U,S,V] = svd(stepproxies,0); % do the decomposition, economy style. U's are the PCs
    pctvar=diag(S.^2)/sum(diag(S.^2));
    ddd=eigenselect(pctvar);
    instpcstokeep(istep)=ddd;
%    if ddd > 1
%        dddd = ddd-1;
%    else
%        dddd = ddd;
%    end
%    instpcstokeep(istep)=dddd;
end
   save('/rootdir/s1/directory1/work1/temanncfr/input/instpcstokeep','instpcstokeep','-ASCII');

   
   
% now get the real reconstructiona..................


clear all; close all;

instpcmeans=load('/rootdir/s1/directory1/work1/temanncfr/input/instpcmeans');
instpcstds=load('/rootdir/s1/directory1/work1/temanncfr/input/instpcstds');
u=load('/rootdir/s1/directory1/work1/temanncfr/input/u');
v=load('/rootdir/s1/directory1/work1/temanncfr/input/v');
s=load('/rootdir/s1/directory1/work1/temanncfr/input/s');
%t=load('/rootdir/s1/directory1/work1/temanncfr/input/instpcstokeep');
for istep=1:19
    
    zh=load(strcat('/rootdir/s1/directory1/work1/temanncfr/zzveri/allprogllagnolut_3d/highf/recon',num2str(istep)));
    zl=load(strcat('/rootdir/s1/directory1/work1/temanncfr/zzveri/allprogllagnolut_3d/lowf/recon',num2str(istep)));

    zz=zh+zl;
    [nyears,ncols]=size(zz);
    means=instpcmeans(1:ncols);
    stdevs=instpcstds(1:ncols);
    zz=zz.*repmat(stdevs,nyears,1);
    zz=zz+repmat(means,nyears,1);
    
    qqq=zz(:,1:ncols)*s(1:ncols,1:ncols)*v(:,1:ncols)'; % transfer PCs to grids
    save(strcat('/rootdir/s1/directory1/work1/temanncfr/reconveri/recon',num2str(istep)),'qqq','-ASCII');
end


clear all; close all;

instmeans=load('/rootdir/s1/directory1/work1/temanncfr/input/instmeans');
inststds=load('/rootdir/s1/directory1/work1/temanncfr/input/inststds');

for istep=1:19
    
    zzz=load(strcat('/rootdir/s1/directory1/work1/temanncfr/reconveri/recon',num2str(istep)));
    [nyears,ncols]=size(zzz);
    means=instmeans(1:ncols);
    stdevs=inststds(1:ncols);
    zzz=zzz.*repmat(stdevs,nyears,1);
    zzz=zzz+repmat(means,nyears,1); % back to original data scale
    
    save(strcat('/rootdir/s1/directory1/work1/temanncfr/reconveri/recon',num2str(istep)),'zzz','-ASCII');
end




%###################################################################
%////////////// calculate multivariate skill(beta) of reconstruction (annual) ////////////////////////
%###################################################################

clear all; close all;

veridat=[];

load /rootdir/s1/directory1/hadcru3v/infillannt_3d/instfilled_3d
load /rootdir/s1/directory1/hadcru3v/infillannt_3d/nanmatrix_3d

act3=instfilled_3d(1:146,:);
tnnd=nanmatrix_3d(1:146,:);
act3=act3.*tnnd;

for istep=1:19
%istep=1

    instfilled_3d=load(strcat('/rootdir/s1/directory1/work1/temanncfr/reconveri/recon',num2str(istep)));
    [iyy,mi]=size(instfilled_3d);
    inf3=instfilled_3d(iyy-145:iyy,:);
    inf3=inf3.*tnnd;
    
    
    act=act3(1:46,1:mi);
    inf=inf3(1:46,1:mi);

[n1,m1]=size(act)


for i=1:mi
    aare(i)=nanmean(act3(47:146,i));
    aace(i)=nanmean(act3(1:46,i));
end

for i=1:mi
    for j=1:46
        actr(j,i)=act(j,i)-aare(i);
        actc(j,i)=act(j,i)-aace(i);
    end
end

%aar = reshape(actr,n1*m1,1);
%aac = reshape(actc,n1*m1,1);
   
%aa = reshape(act,n1*m1,1);
%bb = reshape(inf,n1*m1,1);

aarjb = reshape(actr,n1*m1,1);
aacjb = reshape(actc,n1*m1,1);
   
aajb = reshape(act,n1*m1,1);
bbjb = reshape(inf,n1*m1,1);

njb=0;
for i=1:n1*m1
    if aajb(i,1) > -999999
        njb=njb+1
        aar(njb,1)=aarjb(i,1);
        aac(njb,1)=aacjb(i,1);
        aa(njb,1)=aajb(i,1);
        bb(njb,1)=bbjb(i,1);
    end
end
        
        



cc=aa-bb;
ccw=cc'*cc;
%---- calculate multivariate RE of reconstruction------
aaw=aar'*aar;
RE=1-ccw/aaw

%---- calculate multivariate CE of reconstruction------
aaw=aac'*aac;
CE=1-ccw/aaw
%--------------------------------
    actn=act3(1:46,748:1730); % for northern hemisphere
    infn=inf3(1:46,748:1730);

for i=1:983
    aaren(i)=nanmean(act3(47:146,i+747));
    aacen(i)=nanmean(act3(1:46,i+747));
end
for i=1:983
    for j=1:46
        actrn(j,i)=actn(j,i)-aaren(i);
        actcn(j,i)=actn(j,i)-aacen(i);
    end
end
%aarn = reshape(actrn,46*983,1);
%aacn = reshape(actcn,46*983,1);
%aan = reshape(actn,46*983,1);
%bbn = reshape(infn,46*983,1);

aarnjb = reshape(actrn,46*983,1);
aacnjb = reshape(actcn,46*983,1);
aanjb = reshape(actn,46*983,1);
bbnjb = reshape(infn,46*983,1);


njb=0;
for i=1:46*983
    if aanjb(i,1) > -999999
        njb=njb+1
        aarn(njb,1)=aarnjb(i,1);
        aacn(njb,1)=aacnjb(i,1);
        aan(njb,1)=aanjb(i,1);
        bbn(njb,1)=bbnjb(i,1);
    end
end



ccn=aan-bbn;
ccwn=ccn'*ccn;
%---- calculate multivariate RE of reconstruction------
aawn=aarn'*aarn;
REn=1-ccwn/aawn
%---- calculate multivariate CE of reconstruction------
aawn=aacn'*aacn;
CEn=1-ccwn/aawn
%--------------------------------
%--------------------------------
    acts=act3(1:46,1:747); % for southern hemisphere
    infs=inf3(1:46,1:747);

for i=1:747
    aares(i)=nanmean(act3(47:146,i));
    aaces(i)=nanmean(act3(1:46,i));
end
for i=1:747
    for j=1:46
        actrs(j,i)=acts(j,i)-aares(i);
        actcs(j,i)=acts(j,i)-aaces(i);
    end
end
%aars = reshape(actrs,46*747,1);
%aacs = reshape(actcs,46*747,1);
%aas = reshape(acts,46*747,1);
%bbs = reshape(infs,46*747,1);

aarsjb = reshape(actrs,46*747,1);
aacsjb = reshape(actcs,46*747,1);
aasjb = reshape(acts,46*747,1);
bbsjb = reshape(infs,46*747,1);

njb=0;
for i=1:46*747
    if aasjb(i,1) > -999999
        njb=njb+1
        aars(njb,1)=aarsjb(i,1);
        aacs(njb,1)=aacsjb(i,1);
        aas(njb,1)=aasjb(i,1);
        bbs(njb,1)=bbsjb(i,1);
    end
end




ccs=aas-bbs;
ccws=ccs'*ccs;
%---- calculate multivariate RE of reconstruction------
aaws=aars'*aars;
REs=1-ccws/aaws
%---- calculate multivariate CE of reconstruction------
aaws=aacs'*aacs;
CEs=1-ccws/aaws
%--------------------------------
%--------------------------------
    actq=act3(1:46,14:1586); % for 60N-60S
    infq=inf3(1:46,14:1586);

for i=1:1573
    aareq(i)=nanmean(act3(47:146,i+13));
    aaceq(i)=nanmean(act3(1:46,i+13));
end
for i=1:1573
    for j=1:46
        actrq(j,i)=actq(j,i)-aareq(i);
        actcq(j,i)=actq(j,i)-aaceq(i);
    end
end
%aarq = reshape(actrq,46*1573,1);
%aacq = reshape(actcq,46*1573,1);
%aaq = reshape(actq,46*1573,1);
%bbq = reshape(infq,46*1573,1);

aarqjb = reshape(actrq,46*1573,1);
aacqjb = reshape(actcq,46*1573,1);
aaqjb = reshape(actq,46*1573,1);
bbqjb = reshape(infq,46*1573,1);

njb=0;
for i=1:46*1573
    if aaqjb(i,1) > -999999
        njb=njb+1
        aarq(njb,1)=aarqjb(i,1);
        aacq(njb,1)=aacqjb(i,1);
        aaq(njb,1)=aaqjb(i,1);
        bbq(njb,1)=bbqjb(i,1);
    end
end





ccq=aaq-bbq;
ccwq=ccq'*ccq;
%---- calculate multivariate RE of reconstruction------
aawq=aarq'*aarq;
REq=1-ccwq/aawq
%---- calculate multivariate CE of reconstruction------
aawq=aacq'*aacq;
CEq=1-ccwq/aawq
%--------------------------------

veridat(istep,1:6)=[REn CEn REs CEs RE CE];
%---- end calculating multivariate skill of reconstruction------

end

save('/rootdir/s1/directory1/work1/temanncfr/data/veridatlagnolut_3d','veridat','-ASCII');

%###################################################################
%////////////// calculate multivariate skill(beta) of reconstruction (decadal) ////////////////////////
%###################################################################

clear all; close all;

veridat=[];

load /rootdir/s1/directory1/hadcru3v/infillannt_3d/instfilled_3d
load /rootdir/s1/directory1/hadcru3v/infillannt_3d/nanmatrix_3d

act3=instfilled_3d(1:146,:);
tnnd=nanmatrix_3d(1:146,:);
%%------------- smooth instrumental data series---------------------
stdkeepers=act3;
    cutoff=10; %years
	lowlim=1./(cutoff./2)
%%%%%%%%%%%%%%
[nrows,ncols]=size(stdkeepers) 
gotstart=zeros(1,ncols);
gotstop=zeros(1,ncols);
start=zeros(1,ncols);
stop=zeros(1,ncols);
stop(:)=nrows;
[lowa,lowb]=butter(5,lowlim,'low');

%%%%%%%%%%%%%%  
%%%% This section finds the starting and stopping year for each column. It assumes that
%%%%  there are no missing values between start and stop, otherwise filtering chokes.

for i=1:nrows
	for j=1:ncols
		if ~isnan(stdkeepers(i,j)) & gotstart(j)~=1
			start(j)=i;
			gotstart(j)=1;
		elseif isnan(stdkeepers(i,j)) & gotstart(j)==1 & gotstop(j)~=1
			stop(j)=i-1;
			gotstop(j)=1;	
		end
	end
end
%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%
%%%%	Apply the zero-phase butterworth filter but only to the proxies
	highf=zeros(nrows,ncols);
	lowf=zeros(nrows,ncols); 
	lowf=lowf./lowf;
	highf=highf./highf;
	for k=1:ncols
		clear y y2
		data=stdkeepers(start(k):stop(k),k);
%		lowf(start(k):stop(k),k)=lowpass(data,0.05,2,2);
        frequency = 0.1;
        [smoothed0,icb,ice,mse0] = lowpassmin(data,frequency);
        lowf(start(k):stop(k),k)= smoothed0;
%		lowf(start(k):stop(k),k)=filtfilt(lowa,lowb,data);
		highf(start(k):stop(k),k)=data-lowf(start(k):stop(k),k);
	end
    %%
  act3=  lowf;
  act3=act3.*tnnd;
%%-----------------------------------------------------------------

for istep=1:19

    instfilled_3d=load(strcat('/rootdir/s1/directory1/work1/temanncfr/reconveri/recon',num2str(istep)));
    [iyy,mi]=size(instfilled_3d);
    inf3=instfilled_3d(iyy-145:iyy,:);
%%------------- smooth instrumental data series---------------------
stdkeepers=inf3;
    cutoff=10; %years
	lowlim=1./(cutoff./2)
%%%%%%%%%%%%%%
[nrows,ncols]=size(stdkeepers) 
gotstart=zeros(1,ncols);
gotstop=zeros(1,ncols);
start=zeros(1,ncols);
stop=zeros(1,ncols);
stop(:)=nrows;
[lowa,lowb]=butter(5,lowlim,'low');

%%%%%%%%%%%%%%  
%%%% This section finds the starting and stopping year for each column. It assumes that
%%%%  there are no missing values between start and stop, otherwise filtering chokes.

for i=1:nrows
	for j=1:ncols
		if ~isnan(stdkeepers(i,j)) & gotstart(j)~=1
			start(j)=i;
			gotstart(j)=1;
		elseif isnan(stdkeepers(i,j)) & gotstart(j)==1 & gotstop(j)~=1
			stop(j)=i-1;
			gotstop(j)=1;	
		end
	end
end
%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%
%%%%	Apply the zero-phase butterworth filter but only to the proxies
	highf=zeros(nrows,ncols);
	lowf=zeros(nrows,ncols); 
	lowf=lowf./lowf;
	highf=highf./highf;
	for k=1:ncols
		clear y y2
		data=stdkeepers(start(k):stop(k),k);
%		lowf(start(k):stop(k),k)=lowpass(data,0.05,2,2);
        frequency = 0.1;
        [smoothed0,icb,ice,mse0] = lowpassmin(data,frequency);
        lowf(start(k):stop(k),k)= smoothed0;
%		lowf(start(k):stop(k),k)=filtfilt(lowa,lowb,data);
		highf(start(k):stop(k),k)=data-lowf(start(k):stop(k),k);
	end
    %%
  inf3=  lowf;
  inf3=inf3.*tnnd;
%%-----------------------------------------------------------------
    
    act=act3(1:46,1:mi);
    inf=inf3(1:46,1:mi);

[n1,m1]=size(act)


for i=1:mi
    aare(i)=nanmean(act3(47:146,i));
    aace(i)=nanmean(act3(1:46,i));
end

for i=1:mi
    for j=1:46
        actr(j,i)=act(j,i)-aare(i);
        actc(j,i)=act(j,i)-aace(i);
    end
end

%aar = reshape(actr,n1*m1,1);
%aac = reshape(actc,n1*m1,1);
   
%aa = reshape(act,n1*m1,1);
%bb = reshape(inf,n1*m1,1);

aarjb = reshape(actr,n1*m1,1);
aacjb = reshape(actc,n1*m1,1);
   
aajb = reshape(act,n1*m1,1);
bbjb = reshape(inf,n1*m1,1);

njb=0;
for i=1:n1*m1
    if aajb(i,1) > -999999
        njb=njb+1
        aar(njb,1)=aarjb(i,1);
        aac(njb,1)=aacjb(i,1);
        aa(njb,1)=aajb(i,1);
        bb(njb,1)=bbjb(i,1);
    end
end
        
        



cc=aa-bb;
ccw=cc'*cc;
%---- calculate multivariate RE of reconstruction------
aaw=aar'*aar;
RE=1-ccw/aaw

%---- calculate multivariate CE of reconstruction------
aaw=aac'*aac;
CE=1-ccw/aaw
%--------------------------------
    actn=act3(1:46,748:1730); % for northern hemisphere
    infn=inf3(1:46,748:1730);

for i=1:983
    aaren(i)=nanmean(act3(47:146,i+747));
    aacen(i)=nanmean(act3(1:46,i+747));
end
for i=1:983
    for j=1:46
        actrn(j,i)=actn(j,i)-aaren(i);
        actcn(j,i)=actn(j,i)-aacen(i);
    end
end
%aarn = reshape(actrn,46*983,1);
%aacn = reshape(actcn,46*983,1);
%aan = reshape(actn,46*983,1);
%bbn = reshape(infn,46*983,1);

aarnjb = reshape(actrn,46*983,1);
aacnjb = reshape(actcn,46*983,1);
aanjb = reshape(actn,46*983,1);
bbnjb = reshape(infn,46*983,1);


njb=0;
for i=1:46*983
    if aanjb(i,1) > -999999
        njb=njb+1
        aarn(njb,1)=aarnjb(i,1);
        aacn(njb,1)=aacnjb(i,1);
        aan(njb,1)=aanjb(i,1);
        bbn(njb,1)=bbnjb(i,1);
    end
end



ccn=aan-bbn;
ccwn=ccn'*ccn;
%---- calculate multivariate RE of reconstruction------
aawn=aarn'*aarn;
REn=1-ccwn/aawn
%---- calculate multivariate CE of reconstruction------
aawn=aacn'*aacn;
CEn=1-ccwn/aawn
%--------------------------------
%--------------------------------
    acts=act3(1:46,1:747); % for southern hemisphere
    infs=inf3(1:46,1:747);

for i=1:747
    aares(i)=nanmean(act3(47:146,i));
    aaces(i)=nanmean(act3(1:46,i));
end
for i=1:747
    for j=1:46
        actrs(j,i)=acts(j,i)-aares(i);
        actcs(j,i)=acts(j,i)-aaces(i);
    end
end
%aars = reshape(actrs,46*747,1);
%aacs = reshape(actcs,46*747,1);
%aas = reshape(acts,46*747,1);
%bbs = reshape(infs,46*747,1);

aarsjb = reshape(actrs,46*747,1);
aacsjb = reshape(actcs,46*747,1);
aasjb = reshape(acts,46*747,1);
bbsjb = reshape(infs,46*747,1);

njb=0;
for i=1:46*747
    if aasjb(i,1) > -999999
        njb=njb+1
        aars(njb,1)=aarsjb(i,1);
        aacs(njb,1)=aacsjb(i,1);
        aas(njb,1)=aasjb(i,1);
        bbs(njb,1)=bbsjb(i,1);
    end
end




ccs=aas-bbs;
ccws=ccs'*ccs;
%---- calculate multivariate RE of reconstruction------
aaws=aars'*aars;
REs=1-ccws/aaws
%---- calculate multivariate CE of reconstruction------
aaws=aacs'*aacs;
CEs=1-ccws/aaws
%--------------------------------
%--------------------------------
    actq=act3(1:46,14:1586); % for 60N-60S
    infq=inf3(1:46,14:1586);

for i=1:1573
    aareq(i)=nanmean(act3(47:146,i+13));
    aaceq(i)=nanmean(act3(1:46,i+13));
end
for i=1:1573
    for j=1:46
        actrq(j,i)=actq(j,i)-aareq(i);
        actcq(j,i)=actq(j,i)-aaceq(i);
    end
end
%aarq = reshape(actrq,46*1573,1);
%aacq = reshape(actcq,46*1573,1);
%aaq = reshape(actq,46*1573,1);
%bbq = reshape(infq,46*1573,1);

aarqjb = reshape(actrq,46*1573,1);
aacqjb = reshape(actcq,46*1573,1);
aaqjb = reshape(actq,46*1573,1);
bbqjb = reshape(infq,46*1573,1);

njb=0;
for i=1:46*1573
    if aaqjb(i,1) > -999999
        njb=njb+1
        aarq(njb,1)=aarqjb(i,1);
        aacq(njb,1)=aacqjb(i,1);
        aaq(njb,1)=aaqjb(i,1);
        bbq(njb,1)=bbqjb(i,1);
    end
end





ccq=aaq-bbq;
ccwq=ccq'*ccq;
%---- calculate multivariate RE of reconstruction------
aawq=aarq'*aarq;
REq=1-ccwq/aawq
%---- calculate multivariate CE of reconstruction------
aawq=aacq'*aacq;
CEq=1-ccwq/aawq
%--------------------------------

veridat(istep,1:6)=[REn CEn REs CEs RE CE];
%---- end calculating multivariate skill of reconstruction------

end

save('/rootdir/s1/directory1/work1/temanncfr/data/veridatsmlagnolut_3d','veridat','-ASCII');


%###################################################################
%---------------------------
%////////////// calculate SPATIAL skill(beta) of reconstruction (annual) ////////////////////////
%---------------------------
%###################################################################

clear all; close all;

veridat=[];

load /rootdir/s1/directory1/hadcru3v/infillannt_3d/instfilled_3d
load /rootdir/s1/directory1/hadcru3v/infillannt_3d/nanmatrix_3d

act3=instfilled_3d(1:146,:);
tnnd=nanmatrix_3d(1:146,:);
act3=act3.*tnnd;

for istep=1:19
% istep=1

    instfilled_3d=load(strcat('/rootdir/s1/directory1/work1/temanncfr/reconveri/recon',num2str(istep)));
    [iyy,mi]=size(instfilled_3d);
    inf3=instfilled_3d(iyy-145:iyy,:);
    inf3=inf3.*tnnd;
    
 
%for i=1:mi
%    aare(i)=nanmean(act3(47:146,i));
%    aace(i)=nanmean(act3(1:46,i));
%end

for i=1:mi
    aarejb(i)=nanmean(act3(47:146,i));
    aacejb(i)=nanmean(act3(1:46,i));
end

nj=0;
for i=1:mi
    if aacejb(i)>-99999999
        nj=nj+1;
        aare(nj)=aarejb(i);
        aace(nj)=aacejb(i);
        aaa(:,nj)=act3(:,i);
        bbb(:,nj)=inf3(:,i);
        tjjd(nj)=i;
    end
end

tddd=tjjd;

for i=1:nj
%    aa=act3(1:46,i); bb=inf3(1:46,i);
%    cc=aa-bb;
    
    aajb=aaa(1:46,i); bbjb=bbb(1:46,i);
    ccjb=aajb-bbjb;
    njj=0;
    for jj=1:46
        if aaa(jj,i)>-9999999
            njj=njj+1;
            aa(njj,1)=aaa(jj,i);
            bb(njj,1)=bbb(jj,i);
        end
    end
    cc=aa-bb;
        
    
    
    ccw=cc'*cc;
    aaw=aa'*aa;
    betamulti=1-ccw/aaw;

    rmulti=corrcoef(aa,bb);
    rm(istep,i)=rmulti(1,2)^2;

    %---- calculate multivariate RE of reconstruction------
    aaw=(aa-aare(i))'*(aa-aare(i));
    re(istep,i)=1-ccw/aaw;

    %---- calculate multivariate CE of reconstruction------
    aaw=(aa-aace(i))'*(aa-aace(i));
    ce(istep,i)=1-ccw/aaw;
end
clear aa bb cc aaa bbb tjjd
end
rrm(1:19,1:mi)=NaN; rre(1:19,1:mi)=NaN; cce(1:19,1:mi)=NaN;
rrm(:,[tddd])=rm; rre(:,[tddd])=re; cce(:,[tddd])=ce;
save('/rootdir/s1/directory1/work1/temanncfr/data/verisparmlagnolut_3d','rrm','-ASCII');
save('/rootdir/s1/directory1/work1/temanncfr/data/verisparelagnolut_3d','rre','-ASCII');
save('/rootdir/s1/directory1/work1/temanncfr/data/verispacelagnolut_3d','cce','-ASCII');
gridcont(1:mi)=NaN; gridcont([tddd])=1;
save('/rootdir/s1/directory1/work1/temanncfr/input/gridcont','gridcont','-ASCII');

%###################################################################
%---------------------------
%////////////// calculate SPATIAL skill(beta) of reconstruction (decadal) ////////////////////////
%---------------------------
%###################################################################

clear all; close all;

veridat=[];

load /rootdir/s1/directory1/hadcru3v/infillannt_3d/instfilled_3d
load /rootdir/s1/directory1/hadcru3v/infillannt_3d/nanmatrix_3d

act3=instfilled_3d(1:146,:);
tnnd=nanmatrix_3d(1:146,:);
%%------------- smooth instrumental data series---------------------
stdkeepers=act3;
    cutoff=10; %years
	lowlim=1./(cutoff./2)
%%%%%%%%%%%%%%
[nrows,ncols]=size(stdkeepers) 
gotstart=zeros(1,ncols);
gotstop=zeros(1,ncols);
start=zeros(1,ncols);
stop=zeros(1,ncols);
stop(:)=nrows;
[lowa,lowb]=butter(5,lowlim,'low');

%%%%%%%%%%%%%%  
%%%% This section finds the starting and stopping year for each column. It assumes that
%%%%  there are no missing values between start and stop, otherwise filtering chokes.

for i=1:nrows
	for j=1:ncols
		if ~isnan(stdkeepers(i,j)) & gotstart(j)~=1
			start(j)=i;
			gotstart(j)=1;
		elseif isnan(stdkeepers(i,j)) & gotstart(j)==1 & gotstop(j)~=1
			stop(j)=i-1;
			gotstop(j)=1;	
		end
	end
end
%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%
%%%%	Apply the zero-phase butterworth filter but only to the proxies
	highf=zeros(nrows,ncols);
	lowf=zeros(nrows,ncols); 
	lowf=lowf./lowf;
	highf=highf./highf;
	for k=1:ncols
		clear y y2
		data=stdkeepers(start(k):stop(k),k);
%		lowf(start(k):stop(k),k)=lowpass(data,0.05,2,2);
        frequency = 0.1;
        [smoothed0,icb,ice,mse0] = lowpassmin(data,frequency);
        lowf(start(k):stop(k),k)= smoothed0;
%		lowf(start(k):stop(k),k)=filtfilt(lowa,lowb,data);
		highf(start(k):stop(k),k)=data-lowf(start(k):stop(k),k);
	end
    %%
  act3=  lowf;
  act3=act3.*tnnd;
%%-----------------------------------------------------------------

for istep=1:19

    instfilled_3d=load(strcat('/rootdir/s1/directory1/work1/temanncfr/reconveri/recon',num2str(istep)));
    [iyy,mi]=size(instfilled_3d);
    inf3=instfilled_3d(iyy-145:iyy,:);
%%------------- smooth instrumental data series---------------------
stdkeepers=inf3;
    cutoff=10; %years
	lowlim=1./(cutoff./2)
%%%%%%%%%%%%%%
[nrows,ncols]=size(stdkeepers) 
gotstart=zeros(1,ncols);
gotstop=zeros(1,ncols);
start=zeros(1,ncols);
stop=zeros(1,ncols);
stop(:)=nrows;
[lowa,lowb]=butter(5,lowlim,'low');

%%%%%%%%%%%%%%  
%%%% This section finds the starting and stopping year for each column. It assumes that
%%%%  there are no missing values between start and stop, otherwise filtering chokes.

for i=1:nrows
	for j=1:ncols
		if ~isnan(stdkeepers(i,j)) & gotstart(j)~=1
			start(j)=i;
			gotstart(j)=1;
		elseif isnan(stdkeepers(i,j)) & gotstart(j)==1 & gotstop(j)~=1
			stop(j)=i-1;
			gotstop(j)=1;	
		end
	end
end
%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%
%%%%	Apply the zero-phase butterworth filter but only to the proxies
	highf=zeros(nrows,ncols);
	lowf=zeros(nrows,ncols); 
	lowf=lowf./lowf;
	highf=highf./highf;
	for k=1:ncols
		clear y y2
		data=stdkeepers(start(k):stop(k),k);
%		lowf(start(k):stop(k),k)=lowpass(data,0.05,2,2);
        frequency = 0.1;
        [smoothed0,icb,ice,mse0] = lowpassmin(data,frequency);
        lowf(start(k):stop(k),k)= smoothed0;
%		lowf(start(k):stop(k),k)=filtfilt(lowa,lowb,data);
		highf(start(k):stop(k),k)=data-lowf(start(k):stop(k),k);
	end
    %%
  inf3=  lowf;
  inf3=inf3.*tnnd;
%%-----------------------------------------------------------------
    
 
%for i=1:mi
%    aare(i)=nanmean(act3(47:146,i));
%    aace(i)=nanmean(act3(1:46,i));
%end

for i=1:mi
    aarejb(i)=nanmean(act3(47:146,i));
    aacejb(i)=nanmean(act3(1:46,i));
end

nj=0;
for i=1:mi
    if aacejb(i)>-99999999
        nj=nj+1;
        aare(nj)=aarejb(i);
        aace(nj)=aacejb(i);
        aaa(:,nj)=act3(:,i);
        bbb(:,nj)=inf3(:,i);
        tjjd(nj)=i;
    end
end

tddd=tjjd;

for i=1:nj
%    aa=act3(1:46,i); bb=inf3(1:46,i);
%    cc=aa-bb;
    
    aajb=aaa(1:46,i); bbjb=bbb(1:46,i);
    ccjb=aajb-bbjb;
    njj=0;
    for jj=1:46
        if aaa(jj,i)>-9999999
            njj=njj+1;
            aa(njj,1)=aaa(jj,i);
            bb(njj,1)=bbb(jj,i);
        end
    end
    cc=aa-bb;
        
    
    
    ccw=cc'*cc;
    aaw=aa'*aa;
    betamulti=1-ccw/aaw;

    rmulti=corrcoef(aa,bb);
    rm(istep,i)=rmulti(1,2)^2;

    %---- calculate multivariate RE of reconstruction------
    aaw=(aa-aare(i))'*(aa-aare(i));
    re(istep,i)=1-ccw/aaw;

    %---- calculate multivariate CE of reconstruction------
    aaw=(aa-aace(i))'*(aa-aace(i));
    ce(istep,i)=1-ccw/aaw;
end
clear aa bb cc aaa bbb tjjd
end
rrm(1:19,1:mi)=NaN; rre(1:19,1:mi)=NaN; cce(1:19,1:mi)=NaN;
rrm(:,[tddd])=rm; rre(:,[tddd])=re; cce(:,[tddd])=ce;
save('/rootdir/s1/directory1/work1/temanncfr/data/verisparmsmlagnolut_3d','rrm','-ASCII');
save('/rootdir/s1/directory1/work1/temanncfr/data/verisparesmlagnolut_3d','rre','-ASCII');
save('/rootdir/s1/directory1/work1/temanncfr/data/verispacesmlagnolut_3d','cce','-ASCII');


%###################################################################
%###################################################################
%###################################################################
%---------------------------
%////////////// do regional and global mean series ////////////////////////
%---------------------------
%###################################################################
%###################################################################
%###################################################################


%--------- generate mask files ----------

% (1) mult (NH)  
% (2) mult (SH)   
% (3) mult (globe)
% (4) NH mean 
% (5) SH mean 
% (6) global mean
% (7) extratropical mean  
% (8) tropical mean
% (9) Nino3 mean
% (10) N. Pacific  mean
% (11)  N. Atlantic mean
% (12) N. Land mean
% (13) tropical Atlantic Main Development Region



clear all; close all;
load /rootdir/s1/directory1/work1/temann/data/lsmask_jones.dat
landmask=flipud(lsmask_jones);  % 
landmask=landmask./landmask; % this makes the ocean NaN and all others 1
landmask=reshape(landmask',1,36*72); % now string out the land mask 1:36*72

% global gridboxes
tt(1:36*72)=1;

% cosin weight
for k=1:36
    for j=1:72
        tcoswt((k-1)*72+j)= sin(((k-1)*5+2.5)*pi/180);
    end
end

% weight gridboxes
%tt=tt.*tcoswt;

%################################################
% mask NH, SH, PDO, AMO, Nino3, Extrotripical......
%################################################

masknh(1:36*72)=NaN; masknh(18*72+1:36*72)=1;
masksh(1:36*72)=NaN; masksh(1:18*72)=1;

% PDO region mask 20N-60N, 150E-130W (57.5(30)-22.5(23)), 67-72, 1-10,
maskpdo(1:36*72)=NaN;  % for pdo
    for k=23:30
        for j=1:10
            maskpdo((k-1)*72+j)= 1;
        end
    end
    for k=23:30
        for j=67:72
            maskpdo((k-1)*72+j)= 1;
        end
    end

% North Atlantic region (AMO) mask axis([-80 10 0 70])
maskamo(1:36*72)=NaN;  % for amo
    for k=31:32
        for j=29:38
            maskamo((k-1)*72+j)= 1;
        end
    end
    for k=27:30
        for j=26:34
            maskamo((k-1)*72+j)= 1;
        end
    end
    for k=21:26
        for j=22:33
            maskamo((k-1)*72+j)= 1;
        end
    end
    for k=19:20
        for j=26:34
            maskamo((k-1)*72+j)= 1;
        end
    end

% Nino3 mask 7-18(-147.5~-92.5), 18-19(2.5~-2.5)(150W-90W,5N-5S)
masknino(1:36*72)=NaN;  % for Nino3
    for k=18:19
        for j=7:18
            masknino((k-1)*72+j)= 1;
        end
    end
% mask for tropical (30S-30N), extrotropical (30N-70N)
masktro(1:36*72)=NaN; masktro(12*72+1:24*72)=1;
maskextro(1:36*72)=NaN; maskextro(23*72+1:32*72)=1;

% mask NH land 
masknhland=landmask; masknhland(1:18*72)=NaN;

% mask tropical Atlantic Main Development Region(tamdr)
%(6-18N,20-60W)--[20(7.5N)-22(17.5N),24(-62.5W)-33(-17.5W)]
masktamdr(1:36*72)=NaN;  % for tamdr
    for k=20:22
        for j=24:33
            masktamdr((k-1)*72+j)= 1;
        end
    end


save('/rootdir/s1/directory1/work1/temanncfr/input/tcoswt','tcoswt','-ASCII');
save('/rootdir/s1/directory1/work1/temanncfr/input/masknh','masknh','-ASCII');
save('/rootdir/s1/directory1/work1/temanncfr/input/masksh','masksh','-ASCII');
save('/rootdir/s1/directory1/work1/temanncfr/input/maskgl','tt','-ASCII');
save('/rootdir/s1/directory1/work1/temanncfr/input/maskextro','maskextro','-ASCII');
save('/rootdir/s1/directory1/work1/temanncfr/input/masktro','masktro','-ASCII');
save('/rootdir/s1/directory1/work1/temanncfr/input/masknino','masknino','-ASCII');
save('/rootdir/s1/directory1/work1/temanncfr/input/maskpdo','maskpdo','-ASCII');
save('/rootdir/s1/directory1/work1/temanncfr/input/maskamo','maskamo','-ASCII');
save('/rootdir/s1/directory1/work1/temanncfr/input/masknhland','masknhland','-ASCII');
save('/rootdir/s1/directory1/work1/temanncfr/input/masktamdr','masktamdr','-ASCII');



%-------------- load these mask files --------------------

clear all; close all;

load /rootdir/s1/directory1/work1/temanncfr/input/tcoswt
load /rootdir/s1/directory1/work1/temanncfr/input/masknh
load /rootdir/s1/directory1/work1/temanncfr/input/masksh
load /rootdir/s1/directory1/work1/temanncfr/input/maskgl
load /rootdir/s1/directory1/work1/temanncfr/input/maskextro
load /rootdir/s1/directory1/work1/temanncfr/input/masktro
load /rootdir/s1/directory1/work1/temanncfr/input/masknino
load /rootdir/s1/directory1/work1/temanncfr/input/maskpdo
load /rootdir/s1/directory1/work1/temanncfr/input/maskamo
load /rootdir/s1/directory1/work1/temanncfr/input/masknhland
load /rootdir/s1/directory1/work1/temanncfr/input/masktamdr

% gridboxes that have data to do RE/CE
load /rootdir/s1/directory1/work1/temanncfr/input/gridcont
% gridboxes that have data to do RE/CE---- back to 36*72 gridpoints
load /rootdir/s1/directory1/hadcru3v/enough_3d
load /rootdir/s1/directory1/hadcru3v/toofew_3d
  gridd([enough_3d])=gridcont;
  gridd([toofew_3d])=NaN;

%---- now spatialtemporal mask ----

tcoswt=repmat(tcoswt,146,1);

masknh=masknh.*gridd;
masknh=repmat(masknh,146,1);
masksh=masksh.*gridd;
masksh=repmat(masksh,146,1);
maskgl=maskgl.*gridd;
maskgl=repmat(maskgl,146,1);
maskextro=maskextro.*gridd;
maskextro=repmat(maskextro,146,1);
masktro=masktro.*gridd;
masktro=repmat(masktro,146,1);
masknino=masknino.*gridd;
masknino=repmat(masknino,146,1);
maskpdo=maskpdo.*gridd;
maskpdo=repmat(maskpdo,146,1);
maskamo=maskamo.*gridd;
maskamo=repmat(maskamo,146,1);
masknhland=masknhland.*gridd;
masknhland=repmat(masknhland,146,1);
masktamdr=masktamdr.*gridd;
masktamdr=repmat(masktamdr,146,1);



%###################################################################
%---------------------------
%////////////// calculate RE/CE for reginal/global mean of reconstructions (annual) ////////////////////////
%---------------------------
%###################################################################

veridat=[];

load /rootdir/s1/directory1/hadcru3v/infillannt_3d/instfilled_3d
load /rootdir/s1/directory1/hadcru3v/infillannt_3d/nanmatrix_3d

act3=instfilled_3d(1:146,:);
tnnd=nanmatrix_3d(1:146,:);
act3=act3.*tnnd;
act3(:,[enough_3d])=act3; % back to 36*72 gridpoints
act3(:,[toofew_3d])=NaN;

for istep=1:19
% istep=1

    instfilled_3d=load(strcat('/rootdir/s1/directory1/work1/temanncfr/reconveri/recon',num2str(istep)));
    [iyy,mi]=size(instfilled_3d);
    inf3=instfilled_3d(iyy-145:iyy,:);
    inf3=inf3.*tnnd;
    inf3(:,[enough_3d])=inf3; % back to 36*72 gridpoints
    inf3(:,[toofew_3d])=NaN;
    
    %--- NH mean ---
    www=act3./act3;
    www=(www.*masknh.*tcoswt)';
    www=nanmean(www);
    www=www';
    aaa=(act3.*masknh.*tcoswt)';
    aaa=nanmean(aaa);
    aaa=aaa';
    aaa=aaa./www;
    bbb=(inf3.*masknh.*tcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    bbb=bbb./www;

    aare=nanmean(nanmean(aaa(47:146,1)));
    aace=nanmean(nanmean(aaa(1:46,1)));
    aa=aaa(1:46,1); bb=bbb(1:46,1);
    cc=aa-bb;
    ccw=cc'*cc;
    aaw=(aa-aare)'*(aa-aare);
    R1=1-ccw/aaw
    aaw=(aa-aace)'*(aa-aace);
    C1=1-ccw/aaw
    %aan=aa; bbn=bb;
    rm=corrcoef(aa,bb);
    rr1=rm(1,2)^2;

    %--- SH mean ---
    www=act3./act3;
    www=(www.*masksh.*tcoswt)';
    www=nanmean(www);
    www=www';
    aaa=(act3.*masksh.*tcoswt)';
    aaa=nanmean(aaa);
    aaa=aaa';
    aaa=aaa./www;
    bbb=(inf3.*masksh.*tcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    bbb=bbb./www;

    aare=nanmean(nanmean(aaa(47:146,1)));
    aace=nanmean(nanmean(aaa(1:46,1)));
    aa=aaa(1:46,1); bb=bbb(1:46,1);
    cc=aa-bb;
    ccw=cc'*cc;
    aaw=(aa-aare)'*(aa-aare);
    R2=1-ccw/aaw
    aaw=(aa-aace)'*(aa-aace);
    C2=1-ccw/aaw
    %aan=aa; bbn=bb;
    rm=corrcoef(aa,bb);
    rr2=rm(1,2)^2;
    
    %--- GL mean ---
    www=act3./act3;
    www=(www.*maskgl.*tcoswt)';
    www=nanmean(www);
    www=www';
    aaa=(act3.*maskgl.*tcoswt)';
    aaa=nanmean(aaa);
    aaa=aaa';
    aaa=aaa./www;
    bbb=(inf3.*maskgl.*tcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    bbb=bbb./www;

    aare=nanmean(nanmean(aaa(47:146,1)));
    aace=nanmean(nanmean(aaa(1:46,1)));
    aa=aaa(1:46,1); bb=bbb(1:46,1);
    cc=aa-bb;
    ccw=cc'*cc;
    aaw=(aa-aare)'*(aa-aare);
    R3=1-ccw/aaw
    aaw=(aa-aace)'*(aa-aace);
    C3=1-ccw/aaw
    %aan=aa; bbn=bb;
    rm=corrcoef(aa,bb);
    rr3=rm(1,2)^2;
    
    
    %--- Extropical mean ---
    www=act3./act3;
    www=(www.*maskextro.*tcoswt)';
    www=nanmean(www);
    www=www';
    aaa=(act3.*maskextro.*tcoswt)';
    aaa=nanmean(aaa);
    aaa=aaa';
    aaa=aaa./www;
    bbb=(inf3.*maskextro.*tcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    bbb=bbb./www;

    aare=nanmean(nanmean(aaa(47:146,1)));
    aace=nanmean(nanmean(aaa(1:46,1)));
    aa=aaa(1:46,1); bb=bbb(1:46,1);
    cc=aa-bb;
    ccw=cc'*cc;
    aaw=(aa-aare)'*(aa-aare);
    R4=1-ccw/aaw
    aaw=(aa-aace)'*(aa-aace);
    C4=1-ccw/aaw
    %aan=aa; bbn=bb;
    rm=corrcoef(aa,bb);
    rr4=rm(1,2)^2;
    
    
    %--- Tropical mean ---
    www=act3./act3;
    www=(www.*masktro.*tcoswt)';
    www=nanmean(www);
    www=www';
    aaa=(act3.*masktro.*tcoswt)';
    aaa=nanmean(aaa);
    aaa=aaa';
    aaa=aaa./www;
    bbb=(inf3.*masktro.*tcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    bbb=bbb./www;

    aare=nanmean(nanmean(aaa(47:146,1)));
    aace=nanmean(nanmean(aaa(1:46,1)));
    aa=aaa(1:46,1); bb=bbb(1:46,1);
    cc=aa-bb;
    ccw=cc'*cc;
    aaw=(aa-aare)'*(aa-aare);
    R5=1-ccw/aaw
    aaw=(aa-aace)'*(aa-aace);
    C5=1-ccw/aaw
    %aan=aa; bbn=bb;
    rm=corrcoef(aa,bb);
    rr5=rm(1,2)^2;
    
    
    %--- Nino3 mean ---
    www=act3./act3;
    www=(www.*masknino.*tcoswt)';
    www=nanmean(www);
    www=www';
    aaa=(act3.*masknino.*tcoswt)';
    aaa=nanmean(aaa);
    aaa=aaa';
    aaa=aaa./www;
    bbb=(inf3.*masknino.*tcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    bbb=bbb./www;

    aare=nanmean(nanmean(aaa(47:146,1)));
    aace=nanmean(nanmean(aaa(1:46,1)));
    aa=aaa(1:46,1); bb=bbb(1:46,1);
    
    aajb=aaa(1:46,1); bbjb=bbb(1:46,1);
    ccjb=aajb-bbjb;
    njj=0;
    for jj=1:46
        if aaa(jj,1)>-9999999
            njj=njj+1;
            aa1(njj,1)=aaa(jj,1);   %remove NaN gridpoints in Nino3 region
            bb1(njj,1)=bbb(jj,1);
        end
    end
   
    
    cc=aa1-bb1;
    ccw=cc'*cc;
    aaw=(aa1-aare)'*(aa1-aare);
    R6=1-ccw/aaw
    aaw=(aa1-aace)'*(aa1-aace);
    C6=1-ccw/aaw
    %aan=aa; bbn=bb;
    rm=corrcoef(aa1,bb1);
    rr6=rm(1,2)^2;
    
    
    %--- PDO region mean ---
    www=act3./act3;
    www=(www.*maskpdo.*tcoswt)';
    www=nanmean(www);
    www=www';
    aaa=(act3.*maskpdo.*tcoswt)';
    aaa=nanmean(aaa);
    aaa=aaa';
    aaa=aaa./www;
    bbb=(inf3.*maskpdo.*tcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    bbb=bbb./www;

    aare=nanmean(nanmean(aaa(47:146,1)));
    aace=nanmean(nanmean(aaa(1:46,1)));
    aa=aaa(1:46,1); bb=bbb(1:46,1);
    cc=aa-bb;
    ccw=cc'*cc;
    aaw=(aa-aare)'*(aa-aare);
    R7=1-ccw/aaw
    aaw=(aa-aace)'*(aa-aace);
    C7=1-ccw/aaw
    %aan=aa; bbn=bb;
    rm=corrcoef(aa,bb);
    rr7=rm(1,2)^2;
    
    
    %--- AMO region mean ---
    www=act3./act3;
    www=(www.*maskamo.*tcoswt)';
    www=nanmean(www);
    www=www';
    aaa=(act3.*maskamo.*tcoswt)';
    aaa=nanmean(aaa);
    aaa=aaa';
    aaa=aaa./www;
    bbb=(inf3.*maskamo.*tcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    bbb=bbb./www;

    aare=nanmean(nanmean(aaa(47:146,1)));
    aace=nanmean(nanmean(aaa(1:46,1)));
    aa=aaa(1:46,1); bb=bbb(1:46,1);
    cc=aa-bb;
    ccw=cc'*cc;
    aaw=(aa-aare)'*(aa-aare);
    R8=1-ccw/aaw
    aaw=(aa-aace)'*(aa-aace);
    C8=1-ccw/aaw
    %aan=aa; bbn=bb;
    rm=corrcoef(aa,bb);
    rr8=rm(1,2)^2;
    
    
    %--- NH land mean ---
    www=act3./act3;
    www=(www.*masknhland.*tcoswt)';
    www=nanmean(www);
    www=www';
    aaa=(act3.*masknhland.*tcoswt)';
    aaa=nanmean(aaa);
    aaa=aaa';
    aaa=aaa./www;
    bbb=(inf3.*masknhland.*tcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    bbb=bbb./www;

    aare=nanmean(nanmean(aaa(47:146,1)));
    aace=nanmean(nanmean(aaa(1:46,1)));
    aa=aaa(1:46,1); bb=bbb(1:46,1);
    cc=aa-bb;
    ccw=cc'*cc;
    aaw=(aa-aare)'*(aa-aare);
    R9=1-ccw/aaw
    aaw=(aa-aace)'*(aa-aace);
    C9=1-ccw/aaw
    %aan=aa; bbn=bb;
    rm=corrcoef(aa,bb);
    rr9=rm(1,2)^2;
    
    
    %--- tropical Atlantic main development region(tamdr) mean ---
    www=act3./act3;
    www=(www.*masktamdr.*tcoswt)';
    www=nanmean(www);
    www=www';
    aaa=(act3.*masktamdr.*tcoswt)';
    aaa=nanmean(aaa);
    aaa=aaa';
    aaa=aaa./www;
    bbb=(inf3.*masktamdr.*tcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    bbb=bbb./www;

    aare=nanmean(nanmean(aaa(47:146,1)));
    aace=nanmean(nanmean(aaa(1:46,1)));
    aa=aaa(1:46,1); bb=bbb(1:46,1);
    cc=aa-bb;
    ccw=cc'*cc;
    aaw=(aa-aare)'*(aa-aare);
    R10=1-ccw/aaw
    aaw=(aa-aace)'*(aa-aace);
    C10=1-ccw/aaw
    %aan=aa; bbn=bb;
    rm=corrcoef(aa,bb);
    rr10=rm(1,2)^2;
    
    
    veridat(istep,1:10)=[rr1 rr2 rr3 rr4 rr5 rr6 rr7 rr8 rr9 rr10];
%    veridat(istep,1:30)=[R1 C1 rr1 R2 C2 rr2 R3 C3 rr3 R4 C4 rr4 R5 C5 rr5 R6 C6 rr6 R7 C7 rr7 R8 C8 rr8 R9 C9 rr9 R10 C10 rr10];
end   
    

save('/rootdir/s1/directory1/work1/temanncfr/data/veridatearlylagnolut_3drr','veridat','-ASCII');



%###################################################################
%---------------------------
%////////////// calculate RE/CE for reginal/global mean of reconstructions (decadal) ////////////////////////
%---------------------------
%###################################################################


veridat=[];

load /rootdir/s1/directory1/hadcru3v/infillannt_3d/instfilled_3d
load /rootdir/s1/directory1/hadcru3v/infillannt_3d/nanmatrix_3d

act3=instfilled_3d(1:146,:);
tnnd=nanmatrix_3d(1:146,:);
%%------------- smooth instrumental data series---------------------
stdkeepers=act3;
    cutoff=10; %years
	lowlim=1./(cutoff./2)
%%%%%%%%%%%%%%
[nrows,ncols]=size(stdkeepers) 
gotstart=zeros(1,ncols);
gotstop=zeros(1,ncols);
start=zeros(1,ncols);
stop=zeros(1,ncols);
stop(:)=nrows;
[lowa,lowb]=butter(5,lowlim,'low');

%%%%%%%%%%%%%%  
%%%% This section finds the starting and stopping year for each column. It assumes that
%%%%  there are no missing values between start and stop, otherwise filtering chokes.

for i=1:nrows
	for j=1:ncols
		if ~isnan(stdkeepers(i,j)) & gotstart(j)~=1
			start(j)=i;
			gotstart(j)=1;
		elseif isnan(stdkeepers(i,j)) & gotstart(j)==1 & gotstop(j)~=1
			stop(j)=i-1;
			gotstop(j)=1;	
		end
	end
end
%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%
%%%%	Apply the zero-phase butterworth filter but only to the proxies
	highf=zeros(nrows,ncols);
	lowf=zeros(nrows,ncols); 
	lowf=lowf./lowf;
	highf=highf./highf;
	for k=1:ncols
		clear y y2
		data=stdkeepers(start(k):stop(k),k);
%		lowf(start(k):stop(k),k)=lowpass(data,0.05,2,2);
        frequency = 0.1;
        [smoothed0,icb,ice,mse0] = lowpassmin(data,frequency);
        lowf(start(k):stop(k),k)= smoothed0;
%		lowf(start(k):stop(k),k)=filtfilt(lowa,lowb,data);
		highf(start(k):stop(k),k)=data-lowf(start(k):stop(k),k);
	end
    %%
  act3=  lowf;
  act3=act3.*tnnd;
  act3(:,[enough_3d])=act3; % back to 36*72 gridpoints
  act3(:,[toofew_3d])=NaN;

%%-----------------------------------------------------------------

for istep=1:19

    instfilled_3d=load(strcat('/rootdir/s1/directory1/work1/temanncfr/reconveri/recon',num2str(istep)));
    [iyy,mi]=size(instfilled_3d);
    inf3=instfilled_3d(iyy-145:iyy,:);
%%------------- smooth instrumental data series---------------------
stdkeepers=inf3;
    cutoff=10; %years
	lowlim=1./(cutoff./2)
%%%%%%%%%%%%%%
[nrows,ncols]=size(stdkeepers) 
gotstart=zeros(1,ncols);
gotstop=zeros(1,ncols);
start=zeros(1,ncols);
stop=zeros(1,ncols);
stop(:)=nrows;
[lowa,lowb]=butter(5,lowlim,'low');

%%%%%%%%%%%%%%  
%%%% This section finds the starting and stopping year for each column. It assumes that
%%%%  there are no missing values between start and stop, otherwise filtering chokes.

for i=1:nrows
	for j=1:ncols
		if ~isnan(stdkeepers(i,j)) & gotstart(j)~=1
			start(j)=i;
			gotstart(j)=1;
		elseif isnan(stdkeepers(i,j)) & gotstart(j)==1 & gotstop(j)~=1
			stop(j)=i-1;
			gotstop(j)=1;	
		end
	end
end
%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%
%%%%	Apply the zero-phase butterworth filter but only to the proxies
	highf=zeros(nrows,ncols);
	lowf=zeros(nrows,ncols); 
	lowf=lowf./lowf;
	highf=highf./highf;
	for k=1:ncols
		clear y y2
		data=stdkeepers(start(k):stop(k),k);
%		lowf(start(k):stop(k),k)=lowpass(data,0.05,2,2);
        frequency = 0.1;
        [smoothed0,icb,ice,mse0] = lowpassmin(data,frequency);
        lowf(start(k):stop(k),k)= smoothed0;
%		lowf(start(k):stop(k),k)=filtfilt(lowa,lowb,data);
		highf(start(k):stop(k),k)=data-lowf(start(k):stop(k),k);
	end
    %%
  inf3=  lowf;
  inf3=inf3.*tnnd;
  inf3(:,[enough_3d])=inf3; % back to 36*72 gridpoints
  inf3(:,[toofew_3d])=NaN;
%%-----------------------------------------------------------------
    
    %--- NH mean ---
    www=act3./act3;
    www=(www.*masknh.*tcoswt)';
    www=nanmean(www);
    www=www';
    aaa=(act3.*masknh.*tcoswt)';
    aaa=nanmean(aaa);
    aaa=aaa';
    aaa=aaa./www;
    bbb=(inf3.*masknh.*tcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    bbb=bbb./www;

    aare=nanmean(nanmean(aaa(47:146,1)));
    aace=nanmean(nanmean(aaa(1:46,1)));
    aa=aaa(1:46,1); bb=bbb(1:46,1);
    cc=aa-bb;
    ccw=cc'*cc;
    aaw=(aa-aare)'*(aa-aare);
    R1=1-ccw/aaw
    aaw=(aa-aace)'*(aa-aace);
    C1=1-ccw/aaw
    %aan=aa; bbn=bb;
    rm=corrcoef(aa,bb);
    rr1=rm(1,2)^2;

    %--- SH mean ---
    www=act3./act3;
    www=(www.*masksh.*tcoswt)';
    www=nanmean(www);
    www=www';
    aaa=(act3.*masksh.*tcoswt)';
    aaa=nanmean(aaa);
    aaa=aaa';
    aaa=aaa./www;
    bbb=(inf3.*masksh.*tcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    bbb=bbb./www;

    aare=nanmean(nanmean(aaa(47:146,1)));
    aace=nanmean(nanmean(aaa(1:46,1)));
    aa=aaa(1:46,1); bb=bbb(1:46,1);
    cc=aa-bb;
    ccw=cc'*cc;
    aaw=(aa-aare)'*(aa-aare);
    R2=1-ccw/aaw
    aaw=(aa-aace)'*(aa-aace);
    C2=1-ccw/aaw
    %aan=aa; bbn=bb;
    rm=corrcoef(aa,bb);
    rr2=rm(1,2)^2;
    
    %--- GL mean ---
    www=act3./act3;
    www=(www.*maskgl.*tcoswt)';
    www=nanmean(www);
    www=www';
    aaa=(act3.*maskgl.*tcoswt)';
    aaa=nanmean(aaa);
    aaa=aaa';
    aaa=aaa./www;
    bbb=(inf3.*maskgl.*tcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    bbb=bbb./www;

    aare=nanmean(nanmean(aaa(47:146,1)));
    aace=nanmean(nanmean(aaa(1:46,1)));
    aa=aaa(1:46,1); bb=bbb(1:46,1);
    cc=aa-bb;
    ccw=cc'*cc;
    aaw=(aa-aare)'*(aa-aare);
    R3=1-ccw/aaw
    aaw=(aa-aace)'*(aa-aace);
    C3=1-ccw/aaw
    %aan=aa; bbn=bb;
    rm=corrcoef(aa,bb);
    rr3=rm(1,2)^2;
    
    
    %--- Extropical mean ---
    www=act3./act3;
    www=(www.*maskextro.*tcoswt)';
    www=nanmean(www);
    www=www';
    aaa=(act3.*maskextro.*tcoswt)';
    aaa=nanmean(aaa);
    aaa=aaa';
    aaa=aaa./www;
    bbb=(inf3.*maskextro.*tcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    bbb=bbb./www;

    aare=nanmean(nanmean(aaa(47:146,1)));
    aace=nanmean(nanmean(aaa(1:46,1)));
    aa=aaa(1:46,1); bb=bbb(1:46,1);
    cc=aa-bb;
    ccw=cc'*cc;
    aaw=(aa-aare)'*(aa-aare);
    R4=1-ccw/aaw
    aaw=(aa-aace)'*(aa-aace);
    C4=1-ccw/aaw
    %aan=aa; bbn=bb;
    rm=corrcoef(aa,bb);
    rr4=rm(1,2)^2;
    
    
    %--- Tropical mean ---
    www=act3./act3;
    www=(www.*masktro.*tcoswt)';
    www=nanmean(www);
    www=www';
    aaa=(act3.*masktro.*tcoswt)';
    aaa=nanmean(aaa);
    aaa=aaa';
    aaa=aaa./www;
    bbb=(inf3.*masktro.*tcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    bbb=bbb./www;

    aare=nanmean(nanmean(aaa(47:146,1)));
    aace=nanmean(nanmean(aaa(1:46,1)));
    aa=aaa(1:46,1); bb=bbb(1:46,1);
    cc=aa-bb;
    ccw=cc'*cc;
    aaw=(aa-aare)'*(aa-aare);
    R5=1-ccw/aaw
    aaw=(aa-aace)'*(aa-aace);
    C5=1-ccw/aaw
    %aan=aa; bbn=bb;
    rm=corrcoef(aa,bb);
    rr5=rm(1,2)^2;
    
    
    %--- Nino3 mean ---
    www=act3./act3;
    www=(www.*masknino.*tcoswt)';
    www=nanmean(www);
    www=www';
    aaa=(act3.*masknino.*tcoswt)';
    aaa=nanmean(aaa);
    aaa=aaa';
    aaa=aaa./www;
    bbb=(inf3.*masknino.*tcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    bbb=bbb./www;

    aare=nanmean(nanmean(aaa(47:146,1)));
    aace=nanmean(nanmean(aaa(1:46,1)));
    aa=aaa(1:46,1); bb=bbb(1:46,1);
    
    aajb=aaa(1:46,1); bbjb=bbb(1:46,1);
    ccjb=aajb-bbjb;
    njj=0;
    for jj=1:46
        if aaa(jj,1)>-9999999
            njj=njj+1;
            aa1(njj,1)=aaa(jj,1);   %remove NaN gridpoints in Nino3 region
            bb1(njj,1)=bbb(jj,1);
        end
    end
   
    
    cc=aa1-bb1;
    ccw=cc'*cc;
    aaw=(aa1-aare)'*(aa1-aare);
    R6=1-ccw/aaw
    aaw=(aa1-aace)'*(aa1-aace);
    C6=1-ccw/aaw
    %aan=aa; bbn=bb;
    rm=corrcoef(aa1,bb1);
    rr6=rm(1,2)^2;
    
    
    %--- PDO region mean ---
    www=act3./act3;
    www=(www.*maskpdo.*tcoswt)';
    www=nanmean(www);
    www=www';
    aaa=(act3.*maskpdo.*tcoswt)';
    aaa=nanmean(aaa);
    aaa=aaa';
    aaa=aaa./www;
    bbb=(inf3.*maskpdo.*tcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    bbb=bbb./www;

    aare=nanmean(nanmean(aaa(47:146,1)));
    aace=nanmean(nanmean(aaa(1:46,1)));
    aa=aaa(1:46,1); bb=bbb(1:46,1);
    cc=aa-bb;
    ccw=cc'*cc;
    aaw=(aa-aare)'*(aa-aare);
    R7=1-ccw/aaw
    aaw=(aa-aace)'*(aa-aace);
    C7=1-ccw/aaw
    %aan=aa; bbn=bb;
    rm=corrcoef(aa,bb);
    rr7=rm(1,2)^2;
    
    
    %--- AMO region mean ---
    www=act3./act3;
    www=(www.*maskamo.*tcoswt)';
    www=nanmean(www);
    www=www';
    aaa=(act3.*maskamo.*tcoswt)';
    aaa=nanmean(aaa);
    aaa=aaa';
    aaa=aaa./www;
    bbb=(inf3.*maskamo.*tcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    bbb=bbb./www;

    aare=nanmean(nanmean(aaa(47:146,1)));
    aace=nanmean(nanmean(aaa(1:46,1)));
    aa=aaa(1:46,1); bb=bbb(1:46,1);
    cc=aa-bb;
    ccw=cc'*cc;
    aaw=(aa-aare)'*(aa-aare);
    R8=1-ccw/aaw
    aaw=(aa-aace)'*(aa-aace);
    C8=1-ccw/aaw
    %aan=aa; bbn=bb;
    rm=corrcoef(aa,bb);
    rr8=rm(1,2)^2;
    
    
    %--- NH land mean ---
    www=act3./act3;
    www=(www.*masknhland.*tcoswt)';
    www=nanmean(www);
    www=www';
    aaa=(act3.*masknhland.*tcoswt)';
    aaa=nanmean(aaa);
    aaa=aaa';
    aaa=aaa./www;
    bbb=(inf3.*masknhland.*tcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    bbb=bbb./www;

    aare=nanmean(nanmean(aaa(47:146,1)));
    aace=nanmean(nanmean(aaa(1:46,1)));
    aa=aaa(1:46,1); bb=bbb(1:46,1);
    cc=aa-bb;
    ccw=cc'*cc;
    aaw=(aa-aare)'*(aa-aare);
    R9=1-ccw/aaw
    aaw=(aa-aace)'*(aa-aace);
    C9=1-ccw/aaw
    %aan=aa; bbn=bb;
    rm=corrcoef(aa,bb);
    rr9=rm(1,2)^2;
    
    
    %--- tropical Atlantic main development region(tamdr) mean ---
    www=act3./act3;
    www=(www.*masktamdr.*tcoswt)';
    www=nanmean(www);
    www=www';
    aaa=(act3.*masktamdr.*tcoswt)';
    aaa=nanmean(aaa);
    aaa=aaa';
    aaa=aaa./www;
    bbb=(inf3.*masktamdr.*tcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    bbb=bbb./www;

    aare=nanmean(nanmean(aaa(47:146,1)));
    aace=nanmean(nanmean(aaa(1:46,1)));
    aa=aaa(1:46,1); bb=bbb(1:46,1);
    cc=aa-bb;
    ccw=cc'*cc;
    aaw=(aa-aare)'*(aa-aare);
    R10=1-ccw/aaw
    aaw=(aa-aace)'*(aa-aace);
    C10=1-ccw/aaw
    %aan=aa; bbn=bb;
    rm=corrcoef(aa,bb);
    rr10=rm(1,2)^2;
    
    
    
    veridat(istep,1:10)=[rr1 rr2 rr3 rr4 rr5 rr6 rr7 rr8 rr9 rr10];
%    veridat(istep,1:30)=[R1 C1 rr1 R2 C2 rr2 R3 C3 rr3 R4 C4 rr4 R5 C5 rr5 R6 C6 rr6 R7 C7 rr7 R8 C8 rr8 R9 C9 rr9 R10 C10 rr10];
end   
    

save('/rootdir/s1/directory1/work1/temanncfr/data/veridatearlysmlagnolut_3drr','veridat','-ASCII');