clc
clear all
close all

cn45_280_path = '/discover/nobackup/jkolassa/test_new/';
cn45_280_exp_name = 'GEOSldas_CN45_local_280';

cn45_transient_path = '/discover/nobackup/jkolassa/test_new/';
cn45_transient_exp_base_name = 'GEOSldas_CN45_local_transient';

cn51_280_path = '/discover/nobackup/jkolassa/new/';
cn51_280_exp_name = 'GEOSldas_CN51_local_280';

cn51_transient_path = '/discover/nobackup/jkolassa/new/';
cn51_transient_exp_base_name = 'GEOSldas_CN51_local_transient';

nr_of_transient_cycles = 1;

transient_co2_start_year = [1851];
transient_co2_stop_year = [1891];

start_year = 1980;
stop_year = 2019; 

leapYears = 1980:4:2022;
monthLength_nly = [31 28 31 30 31 30 31 31 30 31 30 31];
monthLength_ly = [31 29 31 30 31 30 31 31 30 31 30 31];

cn51_read_file = ['/discover/nobackup/jkolassa/new/GEOSldas_CN51_local_280/output/SMAP_EASEv2_M36/cat/ens0000/Y1980/M01/GEOSldas_CN51_local_280.tavg24_1d_lnd_Nt.monthly.198001.nc4'];
lat_vec = double(ncread(cn51_read_file,'lat'));
lon_vec = double(ncread(cn51_read_file,'lon'));

for i = 1:nr_of_transient_cycles

    cyc_str = ['c' num2str(i,'%0i')];
    cn51_transient_exp_name = [cn51_transient_exp_base_name '_' cyc_str];
    cn45_transient_exp_name = [cn45_transient_exp_base_name '_' cyc_str];

    mat_index = 0;
    y_index = 0;

    for y = start_year:stop_year

        y_index = y_index + 1;
        y_str = num2str(y,'%04i');
        y_index_str = num2str(y_index,'%04i');

        % define month length based on leap year status
        if sum(y==leapYears)>0
            monthLength = monthLength_ly;
        else
            monthLength = monthLength_nly;
        end


        % loop over months
        for m = 1:12

            mat_index = mat_index + 1;

            m_str = num2str(m,'%02i');

            % define file names

            cn51_transient_read_file = [cn51_transient_path cn51_transient_exp_name '/output/SMAP_EASEv2_M36/cat/ens0000/Y' y_str ...
                           '/M' m_str '/' cn51_transient_exp_name '.tavg24_1d_lnd_Nt.monthly.' y_str m_str '.nc4'];

            cn45_transient_read_file = [cn45_transient_path cn45_transient_exp_name '/output/SMAP_EASEv2_M36/cat/ens0000/Y' y_str ...
                           '/M' m_str '/' cn45_transient_exp_name '.tavg24_1d_lnd_Nt.monthly.' y_str m_str '.nc4'];


            if (i==1)
               cn51_280_read_file = [cn51_280_path cn51_280_exp_name '/output/SMAP_EASEv2_M36/cat/ens0000/Y' y_str ...
                           '/M' m_str '/' cn51_280_exp_name '.tavg24_1d_lnd_Nt.monthly.' y_str m_str '.nc4'];
               cn45_280_read_file = [cn45_280_path cn45_280_exp_name '/output/SMAP_EASEv2_M36/cat/ens0000/Y' y_str ...
                           '/M' m_str '/' cn45_280_exp_name '.tavg24_1d_lnd_Nt.monthly.' y_str m_str '.nc4'];
            end

            % initialize arrays

            if (mat_index == 1)

                cn51_cngpp_transient_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                cn51_cnnee_transient_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                cn51_cncloss_transient_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                cn51_cntlai_transient_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                cn51_cntotc_transient_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                cn51_et_transient_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                cn51_evpsoil_transient_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                cn51_evptrns_transient_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                cn51_evpintr_transient_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                cn51_evpsbln_transient_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);

                cn45_cngpp_transient_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                cn45_cnnee_transient_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                cn45_cncloss_transient_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                cn45_cntlai_transient_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                cn45_cntotc_transient_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                cn45_et_transient_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                cn45_evpsoil_transient_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                cn45_evptrns_transient_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                cn45_evpintr_transient_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                cn45_evpsbln_transient_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);

                if (i ==1)
                   cn51_cngpp_280_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cn51_cnnee_280_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cn51_cncloss_280_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cn51_cntlai_280_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cn51_cntotc_280_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cn51_et_280_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cn51_evpsoil_280_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cn51_evptrns_280_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cn51_evpintr_280_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cn51_evpsbln_280_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);

                   cn45_cngpp_280_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cn45_cnnee_280_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cn45_cncloss_280_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cn45_cntlai_280_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cn45_cntotc_280_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cn45_et_280_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cn45_evpsoil_280_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cn45_evptrns_280_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cn45_evpintr_280_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cn45_evpsbln_280_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                 end
            end % mat_index

            % read data

            if exist(cn51_transient_read_file)
               cn51_cngpp_transient_mat{i}(:,mat_index) = double(ncread(cn51_transient_read_file,'CNGPP'));
               cn51_cnnee_transient_mat{i}(:,mat_index) = double(ncread(cn51_transient_read_file,'CNNEE'));
               cn51_cncloss_transient_mat{i}(:,mat_index) = double(ncread(cn51_transient_read_file,'CLOSS'));
               cn51_cntlai_transient_mat{i}(:,mat_index) = double(ncread(cn51_transient_read_file,'CNTLAI'));
               cn51_cntotc_transient_mat{i}(:,mat_index) = double(ncread(cn51_transient_read_file,'CNTOTC'));
               cn51_evpsoil_transient_mat{i}(:,mat_index) = double(ncread(cn51_transient_read_file,'EVPSOIL'));
               cn51_evptrns_transient_mat{i}(:,mat_index) = double(ncread(cn51_transient_read_file,'EVPTRNS'));
               cn51_evpintr_transient_mat{i}(:,mat_index) = double(ncread(cn51_transient_read_file,'EVPINTR'));
               cn51_evpsbln_transient_mat{i}(:,mat_index) = double(ncread(cn51_transient_read_file,'EVPSBLN'));
               cn51_et_transient_mat{i}(:,mat_index) = cn51_evpsoil_transient_mat{i}(:,mat_index) + cn51_evptrns_transient_mat{i}(:,mat_index) ...
                                      + cn51_evpintr_transient_mat{i}(:,mat_index) + cn51_evpsbln_transient_mat{i}(:,mat_index);
            else
               cn51_cngpp_transient_mat{i}(:,mat_index) = NaN*ones(size(cn51_cngpp_transient_mat{i},1),1);
               cn51_cnnee_transient_mat{i}(:,mat_index) = NaN*ones(size(cn51_cnnee_transient_mat{i},1),1);
               cn51_cncloss_transient_mat{i}(:,mat_index) = NaN*ones(size(cn51_cncloss_transient_mat{i},1),1);
               cn51_cntlai_transient_mat{i}(:,mat_index) = NaN*ones(size(cn51_cntlai_transient_mat{i},1),1);
               cn51_cntotc_transient_mat{i}(:,mat_index) = NaN*ones(size(cn51_cntotc_transient_mat{i},1),1);
               cn51_evpsoil_transient_mat{i}(:,mat_index) = NaN*ones(size(cn51_evpsoil_transient_mat{i},1),1);
               cn51_evptrns_transient_mat{i}(:,mat_index) = NaN*ones(size(cn51_evptrns_transient_mat{i},1),1);
               cn51_evpintr_transient_mat{i}(:,mat_index) = NaN*ones(size(cn51_evpintr_transient_mat{i},1),1);
               cn51_evpsbln_transient_mat{i}(:,mat_index) = NaN*ones(size(cn51_evpsbln_transient_mat{i},1),1);
               cn51_et_transient_mat{i}(:,mat_index) = NaN*ones(size(cn51_et_transient_mat{i},1),1);
            end

            if exist(cn45_transient_read_file)
               cn45_cngpp_transient_mat{i}(:,mat_index) = double(ncread(cn45_transient_read_file,'CNGPP'));
               cn45_cnnee_transient_mat{i}(:,mat_index) = double(ncread(cn45_transient_read_file,'CNNEE'));
               cn45_cncloss_transient_mat{i}(:,mat_index) = double(ncread(cn45_transient_read_file,'CLOSS'));
               cn45_cntlai_transient_mat{i}(:,mat_index) = double(ncread(cn45_transient_read_file,'CNTLAI'));
               cn45_cntotc_transient_mat{i}(:,mat_index) = double(ncread(cn45_transient_read_file,'CNTOTC'));
               cn45_evpsoil_transient_mat{i}(:,mat_index) = double(ncread(cn45_transient_read_file,'EVPSOIL'));
               cn45_evptrns_transient_mat{i}(:,mat_index) = double(ncread(cn45_transient_read_file,'EVPTRNS'));
               cn45_evpintr_transient_mat{i}(:,mat_index) = double(ncread(cn45_transient_read_file,'EVPINTR'));
               cn45_evpsbln_transient_mat{i}(:,mat_index) = double(ncread(cn45_transient_read_file,'EVPSBLN'));
               cn45_et_transient_mat{i}(:,mat_index) = cn45_evpsoil_transient_mat{i}(:,mat_index) + cn45_evptrns_transient_mat{i}(:,mat_index) ... 
                                      + cn45_evpintr_transient_mat{i}(:,mat_index) + cn45_evpsbln_transient_mat{i}(:,mat_index);
            else
               cn45_cngpp_transient_mat{i}(:,mat_index) = NaN*ones(size(cn45_cngpp_transient_mat{i},1),1);
               cn45_cnnee_transient_mat{i}(:,mat_index) = NaN*ones(size(cn45_cnnee_transient_mat{i},1),1);
               cn45_cncloss_transient_mat{i}(:,mat_index) = NaN*ones(size(cn45_cncloss_transient_mat{i},1),1);
               cn45_cntlai_transient_mat{i}(:,mat_index) = NaN*ones(size(cn45_cntlai_transient_mat{i},1),1);
               cn45_cntotc_transient_mat{i}(:,mat_index) = NaN*ones(size(cn45_cntotc_transient_mat{i},1),1);
               cn45_evpsoil_transient_mat{i}(:,mat_index) = NaN*ones(size(cn45_evpsoil_transient_mat{i},1),1);
               cn45_evptrns_transient_mat{i}(:,mat_index) = NaN*ones(size(cn45_evptrns_transient_mat{i},1),1);
               cn45_evpintr_transient_mat{i}(:,mat_index) = NaN*ones(size(cn45_evpintr_transient_mat{i},1),1);
               cn45_evpsbln_transient_mat{i}(:,mat_index) = NaN*ones(size(cn45_evpsbln_transient_mat{i},1),1);
               cn45_et_transient_mat{i}(:,mat_index) = NaN*ones(size(cn45_et_transient_mat{i},1),1);
            end

            if (i == 1)
               if exist(cn51_280_read_file)
                  cn51_cngpp_280_mat(:,mat_index) = double(ncread(cn51_280_read_file,'CNGPP'));
                  cn51_cnnee_280_mat(:,mat_index) = double(ncread(cn51_280_read_file,'CNNEE'));
                  cn51_cncloss_280_mat(:,mat_index) = double(ncread(cn51_280_read_file,'CLOSS'));
                  cn51_cntlai_280_mat(:,mat_index) = double(ncread(cn51_280_read_file,'CNTLAI'));
                  cn51_cntotc_280_mat(:,mat_index) = double(ncread(cn51_280_read_file,'CNTOTC'));
                  cn51_evpsoil_280_mat(:,mat_index) = double(ncread(cn51_280_read_file,'EVPSOIL'));
                  cn51_evptrns_280_mat(:,mat_index) = double(ncread(cn51_280_read_file,'EVPTRNS'));
                  cn51_evpintr_280_mat(:,mat_index) = double(ncread(cn51_280_read_file,'EVPINTR'));
                  cn51_evpsbln_280_mat(:,mat_index) = double(ncread(cn51_280_read_file,'EVPSBLN'));
                  cn51_et_280_mat(:,mat_index) = cn51_evpsoil_280_mat(:,mat_index) + cn51_evptrns_280_mat(:,mat_index) ...
                                              + cn51_evpintr_280_mat(:,mat_index) + cn51_evpsbln_280_mat(:,mat_index);
                else
                  cn51_cngpp_280_mat(:,mat_index) = NaN*ones(size(cn51_cngpp_280_mat,1),1);
                  cn51_cnnee_280_mat(:,mat_index) = NaN*ones(size(cn51_cnnee_280_mat,1),1);
                  cn51_cncloss_280_mat(:,mat_index) = NaN*ones(size(cn51_cncloss_280_mat,1),1);
                  cn51_cntlai_280_mat(:,mat_index) = NaN*ones(size(cn51_cntlai_280_mat,1),1);
                  cn51_cntotc_280_mat(:,mat_index) = NaN*ones(size(cn51_cntotc_280_mat,1),1);
                  cn51_evpsoil_280_mat(:,mat_index) = NaN*ones(size(cn51_evpsoil_280_mat,1),1);
                  cn51_evptrns_280_mat(:,mat_index) = NaN*ones(size(cn51_evptrns_280_mat,1),1);
                  cn51_evpintr_280_mat(:,mat_index) = NaN*ones(size(cn51_evpintr_280_mat,1),1);
                  cn51_evpsbln_280_mat(:,mat_index) = NaN*ones(size(cn51_evpsbln_280_mat,1),1);
                  cn51_et_280_mat(:,mat_index) = NaN*ones(size(cn51_et_280_mat,1),1);
                end

               if exist(cn45_280_read_file)
                  cn45_cngpp_280_mat(:,mat_index) = double(ncread(cn45_280_read_file,'CNGPP'));
                  cn45_cnnee_280_mat(:,mat_index) = double(ncread(cn45_280_read_file,'CNNEE'));
                  cn45_cncloss_280_mat(:,mat_index) = double(ncread(cn45_280_read_file,'CLOSS'));
                  cn45_cntlai_280_mat(:,mat_index) = double(ncread(cn45_280_read_file,'CNTLAI'));
                  cn45_cntotc_280_mat(:,mat_index) = double(ncread(cn45_280_read_file,'CNTOTC'));
                  cn45_evpsoil_280_mat(:,mat_index) = double(ncread(cn45_280_read_file,'EVPSOIL'));
                  cn45_evptrns_280_mat(:,mat_index) = double(ncread(cn45_280_read_file,'EVPTRNS'));
                  cn45_evpintr_280_mat(:,mat_index) = double(ncread(cn45_280_read_file,'EVPINTR'));
                  cn45_evpsbln_280_mat(:,mat_index) = double(ncread(cn45_280_read_file,'EVPSBLN'));
                  cn45_et_280_mat(:,mat_index) = cn45_evpsoil_280_mat(:,mat_index) + cn45_evptrns_280_mat(:,mat_index) ... 
                                              + cn45_evpintr_280_mat(:,mat_index) + cn45_evpsbln_280_mat(:,mat_index);
                else
                  cn45_cngpp_280_mat(:,mat_index) = NaN*ones(size(cn45_cngpp_280_mat,1),1);
                  cn45_cnnee_280_mat(:,mat_index) = NaN*ones(size(cn45_cnnee_280_mat,1),1);
                  cn45_cncloss_280_mat(:,mat_index) = NaN*ones(size(cn45_cncloss_280_mat,1),1);
                  cn45_cntlai_280_mat(:,mat_index) = NaN*ones(size(cn45_cntlai_280_mat,1),1);
                  cn45_cntotc_280_mat(:,mat_index) = NaN*ones(size(cn45_cntotc_280_mat,1),1);
                  cn45_evpsoil_280_mat(:,mat_index) = NaN*ones(size(cn45_evpsoil_280_mat,1),1);
                  cn45_evptrns_280_mat(:,mat_index) = NaN*ones(size(cn45_evptrns_280_mat,1),1);
                  cn45_evpintr_280_mat(:,mat_index) = NaN*ones(size(cn45_evpintr_280_mat,1),1);
                  cn45_evpsbln_280_mat(:,mat_index) = NaN*ones(size(cn45_evpsbln_280_mat,1),1);
                  cn45_et_280_mat(:,mat_index) = NaN*ones(size(cn45_et_280_mat,1),1);
                end

                datenum_vec(mat_index) = datenum(y,m,1);
                year_vec(mat_index) = y;
                month_vec(mat_index) = m;

            end % i == 1

        end % m
     end % y
end % i


for i = 1:nr_of_transient_cycles

    ann_index = 0;

    for y = start_year:stop_year

        ann_index = ann_index + 1;
        year_index = (year_vec==y);

        tmp = cn51_cngpp_transient_mat{i}(:,year_index);
        mean_global_ann_transient_cngpp_cn51{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_transient_cngpp_cn51{i}(:,ann_index) = nanmean(tmp,2);

        tmp = cn51_cnnee_transient_mat{i}(:,year_index);
        mean_global_ann_transient_cnnee_cn51{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_transient_cnnee_cn51{i}(:,ann_index) = nanmean(tmp,2);

        tmp = cn51_cncloss_transient_mat{i}(:,year_index);
        mean_global_ann_transient_cncloss_cn51{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_transient_cncloss_cn51{i}(:,ann_index) = nanmean(tmp,2);

        tmp = cn51_cntlai_transient_mat{i}(:,year_index);
        mean_global_ann_transient_cntlai_cn51{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_transient_cntlai_cn51{i}(:,ann_index) = nanmean(tmp,2);

        tmp = cn51_cntotc_transient_mat{i}(:,year_index);
        mean_global_ann_transient_cntotc_cn51{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_transient_cntotc_cn51{i}(:,ann_index) = nanmean(tmp,2);

        tmp = cn51_evpsoil_transient_mat{i}(:,year_index);
        mean_global_ann_transient_evpsoil_cn51{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_transient_evpsoil_cn51{i}(:,ann_index) = nanmean(tmp,2);

        tmp = cn51_evptrns_transient_mat{i}(:,year_index);
        mean_global_ann_transient_evptrns_cn51{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_transient_evptrns_cn51{i}(:,ann_index) = nanmean(tmp,2);

        tmp = cn51_evpintr_transient_mat{i}(:,year_index);
        mean_global_ann_transient_evpintr_cn51{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_transient_evpintr_cn51{i}(:,ann_index) = nanmean(tmp,2);

        tmp = cn51_evpsbln_transient_mat{i}(:,year_index);
        mean_global_ann_transient_evpsbln_cn51{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_transient_evpsbln_cn51{i}(:,ann_index) = nanmean(tmp,2);

        tmp = cn51_et_transient_mat{i}(:,year_index);
        mean_global_ann_transient_et_cn51{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_transient_et_cn51{i}(:,ann_index) = nanmean(tmp,2);

        tmp = cn45_cngpp_transient_mat{i}(:,year_index);
        mean_global_ann_transient_cngpp_cn45{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_transient_cngpp_cn45{i}(:,ann_index) = nanmean(tmp,2);

        tmp = cn45_cnnee_transient_mat{i}(:,year_index);
        mean_global_ann_transient_cnnee_cn45{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_transient_cnnee_cn45{i}(:,ann_index) = nanmean(tmp,2);

        tmp = cn45_cncloss_transient_mat{i}(:,year_index);
        mean_global_ann_transient_cncloss_cn45{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_transient_cncloss_cn45{i}(:,ann_index) = nanmean(tmp,2);

        tmp = cn45_cntlai_transient_mat{i}(:,year_index);
        mean_global_ann_transient_cntlai_cn45{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_transient_cntlai_cn45{i}(:,ann_index) = nanmean(tmp,2);

        tmp = cn45_cntotc_transient_mat{i}(:,year_index);
        mean_global_ann_transient_cntotc_cn45{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_transient_cntotc_cn45{i}(:,ann_index) = nanmean(tmp,2);

        tmp = cn45_evpsoil_transient_mat{i}(:,year_index);
        mean_global_ann_transient_evpsoil_cn45{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_transient_evpsoil_cn45{i}(:,ann_index) = nanmean(tmp,2);

        tmp = cn45_evptrns_transient_mat{i}(:,year_index);
        mean_global_ann_transient_evptrns_cn45{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_transient_evptrns_cn45{i}(:,ann_index) = nanmean(tmp,2);

        tmp = cn45_evpintr_transient_mat{i}(:,year_index);
        mean_global_ann_transient_evpintr_cn45{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_transient_evpintr_cn45{i}(:,ann_index) = nanmean(tmp,2);

        tmp = cn45_evpsbln_transient_mat{i}(:,year_index);
        mean_global_ann_transient_evpsbln_cn45{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_transient_evpsbln_cn45{i}(:,ann_index) = nanmean(tmp,2);

        tmp = cn45_et_transient_mat{i}(:,year_index);
        mean_global_ann_transient_et_cn45{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_transient_et_cn45{i}(:,ann_index) = nanmean(tmp,2);

        if (i == 1)

           tmp = cn51_cngpp_280_mat(:,year_index);
           mean_global_ann_280_cngpp_cn51(ann_index) = nanmean(tmp(:));
           mean_local_ann_280_cngpp_cn51(:,ann_index) = nanmean(tmp,2);

           tmp = cn51_cnnee_280_mat(:,year_index);
           mean_global_ann_280_cnnee_cn51(ann_index) = nanmean(tmp(:));
           mean_local_ann_280_cnnee_cn51(:,ann_index) = nanmean(tmp,2);

           tmp = cn51_cncloss_280_mat(:,year_index);
           mean_global_ann_280_cncloss_cn51(ann_index) = nanmean(tmp(:));
           mean_local_ann_280_cncloss_cn51(:,ann_index) = nanmean(tmp,2);

           tmp = cn51_cntlai_280_mat(:,year_index);
           mean_global_ann_280_cntlai_cn51(ann_index) = nanmean(tmp(:));
           mean_local_ann_280_cntlai_cn51(:,ann_index) = nanmean(tmp,2);

           tmp = cn51_cntotc_280_mat(:,year_index);
           mean_global_ann_280_cntotc_cn51(ann_index) = nanmean(tmp(:));
           mean_local_ann_280_cntotc_cn51(:,ann_index) = nanmean(tmp,2);

           tmp = cn51_evpsoil_280_mat(:,year_index);
           mean_global_ann_280_evpsoil_cn51(ann_index) = nanmean(tmp(:));
           mean_local_ann_280_evpsoil_cn51(:,ann_index) = nanmean(tmp,2);

           tmp = cn51_evptrns_280_mat(:,year_index);
           mean_global_ann_280_evptrns_cn51(ann_index) = nanmean(tmp(:));
           mean_local_ann_280_evptrns_cn51(:,ann_index) = nanmean(tmp,2);

           tmp = cn51_evpintr_280_mat(:,year_index);
           mean_global_ann_280_evpintr_cn51(ann_index) = nanmean(tmp(:));
           mean_local_ann_280_evpintr_cn51(:,ann_index) = nanmean(tmp,2);

           tmp = cn51_evpsbln_280_mat(:,year_index);
           mean_global_ann_280_evpsbln_cn51(ann_index) = nanmean(tmp(:));
           mean_local_ann_280_evpsbln_cn51(:,ann_index) = nanmean(tmp,2);

           tmp = cn51_et_280_mat(:,year_index);
           mean_global_ann_280_et_cn51(ann_index) = nanmean(tmp(:));
           mean_local_ann_280_et_cn51(:,ann_index) = nanmean(tmp,2);

           tmp = cn45_cngpp_280_mat(:,year_index);
           mean_global_ann_280_cngpp_cn45(ann_index) = nanmean(tmp(:));
           mean_local_ann_280_cngpp_cn45(:,ann_index) = nanmean(tmp,2);

           tmp = cn45_cnnee_280_mat(:,year_index);
           mean_global_ann_280_cnnee_cn45(ann_index) = nanmean(tmp(:));
           mean_local_ann_280_cnnee_cn45(:,ann_index) = nanmean(tmp,2);

           tmp = cn45_cncloss_280_mat(:,year_index);
           mean_global_ann_280_cncloss_cn45(ann_index) = nanmean(tmp(:));
           mean_local_ann_280_cncloss_cn45(:,ann_index) = nanmean(tmp,2);

           tmp = cn45_cntlai_280_mat(:,year_index);
           mean_global_ann_280_cntlai_cn45(ann_index) = nanmean(tmp(:));
           mean_local_ann_280_cntlai_cn45(:,ann_index) = nanmean(tmp,2);

           tmp = cn45_cntotc_280_mat(:,year_index);
           mean_global_ann_280_cntotc_cn45(ann_index) = nanmean(tmp(:));
           mean_local_ann_280_cntotc_cn45(:,ann_index) = nanmean(tmp,2);

           tmp = cn45_evpsoil_280_mat(:,year_index);
           mean_global_ann_280_evpsoil_cn45(ann_index) = nanmean(tmp(:));
           mean_local_ann_280_evpsoil_cn45(:,ann_index) = nanmean(tmp,2);

           tmp = cn45_evptrns_280_mat(:,year_index);
           mean_global_ann_280_evptrns_cn45(ann_index) = nanmean(tmp(:));
           mean_local_ann_280_evptrns_cn45(:,ann_index) = nanmean(tmp,2);

           tmp = cn45_evpintr_280_mat(:,year_index);
           mean_global_ann_280_evpintr_cn45(ann_index) = nanmean(tmp(:));
           mean_local_ann_280_evpintr_cn45(:,ann_index) = nanmean(tmp,2);

           tmp = cn45_evpsbln_280_mat(:,year_index);
           mean_global_ann_280_evpsbln_cn45(ann_index) = nanmean(tmp(:));
           mean_local_ann_280_evpsbln_cn45(:,ann_index) = nanmean(tmp,2);

           tmp = cn45_et_280_mat(:,year_index);
           mean_global_ann_280_et_cn45(ann_index) = nanmean(tmp(:));
           mean_local_ann_280_et_cn45(:,ann_index) = nanmean(tmp,2);

        end

     end
end

legend_str = {};

color_mat = lines(nr_of_transient_cycles);

figure
for i = 1:nr_of_transient_cycles
    plot(start_year:stop_year,mean_global_ann_transient_cngpp_cn51{i},'LineWidth',2,'LineStyle','-','Color',color_mat(i,:));
    hold on
    legend_str{i} =  ['CN51 cycle ' num2str(i) ' (transient CO2 ' num2str(transient_co2_start_year(i)) ' - ' num2str(transient_co2_stop_year(i)) ')'];
end
for i = 1:nr_of_transient_cycles
    plot(start_year:stop_year,mean_global_ann_transient_cngpp_cn45{i},'LineWidth',2,'LineStyle','--','Color',color_mat(i,:));
    hold on
    legend_str{i+nr_of_transient_cycles} =  ['CN45 cycle ' num2str(i) ' (transient CO2 ' num2str(transient_co2_start_year(i)) ' - ' num2str(transient_co2_stop_year(i)) ')'];
end
plot(start_year:stop_year,mean_global_ann_280_cngpp_cn51,'LineWidth',2,'LineStyle','-','Color','k');
hold on
plot(start_year:stop_year,mean_global_ann_280_cngpp_cn45,'LineWidth',2,'LineStyle','-','Color','m');
hold off
grid on
xlabel('year')
ylabel('GPP [kg m^{-2} s^{-1}]')
title('annual average GPP [kg m^{-2} s^{-1}]')
legend_str{2*nr_of_transient_cycles+1} = 'CN51 cycle 1 (@280ppm)';
legend_str{2*nr_of_transient_cycles+2} = 'CN45 cycle 1 (@280ppm)';
legend(legend_str,'Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/GPP_ann_CN51_CN45_cycles_280_transient'],'png')

% compute and plot monthly average data

mean_global_mon_280_cngpp_cn51 = nanmean(cn51_cngpp_280_mat,1);
mean_global_mon_280_cngpp_cn45 = nanmean(cn45_cngpp_280_mat,1);
mean_global_mon_transient_cngpp_cn51{1} = nanmean(cn51_cngpp_transient_mat{1},1);
mean_global_mon_transient_cngpp_cn45{1} = nanmean(cn45_cngpp_transient_mat{1},1);

figure
plot(datenum_vec,mean_global_mon_transient_cngpp_cn51{1},'LineWidth',2,'LineStyle','-','Color',color_mat(1,:));
hold on
plot(datenum_vec,mean_global_mon_transient_cngpp_cn45{1},'LineWidth',2,'LineStyle','--','Color',color_mat(1,:));
hold on
plot(datenum_vec,mean_global_mon_280_cngpp_cn51,'LineWidth',2,'LineStyle','-','Color','k');
hold on
plot(datenum_vec,mean_global_mon_280_cngpp_cn45,'LineWidth',2,'LineStyle','-','Color','m');
hold off 
grid on
xlabel('year')
ylabel('GPP [kg m^{-2} s^{-1}]')
title('monthly average GPP [kg m^{-2} s^{-1}]')
legend_str{2*nr_of_transient_cycles+1} = 'CN51 cycle 1 (@280ppm)';
legend_str{2*nr_of_transient_cycles+2} = 'CN45 cycle 1 (@280ppm)';
legend(legend_str,'Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/GPP_mon_CN51_CN45_cycles_280_transient'],'png')

