%--------------------------super ensembling prcp

loaddods('http://iridl.ldeo.columbia.edu/expert/SOURCES/.UEA/.CRU/.TS2p1/.monthly/.prcp/X/64.6875/92.8125/RANGEEDGES/Y/1.904735/32.38044/RANGEEDGES%5BX/Y%5Daverage/T/%28Jan%201980%29%28Dec%202001%29RANGE/T/3/runningAverage/T/%28Jun-Aug%29VALUES/dods')

obs=prcp.prcp;

clear model
loaddods('http://iridl.ldeo.columbia.edu/expert/IRIONLY/SOURCES/.ECMWF/.DEMETER/.CERFACS/.tp/S/%28Jan%201980%29%28Dec%202001%29RANGE/X/64.6875/92.8125/RANGEEDGES/Y/1.904735/32.38044/RANGEEDGES%5BX/Y%5Daverage/S/%281%20May%29VALUES/L/1.5/3.5/RANGEEDGES%5BL%5Daverage/dods')

model(:,1)=mean(tp.tp*1000*60*60*24*30,2); 

loaddods('http://iridl.ldeo.columbia.edu/expert/IRIONLY/SOURCES/.ECMWF/.DEMETER/.ECMWF/.tp/S/%28Jan%201980%29%28Dec%202001%29RANGE/X/64.6875/92.8125/RANGEEDGES/Y/1.904735/32.38044/RANGEEDGES%5BX/Y%5Daverage/S/%281%20May%29VALUES/L/1.5/3.5/RANGEEDGES%5BL%5Daverage/dods')

model(:,2)=mean(tp.tp*1000*60*60*24*30,2); 

loaddods('http://iridl.ldeo.columbia.edu/expert/IRIONLY/SOURCES/.ECMWF/.DEMETER/.INGV/.tp/S/%28Jan%201980%29%28Dec%202001%29RANGE/X/64.6875/92.8125/RANGEEDGES/Y/1.904735/32.38044/RANGEEDGES%5BX/Y%5Daverage/S/%281%20May%29VALUES/L/1.5/3.5/RANGEEDGES%5BL%5Daverage/dods')

model(:,3)=mean(tp.tp*1000*60*60*24*30,2); 

loaddods('http://iridl.ldeo.columbia.edu/expert/IRIONLY/SOURCES/.ECMWF/.DEMETER/.LODYC/.tp/S/%28Jan%201980%29%28Dec%202001%29RANGE/X/64.6875/92.8125/RANGEEDGES/Y/1.904735/32.38044/RANGEEDGES%5BX/Y%5Daverage/S/%281%20May%29VALUES/L/1.5/3.5/RANGEEDGES%5BL%5Daverage/dods')

model(:,4)=mean(tp.tp*1000*60*60*24*30,2); 

loaddods('http://iridl.ldeo.columbia.edu/expert/IRIONLY/SOURCES/.ECMWF/.DEMETER/.MeteoFrance/.tp/S/%28Jan%201980%29%28Dec%202001%29RANGE/X/64.6875/92.8125/RANGEEDGES/Y/1.904735/32.38044/RANGEEDGES%5BX/Y%5Daverage/S/%281%20May%29VALUES/L/1.5/3.5/RANGEEDGES%5BL%5Daverage/dods')

model(:,5)=mean(tp.tp*1000*60*60*24*30,2); 

loaddods('http://iridl.ldeo.columbia.edu/expert/IRIONLY/SOURCES/.ECMWF/.DEMETER/.MPI/.tp/S/%28Jan%201980%29%28Dec%202001%29RANGE/X/64.6875/92.8125/RANGEEDGES/Y/1.904735/32.38044/RANGEEDGES%5BX/Y%5Daverage/S/%281%20May%29VALUES/L/1.5/3.5/RANGEEDGES%5BL%5Daverage/dods')

model(:,6)=mean(tp.tp*1000*60*60*24*30,2); 

loaddods('http://iridl.ldeo.columbia.edu/expert/IRIONLY/SOURCES/.ECMWF/.DEMETER/.UKMO/.tp/S/%28Jan%201980%29%28Dec%202001%29RANGE/X/64.6875/92.8125/RANGEEDGES/Y/1.904735/32.38044/RANGEEDGES%5BX/Y%5Daverage/S/%281%20May%29VALUES/L/1.5/3.5/RANGEEDGES%5BL%5Daverage/dods')

model(:,7)=mean(tp.tp*1000*60*60*24*30,2); 



for ii=1:22
  ind = setdiff(1:22,ii);
  b=regress(obs(ind),[model(ind,:),ones(21,1)]);
  y(ii) = [model(ii,:),1]*b;
  obsm = mean(obs(ind)); ya = obs(ind)-obsm;
  xm = mean(model(ind,:)); xa = model(ind,:)-xm(ones(21,1),:);
  b = (ya'*xa)./diag(xa'*xa)';
  ysep(ii) = (model(ii,:)-xm)*b'/7+obsm;
  [u,s,v] = svd(xa','econ');
  b=regress(obs(ind),[v(:,1),ones(21,1)]);
  ypcr(ii) = [(inv(s(1,1))*u(:,1)'*(model(ii,:)'-xm'))',1]*b;
  ymm(ii) = mean(model(ii,:),2)-mean(mean(model(ind,:),2))+obsm;
end

h=plot(T/12+1960,[obs,ymm']);
set(h(1),'color','k','linewidth',2)
set(h(2:end),'color',[0,0,0]+0.5)
set(h(end),'linestyle','--')
set(gca,'ylim',[150,240])
legend('obs','model')
title(['JJA prcp correlation = ',num2str(corr2(obs,ymm'),2)])

print -dpdf super1p
print -depsc super1p
pause

b=regress(obs,[model,ones(22,1)]);
yin = [model,ones(22,1)]*b;

h=plot(T/12+1960,[obs,yin,ymm']);
set(h(1),'color','k','linewidth',2)
set(h(2),'color','r','linewidth',2)
set(h(3:end),'color',[0,0,0]+0.5)
set(gca,'ylim',[150,240])
legend('obs','regression','model')

title(['JJA prcp correlation = ',num2str(corr2(obs,yin),2)])
print -dpdf super2p
print -depsc super2p
pause

h=plot(T/12+1960,[obs,yin,y',ymm']);
set(h(1),'color','k','linewidth',2)
set(h(2),'color','r','linewidth',2)
set(h(3),'color','c','linewidth',2)
set(h(4:end),'color',[0,0,0]+0.5)
set(gca,'ylim',[150,240])
legend('obs','regression','cv regression', 'model')

title(['JJA prcp correlation = ',num2str(corr2(obs,y'),2)])
print -dpdf super3p
print -depsc super3p
pause


h=plot(T/12+1960,[obs,yin,y',ysep',ymm']);
set(h(1),'color','k','linewidth',2)
set(h(2),'color','r','linewidth',2)
set(h(3),'color','c','linewidth',2)
set(h(4),'color','g','linewidth',2)
set(h(5:end),'color',[0,0,0]+0.5)
set(gca,'ylim',[150,240])
legend('obs','regression','cv regression', 'sep regres','model')


title(['JJA prcp correlation = ',num2str(corr2(obs,y'),2)])
print -dpdf super4p
print -depsc super4p
pause


h=plot(T/12+1960,[obs,yin,y',ysep',ypcr',ymm']);
set(h(1),'color','k','linewidth',2)
set(h(2),'color','r','linewidth',2)
set(h(3),'color','c','linewidth',2)
set(h(4),'color','g','linewidth',2)
set(h(5),'color','m','linewidth',2)
set(h(6:end),'color',[0,0,0]+0.5)
set(gca,'ylim',[150,240])
legend('obs','regression','cv regression', 'sep regres','pcr','model')

title(['JJA prcp correlation = ',num2str(corr2(obs,ypcr'),2)])
print -dpdf super5p
print -depsc super5p
pause

%--------------------------super ensembling temp

loaddods('http://iridl.ldeo.columbia.edu/expert/SOURCES/.UEA/.CRU/.TS2p1/.monthly/.mean/.temp/X/64.6875/92.8125/RANGEEDGES/Y/1.904735/32.38044/RANGEEDGES%5BX/Y%5Daverage/T/%28Jan%201980%29%28Dec%202001%29RANGE/T/3/runningAverage/T/%28Jun-Aug%29VALUES/dods')

obs=temp.temp;

clear model
loaddods('http://iridl.ldeo.columbia.edu/expert/IRIONLY/SOURCES/.ECMWF/.DEMETER/.CERFACS/.t2m/S/%28Jan%201980%29%28Dec%202001%29RANGE/X/64.6875/92.8125/RANGEEDGES/Y/1.904735/32.38044/RANGEEDGES%5BX/Y%5Daverage/S/%281%20May%29VALUES/L/1.5/3.5/RANGEEDGES%5BL%5Daverage/dods')

model(:,1)=mean(t2m.t2m-273.15,2); 

loaddods('http://iridl.ldeo.columbia.edu/expert/IRIONLY/SOURCES/.ECMWF/.DEMETER/.ECMWF/.t2m/S/%28Jan%201980%29%28Dec%202001%29RANGE/X/64.6875/92.8125/RANGEEDGES/Y/1.904735/32.38044/RANGEEDGES%5BX/Y%5Daverage/S/%281%20May%29VALUES/L/1.5/3.5/RANGEEDGES%5BL%5Daverage/dods')

model(:,2)=mean(t2m.t2m-273.15,2); 

loaddods('http://iridl.ldeo.columbia.edu/expert/IRIONLY/SOURCES/.ECMWF/.DEMETER/.INGV/.t2m/S/%28Jan%201980%29%28Dec%202001%29RANGE/X/64.6875/92.8125/RANGEEDGES/Y/1.904735/32.38044/RANGEEDGES%5BX/Y%5Daverage/S/%281%20May%29VALUES/L/1.5/3.5/RANGEEDGES%5BL%5Daverage/dods')

model(:,3)=mean(t2m.t2m-273.15,2); 

loaddods('http://iridl.ldeo.columbia.edu/expert/IRIONLY/SOURCES/.ECMWF/.DEMETER/.LODYC/.t2m/S/%28Jan%201980%29%28Dec%202001%29RANGE/X/64.6875/92.8125/RANGEEDGES/Y/1.904735/32.38044/RANGEEDGES%5BX/Y%5Daverage/S/%281%20May%29VALUES/L/1.5/3.5/RANGEEDGES%5BL%5Daverage/dods')

model(:,4)=mean(t2m.t2m-273.15,2); 

loaddods('http://iridl.ldeo.columbia.edu/expert/IRIONLY/SOURCES/.ECMWF/.DEMETER/.MeteoFrance/.t2m/S/%28Jan%201980%29%28Dec%202001%29RANGE/X/64.6875/92.8125/RANGEEDGES/Y/1.904735/32.38044/RANGEEDGES%5BX/Y%5Daverage/S/%281%20May%29VALUES/L/1.5/3.5/RANGEEDGES%5BL%5Daverage/dods')

model(:,5)=mean(t2m.t2m-273.15,2); 

loaddods('http://iridl.ldeo.columbia.edu/expert/IRIONLY/SOURCES/.ECMWF/.DEMETER/.MPI/.t2m/S/%28Jan%201980%29%28Dec%202001%29RANGE/X/64.6875/92.8125/RANGEEDGES/Y/1.904735/32.38044/RANGEEDGES%5BX/Y%5Daverage/S/%281%20May%29VALUES/L/1.5/3.5/RANGEEDGES%5BL%5Daverage/dods')

model(:,6)=mean(t2m.t2m-273.15,2); 

loaddods('http://iridl.ldeo.columbia.edu/expert/IRIONLY/SOURCES/.ECMWF/.DEMETER/.UKMO/.t2m/S/%28Jan%201980%29%28Dec%202001%29RANGE/X/64.6875/92.8125/RANGEEDGES/Y/1.904735/32.38044/RANGEEDGES%5BX/Y%5Daverage/S/%281%20May%29VALUES/L/1.5/3.5/RANGEEDGES%5BL%5Daverage/dods')

model(:,7)=mean(t2m.t2m-273.15,2); 

for ii=1:22
  ind = setdiff(1:22,ii);
  b=regress(obs(ind),[model(ind,:),ones(21,1)]);
  y(ii) = [model(ii,:),1]*b;
  obsm = mean(obs(ind)); ya = obs(ind)-obsm;
  xm = mean(model(ind,:)); xa = model(ind,:)-xm(ones(21,1),:);
  b = (ya'*xa)./diag(xa'*xa)';
  ysep(ii) = (model(ii,:)-xm)*b'/7+obsm;
  [u,s,v] = svd(xa','econ');
  b=regress(obs(ind),[v(:,1),ones(21,1)]);
  ypcr(ii) = [(inv(s(1,1))*u(:,1)'*(model(ii,:)'-xm'))',1]*b;
  ymm(ii) = mean(model(ii,:),2)-mean(mean(model(ind,:),2))+obsm;
end

h=plot(T/12+1960,[obs,ymm']);
set(h(1),'color','k','linewidth',2)
set(h(2:end),'color',[0,0,0]+0.5)
set(h(end),'linestyle','--')
legend('obs','model')
title(['JJA t2m correlation = ',num2str(corr2(obs,ymm'),2)])


print -dpdf super1t
print -depsc super1t
pause

b=regress(obs,[model,ones(22,1)]);
yin = [model,ones(22,1)]*b;

h=plot(T/12+1960,[obs,yin,ymm']);
set(h(1),'color','k','linewidth',2)
set(h(2),'color','r','linewidth',2)
set(h(3:end),'color',[0,0,0]+0.5)
legend('obs','regression','model')

title(['JJA t2m correlation = ',num2str(corr2(obs,yin),2)])

print -dpdf super2t
print -depsc super2t
pause

for ii=1:22
  ind = setdiff(1:22,ii);
  b=regress(obs(ind),[model(ind,:),ones(21,1)]);
  y(ii) = [model(ii,:),1]*b;
  obsm = mean(obs(ind)); ya = obs(ind)-obsm;
  xm = mean(model(ind,:)); xa = model(ind,:)-xm(ones(21,1),:);
  b = (ya'*xa)./diag(xa'*xa)';
  ysep(ii) = (model(ii,:)-xm)*b'/7+obsm;
  [u,s,v] = svd(xa','econ');
  b=regress(obs(ind),[v(:,1:3),ones(21,1)]);
  ypcr(ii) = [(inv(s(1:2,1:2))*u(:,1:2)'*(model(ii,:)'-xm'))',1]*b;
  ymm(ii) = mean(model(ii,:),2)-mean(mean(model(ind,:),2))+obsm;
end

h=plot(T/12+1960,[obs,yin,y',ymm']);
set(h(1),'color','k','linewidth',2)
set(h(2),'color','r','linewidth',2)
set(h(3),'color','c','linewidth',2)
set(h(4:end),'color',[0,0,0]+0.5)
legend('obs','regression','cv regression', 'model')

title(['JJA t2m correlation = ',num2str(corr2(obs,y'),2)])
print -dpdf super3t
print -depsc super3t
pause

h=plot(T/12+1960,[obs,yin,y',ysep',ymm']);
set(h(1),'color','k','linewidth',2)
set(h(2),'color','r','linewidth',2)
set(h(3),'color','c','linewidth',2)
set(h(4),'color','g','linewidth',2)
set(h(5:end),'color',[0,0,0]+0.5)
legend('obs','regression','cv regression', 'sep regres','model')


title(['JJA t2m correlation = ',num2str(corr2(obs,y'),2)])
print -dpdf super4t
print -depsc super4t
pause


h=plot(T/12+1960,[obs,yin,y',ysep',ypcr',ymm']);
set(h(1),'color','k','linewidth',2)
set(h(2),'color','r','linewidth',2)
set(h(3),'color','c','linewidth',2)
set(h(4),'color','g','linewidth',2)
set(h(5),'color','m','linewidth',2)
set(h(6:end),'color',[0,0,0]+0.5)
legend('obs','regression','cv regression', 'sep regres','pcr','model')


title(['JJA t2m correlation = ',num2str(corr2(obs,ypcr'),2)])
print -dpdf super5t
print -depsc super5t
pause


