 % preparing high-f/low-f hemispherical land(iCRU) or land+ocean(iHAD) mean surface temperature
 % series based on our infilled global instrumental gridbox.

 clear all; close all;

mask=load('/holocene/s1/zuz10/work1/temann/data/lsmask_jones.dat');
landmask=mask./mask; % this makes the ocean NaN and all others 1
landmask=flipdim(landmask,1); % make initial begins with SH (-87.5) which match to Hadcru structure
landmask=reshape(landmask',1,2592); % now string out the land mask 1:72*36

ncols=length(landmask)
nmis=0;
good=1;
for i=1:ncols
        if landmask(i)> -99999 %(atually land==1)
                lands(good)=i; % land gridboxes
        good=good+1;
        else
        nmis=nmis+1;
        oceans(nmis)=i; % ocean gridboxes
        end
end
    

load /holocene/s1/zuz10/hadcru3v/infillannt/instfilled
load /holocene/s1/zuz10/hadcru3v/enough
load /holocene/s1/zuz10/hadcru3v/toofew

  
  %###########################################
  %////////////////////////////////////////////
  % use % for line 48 and 61, then change land-->laoc to get the hemisphere mean (from land mean)
  %////////////////////////////////////////////
  %###########################################
  
  annfull(:,[enough])=instfilled;
  %-----
  annfull(:,[toofew])=NaN; % for land + ocean----iHAD
  %-----
%  annfull(:,[toofew])=NaN; 
%  annfull(:,[oceans])=NaN;  % only land gridboxes--iCRU
  %------

% notice the follow for estimating area weight parameter < latwt >  
% the first 72 data are SH -87.5 and -177.5 -172.5-------172.5 177.5, then SH -82.5 and -177.5 -172.5-------172.5 177.5
% the last 72 data are NH 87.5 and -177.5 -172.5-------172.5 177.5.

for i=1:36
    for j=1:72
        latwt((i-1)*72+j)=cos(((i-1)*5-87.5)*pi/180);
    end
end
%---------
  latwt([toofew])=NaN; % for land + ocean----iHAD
%---------
%  latwt([toofew])=NaN; 
%  latwt([oceans])=NaN; %% only land gridboxes--iCRU
%---------

% for South Hemisphere Mean
shrecon=annfull(:,1:1296);
shlatwt=latwt(1:1296);
shrecon=shrecon.*repmat(shlatwt,157,1);
shreconmean=nansum(shrecon')./nansum(shlatwt');

% for North Hemisphere Mean
nhrecon=annfull(:,1297:2592);
nhlatwt=latwt(1297:2592);
nhrecon=nhrecon.*repmat(nhlatwt,157,1);
nhreconmean=nansum(nhrecon')./nansum(nhlatwt');

% for Global Mean
glreconmean=(shreconmean+nhreconmean)/2;

% make all ocean gridboxes = NaNs, we get land mean or don't make all ocean gridboxes = NaNs, we get land+ocean (laoc) mean
save('/holocene/s1/zuz10/work1/temann/newtemp/shxxx','shreconmean','-ascii')
save('/holocene/s1/zuz10/work1/temann/newtemp/nhxxx','nhreconmean','-ascii')
save('/holocene/s1/zuz10/work1/temann/newtemp/glxxx','glreconmean','-ascii')
 


% /////////////work place = /holocene/s1/zuz10/work1/temann/dodata/////////////////////////

%----------------------------
clear all; close all

NH=load('/holocene/s1/zuz10/work1/temann/newtemp/nhxxx');
SH=load('/holocene/s1/zuz10/work1/temann/newtemp/shxxx');
GG=load('/holocene/s1/zuz10/work1/temann/newtemp/glxxx');

aaa=(1850:1:2006)';
%iCRU_NH_reform=[aaa NH'];
%iCRU_SH_reform=[aaa SH'];
%save('/holocene/s1/zuz10/work1/temann/newtemp/iCRU_NH_reform','iCRU_NH_reform','-ascii')
%save('/holocene/s1/zuz10/work1/temann/newtemp/iCRU_SH_reform','iCRU_SH_reform','-ascii')
%iHAD_NH_reform=[aaa NH'];
%iHAD_SH_reform=[aaa SH'];
%save('/holocene/s1/zuz10/work1/temann/newtemp/iHAD_NH_reform','iHAD_NH_reform','-ascii')
%save('/holocene/s1/zuz10/work1/temann/newtemp/iHAD_SH_reform','iHAD_SH_reform','-ascii')



nt=NH(1:146)';
st=SH(1:146)';
gl=GG(1:146)';

save('/holocene/s1/zuz10/work1/temann/newtemp/tempn95infxxx','nt','-ascii')
save('/holocene/s1/zuz10/work1/temann/newtemp/temps95infxxx','st','-ascii')
save('/holocene/s1/zuz10/work1/temann/newtemp/tempg95infxxx','gl','-ascii')

%------------ get standardized NH/SM/GL mean (1850-1995)
respliced=nt;
% standardize sst/sat based on interval 1850-1995 for the full series(1850-1998)
[nrows,ncols]=size(respliced)
means=nanmean(respliced(1:146,:));
stdevs=nanstd(respliced(1:146,:));
respliced=respliced-repmat(means,nrows,1);
respliced=respliced./repmat(stdevs,nrows,1);
stdtempn95=respliced;
save('/holocene/s1/zuz10/work1/temann/newtemp/stdtempn95infxxx','stdtempn95','-ascii')

respliced=st;
% standardize sst/sat based on interval 1850-1995 for the full series(1850-1998)
[nrows,ncols]=size(respliced)
means=nanmean(respliced(1:146,:));
stdevs=nanstd(respliced(1:146,:));
respliced=respliced-repmat(means,nrows,1);
respliced=respliced./repmat(stdevs,nrows,1);
stdtemps95=respliced;
save('/holocene/s1/zuz10/work1/temann/newtemp/stdtemps95infxxx','stdtemps95','-ascii')

respliced=gl;
% standardize sst/sat based on interval 1850-1995 for the full series(1850-1998)
[nrows,ncols]=size(respliced)
means=nanmean(respliced(1:146,:));
stdevs=nanstd(respliced(1:146,:));
respliced=respliced-repmat(means,nrows,1);
respliced=respliced./repmat(stdevs,nrows,1);
stdtempg95=respliced;
save('/holocene/s1/zuz10/work1/temann/newtemp/stdtempg95infxxx','stdtempg95','-ascii')

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

clear all

load /holocene/s1/zuz10/work1/temann/newtemp/stdtempn95infxxx

    cutoff=20; %years
	lowlim=1./(cutoff./2)

stdkeepers=stdtempn95infxxx; % because standardization have done in doannsst.m


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[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.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('/holocene/s1/zuz10/work1/temann/newtemp/nhhinfxxx','highf','-ascii')
save('/holocene/s1/zuz10/work1/temann/newtemp/nhlinfxxx','lowf','-ascii')
%save('/holocene/s1/zuz10/work1/temann/newtemp/nhhinfxxxicru','highf','-ascii')
%save('/holocene/s1/zuz10/work1/temann/newtemp/nhlinfxxxicru','lowf','-ascii')
%save('/holocene/s1/zuz10/work1/temann/newtemp/nhhinfxxxihad','highf','-ascii')
%save('/holocene/s1/zuz10/work1/temann/newtemp/nhlinfxxxihad','lowf','-ascii')

%%%%%%%

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

clear all

load /holocene/s1/zuz10/work1/temann/newtemp/stdtemps95infxxx

    cutoff=20; %years
	lowlim=1./(cutoff./2)

stdkeepers=stdtemps95infxxx; % because standardization have done in doannsst.m


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[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.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('/holocene/s1/zuz10/work1/temann/newtemp/shhinfxxx','highf','-ascii')
save('/holocene/s1/zuz10/work1/temann/newtemp/shlinfxxx','lowf','-ascii')
%save('/holocene/s1/zuz10/work1/temann/newtemp/shhinfxxxicru','highf','-ascii')
%save('/holocene/s1/zuz10/work1/temann/newtemp/shlinfxxxicru','lowf','-ascii')
%save('/holocene/s1/zuz10/work1/temann/newtemp/shhinfxxxihad','highf','-ascii')
%save('/holocene/s1/zuz10/work1/temann/newtemp/shlinfxxxihad','lowf','-ascii')

%%%%%%%

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

clear all

load /holocene/s1/zuz10/work1/temann/newtemp/stdtempg95infxxx

    cutoff=20; %years
	lowlim=1./(cutoff./2)

stdkeepers=stdtempg95infxxx; % because standardization have done in doannsst.m


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[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.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('/holocene/s1/zuz10/work1/temann/newtemp/glhinfxxx','highf','-ascii')
save('/holocene/s1/zuz10/work1/temann/newtemp/gllinfxxx','lowf','-ascii')

%%%%%%%

