clear expr ph

%expr(1)=struct('name','N1850C5OL45OCL32_18mar2016_f19_tn11', ...
%               'name_disp','N1850C5OL45OCL32\_18mar2016\_f19\_tn11', ...
%               'first_year',1, ...
%               'last_year',200, ...
%               'display_first_year',1, ...
%               'add_to_expr',[], ...
%               'color',cmu.colors('forest green (web)'), ...
%               'print',1, ...
%               'grid_file','/hexagon/work/shared/noresm/inputdata/ocn/micom/tnx1v1/20120120/grid.nc', ...
%               'heat_toa_path','../atm', ...
%               'monthly_diagnostics',1,...
%               'path','/fimm/work/matsbn/hexagon-work/oyvinds/archive/N1850C5OL45OCL32_18mar2016_f19_tn11/ocn/hist');
%expr(2)=struct('name','N1850C5OL45OCL32_18mar2016_f09_tn11', ...
%               'name_disp','N1850C5OL45OCL32\_18mar2016\_f09\_tn11', ...
%               'first_year',1, ...
%               'last_year',200, ...
%               'display_first_year',1, ...
%               'add_to_expr',[], ...
%               'color',cmu.colors('red'), ...
%               'print',0, ...
%               'grid_file','/hexagon/work/shared/noresm/inputdata/ocn/micom/tnx1v1/20120120/grid.nc', ...
%               'heat_toa_path','../atm', ...
%               'monthly_diagnostics',1,...
%               'path','/fimm/work/matsbn/hexagon-work/oyvinds/archive/N1850C5OL45OCL32_18mar2016_f09_tn11/ocn/hist');
%expr(3)=struct('name','N1850C5OL45OCL32_f19_tn11_EnergyFix_01', ...
%               'name_disp','N1850C5OL45OCL32\_f19\_tn11\_EnergyFix\_01', ...
%               'first_year',1, ...
%               'last_year',200, ...
%               'display_first_year',1, ...
%               'add_to_expr',[], ...
%               'color',cmu.colors('black'), ...
%               'print',0, ...
%               'grid_file','/hexagon/work/shared/noresm/inputdata/ocn/micom/tnx1v1/20120120/grid.nc', ...
%               'heat_toa_path','../atm', ...
%               'monthly_diagnostics',1,...
%               'path','/fimm/work/matsbn/hexagon-work/matsbn/archive/N1850C5OL45OCL32_f19_tn11_EnergyFix_01/ocn/hist');
%expr(4)=struct('name','N1850C5OL45OCL32_21mar2016_f19_tn11', ...
%               'name_disp','N1850C5OL45OCL32\_21mar2016\_f19\_tn11', ...
%               'first_year',1, ...
%               'last_year',200, ...
%               'display_first_year',1, ...
%               'add_to_expr',[], ...
%               'color',cmu.colors('magenta'), ...
%               'print',0, ...
%               'grid_file','/hexagon/work/shared/noresm/inputdata/ocn/micom/tnx1v1/20120120/grid.nc', ...
%               'heat_toa_path','../atm', ...
%               'monthly_diagnostics',1,...
%               'path','/fimm/work/matsbn/hexagon-work/alfgr/archive/N1850C5OL45OCL32_21mar2016_f19_tn11/ocn/hist');
%expr(1)=struct('name','N1850C5OL45OCL32_30mar2016_f19_tn11', ...
%               'name_disp','N1850C5OL45OCL32\_30mar2016\_f19\_tn11', ...
%               'first_year',1, ...
%               'last_year',200, ...
%               'display_first_year',1, ...
%               'add_to_expr',[], ...
%               'color',cmu.colors('forest green (web)'), ...
%               'print',1, ...
%               'grid_file','/hexagon/work/shared/noresm/inputdata/ocn/micom/tnx1v1/20120120/grid.nc', ...
%               'heat_toa_path','../atm', ...
%               'monthly_diagnostics',1,...
%               'path','/fimm/work/matsbn/hexagon-work/matsbn/archive/N1850C5OL45OCL32_30mar2016_f19_tn11/ocn/hist');
%%              'path','/fimm/work/matsbn/hexagon-work/matsbn/noresm/N1850C5OL45OCL32_30mar2016_f19_tn11/run');
%expr(2)=struct('name','N1850C5OL45OCL32_01apr2016_f09_tn11', ...
%               'name_disp','N1850C5OL45OCL32\_01apr2016\_f09\_tn11', ...
%               'first_year',1, ...
%               'last_year',200, ...
%               'display_first_year',1, ...
%               'add_to_expr',[], ...
%               'color',cmu.colors('red'), ...
%               'print',0, ...
%               'grid_file','/hexagon/work/shared/noresm/inputdata/ocn/micom/tnx1v1/20120120/grid.nc', ...
%               'heat_toa_path','../atm', ...
%               'monthly_diagnostics',1,...
%               'path','/fimm/work/matsbn/hexagon-work/oyvinds/archive/N1850C5OL45OCL32_01apr2016_f09_tn11/ocn/hist');
%%              'path','/fimm/work/matsbn/hexagon-work/oyvinds/noresm/N1850C5OL45OCL32_01apr2016_f09_tn11/run');
%expr(3)=struct('name','N1850C5OL45L32_f09_tn0251_T02', ...
%               'name_disp','N1850C5OL45L32\_f09\_tn0251\_T02', ...
%               'first_year',1, ...
%               'last_year',200, ...
%               'display_first_year',1, ...
%               'add_to_expr',[], ...
%               'color',cmu.colors('blue'), ...
%               'print',0, ...
%               'grid_file','/shared/noresm/inputdata/ocn/micom/tnx0.25v1/20130930/grid.nc', ...
%               'heat_toa_path','../atm', ...
%               'monthly_diagnostics',1,...
%               'path','/fimm/work/matsbn/hexagon-work/agu002/archive/N1850C5OL45L32_f09_tn0251_T02/ocn/hist');
%%              'path','/fimm/work/matsbn/hexagon-work/agu002/noresm/N1850C5OL45L32_f09_tn0251_T02/run');
%%%
expr(1)=struct('name','NAER1850CNOC_f19_g16_03', ...
               'name_disp','CMIP5', ...
               'first_year',1, ...
               'last_year',100, ...
               'display_first_year',1, ...
               'add_to_expr',[], ...
               'color',cmu.colors('black'), ...
               'print',0, ...
               'grid_file','/fimm/work/matsbn/hexagon-work-common/shared/noresm/inputdata/ocn/micom/gx1v6/20100629/grid.nc', ...
               'heat_toa_path','../atm', ...
               'monthly_diagnostics',1,...
               'path','/work/matsbn/archive/NAER1850CNOC_f19_g16_03/ocn/hist');
expr(2)=struct('name','N1850_f09_tn11_NorESM_15082017', ...
               'name_disp','N1850\_f09\_tn11\_NorESM\_15082017', ...
               'first_year',1, ...
               'last_year',100, ...
               'display_first_year',1, ...
               'add_to_expr',[], ...
               'color',cmu.colors('blue'), ...
               'print',1, ...
               'grid_file','/fimm/work/matsbn/hexagon-work-common/shared/noresm/inputdata/ocn/micom/tnx1v1/20120120/grid.nc', ...
               'heat_toa_path','../atm', ...
               'monthly_diagnostics',1,...
               'path','/hexagon/work/agu002/archive/N1850_f09_tn11_NorESM_15082017/ocn/hist');
expr(3)=struct('name','N1850_f09_tn11_NorESM_NorESM_15082017_20years_AM', ...
               'name_disp','N1850\_f09\_tn11\_NorESM\_NorESM\_15082017\_20years\_AM', ...
               'first_year',1, ...
               'last_year',100, ...
               'display_first_year',1, ...
               'add_to_expr',[], ...
               'color',cmu.colors('red'), ...
               'print',0, ...
               'grid_file','/fimm/work/matsbn/hexagon-work-common/shared/noresm/inputdata/ocn/micom/tnx1v1/20120120/grid.nc', ...
               'heat_toa_path','../atm', ...
               'monthly_diagnostics',1,...
               'path','/fimm/work/matsbn/vilje-work/agu002/archive/N1850_f09_tn11_NorESM_NorESM_15082017_20years_AM/ocn/hist');
expr(4)=struct('name','N1850_f09_tn14_NorESM_15082017_20years', ...
               'name_disp','N1850\_f09\_tn14\_NorESM\_15082017\_20years', ...
               'first_year',1, ...
               'last_year',100, ...
               'display_first_year',1, ...
               'add_to_expr',[], ...
               'color',cmu.colors('forest green (web)'), ...
               'print',0, ...
               'grid_file','/fimm/work/matsbn/hexagon-work-common/shared/noresm/inputdata/ocn/micom/tnx1v4/20170601/grid.nc', ...
               'heat_toa_path','../atm', ...
               'monthly_diagnostics',1,...
               'path','/fimm/work/matsbn/vilje-work/agu002/archive/N1850_f09_tn14_NorESM_15082017_20years/ocn/hist');
expr(5)=struct('name','N1850_f19_tn11_230815', ...
               'name_disp','N1850\_f19\_tn11\_230815', ...
               'first_year',1, ...
               'last_year',100, ...
               'display_first_year',1, ...
               'add_to_expr',[], ...
               'color',cmu.colors('pumpkin'), ...
               'print',0, ...
               'grid_file','/fimm/work/matsbn/hexagon-work-common/shared/noresm/inputdata/ocn/micom/tnx1v1/20120120/grid.nc', ...
               'heat_toa_path','../atm', ...
               'monthly_diagnostics',1,...
               'path','/fimm/work/matsbn/vilje-work/oyvinds/archive/N1850_f19_tn11_230815/ocn/hist');
%expr(5)=struct('name','N1850OC_f09_tn11_NorESM_NorESM_15082017_20years_AM', ...
%               'name_disp','N1850OC\_f09\_tn11\_NorESM\_NorESM\_15082017\_20years\_AM', ...
%               'first_year',1, ...
%               'last_year',50, ...
%               'display_first_year',1, ...
%               'add_to_expr',[], ...
%               'color',cmu.colors('pumpkin'), ...
%               'print',0, ...
%               'grid_file','/fimm/work/matsbn/hexagon-work-common/shared/noresm/inputdata/ocn/micom/tnx1v1/20120120/grid.nc', ...
%               'heat_toa_path','../atm', ...
%               'monthly_diagnostics',1,...
%               'path','/fimm/work/matsbn/vilje-work/agu002/archive/N1850OC_f09_tn11_NorESM_NorESM_15082017_20years_AM/ocn/hist');


mw=ones(1,12)/12;
mw=[31 28 31 30 31 30 31 31 30 31 30 31]/365;
g=9.806;
rearth=6.37122e6;
cp=3990;
picpath='/fimm/home/nersc/matsbn/NorClim/diag/pic/';
fontsize=12;
linewidth=2.0;

legend_list=[];

for nexp=1:length(expr)

  % Load stored time series data if available
  if exist([expr(nexp).heat_toa_path '/heat_toa_ts_' expr(nexp).name '.mat'])
    load([expr(nexp).heat_toa_path '/heat_toa_ts_' expr(nexp).name '.mat'])
    expr(nexp).year_list_toa=year_list;
    expr(nexp).FSNT=FSNT;
    expr(nexp).FLNT=FLNT;

    ncid=netcdf.open(expr(nexp).grid_file,'NC_NOWRITE');
    varid=netcdf.inqVarID(ncid,'parea');
    area=netcdf.getVar(ncid,varid);
    varid=netcdf.inqVarID(ncid,'pdepth');
    depth=netcdf.getVar(ncid,varid);
    nreg=netcdf.getAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'nreg');
    netcdf.close(ncid)

    if nreg==2
      area(:,end)=nan;
      depth(:,end)=nan;
    end

    expr(nexp).heat_to_dt=365*86400*nansum(area(:))/(3990*1024*nansum(area(:).*depth(:)));
  else
    expr(nexp).year_list_toa=[];
    expr(nexp).FSNT=[];
    expr(nexp).FLNT=[];
  end
  if exist(['temp_ts_' expr(nexp).name '.mat'])
    load(['temp_ts_' expr(nexp).name '.mat'])
  else
    year_list=[];
    temp=[];
    mass=[];
  end

  % If stored time series data is continous in time and cover the
  % requested time period, skip reading data from disk
  if isempty(year_list)||any(diff(year_list)~=1)|| ...
     year_list(1)>expr(nexp).first_year||year_list(end)<expr(nexp).last_year

    % Try to complement the time series by reading and transforming data
    % from disk
    year_missing=[];
    years_added=0;
    grid_info_read=0;
    for year=expr(nexp).first_year:expr(nexp).last_year
      cyear=sprintf('%4.4d',year);
      if isempty(find(year_list==year))
        if expr(nexp).monthly_diagnostics
          if exist([expr(nexp).path '/' expr(nexp).name '.micom.hm.' ...
                    cyear '-01.nc'])
            complete_year=1;
            disp([expr(nexp).name ' ' cyear])
            for month=1:12
              cmonth=sprintf('%2.2d',month);
              fname=[expr(nexp).path '/' expr(nexp).name '.micom.hm.' ...
                     cyear '-' cmonth '.nc'];
              if exist(fname)

                if ~grid_info_read
                  ncid=netcdf.open(expr(nexp).grid_file,'NC_NOWRITE');
                  varid=netcdf.inqVarID(ncid,'parea');
                  area=netcdf.getVar(ncid,varid);
                  varid=netcdf.inqVarID(ncid,'pmask');
                  mask=netcdf.getVar(ncid,varid);
                  varid=netcdf.inqVarID(ncid,'plat');
                  lat=netcdf.getVar(ncid,varid);
                  varid=netcdf.inqVarID(ncid,'plon');
                  lon=netcdf.getVar(ncid,varid);
                  nreg=netcdf.getAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'nreg');
                  netcdf.close(ncid)

                  area(find(mask==0))=nan;
                  area_sum=nansum(area(:));

                  grid_info_read=1;
                end

                ncid=netcdf.open(fname,'NC_NOWRITE');
                varid=netcdf.inqVarID(ncid,'temp');
                tmp=netcdf.getVar(ncid,varid);
                scale_factor=netcdf.getAtt(ncid,varid,'scale_factor');
                add_offset=netcdf.getAtt(ncid,varid,'add_offset');
                fill_value=netcdf.getAtt(ncid,varid,'_FillValue');
                temp3d=ones(size(tmp))*nan;
                ind=find(tmp~=fill_value);
                temp3d(ind)=double(tmp(ind))*scale_factor+add_offset;
                if nreg==2
                  temp3d(:,end,:)=nan;
                end
                varid=netcdf.inqVarID(ncid,'dp');
                tmp=netcdf.getVar(ncid,varid);
                scale_factor=netcdf.getAtt(ncid,varid,'scale_factor');
                add_offset=netcdf.getAtt(ncid,varid,'add_offset');
                fill_value=netcdf.getAtt(ncid,varid,'_FillValue');
                dp3d=ones(size(tmp))*nan;
                ind=find(tmp~=fill_value);
                dp3d(ind)=double(tmp(ind))*scale_factor+add_offset;
                if nreg==2
                  dp3d(:,end,:)=nan;
                end
                netcdf.close(ncid)

                temp_tmp(month)=nansum(reshape(nansum(temp3d.*dp3d,3).*area,1,[]))/nansum(reshape(nansum(dp3d,3).*area,1,[]));
                mass_tmp(month)=nansum(reshape(nansum(dp3d,3).*area,1,[]))/g;

              else
                complete_year=0;
                break
              end
            end
            if complete_year
              years_added=1;
              year_list(end+1)=year;
              temp=[temp temp_tmp];
              mass=[mass mass_tmp];
            else
              year_missing=[year_missing year];
            end
          else
            year_missing=[year_missing year];
          end
        else
          fname=[expr(nexp).path '/' expr(nexp).name '.micom.hy.' ...
                 cyear '.nc'];
          if exist(fname)
            disp([expr(nexp).name ' ' cyear])

            if ~grid_info_read
              ncid=netcdf.open(expr(nexp).grid_file,'NC_NOWRITE');
              varid=netcdf.inqVarID(ncid,'parea');
              area=netcdf.getVar(ncid,varid);
              varid=netcdf.inqVarID(ncid,'pmask');
              mask=netcdf.getVar(ncid,varid);
              varid=netcdf.inqVarID(ncid,'plat');
              lat=netcdf.getVar(ncid,varid);
              varid=netcdf.inqVarID(ncid,'plon');
              lon=netcdf.getVar(ncid,varid);
              nreg=netcdf.getAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'nreg');
              netcdf.close(ncid)

              area(find(mask==0))=nan;
              area_sum=nansum(area(:));

              grid_info_read=1;
            end

            ncid=netcdf.open(fname,'NC_NOWRITE');
            varid=netcdf.inqVarID(ncid,'temp');
            tmp=netcdf.getVar(ncid,varid);
            scale_factor=netcdf.getAtt(ncid,varid,'scale_factor');
            add_offset=netcdf.getAtt(ncid,varid,'add_offset');
            fill_value=netcdf.getAtt(ncid,varid,'_FillValue');
            temp3d=ones(size(tmp))*nan;
            ind=find(tmp~=fill_value);
            temp3d(ind)=double(tmp(ind))*scale_factor+add_offset;
            if nreg==2
              temp3d(:,end,:)=nan;
            end
            varid=netcdf.inqVarID(ncid,'dp');
            tmp=netcdf.getVar(ncid,varid);
            scale_factor=netcdf.getAtt(ncid,varid,'scale_factor');
            add_offset=netcdf.getAtt(ncid,varid,'add_offset');
            fill_value=netcdf.getAtt(ncid,varid,'_FillValue');
            dp3d=ones(size(tmp))*nan;
            ind=find(tmp~=fill_value);
            dp3d(ind)=double(tmp(ind))*scale_factor+add_offset;
            if nreg==2
              dp3d(:,end,:)=nan;
            end
            netcdf.close(ncid)

            temp_tmp=nansum(reshape(nansum(temp3d.*dp3d,3).*area,1,[]))/nansum(reshape(nansum(dp3d,3).*area,1,[]));
            mass_tmp=nansum(reshape(nansum(dp3d,3).*area,1,[]))/g;

            temp_tmp=temp_tmp*ones(1,12);
            mass_tmp=mass_tmp*ones(1,12);

            years_added=1;
            year_list(end+1)=year;
            temp=[temp temp_tmp];
            mass=[mass mass_tmp];
          else
            year_missing=[year_missing year];
          end
        end
      end
    end

    if ~isempty(year_missing)
      disp([expr(nexp).name ': Could not find data for the years: ' ...
           num2str(year_missing)])
    end

    % If the time series has been extended, sort the series and store
    % the updated time series
    if years_added
      [year_list svind]=sort(year_list);
      smind=(ones(12,1)*(svind-1)*12)+(1:12)'*ones(1,length(svind));
      temp=reshape(temp(smind),1,[]);
      mass=reshape(mass(smind),1,[]);
      save(['temp_ts_' expr(nexp).name '.mat'], ...
           'year_list','temp','mass')
    end

  end

  expr(nexp).year_list=year_list;
  expr(nexp).temp=temp;
  expr(nexp).mass=mass;
  if ~isempty(expr(nexp).name_disp)&&isempty(expr(nexp).add_to_expr)
    legend_list(end+1)=nexp;
  end

end

% Combine time series if requested
for nexp=1:length(expr)
  nexpa=expr(nexp).add_to_expr;
  if ~isempty(nexpa) && ~isempty(expr(nexpa).year_list)
    year_list=expr(nexpa).year_list;
    temp=expr(nexpa).temp;
    mass=expr(nexpa).mass;
    ind=find(year_list>expr(nexpa).last_year,1,'first');
    if isempty(ind)
      year_list=[year_list expr(nexp).year_list];
      temp=[temp expr(nexp).temp];
      mass=[mass expr(nexp).mass];
    else
      year_list=[year_list(1:ind-1) expr(nexp).year_list];
      temp=[temp(1:(ind-1)*12) expr(nexp).temp];
      mass=[mass(1:(ind-1)*12) expr(nexp).mass];
    end
    [year_list svind]=sort(year_list);
    smind=(ones(12,1)*(svind-1)*12)+(1:12)'*ones(1,length(svind));
    temp=reshape(temp(smind),1,[]);
    mass=reshape(mass(smind),1,[]);
    expr(nexpa).year_list=year_list;
    expr(nexpa).temp=temp;
    expr(nexpa).mass=mass;
    year_list_toa=expr(nexpa).year_list_toa;
    FSNT=expr(nexpa).FSNT;
    FLNT=expr(nexpa).FLNT;
    ind=find(year_list_toa>expr(nexpa).last_year,1,'first');
    if isempty(ind)
      year_list_toa=[year_list_toa expr(nexp).year_list_toa];
      FSNT=[FSNT expr(nexp).FSNT];
      FLNT=[FLNT expr(nexp).FLNT];
    else
      year_list_toa=[year_list_toa(1:ind-1) expr(nexp).year_list_toa];
      FSNT=[FSNT(1:(ind-1)*12) expr(nexp).FSNT];
      FLNT=[FLNT(1:(ind-1)*12) expr(nexp).FLNT];
    end
    [year_list_toa svind]=sort(year_list_toa);
    smind=(ones(12,1)*(svind-1)*12)+(1:12)'*ones(1,length(svind));
    FSNT=reshape(FSNT(smind),1,[]);
    FLNT=reshape(FLNT(smind),1,[]);
    expr(nexpa).year_list_toa=year_list_toa;
    expr(nexpa).FSNT=FSNT;
    expr(nexpa).FLNT=FLNT;
    expr(nexpa).first_year=min(expr(nexpa).first_year,expr(nexp).first_year);
    expr(nexpa).last_year=max(expr(nexpa).last_year,expr(nexp).last_year);
  end
end


% Plot the time series

%figure('visible','off')
clf
hold on

temp_rec_shown=0;
for nexp=1:length(expr)
  if isempty(expr(nexp).add_to_expr)
    year_list=expr(nexp).year_list;
    temp=expr(nexp).temp;
    mass=expr(nexp).mass;

%   temp_yr=mon2ann(temp);
    temp_yr=mw*reshape(temp,12,[]);
    mass_yr=mw*reshape(mass,12,[]);

    iyf=find(year_list>=expr(nexp).first_year,1);
    iyl=find(year_list<=expr(nexp).last_year,1,'last');
    year_list=year_list(iyf:iyl);
    temp_yr=temp_yr(iyf:iyl);
    year_discont=[diff(year_list)~=1 1];
    iyl=0;
    while iyl<length(year_list)
      iyf=iyl+1;
      iyl=iyl+find(year_discont(iyf:end)==1,1);
      time=expr(nexp).display_first_year-expr(nexp).first_year+ ...
           year_list(iyf:iyl);
      ph(nexp)=plot(time,temp_yr(iyf:iyl), ...
                    '-','LineWidth',linewidth,'Color',expr(nexp).color);
    end

    if ~isempty(expr(nexp).year_list_toa)&&all(diff(year_list)==1)&& ...
       expr(nexp).year_list_toa(1)<=expr(nexp).first_year
      year_list=expr(nexp).year_list_toa;
      FSNT=expr(nexp).FSNT;
      FLNT=expr(nexp).FLNT;

%     FNET_yr=mon2ann(FSNT-FLNT);
      FNET_yr=mw*reshape(FSNT-FLNT,12,[]);
      iyf=find(year_list>=expr(nexp).first_year,1);
      iyl=find(year_list<=expr(nexp).last_year,1,'last');
      year_list=year_list(iyf:iyl);
      FNET_yr=FNET_yr(iyf:iyl);
      temp_rec_yr=zeros(size(year_list));
      temp_rec_yr=temp_yr(1);
      q=365*86400*4*pi*rearth^2/cp;
%     mean(FNET_yr(151:200))
      for n=1:length(year_list)-1
        temp_rec_yr(n+1)=temp_rec_yr(n) ...
                        +0.5*q*(FNET_yr(n  )/mass_yr(n  ) ...
                               +FNET_yr(n+1)/mass_yr(n+1));
      end
      time=expr(nexp).display_first_year-expr(nexp).first_year+year_list;
      plot(time,temp_rec_yr, ...
           '--','LineWidth',linewidth,'Color',expr(nexp).color);
      temp_rec_shown=1;
    end

  end
end


if ~isempty(legend_list)
  h=legend(ph(legend_list),expr(legend_list).name_disp,'location','Best');
  set(h,'FontSize',fontsize)
  legend boxoff
end
grid on
box on
axis tight
set(gca,'FontSize',fontsize) 
xlabel('year','FontSize',fontsize) 
ylabel('\circC','FontSize',fontsize) 
if ~temp_rec_shown
  title('Global mean temperature','FontSize',fontsize) 
else
  title('Global mean temperature (solid line).\newlineTemperature evolution reconstructed from TOA heat balance (dashed line)','FontSize',fontsize) 
end

% Print the time series

for nexp=length(expr):-1:1
  if expr(nexp).print
    if ~exist([picpath '/' expr(nexp).name])
      mkdir([picpath '/' expr(nexp).name])
    end
    print (gcf,'-dpng', '-r300', [picpath '/' expr(nexp).name '/temp_ts.png'])
    break
  end
end
