

% LOAD PARAMETERS FROM REAL CORAL DATA
load('pages2kv2_corals_b.mat')

mon_avg_o=4;
mon_avg_f=3;
%mdl_prior='cesm_lme010';
%mdl_prior='cesm_lme009';
%p_yrs=[850:1848]-850+1; 
%c_yrs=1900:2000;

%mdl_prior='gfdlesm2g_pi';
mdl_prior='gfdlesm2m_pi';
p_yrs=[1:499]; 
c_yrs=1900:2000;


[~,ip] = intersect(1:2012,c_yrs);

%-------------------------
% LINEAR STATISTICAL PSM
%-------------------------

% LOAD OBSERVATIONAL DATA
obs_dtp1='soda_sst';% sea surface temperature
[O_sst,xlat,xlon,x_yrs] = load_obs_S(obs_dtp1,mon_avg_o,mon_avg_f);

obs_dtp1='soda_sss';% sea surface salinity
[O_sss] = load_obs_S(obs_dtp1,mon_avg_o,mon_avg_f);

% Get locations of proxies in obs
mask=ones(size(O_sst(:,1))); % mask out land lat/lon values
mask(isnan(O_sst(:,1)))=NaN;
[pl_oi] = nearest_latlon(xlat,xlon,plat,plon,mask);

% Just pick unique locations
%----- Just take unique locations -------
[pl_o,id_u,id_c]=unique(pl_oi,'stable');

% Average the PSM parameters
b1=zeros(3,length(pl_o));
for i=1:length(pl_o)
   b1(:,i)=mean(b(:,find(id_c==i)),2);
end
b=b1;

pNum=length(pl_o);
plat=plat(id_u);plon=plon(id_u);

% LOAD MODEL DATA
[M_sst,~,~,~,olat,olon]=load_S(mdl_prior,'s',mon_avg_o,mon_avg_f,p_yrs);
[M_sss]=load_S(mdl_prior,'y',mon_avg_o,mon_avg_f,p_yrs);

% Locations of proxies in the model
% mask out land lat/lon values
mask=ones(size(M_sst(:,1)));
mask(isnan(M_sst(:,1)))=NaN;
[pl_m] = nearest_latlon(olat,olon,plat,plon,mask);

% Bias-correct model output at proxy locations only
[M_sstbc] = bias_correct(O_sst(pl_o,:),M_sst(pl_m,:),'n');
[M_sssbc] = bias_correct(O_sss(pl_o,:),M_sss(pl_m,:),'n');

% Free up space
clear M_sst M_sss

% Years of overlap
[~,io] = intersect(x_yrs,c_yrs);


% Bi-variate linear fit with 'observed' SST and SSS and
% transfer those fits to the model data
% b=zeros(3,pNum);
% r=zeros(pNum,length(io));
% for k=1:pNum
%     o_sst=O_sst(pl_o(k),io);
%     o_sss=O_sss(pl_o(k),io);
%     pd=prxs(ip,k);
%     [b(:,k),~,r(k,:)]=regress(pd(:),[ones(size(o_sst(:))) o_sst(:) o_sss(:)]);
% end

%save('pages2kv2_corals_b.mat','b','r','plat','plon')

%-------------------------------------------
% APPLY LINEAR FIT TO MODEL DATA AND OUTPUT
%-------------------------------------------

% ps_crl=zeros(pNum,length(p_yrs));
% for jj=1:pNum
%     ps_crl(jj,:)=b(1,jj) + b(2,jj)*M_sstbc(jj,:) + b(3,jj)*M_sssbc(jj,:);
% end

ps_crl=zeros(pNum,length(p_yrs));
for jj=1:pNum
    ps_crl(jj,:)=zscore(b(1,jj) + b(2,jj)*M_sstbc(jj,:) + b(3,jj)*M_sssbc(jj,:));
end

save(['/d1/nsteiger/proxy-data/psm-pseudo-proxies/pseudo_corals_' mdl_prior '_monavg_' num2str(mon_avg_o) '_'...
    num2str(mon_avg_f) '_pyrs_' num2str(p_yrs(1)) '_' ...
    num2str(p_yrs(end))],'ps_crl','plat','plon'); 


%mdl_prior='cesm_lme010';
%mdl_prior='cesm_lme009';
%save(['/d1/nsteiger/proxy-data/psm-pseudo-proxies/pseudo_vsltrees_' mdl_prior '_pyrs_' num2str(850) '_' ...
%    num2str(1849)],'p_d','p_lat','p_lon'); 

% for jj=1:pNum
%     
%     % Obs error estimate (obs error variance)
%     proxy{ii}.R=nanmean(r(jj,:).^2).*ones(1,itr); % MSE = evar + bias^2
%     
%     proxy{ii}.HXb=b(1,jj) + b(2,jj)*M_sstbc(jj,:) + b(3,jj)*M_sssbc(jj,:);
%     proxy{ii}.ptype=archv_tp{jj}; % save type
%     proxy{ii}.b=b(:,jj); % save type
%     
%     % Make each age ensemble the same (will be modified later?)
%     age_yng=repmat(year',[1 itr]);
%     age_old=repmat(year',[1 itr]);
%     % no age ensemble
%     %age_yng=year';
%     %age_old=year';
%     
%     % Put proxies into structure
%     proxy{ii}.age_yng=age_yng;
%     proxy{ii}.age_old=age_old;
%     proxy{ii}.data=repmat(prxs(:,jj)',[itr,1]); % row vector ensemble
%     proxy{ii}.location=[plat(jj) plon(jj)];
%     
%     ii=ii+1;
% end




