MAKEPROXY;
clear all; close all;
load ./instdata/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(1:146,:);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
means=nanmean(keepers(1:146,:));
stdevs=nanstd(keepers(1:146,:));
stdkeepers=keepers-repmat(means,146,1);
stdkeepers=stdkeepers./repmat(stdevs,146,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,146,1);
instpcstds=nanstd(instpcs(:,:));
instpcs=instpcs./repmat(instpcstds,146,1);

save('./recon/instmeans','means','-ASCII');
save('./recon/inststds','stdevs','-ASCII');

save('./recon/instpcmeans','instpcmeans','-ASCII');
save('./recon/instpcstds','instpcstds','-ASCII');
save('./recon/u','U','-ASCII');
save('./recon/v','V','-ASCII');
save('./recon/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('./recon/instpchigh','highf','-ASCII');
   save('./recon/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('./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 ./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(4,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(5,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('./recon/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('./recon/proxy',num2str(istep)));
    keepers=z;
    [nrows,ncols]=size(keepers)

    means=nanmean(keepers(nrows-145:nrows,:)); % normalization over 1850-1995/ or whatever intervals
    stdevs=nanstd(keepers(nrows-145:nrows,:));
    stdkeepers=keepers-repmat(means,nrows,1);
    stdkeepers=stdkeepers./repmat(stdevs,nrows,1);
    
    save(strcat('./recon/proxy',num2str(istep)),'stdkeepers','-ASCII');
end


for istep=1:19
    
    z=load(strcat('./recon/proxy',num2str(istep)));
    keepers=z;
    [nrows,ncols]=size(keepers)

    stepproxies=keepers(nrows-145: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;
end
   save('./recon/instpcstokeep','instpcstokeep','-ASCII');

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

clear all; close all;

instpcmeans=load('./recon/instpcmeans');
instpcstds=load('./recon/instpcstds');
u=load('./recon/u');
v=load('./recon/v');
s=load('./recon/s');
%t=load('/rootdir/s1/directory1/work1/temanncfr/input/instpcstokeep');
for istep=1:19
    
    zh=load(strcat('./highf/recon',num2str(istep)));
    zl=load(strcat('./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('./recon/recon',num2str(istep)),'qqq','-ASCII');
end


clear all; close all;

instmeans=load('./recon/instmeans');
inststds=load('./recon/inststds');

for istep=1:19
    
    zzz=load(strcat('./recon/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('./recon/recon',num2str(istep)),'zzz','-ASCII');
end




%=====================================================

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


%A)	Glb---(87.5S-87.5N, 177.5W-177.5E)
%B)	NH---(2.5N-87.5N, 177.5W-177.5E)
%C)	SH—(2.5S-87.5S, 177.5W-177.5E)
%D)	Trop---(27.5S-27.5N, 177.5W-177.5E)
%E)	Xtrop---(27.5N-67.5N, 177.5W-177.5E)
%F)	NH land—I used Scott’s landmask data to find the gridboxes over NH land.
%G)	Nino3---(2.5S-2.5N, 92.5W-147.5W)
%H)	PDO---(22.5N-57.5N, 152.5E-132.5W)
%I)	AMO---[(2.5N-7.5N, 12.5W-52.5W), (12.5N-37.5N, 17.5W-72.5W), (42.5N-57.5N, 12.5W-52.5W), (62.5N-67.5N, 37.5W-7.5E)]
%J)	Atlantic main development region---(7.5N-17.5N, 17.5W-62.5W)




%--------- 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 ./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('./data/tcoswt','tcoswt','-ASCII');
save('./data/masknh','masknh','-ASCII');
save('./data/masksh','masksh','-ASCII');
save('./data/maskgl','tt','-ASCII');
save('./data/maskextro','maskextro','-ASCII');
save('./data/masktro','masktro','-ASCII');
save('./data/masknino','masknino','-ASCII');
save('./data/maskpdo','maskpdo','-ASCII');
save('./data/maskamo','maskamo','-ASCII');
save('./data/masknhland','masknhland','-ASCII');
save('./data/masktamdr','masktamdr','-ASCII');



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

clear all; close all;

load ./data/tcoswt
load ./data/masknh
load ./data/masksh
load ./data/maskgl
load ./data/maskextro
load ./data/masktro
load ./data/masknino
load ./data/maskpdo
load ./data/maskamo
load ./data/masknhland
load ./data/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 ./instdata/enough_3d
load ./instdata/toofew_3d
%  gridd([enough_3d])=gridcont;
  gridd=(1:1:2592);
  gridd([enough_3d])=1;
  gridd([toofew_3d])=NaN;

%###################################################################



t=load('./instdata/instfilled_3d');
tt=t(147:157,:); % for adding inst data from 1996-2006
ttt=t(1:157,:); % to replace EOF instrumental data with infilled climate (1850-2006)

for istep=1:19
% istep=5

recon=[];

    instfilled_3d1=load(strcat('./recon/recon',num2str(istep)));
    [iyy,jjj]=size(instfilled_3d1);
    instfilled_3d(1:iyy+11,1:jjj)=NaN;
    instfilled_3d(1:iyy,1:jjj)=instfilled_3d1;
    instfilled_3d(iyy+1:iyy+11,1:jjj)=tt;
    instfilled_3d(iyy-145:iyy+11,1:jjj)=ttt;% to replace EOF instrumental data with infilled climate (1850-2006
    inf3=instfilled_3d;
    inf3(:,[enough_3d])=inf3; % back to 36*72 gridpoints
    inf3(:,[toofew_3d])=NaN;
    
%---- now spatialtemporal mask for each step ----

ztcoswt=repmat(tcoswt,207+(istep-1)*100,1);

zmasknh=masknh.*gridd;
zmasknh=repmat(zmasknh,207+(istep-1)*100,1);
zmasksh=masksh.*gridd;
zmasksh=repmat(zmasksh,207+(istep-1)*100,1);
zmaskgl=maskgl.*gridd;
zmaskgl=repmat(zmaskgl,207+(istep-1)*100,1);
zmaskextro=maskextro.*gridd;
zmaskextro=repmat(zmaskextro,207+(istep-1)*100,1);
zmasktro=masktro.*gridd;
zmasktro=repmat(zmasktro,207+(istep-1)*100,1);
zmasknino=masknino.*gridd;
zmasknino=repmat(zmasknino,207+(istep-1)*100,1);
zmaskpdo=maskpdo.*gridd;
zmaskpdo=repmat(zmaskpdo,207+(istep-1)*100,1);
zmaskamo=maskamo.*gridd;
zmaskamo=repmat(zmaskamo,207+(istep-1)*100,1);
zmasknhland=masknhland.*gridd;
zmasknhland=repmat(zmasknhland,207+(istep-1)*100,1);
zmasktamdr=masktamdr.*gridd;
zmasktamdr=repmat(zmasktamdr,207+(istep-1)*100,1);
%--------------------------------------------    
    %--- NH mean ---
    www=inf3./inf3;
    www=(www.*zmasknh.*ztcoswt)';
    www=nanmean(www);
    www=www';
    bbb=(inf3.*zmasknh.*ztcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    a1=bbb./www;

    
    %--- SH mean ---
    www=inf3./inf3;
    www=(www.*zmasksh.*ztcoswt)';
    www=nanmean(www);
    www=www';
    bbb=(inf3.*zmasksh.*ztcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    a2=bbb./www;
    
    %--- GL mean ---
    www=inf3./inf3;
    www=(www.*zmaskgl.*ztcoswt)';
    www=nanmean(www);
    www=www';
    bbb=(inf3.*zmaskgl.*ztcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    a3=bbb./www;
    
    %--- Extropical mean ---
    www=inf3./inf3;
    www=(www.*zmaskextro.*ztcoswt)';
    www=nanmean(www);
    www=www';
    bbb=(inf3.*zmaskextro.*ztcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    a4=bbb./www;
    
    %--- Tropical mean ---
    www=inf3./inf3;
    www=(www.*zmasktro.*ztcoswt)';
    www=nanmean(www);
    www=www';
    bbb=(inf3.*zmasktro.*ztcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    a5=bbb./www;
    
    %--- Nino3 mean ---
    www=inf3./inf3;
    www=(www.*zmasknino.*ztcoswt)';
    www=nanmean(www);
    www=www';
    bbb=(inf3.*zmasknino.*ztcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    a6=bbb./www;
    
    %--- PDO region mean ---
    www=inf3./inf3;
    www=(www.*zmaskpdo.*ztcoswt)';
    www=nanmean(www);
    www=www';
    bbb=(inf3.*zmaskpdo.*ztcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    a7=bbb./www;
    
    %--- AMO region mean ---
    www=inf3./inf3;
    www=(www.*zmaskamo.*ztcoswt)';
    www=nanmean(www);
    www=www';
    bbb=(inf3.*zmaskamo.*ztcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    a8=bbb./www;
    
    %--- NH land mean ---
    www=inf3./inf3;
    www=(www.*zmasknhland.*ztcoswt)';
    www=nanmean(www);
    www=www';
    bbb=(inf3.*zmasknhland.*ztcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    a9=bbb./www;
    
    %--- tropical Atlantic Main development region mean ---
    www=inf3./inf3;
    www=(www.*zmasktamdr.*ztcoswt)';
    www=nanmean(www);
    www=www';
    bbb=(inf3.*zmasktamdr.*ztcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    a10=bbb./www;
    
    recon=[a1 a2 a3 a4 a5 a6 a7 a8 a9 a10]; 
    save(strcat('./recon/reconserieslagnolut_3d',num2str(istep)),'recon','-ASCII');
    save(strcat('./recon/reconserieslagnolut_3dgrid',num2str(istep)),'inf3','-ASCII');
    
    clear instfilled_3d1 instfilled_3d www bbb recon ztcoswt zmasknh zmasksh zmaskgl zmaskextro zmasktro zmasknino zmaskpdo zmaskamo zmasknhland zmasktamdr
end   
    



%###################################################################
%---------------------------
%////////////// reginal/global series of reconstructions (decadal) ////////////////////////
%---------------------------
%###################################################################

for istep=1:19

recon=[];

    instfilled_3d1=load(strcat('./recon/recon',num2str(istep)));
    [iyy,jjj]=size(instfilled_3d1);
    instfilled_3d(1:iyy+11,1:jjj)=NaN;
    instfilled_3d(1:iyy,1:jjj)=instfilled_3d1;
    instfilled_3d(iyy+1:iyy+11,1:jjj)=tt;
    inf3=instfilled_3d;
%%------------- 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(:,[enough_3d])=inf3; % back to 36*72 gridpoints
  inf3(:,[toofew_3d])=NaN;
%%-----------------------------------------------------------------
 %---- now spatialtemporal mask for each step ----

ztcoswt=repmat(tcoswt,207+(istep-1)*100,1);

zmasknh=masknh.*gridd;
zmasknh=repmat(zmasknh,207+(istep-1)*100,1);
zmasksh=masksh.*gridd;
zmasksh=repmat(zmasksh,207+(istep-1)*100,1);
zmaskgl=maskgl.*gridd;
zmaskgl=repmat(zmaskgl,207+(istep-1)*100,1);
zmaskextro=maskextro.*gridd;
zmaskextro=repmat(zmaskextro,207+(istep-1)*100,1);
zmasktro=masktro.*gridd;
zmasktro=repmat(zmasktro,207+(istep-1)*100,1);
zmasknino=masknino.*gridd;
zmasknino=repmat(zmasknino,207+(istep-1)*100,1);
zmaskpdo=maskpdo.*gridd;
zmaskpdo=repmat(zmaskpdo,207+(istep-1)*100,1);
zmaskamo=maskamo.*gridd;
zmaskamo=repmat(zmaskamo,207+(istep-1)*100,1);
zmasknhland=masknhland.*gridd;
zmasknhland=repmat(zmasknhland,207+(istep-1)*100,1);
zmasktamdr=masktamdr.*gridd;
zmasktamdr=repmat(zmasktamdr,207+(istep-1)*100,1);
%--------------------------------------------    
    %--- NH mean ---
    www=inf3./inf3;
    www=(www.*zmasknh.*ztcoswt)';
    www=nanmean(www);
    www=www';
    bbb=(inf3.*zmasknh.*ztcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    a1=bbb./www;

    
    %--- SH mean ---
    www=inf3./inf3;
    www=(www.*zmasksh.*ztcoswt)';
    www=nanmean(www);
    www=www';
    bbb=(inf3.*zmasksh.*ztcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    a2=bbb./www;
    
    %--- GL mean ---
    www=inf3./inf3;
    www=(www.*zmaskgl.*ztcoswt)';
    www=nanmean(www);
    www=www';
    bbb=(inf3.*zmaskgl.*ztcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    a3=bbb./www;
    
    %--- Extropical mean ---
    www=inf3./inf3;
    www=(www.*zmaskextro.*ztcoswt)';
    www=nanmean(www);
    www=www';
    bbb=(inf3.*zmaskextro.*ztcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    a4=bbb./www;
    
    %--- Tropical mean ---
    www=inf3./inf3;
    www=(www.*zmasktro.*ztcoswt)';
    www=nanmean(www);
    www=www';
    bbb=(inf3.*zmasktro.*ztcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    a5=bbb./www;
    
    %--- Nino3 mean ---
    www=inf3./inf3;
    www=(www.*zmasknino.*ztcoswt)';
    www=nanmean(www);
    www=www';
    bbb=(inf3.*zmasknino.*ztcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    a6=bbb./www;
    
    %--- PDO region mean ---
    www=inf3./inf3;
    www=(www.*zmaskpdo.*ztcoswt)';
    www=nanmean(www);
    www=www';
    bbb=(inf3.*zmaskpdo.*ztcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    a7=bbb./www;
    
    %--- AMO region mean ---
    www=inf3./inf3;
    www=(www.*zmaskamo.*ztcoswt)';
    www=nanmean(www);
    www=www';
    bbb=(inf3.*zmaskamo.*ztcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    a8=bbb./www;
    
    %--- NH land mean ---
    www=inf3./inf3;
    www=(www.*zmasknhland.*ztcoswt)';
    www=nanmean(www);
    www=www';
    bbb=(inf3.*zmasknhland.*ztcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    a9=bbb./www;
    
    %--- tropical Atlantic Main development region mean ---
    www=inf3./inf3;
    www=(www.*zmasktamdr.*ztcoswt)';
    www=nanmean(www);
    www=www';
    bbb=(inf3.*zmasktamdr.*ztcoswt)';
    bbb=nanmean(bbb);
    bbb=bbb';
    a10=bbb./www;
    
    recon=[a1 a2 a3 a4 a5 a6 a7 a8 a9 a10]; 
    save(strcat('./recon/reconseriessmlagnolut_3d',num2str(istep)),'recon','-ASCII');
    save(strcat('./recon/reconseriessmlagnolut_3dgrid',num2str(istep)),'inf3','-ASCII');
    
    clear www bbb recon ztcoswt zmasknh zmasksh zmaskgl zmaskextro zmasktro zmasknino zmaskpdo zmaskamo zmasknhland zmasktamdr
end   
   
