
clear all;

% program to do proxy gridbox CPS against instrumental data 
% Created December 2006 by skm (from earlier work)
% Last modified October 2007 by skm

%%% ****run gridproxy.m first to get screened proxies into directory './data/cps/...'

period='la100'; %% choices: 'whole', 'ea100', 'la100'
expe='_noluter'; %% choices: '', '_nodendro', '_noluter'

sinst=1850;
einst=2006;
cutfreq=0.1; %cut off frequency for lowpass

switch lower(period)
  case 'whole'
    scalib=1850; %start of calibration period
    ecalib=1995; %end of calibration period
    svalid=1850; %start of validation period
    evalid=1995; %end of validation period
  case 'ea100'
    scalib=1850; %start of calibration period
    ecalib=1949; %end of calibration period
    svalid=1950; %start of validation period
    evalid=1995; %end of validation period
  case 'la100'
    scalib=1896; %start of calibration period
    ecalib=1995; %end of calibration period
    svalid=1850; %start of validation period
    evalid=1895; %end of validation period
end

pdir=strcat('data/cps/',period,expe,'/')

load('data/cps/infilled2poles_1850');
%% loads anngridinst and keeplonglat

expe=strcat(expe,'_newgrid');

% save dimensions of instrument data
[nyears,ngrid]=size(anngridinst)
smoothinst=NaN([nyears, ngrid]);

% for each gridbox, filter the instrument data
for i=1:ngrid
  [smoothed, icb, ice, mse0]=lowpassmin(anngridinst(:,i), cutfreq); 
  smoothinst(:,i)=smoothed;
  %lowpass filter with a cutoff frequency of 0.05 -- 20 years
end

clear anngridinst;

nyears_ptot=1996; % years 0 - 1995
reg_ts=NaN([nyears_ptot,ngrid,19]);

for n=1:19 %loop over each proxy set (back to year n*100)
  spyr=1900-(n*100);
  epyr=1995;
   load([pdir,'grid_GL_proxy',int2str(n),'.mat']) %gets proxy data
   rebuiltproxymatrix=y;
   usedproxylonglat=loc'; 
   [nyears_p,nproxy]=size(rebuiltproxymatrix);
   smoothproxy=NaN(nyears_p,nproxy);
   gridproxy=zeros(nyears_p,ngrid);
   cntproxy=zeros([nyears_p,ngrid]);
   
   for i=1:nproxy % normalize
     tmp_proxy=rebuiltproxymatrix(:,i);
     meancal=nanmean(tmp_proxy(scalib-spyr+1:ecalib-spyr+1));
     stdcal=nanstd(tmp_proxy(scalib-spyr+1:ecalib-spyr+1));
     tmp_proxy=(tmp_proxy-meancal)/stdcal;
     smoothproxy(1:nyears_p,i)=tmp_proxy;
   end
   
   clear rebuiltproxymatrix;

   % average all proxies for that record into 5x5 grid boxes.
    glat=keeplonglat(:,2)';
    glon=keeplonglat(:,1)';
    for i=1:ngrid
        % grid box corners
        gwest=glon(i)-2.5;
        geast=glon(i)+2.5;
        gsouth=glat(i)-2.5;
        gnorth=glat(i)+2.5;
        proxybox=find(usedproxylonglat(:,1) >= gwest & usedproxylonglat(:,1) <= geast & usedproxylonglat(:,2) >= gsouth & usedproxylonglat(:,2) <= gnorth);            
        if (~isempty(proxybox))
            for yr=1:nyears_p
                tmp_proxy=smoothproxy(yr,proxybox);
                %inds=find(isfinite(tmp_proxy));
                gridproxy(yr,i)=gridproxy(yr,i)+nansum(tmp_proxy);
                cntproxy(yr,i)=cntproxy(yr,i)+numel(proxybox);
            end % for yr
        end % if isempty
    end % for i=1:ngrid

    inds=find(cntproxy > 0);
    ibad=find(cntproxy == 0);
    gridproxy(inds)=gridproxy(inds)./cntproxy(inds);
    gridproxy(ibad)=NaN;
 
    % calibrate proxies against instrument data
    % period
    for i=1:ngrid
      calib_inst=smoothinst(scalib-sinst+1:ecalib-sinst+1,i);
      calib_proxy=gridproxy(scalib-spyr+1:ecalib-spyr+1,i);
      inds=find(isfinite(calib_proxy));
      if (~isempty(inds))
       meaninst=mean(calib_inst);
       stdinst=std(calib_inst);
       stdprox=nanstd(calib_proxy);
       reg_ts(spyr+1:epyr+1,i,n)=(gridproxy(:,i)*stdinst/stdprox)+meaninst;
      end
    end
end %for n=1:19

%%% save 5X5 gridboxes
save(strcat('data/cps/global_regts_',period,expe),'reg_ts', 'glat', 'glon'); %print out 5x5 gridboxes 
awghts=cos(deg2rad(glat));  %% area-weighting for averaging

%%% regrid to 15X15 globally between and 30 and 90
clon1=-180;
clon2=180;
clat1=30;
clat2=90;

igood=find(glon > clon2 | glon < clon1 | ...
     glat < clat1 | glat > clat2);
saveregts=reg_ts(:,igood,:);
saveglat=glat(igood);
saveglon=glon(igood);
ngd=numel(igood)

ngd2=0;
for ila=clat1:15:clat2-15
for ilon=clon1:15:clon2-15
ngd2=ngd2+1;

for ip=1:19
for iyr=1:1996
   ibox=find(glat >= ila & glat < ila+15 & glon>= ilon & glon < ilon+15);
   tmprts=reg_ts(iyr,ibox,ip);
   tmpwts=awghts(ibox);
   igts=find(isfinite(tmprts));
   if (~isempty(igts))  
      newregts(iyr,ngd2,ip)=sum(tmprts(igts).*tmpwts(igts))./sum(tmpwts(igts));
   else
      newregts(iyr,ngd2,ip)=NaN;
   end
   newglat(ngd2)=ila+7.5;
   newglon(ngd2)=ilon+7.5;
end
end

end
end

reg_ts=[saveregts newregts];
glat=[saveglat newglat];
glon=[saveglon newglon];

%% save final grid gridboxes
save(strcat('data/cps/global_regts_15deg_',period,expe),'reg_ts', 'glat', 'glon'); %% save 15X15 gridboxes

awghts=cos(deg2rad(glat));

%%% calculate hemispheric means
for yr=1:nyears_ptot
for n=1:19  
  tmpts=reg_ts(yr,:,n);
  nhinds=find(glat > 0. & isfinite(tmpts));
  shinds=find(glat < 0. & isfinite(tmpts));
  if (~isempty(nhinds))
    NHmean(yr,n)=nansum(tmpts(nhinds).*awghts(nhinds))./sum(awghts(nhinds));
  else
    NHmean(yr,n)=NaN;
  end
  if (~isempty(shinds))
    SHmean(yr,n)=nansum(tmpts(shinds).*awghts(shinds))./sum(awghts(shinds));
  else
    SHmean(yr,n)=NaN;
  end
end
end

%rescale hemispheric means
for icase=1:4 
%% five cases are 
%% 1) CRU
%% 2) iCRU
%% 3) HAD
%% 4) iHad

 switch icase
 case 1
  
  nhfilledmean=load('data/CRU_NH_reform');
  shfilledmean=load('data/CRU_SH_reform');
  expe2=strcat(expe,'_CRU')

 case 2
 
  nhfilledmean=load('data/iCRU_NH_reform');
  shfilledmean=load('data/iCRU_SH_reform');
  expe2=strcat(expe,'_iCRU')
    
 case 3
 
  nhfilledmean=load('data/HAD_NH_reform');
  shfilledmean=load('data/HAD_SH_reform');
  expe2=strcat(expe,'_HAD')

 case 4
  
  nhfilledmean=load('data/iHAD_NH_reform');
  shfilledmean=load('data/iHAD_SH_reform');
  expe2=strcat(expe,'_iHAD')

 end %% switch icase

nhinst=nhfilledmean(:,2);
shinst=shfilledmean(:,2);

[smoothed, icb, ice, mse0]=lowpassmin(nhinst, cutfreq); 
nhinst=smoothed;
[smoothed, icb, ice, mse0]=lowpassmin(shinst, cutfreq); 
shinst=smoothed;

nhinstm=mean(nhinst(scalib-sinst+1:ecalib-sinst+1));
nhinsts=std(nhinst(scalib-sinst+1:ecalib-sinst+1));
shinstm=mean(shinst(scalib-sinst+1:ecalib-sinst+1));
shinsts=std(shinst(scalib-sinst+1:ecalib-sinst+1));

for n=1:19

nhgrmean=mean(NHmean(scalib-1:ecalib-1,n));
nhgrstd=std(NHmean(scalib-1:ecalib-1,n));
shgrmean=mean(SHmean(scalib-1:ecalib-1,n));
shgrstd=std(SHmean(scalib-1:ecalib-1,n));

NHmean(:,n)=((NHmean(:,n)-nhgrmean)*nhinsts/nhgrstd)+nhinstm;
SHmean(:,n)=((SHmean(:,n)-shgrmean)*shinsts/shgrstd)+shinstm;

end

%splice and save

gbcps=NaN([nyears_ptot,1]);
for n=1:19
  gbcps((n-1)*100+1:n*100)=NHmean((n-1)*100+1:n*100,20-n);
end
gbcps(1901:1996)=NHmean(1901:1996,1);

save(strcat('data/cps/results/NH_',period,expe2,'_gridbox.mat'),'gbcps');

gbcps=NaN([nyears_ptot,1]);
for n=1:19
  gbcps((n-1)*100+1:n*100)=SHmean((n-1)*100+1:n*100,20-n);
end
gbcps(1901:1996)=SHmean(1901:1996,1);

save(strcat('data/cps/results/SH_',period,expe2,'_gridbox.mat'),'gbcps');

%% save unspliced for validation and plotting.

gbcpsfull=NHmean;
save(strcat('data/cps/results/NH_',period,expe2,'_gbeachyear.mat'),'gbcpsfull');
gbcpsfull=SHmean;
save(strcat('data/cps/results/SH_',period,expe2,'_gbeachyear.mat'),'gbcpsfull');

end %% for icase
res=0;
