format compact;clear all

en_number = 10;

cases = {'b.e21.BSSP245smbb.f09_g17.001';...
         'b.e21.BSSP245smbb.f09_g17.002';...
         'b.e21.BSSP245smbb.f09_g17.003';...
         'b.e21.BSSP245smbb.f09_g17.004';...
         'b.e21.BSSP245smbb.f09_g17.005';...
         'b.e21.BSSP245smbb.f09_g17.006';...
         'b.e21.BSSP245smbb.f09_g17.007';...
         'b.e21.BSSP245smbb.f09_g17.009';...
         'b.e21.BSSP245smbb.f09_g17.010';...
         'b.e21.BSSP245smbb.f09_g17.011';... 
         'b.e21.BSSP245smbb.f09_g17.MCB-050PCT.001';...
         'b.e21.BSSP245smbb.f09_g17.MCB-050PCT.002';...
         'b.e21.BSSP245smbb.f09_g17.MCB-050PCT.003';...
         'b.e21.BSSP245smbb.f09_g17.MCB-050PCT.004';...
         'b.e21.BSSP245smbb.f09_g17.MCB-050PCT.005';...
         'b.e21.BSSP245smbb.f09_g17.MCB-050PCT.006';...
         'b.e21.BSSP245smbb.f09_g17.MCB-050PCT.007';...
         'b.e21.BSSP245smbb.f09_g17.MCB-050PCT.009';...
         'b.e21.BSSP245smbb.f09_g17.MCB-050PCT.010';...
         'b.e21.BSSP245smbb.f09_g17.MCB-050PCT.011'};

vars = {'LANDFRAC';'TREFHT';'PRECT'};

dates = {'201501-210012';...
         '201501-210012';...
         '201501-210012';...
         '201501-210012';...
         '201501-210012';...
         '201501-210012';...
         '201501-210012';...
         '201501-210012';...
         '201501-210012';...
         '201501-210012';...
         '203501-206912';...
         '203501-206912';...
         '203501-206912';...
         '203501-206912';...
         '203501-206912';...
         '203501-206912';...
         '203501-206912';...
         '203501-206912';...
         '203501-206912';...
         '203501-206912'};

data = struct('casename',{},'varname',{},'lon',{},'lat',{},'date',{},'year',{},'month',{},'gw',{},'value',{},'global_avg',{},'annual_avg',{},'annual_avg_year',{});

vars_name = {'land fraction';'(a) Global Mean Temperature';'(b) Inter-hemispheric Temperature Gradient';...
             '(c) Equator-to-pole Temperature Gradient';'(d) Global mean precipitation (mm/day)'};

for i=1:length(cases);
for j=1:length(vars)+2;
    
    if( j<=length(vars) )
     if( i<= 10)
      f = ['../../' char(cases(i)) '/atm/proc/tseries/month_1/' char(cases(i)) '.cam.h0.' char(vars(j)) '.201501-206412.nc']
      f2 = ['../../' char(cases(i)) '/atm/proc/tseries/month_1/' char(cases(i)) '.cam.h0.' char(vars(j)) '.206501-210012.nc'];
      date1 = ncread(f,'date');
      date2 = ncread(f2,'date');
      data(i,j).date = [ncread(f,'date');ncread(f2,'date')];
      if(j==length(vars)+2)
      data(i,j).value = ncread(f,char(vars(j)))*1e2*86400;
      data(i,j).value(:,:,length(date1)+1:length(date1)+length(date2)) = ncread(f2,char(vars(j)))*1e2*86400;
      else
      data(i,j).value = ncread(f,char(vars(j)));
      data(i,j).value(:,:,length(date1)+1:length(date1)+length(date2)) = ncread(f2,char(vars(j)));
      end 
      clear date1 date2
     else
      f = ['../' char(cases(i)) '/atm/proc/tseries/month_1/' char(cases(i)) '.cam.h0.' char(vars(j)) '.' char(dates(i)) '.nc']
      data(i,j).date = ncread(f,'date');
      if(j==length(vars)+2)
      data(i,j).value = ncread(f,char(vars(j)))*1e2*86400;
      else
      data(i,j).value = ncread(f,char(vars(j)));
      end
     end

    data(i,j).varname = char(vars(j));
    data(i,j).lon = ncread(f,'lon');
    data(i,j).lat = ncread(f,'lat');
    data(i,j).gw = ncread(f,'gw');

    else
      data(i,j).casename = data(i,j-1).casename;
      data(i,j).varname = char(vars_name(j));
      data(i,j).lon = data(i,1).lon;
      data(i,j).lat = data(i,1).lat;
      data(i,j).date = data(i,1).date;
      data(i,j).gw = data(i,1).gw;   
      data(i,j).value = data(i,2).value;
      if( j==length(vars) )
        for jj=1:length(data(1,1).lat)
		 data(i,j).value(:,jj,:) = data(i,j).value(:,jj,:).*sin(data(1,1).lat(jj)/180*pi);
        end
      end
      if( j==length(vars)+1 ) 
        for jj=1:length(data(1,1).lat)
		 data(i,j).value(:,jj,:) = data(i,j).value(:,jj,:).*0.5.*(3*sin(data(1,1).lat(jj)/180*pi).*sin(data(1,1).lat(jj)/180*pi)-1);
        end 
      end
    end
    [data(i,j).year data(i,j).month] = set_yr_mon(data(i,j).date);

    for n=1:length(data(i,j).date)
        data(i,j).global_avg(:,n) =  global_avg(squeeze(data(i,j).value(:,:,n)),data(i,j).gw,squeeze(data(i,1).value(:,:,n)));
    end 
    
    for m=1:3
       [data(i,j).annual_avg(m,:) data(i,j).annual_avg_year] = annual_avg_1d(squeeze(data(i,j).global_avg(m,:)),data(i,j).year,data(i,j).month);
    end
end
end

panels = [.1 .76  .7 .17;...
	  .1 .54 .7 .17;...
	  .1 .32 .7 .17;...
	  .1 .1 .7 .17];


figure(1);orient tall;orient landscape
for j=2:length(vars)+2;
subplot('position',panels(j-1,:))

for i=1:en_number
hold on
a(i) = avg(data(i,j).annual_avg(1,:),data(i,j).annual_avg_year,2020,2039);
for year = 2020:2069
aa(i,year-2020+1) = avg(data(i,j).annual_avg(1,:),data(i,j).annual_avg_year,year,year);   
end
end

for i=en_number+1:length(cases)
hold on
b(i-en_number) = avg(data(i,j).annual_avg(1,:),data(i,j).annual_avg_year,2050,2069);
for year = 2035:2069
bb(i-en_number,year-2035+1) = avg(data(i,j).annual_avg(1,:),data(i,j).annual_avg_year,year,year);   
end
end

yy1 = mean(aa,1)-2*std(aa,1);
yy2 = mean(aa,1)+2*std(aa,1);
hhh = patch([[2020:1:2069] fliplr([2020:1:2069])],[yy1 fliplr(yy2)],[0.6 0.6 0.6]);
hold on;clear yy1 yy2 hhh
yy1 = mean(bb,1)-2*std(bb,1);
yy2 = mean(bb,1)+2*std(bb,1);
hhh = patch([[2035:1:2069] fliplr([2035:1:2069])],[yy1 fliplr(yy2)],[0.6 0 0]);
hold on;clear yy1 yy2 hhh
alpha(0.7)

l1=plot([2020:1:2069],mean(aa,1),'k-');
l2=plot([2035:1:2069],mean(bb,1),'r-');
set(l1,'linewidth',1)
set(l2,'linewidth',1)

set(title(char(vars_name(j))),'fontsize',10)
set(gca,'xlim',[2020 2070])
set(plot([2020 2039],[mean(a) mean(a)],'k--'),'linewidt',1)
set(plot([2050 2069],[mean(b) mean(b)],'r--'),'linewidt',1)


%clear a aa bb 
box on
set(gca,'fontsize',8)
if(j<length(vars)+2)
ylabel('Temperature (K)')
end

legend([l1 l2],{'CESM SSP2.4-5','MCB 5%'},'Location','Northwest')
legend('boxoff')

end