% preparing input data of reconstruction based on screened proxy network

%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
% high-f proxy for full/subproxy network
%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

clear all; close all;

z=load('/holocene/s1/zuz10/work1/temann/data/rtable1209');
%---------------------------------------------------
load /holocene/s1/zuz10/work1/temann/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;   
%%--------------------------
%------------------------------------------------------------------

%----------------------------------------------------------------
%                                                               %
%  If you want a non-dendro dataset, then you just need to      %
%  delete "z(3,i)==9000/1 (9000/1 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); 

%--------------------------------------
nc1=0; yc1=[]; zzc1=[];
for i=1:m1-1 % 
    if (z(3,i)==9000 |  z(3,i)==8000 |  z(3,i)==7500 | 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 & (z(4,i))>=0.106
        nc1=nc1+1;
        yc1(1:kkk-kk+1,nc1)=x(kk:kkk,i+1);
        zzc1(:,nc1)=z(:,i);
    end
end
nc2=0; yc2=[]; zzc2=[];
for i=1:m1-1 % 
    if (z(3,i)==7000)  &...
         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 & (z(4,i))<=-0.106
        nc2=nc2+1;
        yc2(1:kkk-kk+1,nc2)=x(kk:kkk,i+1);
        zzc2(:,nc2)=z(:,i);
    end
end
nc3=0; yc3=[]; zzc3=[];
for i=1:m1-1 % 
    if (z(3,i)==6000 |  z(3,i)==5000)  &...
         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(4,i))>=0.136
        nc3=nc3+1;
        yc3(1:kkk-kk+1,nc3)=x(kk:kkk,i+1);
        zzc3(:,nc3)=z(:,i);
    end
end

y=[yc1 yc2 yc3];
zz=[zzc1 zzc2 zzc3];

%num(istep,1)=n;

%save(strcat('/holocene/s1/zuz10/work1/temann/info/info',num2str(istep)),'zz','-ASCII');
save(strcat('/holocene/s1/zuz10/work1/temann/zzrecon1209/nhnhscrihad/highf/proxy',num2str(istep)),'y','-ASCII'); % no normalize again this time
clear y zz yc1 yc2 yc3 zzc1 zzc2 zzc3
end
%save('/holocene/s1/zuz10/work1/temann/info/numh','num','-ASCII');

%----- do the climate data and ttls parameters ----
clear all; close all;

qqq=load('/holocene/s1/zuz10/work1/temann/newtemp/nhhinfxxx'); % ---nh temp


climate=qqq; % do not need to normalize again this time
save('/holocene/s1/zuz10/work1/temann/zzrecon1209/nhnhscrihad/highf/climate','climate','-ascii')


%-- estimate ttls values -----
for istep=1:19
    
    z=load(strcat('/holocene/s1/zuz10/work1/temann/zzrecon1209/nhnhscrihad/highf/proxy',num2str(istep)));
    [nrows,ncols]=size(z)
    
    ppp=z(nrows-145:nrows,:);%over 1850-1995
    cli=qqq
    
    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('/holocene/s1/zuz10/work1/temann/zzrecon1209/nhnhscrihad/highf/eofnumb','eofnumb','-ascii')
%------------------------------------------------------------------





%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
% low-f proxy for full/subproxy network
%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

clear all; close all;

z=load('/holocene/s1/zuz10/work1/temann/data/rtable1209');

%---------------------------------------------------
load /holocene/s1/zuz10/work1/temann/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;   
%%--------------------------
%------------------------------------------------------------------


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)


nc1=0; yc1=[]; zzc1=[];
for i=1:m1-1 % 
    if (z(3,i)==9000 |  z(3,i)==8000 |  z(3,i)==7500 | 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 & (z(4,i))>=0.106
        nc1=nc1+1;
        yc1(1:kkk-kk+1,nc1)=x(kk:kkk,i+1);
        zzc1(:,nc1)=z(:,i);
    end
end
nc2=0; yc2=[]; zzc2=[];
for i=1:m1-1 % 
    if (z(3,i)==7000)  &...
         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 & (z(4,i))<=-0.106
        nc2=nc2+1;
        yc2(1:kkk-kk+1,nc2)=x(kk:kkk,i+1);
        zzc2(:,nc2)=z(:,i);
    end
end
nc3=0; yc3=[]; zzc3=[];
for i=1:m1-1 % 
    if (z(3,i)==6000 |  z(3,i)==5000)  &...
         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(4,i))>=0.136
        nc3=nc3+1;
        yc3(1:kkk-kk+1,nc3)=x(kk:kkk,i+1);
        zzc3(:,nc3)=z(:,i);
    end
end


nd1=0; yd1=[]; zzd1=[];
for i=1:m1-1 % 
    if (z(3,i)==9001 |  z(3,i)==8001 |  z(3,i)==7501 | z(3,i)==4001 |  z(3,i)==3001 |  z(3,i)==2001)  &...
         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 & (z(5,i))>=0.337
        nd1=nd1+1;
        yd1(1:kkk-kk+1,nd1)=x(kk:kkk,i+1);
        zzd1(:,nd1)=z(:,i);
    end
end
nd2=0; yd2=[]; zzd2=[];
for i=1:m1-1 % 
    if (z(3,i)==7001)  &...
         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 & (z(5,i))<=-0.337
        nd2=nd2+1;
        yd2(1:kkk-kk+1,nd2)=x(kk:kkk,i+1);
        zzd2(:,nd2)=z(:,i);
    end
end
nd3=0; yd3=[]; zzd3=[];
for i=1:m1-1 % 
    if (z(3,i)==6001 |  z(3,i)==5001)  &...
         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(5,i))>=0.425
        nd3=nd3+1;
        yd3(1:kkk-kk+1,nd3)=x(kk:kkk,i+1);
        zzd3(:,nd3)=z(:,i);
    end
end

%num(istep,1)=nc+nd;
num(istep,1)=nc1+nc2+nc3;
num(istep,2)=nd1+nd2+nd3;
y=[yc1 yc2 yc3 yd1 yd2 yd3];
zz=[zzc1 zzc2 zzc3 zzd1 zzd2 zzd3];

%save(strcat('/holocene/s1/zuz10/work1/temann/info/info',num2str(istep)),'zz','-ASCII');
save(strcat('/holocene/s1/zuz10/work1/temann/zzrecon1209/nhnhscrihad/lowf/proxy',num2str(istep)),'y','-ASCII'); % no normalize again this time
clear y zz yc1 yc2 yc3 yd1 yd2 yd3 zzc1 zzc2 zzc3 zzd1 zzd2 zzd3
end
save('/holocene/s1/zuz10/work1/temann/info/numl','num','-ASCII');


%----- do the climate data and ttls parameters ----
clear all; close all;

qqq=load('/holocene/s1/zuz10/work1/temann/newtemp/nhlinfxxx');  % ---nh temp


climate=qqq; % do not need to normalize again this time
save('/holocene/s1/zuz10/work1/temann/zzrecon1209/nhnhscrihad/lowf/climate','climate','-ascii')

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

%-- estimate ttls values -----
for istep=1:19
    
    z=load(strcat('/holocene/s1/zuz10/work1/temann/zzrecon1209/nhnhscrihad/lowf/proxy',num2str(istep)));
    [nrows,ncols]=size(z)
    
    ppp=z(nrows-145:nrows,:);%over 1850-1995
    cli=qqq
    
    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.5))+1;
    if ttlslowfpar==0
    	ttlslowfpar=1;
    end
    eofnumb(istep)=ttlslowfpar;
end
save('/holocene/s1/zuz10/work1/temann/zzrecon1209/nhnhscrihad/lowf/eofnumb','eofnumb','-ascii')
%------------------------------------------------------------------
%------------------------------------------------------------------


