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%201981%29%28Dec%202002%29RANGE/T/3/runningAverage/T/%28Jun-Aug%29VALUES/dods')

obs=temp;
loaddods('http://iridl.ldeo.columbia.edu/expert/SOURCES/.NOAA/.NCEP/.EMC/.CFS/.MONTHLY/.surface/.temp/X/64.6875/92.8125/RANGEEDGES/Y/1.904735/32.38044/RANGEEDGES%5BX/Y%5Daverage/S/%28Jan%201981%29%28Dec%202002%29RANGE/S/%281%20May%29VALUES/L/1.5/3.5/RANGEEDGES%5BL%5Daverage/dods')

temp=temp-273.15; 
[nt,m]=size(temp)

h=plot(T/12+1960,obs,'ro-');
set(h,'MarkerFaceColor','r')
hold on
h=plot(T/12+1960,mean(temp,2),'ko-');
set(h,'MarkerFaceColor','k')
h=plot(T/12+1960,temp,'+');   
set(h,'color',[1,1,1]*0.5)
plot(T/12+1960,T*0+prctile(obs,85),'--k');   
hold off
legend('obs','CFS')
title('JJA t2m')

%cfs bias
print -dpdf figure1
print -depsc figure1
pause
%-------------------------------------------------------------------

cm = mean(temp(:));
co = mean(obs);

h=plot(T/12+1960,obs-co,'ro-');
set(h,'MarkerFaceColor','r')
hold on
h=plot(T/12+1960,mean(temp,2)-cm,'ko-');
set(h,'MarkerFaceColor','k')
h=plot(T/12+1960,temp-cm,'+');   
set(h,'color',[1,1,1]*0.5)
plot(T/12+1960,T*0+prctile(obs-co,85),'--k');   
hold off
legend('obs','CFS')
title('JJA t2m')

%cfs bias
print -dpdf figure2a
print -depsc figure2a
pause

%-------------------------------------------------------------------
obs_rank = ranks(obs);
obs_prctile = obs_rank/length(obs_rank)*100;

model_rank = ranks(temp(:));
model_prctile = model_rank/length(model_rank)*100;
model_prctile = reshape(model_prctile,nt,m);

mean_rank = ranks([mean(temp,2);temp(:)]); 
mean_prctile=mean_rank(1:nt)/(nt+m*nt)*100;

h=plot(T/12+1960,obs_prctile,'ro-');
set(h,'MarkerFaceColor','r')
hold on
h=plot(T/12+1960,mean_prctile,'ko-');
set(h,'MarkerFaceColor','k')
h=plot(T/12+1960,model_prctile,'+');   
set(h,'color',[1,1,1]*0.5)
plot(T/12+1960,T*0+85,'--k');   
hold off
legend('obs','CFS')
title('JJA t2m percentiles')

%cfs bias fixed by percentiles
print -dpdf figure2
print -depsc figure2
pause
%-------------------------------------------------------------------

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%201960%29%28Dec%202001%29RANGE/T/3/runningAverage/T/%28Jun-Aug%29VALUES/dods')

obs=temp;

loaddods('http://iridl.ldeo.columbia.edu/expert/IRIONLY/SOURCES/.ECMWF/.DEMETER/.ECMWF/.t2m/S/%28Jan%201960%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')

mod1=t2m-273.15; 

loaddods('http://iridl.ldeo.columbia.edu/expert/IRIONLY/SOURCES/.ECMWF/.DEMETER/.UKMO/.t2m/S/%28Jan%201960%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')

mod2=t2m-273.15; 

h=plot(T/12+1960,obs,'ro-');
set(h,'MarkerFaceColor','r')
hold on

h=plot(T/12+1960,mean(mod1,2),'ko-');
set(h,'MarkerFaceColor',[1,1,1]*0.5)

h=plot(T/12+1960,mean(mod2,2),'bo-');
set(h,'MarkerFaceColor',[135 206 250]/255)

h=plot(T/12+1960,mod1,'+');   
set(h,'color',[1,1,1]*0.5)

h=plot(T/12+1960,mod2,'+');   
set(h,'color',[135 206 250]/255)

plot(T/12+1960,T*0+prctile(obs,85),'--k');   
hold off
legend('obs','model1','model2','location','northwest')
title('JJA t2m')

print -dpdf figure3
print -depsc figure3
pause
%----------------------------------------------------------------------------
mod1a = mod1-mean([mod1(:);mod2(:)]);
mod2a = mod2-mean([mod2(:):mod2(:)]);
obsa = obs-mean(obs(:));

h=plot(T/12+1960,obsa,'ro-');
set(h,'MarkerFaceColor','r')
hold on

h=plot(T/12+1960,mean(mod1a,2),'ko-');
set(h,'MarkerFaceColor',[1,1,1]*0.5)

h=plot(T/12+1960,mean(mod2a,2),'bo-');
set(h,'MarkerFaceColor',[135 206 250]/255)

h=plot(T/12+1960,mod1a,'+');   
set(h,'color',[1,1,1]*0.5)

h=plot(T/12+1960,mod2a,'+');   
set(h,'color',[135 206 250]/255)

plot(T/12+1960,T*0+prctile(obsa,85),'--k');   
hold off
legend('obs','model1','model2','location','northwest')
title('JJA t2m')

print -dpdf figure3a
print -depsc figure3a
pause
%----------------------------------------------------------------------------

mod1a = mod1-mean(mod1(:));
mod2a = mod2-mean(mod2(:));
obsa = obs-mean(obs(:));
h=plot(T/12+1960,obsa,'ro-');
set(h,'MarkerFaceColor','r')
hold on

h=plot(T/12+1960,mean(mod1a,2),'ko-');
set(h,'MarkerFaceColor',[1,1,1]*0.5)

h=plot(T/12+1960,mean(mod2a,2),'bo-');
set(h,'MarkerFaceColor',[135 206 250]/255)

h=plot(T/12+1960,mod1a,'+');   
set(h,'color',[1,1,1]*0.5)

h=plot(T/12+1960,mod2a,'+');   
set(h,'color',[135 206 250]/255)

plot(T/12+1960,T*0+prctile(obsa,85),'--k');   
hold off
legend('obs','model1','model2','location','northwest')
title('JJA t2m')

print -dpdf figure4
print -depsc figure4
pause

%----------------------------------------------------------------------------
mod2a = 2*mod2a;
h=plot(T/12+1960,obsa,'ro-');
set(h,'MarkerFaceColor','r')
hold on

h=plot(T/12+1960,mean(mod1a,2),'ko-');
set(h,'MarkerFaceColor',[1,1,1]*0.5)

h=plot(T/12+1960,mean(mod2a,2),'bo-');
set(h,'MarkerFaceColor',[135 206 250]/255)

h=plot(T/12+1960,mod1a,'+');   
set(h,'color',[1,1,1]*0.5)

h=plot(T/12+1960,mod2a,'+');   
set(h,'color',[135 206 250]/255)

plot(T/12+1960,T*0+prctile(obsa,85),'--k');   
hold off
legend('obs','model1','model2','location','northwest')
title('JJA t2m')

print -dpdf figure5
print -depsc figure5
pause

%----------------------------------------------------------------------------
mod1a = mod1a/std(mod1a(:));
mod2a = mod2a/std(mod2a(:));
obsa = obsa/std(obsa);
h=plot(T/12+1960,obsa,'ro-');
set(h,'MarkerFaceColor','r')
hold on

h=plot(T/12+1960,mean(mod1a,2),'ko-');
set(h,'MarkerFaceColor',[1,1,1]*0.5)

h=plot(T/12+1960,mean(mod2a,2),'bo-');
set(h,'MarkerFaceColor',[135 206 250]/255)

h=plot(T/12+1960,mod1a,'+');   
set(h,'color',[1,1,1]*0.5)

h=plot(T/12+1960,mod2a,'+');   
set(h,'color',[135 206 250]/255)

plot(T/12+1960,T*0+prctile(obsa,85),'--k');   
hold off
legend('obs','model1','model2','location','northwest')
title('JJA t2m')

print -dpdf figure6
print -depsc figure6
pause


return


mc=10000;
m=zeros(1,mc);
wrk = temp(~isnan(temp));
for imc = 1:mc
  ind=randperm(length(wrk));
  m(imc)=mean(wrk(ind(1:15)));
end
ml = prctile(m,2.5);
mu = prctile(m,97.5);


plot(S/12+1960,temp,'+');   
hold on
h=plot(S/12+1960,mean(temp,2),'o');
set(h,'MarkerFaceColor','k')
plot(S/12+1960,S*0+ml,'--')
plot(S/12+1960,S*0+mu,'--')
h=plot(T/12+1960,obs,'ro-');
set(h,'MarkerFaceColor','r')
h=plot(S/12+1960,mean(temp,2)-mean(wrk)+mean(obs),'go-');
set(h,'MarkerFaceColor','g')
hold off

corr2(obs,mean(temp(1:22,:),2))

  