
%%%%%%%%%%%%%% highf %%%%%%%%%%%%%%%


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 over 1896-1995 (for early-miss 1850-1895)
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);

% 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/directory2/work1/temanncfr/input/instpchigh','highf','-ASCII');
   save('/rootdir/s1/directory2/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/directory2/work1/temanncfr/input/proxy',num2str(istep)),'yyy','-ASCII'); % no normalize again this time
clear y zz yc yd zzc zzd yyy mmm nnnd nnnc
end

clear all;

for istep=1:19
    
    z=load(strcat('/rootdir/s1/directory2/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/directory2/work1/temanncfr/input/proxy',num2str(istep)),'stdkeepers','-ASCII');
end


for istep=1:19
    
    z=load(strcat('/rootdir/s1/directory2/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);
    %if ddd > 1
    %    dddd = ddd-1;
    %else
    %    dddd = ddd;
    %end
    %instpcstokeep(istep)=dddd;
    instpcstokeep(istep)=ddd;
end
   save('/rootdir/s1/directory2/work1/temanncfr/input/instpcstokeep','instpcstokeep','-ASCII');
   nnd=instpcstokeep';
   save('/rootdir/s1/directory2/work1/temanncfr/input/pcnnd1','nnd','-ASCII');


% now get input climate data.....
clear all; close all;
t=load('/rootdir/s1/directory2/work1/temanncfr/input/instpcstokeep');
cc=load('/rootdir/s1/directory2/work1/temanncfr/input/instpchigh'); % for high inst pcs
%cc=load('d:\scozztemann3\input\instpclow'); % for low inst pcs
for istep=1:19
    ddd=cc(1:100,1:t(istep));
   
    save(strcat('/rootdir/s1/directory2/work1/temanncfr/input/climate',num2str(istep)),'ddd','-ASCII');
end

%###########################################################################################
%%%%%%%% NOW YOU NEED TO RUN INPUT1.m (below) to get input data of high-f  %%%%%%%%%%%%
%###########################################################################################

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/0.136(90%)
%    row 5:  r -value for low-pass data over 1850-1995,         n=14 dof... |r|=0.496/0.425(90%)
%    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/itrdbmatrixh
x=itrdbmatrixh;
[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)
%------------------------------------------------------------------

%----------------------------------------------------------------
%                                                               %
%  If you want a non-dendro dataset, then you just need to      %
%  delete "z(3,i)==9000 | " from line 73 (9000 is tree-ring ID) %
%                                                               %
%----------------------------------------------------------------

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)

n=0;
for i=1:m1-1 % 
%    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)  &...
    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
        n=n+1;
        y(1:kkk-kk+1,n)=x(kk:kkk,i+1);
        zz(:,n)=z(:,i);
    end
end
num(istep,1)=n;

save(strcat('/rootdir/s1/directory2/work1/temanncfr/input/proxy',num2str(istep)),'y','-ASCII'); % no normalize again this time
clear y zz 
end

%----- estimate ttls values ----
clear all; close all;

for istep=1:19
    
    z=load(strcat('/rootdir/s1/directory2/work1/temanncfr/input/proxy',num2str(istep)));
        [nrows,ncols]=size(z);
    qqq=load(strcat('/rootdir/s1/directory2/work1/temanncfr/input/climate',num2str(istep)));
        
    ppp=z(nrows-99:nrows,:);%over 1850-1995
    cli=qqq;
%    ppp=z(nrows-99:nrows,:);%over 1896-1995 (for early-miss 1856-1895)
%    cli=qqq(41:140,:)
%    ppp=z(nrows-139:nrows-40,:);%over 1856-1955   (for late-miss 1956-1995)
%    cli=qqq(1:100,:)
    
    keepers=[cli ppp];
    [nrows,ncols]=size(keepers);
    means=nanmean(keepers); % 
    stdevs=nanstd(keepers);
    stdkeepers=keepers-repmat(means,nrows,1);
    stdkeepers=stdkeepers./repmat(stdevs,nrows,1);
    
    [U,S,V] = svd(stdkeepers,0); % do the decomposition, economy style. U's are the PCs
    highfevalues=diag(S.^2);
    highfpctvar=diag(S).^2./sum(diag(S).^2); % that gives the pct for each pc
    mbest=eigenselect(highfpctvar);
    eofnumb(istep)=mbest;
end
save('/rootdir/s1/directory2/work1/temanncfr/zzveri/allprogllagnolut_3d/highf/eofnumb','eofnumb','-ascii')
nnd=eofnumb';
save('/rootdir/s1/directory2/work1/temanncfr/zzveri/allprogllagnolut_3d/highf/eofnumbnnd','nnd','-ascii')

clear all; close all;

for istep=1:19
    
    z=load(strcat('/rootdir/s1/directory2/work1/temanncfr/input/proxy',num2str(istep)));
    qqq=load(strcat('/rootdir/s1/directory2/work1/temanncfr/input/climate',num2str(istep)));
    save(strcat('/rootdir/s1/directory2/work1/temanncfr/zzveri/allprogllagnolut_3d/highf/proxy',num2str(istep)),'z','-ASCII');
    save(strcat('/rootdir/s1/directory2/work1/temanncfr/zzveri/allprogllagnolut_3d/highf/climate',num2str(istep)),'qqq','-ASCII');
end



%%%%%%%%%%%%%% lowf %%%%%%%%%%%%%%%


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);

% 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/directory2/work1/temanncfr/input/instpchigh','highf','-ASCII');
   save('/rootdir/s1/directory2/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/directory2/work1/temanncfr/input/proxy',num2str(istep)),'yyy','-ASCII'); % no normalize again this time
clear y zz yc yd zzc zzd yyy mmm nnnd nnnc
end

clear all;

for istep=1:19
    
    z=load(strcat('/rootdir/s1/directory2/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/directory2/work1/temanncfr/input/proxy',num2str(istep)),'stdkeepers','-ASCII');
end
clear all;

for istep=1:19
    
    z=load(strcat('/rootdir/s1/directory2/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);
    %if ddd > 1
    %    dddd = ddd-1;
    %else
    %    dddd = ddd;
    %end
    %instpcstokeep(istep)=dddd;
    instpcstokeep(istep)=ddd;
end
   save('/rootdir/s1/directory2/work1/temanncfr/input/instpcstokeep','instpcstokeep','-ASCII');
   nnd=instpcstokeep';
   save('/rootdir/s1/directory2/work1/temanncfr/input/pcnnd2','nnd','-ASCII');

% now get input climate data.....
clear all; close all;
t=load('/rootdir/s1/directory2/work1/temanncfr/input/instpcstokeep');
cc=load('/rootdir/s1/directory2/work1/temanncfr/input/instpclow'); % for high inst pcs
%cc=load('d:\scozztemann3\input\instpclow'); % for low inst pcs
for istep=1:19
    ddd=cc(1:100,1:t(istep));
   
    save(strcat('/rootdir/s1/directory2/work1/temanncfr/input/climate',num2str(istep)),'ddd','-ASCII');
end

%###########################################################################################
%%%%%%%% NOW YOU NEED TO RUN INPUT1.m (below) to get input data of high-f  %%%%%%%%%%%%
%###########################################################################################

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/itrdbmatrixl
x=itrdbmatrixl;
[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/directory2/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

%----- estimate ttls values ----
clear all; close all;

for istep=1:19
    
    z=load(strcat('/rootdir/s1/directory2/work1/temanncfr/input/proxy',num2str(istep)));
        [nrows,ncols]=size(z);
    qqq=load(strcat('/rootdir/s1/directory2/work1/temanncfr/input/climate',num2str(istep)));
        
    ppp=z(nrows-99:nrows,:);%over 1850-1995
    cli=qqq;
%    ppp=z(nrows-99:nrows,:);%over 1896-1995 (for early-miss 1856-1895)
%    cli=qqq(41:140,:)
%    ppp=z(nrows-139:nrows-40,:);%over 1856-1955   (for late-miss 1956-1995)
%    cli=qqq(1:100,:)
    
    keepers=[cli ppp];
    [nrows,ncols]=size(keepers);
    means=nanmean(keepers); % 
    stdevs=nanstd(keepers);
    stdkeepers=keepers-repmat(means,nrows,1);
    stdkeepers=stdkeepers./repmat(stdevs,nrows,1);
    
    [U,S,V] = svd(stdkeepers,0); % do the decomposition, economy style. U's are the PCs
    lowfevalues=diag(S.^2);
    lowfpctvar=diag(S).^2./sum(diag(S).^2); % that gives the pct for each pc
    cumlowfpctvar=cumsum(lowfpctvar);
    ttlslowfpar=length(find(cumlowfpctvar<0.33))+1;
    if ttlslowfpar==0
    	ttlslowfpar=1;
    end
    eofnumb(istep)=ttlslowfpar;
end
save('/rootdir/s1/directory2/work1/temanncfr/zzveri/allprogllagnolut_3d/lowf/eofnumb','eofnumb','-ascii')
nnd=eofnumb';
save('/rootdir/s1/directory2/work1/temanncfr/zzveri/allprogllagnolut_3d/lowf/eofnumbnnd','nnd','-ascii')

clear all; close all;

for istep=1:19
    
    z=load(strcat('/rootdir/s1/directory2/work1/temanncfr/input/proxy',num2str(istep)));
    qqq=load(strcat('/rootdir/s1/directory2/work1/temanncfr/input/climate',num2str(istep)));
    save(strcat('/rootdir/s1/directory2/work1/temanncfr/zzveri/allprogllagnolut_3d/lowf/proxy',num2str(istep)),'z','-ASCII');
    save(strcat('/rootdir/s1/directory2/work1/temanncfr/zzveri/allprogllagnolut_3d/lowf/climate',num2str(istep)),'qqq','-ASCII');
end


