clc
clear all
close all

base_path = '/discover/nobackup/projects/geoscm/fzeng/Catchment-CN5.1/local_runs/';
base_exp_name = 'GEOSldas_CN51_local_280';
%base_path = '/discover/nobackup/jkolassa/CN51_test_local/';
%base_exp_name = 'GEOSldas_CN51_local_280_test';
base_path = '/discover/nobackup/jkolassa/new/';
base_exp_name = 'GEOSldas_CN51_local_280';
base_path = '/discover/nobackup/jkolassa/CN51_config_change/';
base_exp_name = 'GEOSldas_CN51_local_280';
base_path = '/discover/nobackup/jkolassa/CN51_config_test/';
base_exp_name = 'GEOSldas_CN51_local_280';
%base_path = '/discover/nobackup/projects/geoscm/fzeng/Catchment-CN5.1/local_runs/';
%base_exp_name = 'GEOSldas_CN51_local_280';
transient_path = '/discover/nobackup/jkolassa/new/old_cycle_exps/';
transient_exp_name = 'GEOSldas_CN51_local_transient';
nr_of_cycles = 6;
nr_of_transient_cycles = 4;

transient_co2_start_year = [1851; 1891; 1931; 1971];
transient_co2_stop_year = [1891; 1931; 1971; 2011];

cn40_read_path = '/discover/nobackup/projects/geoscm/fzeng/Catchment-CN5.1/reference_runs/GEOSldas_CN40_local/output/SMAP_EASEv2_M09_GLOBAL/cat/ens0000/';
cn45_read_path = '/discover/nobackup/projects/geoscm/fzeng/GEOSldas45_M36_rst_ldas06_15/output/SMAP_EASEv2_M36/cat/ens0000/';
clm5_read_path = ['/discover/nobackup/projects/geoscm/fzeng/CLM/archive/clm5_test/lnd/hist/'];
fluxcom_read_path = '/discover/nobackup/elee15/data/gpp/FluxCom/';
fluxsat_read_path = '/discover/nobackup/elee15/data/gpp/FluxSat/monthly/';
maiac_read_path = '/discover/nobackup/elee15/data/gpp/maiac/';
vpm_read_path = '/discover/nobackup/elee15/data/gpp/vpm/';
clm5_of_path = '/discover/nobackup/jkolassa/CLM5_sim/';
cn51_noStress_read_path = '/discover/nobackup/jkolassa/CN51_noStress/GEOSldas_CN51_local_280/output/SMAP_EASEv2_M36/cat/ens0000/';
cn51_noStress_efpot_read_path = '/discover/nobackup/jkolassa/CN51_noStress_efpot/GEOSldas_CN51_local_280/output/SMAP_EASEv2_M36/cat/ens0000/';
cn51_sat_read_path = '/discover/nobackup/jkolassa/CN51_local_sat/GEOSldas_CN51_local_280/output/SMAP_EASEv2_M36/cat/ens0000/';
cn51_xroot_read_path = '/discover/nobackup/jkolassa/CN51_xroot/GEOSldas_CN51_local_280/output/SMAP_EASEv2_M36/cat/ens0000/';
cn51_wetin_read_path = '/discover/nobackup/jkolassa/CN51_wetin/GEOSldas_CN51_local_280/output/SMAP_EASEv2_M36/cat/ens0000/';

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];

% read PFT distribution data

cn51_rst_file = '/discover/nobackup/projects/geoscm/fzeng/Catchment-CN5.1/local_runs/GEOSldas_CN51_local_280/output/SMAP_EASEv2_M36/rs/ens0000/Y1980/M01/GEOSldas_CN51_local_280.catchcnclm51_internal_rst.19800101_0000';
cn40_rst_file = '/discover/nobackup/projects/geoscm/fzeng/Catchment-CN5.1/reference_runs/GEOSldas_CN40_local/output/SMAP_EASEv2_M09_GLOBAL/rs/ens0000/Y1981/M01/GEOSldas_CN40_local.catchcnclm40_internal_rst.19810101_0000';

cn51_ity = double(ncread(cn51_rst_file,'ITY'));
cn51_fvg = double(ncread(cn51_rst_file,'FVG'));
cn40_ity = double(ncread(cn40_rst_file,'ITY'));
cn40_fvg = double(ncread(cn40_rst_file,'FVG'));

cn51_dom_pft = NaN*ones(length(cn51_ity),1);
cn51_hom_index = zeros(length(cn51_ity),1);
cn51_max_frac = NaN*ones(length(cn51_ity),1);

cn40_dom_pft = NaN*ones(length(cn40_ity),1);
cn40_hom_index = zeros(length(cn40_ity),1);
cn40_max_frac = NaN*ones(length(cn40_ity),1);

for i = 1:length(cn51_ity)

    % CN51
    [max_frac,max_index] = max(cn51_fvg(i,:));

    cn51_dom_pft(i) = cn51_ity(i,max_index);
    cn51_max_frac(i) = max_frac;    

    if (max_frac>0.7)
       cn51_hom_index(i) = 1;
    end

    % CN40
    [max_frac,max_index] = max(cn40_fvg(i,:));

    cn40_dom_pft(i) = cn40_ity(i,max_index);
    cn40_max_frac(i) = max_frac;

    if (max_frac>0.7)
       cn40_hom_index(i) = 1;
    end

end

cn51_pft_name = {...
           'Needleleaf evergreen temperate tree',...
           'Needleleaf evergreen boreal tree',...
           'Needleleaf deciduous boreal tree',...
           'Broadleaf evergreen tropical tree',...
           'Broadleaf evergreen temperate tree',...
           'Broadleaf deciduous tropical tree',...
           'Broadleaf deciduous temperate tree',...
           'Broadleaf deciduous boreal tree',...
           'Broadleaf evergreen temperate shrub',...
           'Broadleaf deciduous temperate shrub',...
           'Broadleaf deciduous boreal shrub',...
           'Arctic C3 grass',...
           'Non-arctic C3 grass',...
           'C4 grass',...
           'C3 crop'...
           };

cn51_pft_name_short = {...
           'NETT',...
           'NEBT',...
           'NDBT',...
           'BETrT',...
           'BETT',...
           'BDTrT',...
           'BDTT',...
           'BDBT',...
           'BETS',...
           'BDTS',...
           'BDBS',...
           'AC3',...
           'NAC3',...
           'C4',...
           'C3Cr'...
           };


% compute mapping indices

cn51_read_file = ['/discover/nobackup/projects/geoscm/fzeng/Catchment-CN5.1/local_runs/GEOSldas_CN51_local_280/output/SMAP_EASEv2_M36/cat/ens0000/Y1980/M01/GEOSldas_CN51_local_280.tavg24_1d_lnd_Nt.monthly.198001.nc4'];
cn45_read_file = [cn45_read_path 'Y2001/M01/GEOSldas45_M36_rst_ldas06_15.tavg24_1d_lnd_Nt.monthly.200101.nc4'];
cn40_read_file = [cn40_read_path 'Y1981/M01/GEOSldas_CN40_local.tavg24_1d_lnd_Nt.monthly.198101.nc4'];
fluxcom_gpp_read_file = [fluxcom_read_path 'GPP.RS.FP-ALL.MLM-ALL.METEO-NONE.720_360.monthly.2001.nc'];
fluxsat_gpp_read_file = [fluxsat_read_path 'SatFlux_GPP_200301.nc'];
maiac_gpp_read_file = [maiac_read_path 'GPP_MAIAC_200301.nc'];
vpm_gpp_read_file = [vpm_read_path 'VPM_GPP_200301.nc'];
clm5_gpp_read_file = [clm5_of_path 'clm50_release-clm5.0.20_1deg_GSWP3V1_isofix2_hist.clm2.h0.GPP.185001-201412.nc'];
clm5_lai_read_file = [clm5_of_path 'clm50_release-clm5.0.20_1deg_GSWP3V1_isofix2_hist.clm2.h0.TLAI.185001-201412.nc'];
clm5_nee_read_file = [clm5_of_path 'clm50_release-clm5.0.20_1deg_GSWP3V1_isofix2_hist.clm2.h0.NEE.185001-201412.nc'];

lat_vec = double(ncread(cn51_read_file,'lat'));
lon_vec = double(ncread(cn51_read_file,'lon'));

cn45_lat_vec = double(ncread(cn45_read_file,'lat'));
cn45_lon_vec = double(ncread(cn45_read_file,'lon'));

cn40_lat_vec = double(ncread(cn40_read_file,'lat'));
cn40_lon_vec = double(ncread(cn40_read_file,'lon'));

fc_lat_vec = ncread(fluxcom_gpp_read_file,'latitude');
fc_lon_vec = ncread(fluxcom_gpp_read_file,'longitude');

fs_lat_vec = ncread(fluxsat_gpp_read_file,'latitude');
fs_lon_vec = ncread(fluxsat_gpp_read_file,'longitude');

maiac_lat_vec = ncread(maiac_gpp_read_file,'latitude');
maiac_lon_vec = ncread(maiac_gpp_read_file,'longitude');

vpm_lat_vec = ncread(vpm_gpp_read_file,'latitude');
vpm_lon_vec = ncread(vpm_gpp_read_file,'longitude');

clm5_read_file = [clm5_read_path 'clm5_test.clm2.h0.0001-01.nc'];
clm5_lat_vec = ncread(clm5_read_file,'lat');
clm5_lon_vec = ncread(clm5_read_file,'lon');
%convert CLM lons from [0 360] to [-180 180]
clm5_lon_vec(clm5_lon_vec>180) = clm5_lon_vec(clm5_lon_vec>180)-360;

% read released CLM simulation data

clm5_of_lat_vec = ncread(clm5_gpp_read_file,'lat');
clm5_of_lon_vec = ncread(clm5_gpp_read_file,'lon');
clm5_of_time_vec = ncread(clm5_gpp_read_file,'time');
clm5_of_gpp_vec = ncread(clm5_gpp_read_file,'GPP');
clm5_of_lai_vec = ncread(clm5_lai_read_file,'TLAI');
clm5_of_nee_vec = ncread(clm5_nee_read_file,'NEE');
%convert CLM lons from [0 360] to [-180 180]
clm5_of_lon_vec(clm5_of_lon_vec>180) = clm5_of_lon_vec(clm5_of_lon_vec>180)-360;

for i = 1:length(lat_vec)
    lat_dist = sqrt((fc_lat_vec-lat_vec(i)).^2);
    fc_lat_index(i) = find(lat_dist==min(lat_dist));

    lon_dist = sqrt((fc_lon_vec-lon_vec(i)).^2);
    fc_lon_index(i) = find(lon_dist==min(lon_dist));

    dist_vec = sqrt((cn45_lat_vec-lat_vec(i)).^2 + (cn45_lon_vec-lon_vec(i)).^2);
    min_index = find(dist_vec==min(dist_vec));
    data_index(i) = min_index(1);

%    dist_vec = sqrt((cn40_lat_vec-lat_vec(i)).^2 + (cn40_lon_vec-lon_vec(i)).^2);
%    min_index = find(dist_vec==min(dist_vec));
%    cn40_data_index(i) = min_index(1);

    lat_dist = sqrt((fs_lat_vec-lat_vec(i)).^2);
    fs_lat_index(i) = find(lat_dist==min(lat_dist));

    lon_dist = sqrt((fs_lon_vec-lon_vec(i)).^2);
    fs_lon_index(i) = find(lon_dist==min(lon_dist));

    lat_dist = sqrt((maiac_lat_vec-lat_vec(i)).^2);
    maiac_lat_index(i) = find(lat_dist==min(lat_dist));

    lon_dist = sqrt((maiac_lon_vec-lon_vec(i)).^2);
    maiac_lon_index(i) = find(lon_dist==min(lon_dist));

    lat_dist = sqrt((vpm_lat_vec-lat_vec(i)).^2);
    vpm_lat_index(i) = find(lat_dist==min(lat_dist));

    lon_dist = sqrt((vpm_lon_vec-lon_vec(i)).^2);
    vpm_lon_index(i) = find(lon_dist==min(lon_dist));

    dist_vec = sqrt((clm5_lat_vec-lat_vec(i)).^2);
    min_index = find(dist_vec==min(dist_vec));
    clm5_lat_index(i) = min_index(1);

    dist_vec = sqrt((clm5_lon_vec-lon_vec(i)).^2);
    min_index = find(dist_vec==min(dist_vec));
    clm5_lon_index(i) = min_index(1);

    dist_vec = sqrt((clm5_of_lat_vec-lat_vec(i)).^2);
    min_index = find(dist_vec==min(dist_vec));
    clm5_of_lat_index(i) = min_index(1);

    dist_vec = sqrt((clm5_of_lon_vec-lon_vec(i)).^2);
    min_index = find(dist_vec==min(dist_vec));
    clm5_of_lon_index(i) = min_index(1);
end

for i = 1:nr_of_cycles

    if (i>1)
       cyc_str = ['c' num2str(i,'%0i')];
       exp_name = [base_exp_name '_' cyc_str];
       transient_name = [transient_exp_name '_' cyc_str];
    else 
       exp_name = base_exp_name;
       cyc_str = ['c' num2str(i,'%0i')];
       transient_name = [transient_exp_name '_' cyc_str];
    end 

    mat_index = 0;
    y_index = 0;
    mat_index_clm5_mod = (start_year-1850)*12;

    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


        cn51_rst_file = [base_path exp_name '/output/SMAP_EASEv2_M36/rs/ens0000/Y' y_str ... 
                        '/M01/' exp_name '.catchcnclm51_internal_rst.' y_str '0101_0000'];

        if (y_index==1)
           srfexc_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1));
           rzexc_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1));
           catdef_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1));
        end

        if exist(cn51_rst_file)
           srfexc_mat{i}(:,y_index) = double(ncread(cn51_rst_file,'SRFEXC'));
           rzexc_mat{i}(:,y_index) = double(ncread(cn51_rst_file,'RZEXC'));
           catdef_mat{i}(:,y_index) = double(ncread(cn51_rst_file,'CATDEF'));
        else
           srfexc_mat{i}(:,y_index) = NaN*ones(size(srfexc_mat{i},1),1);
           rzexc_mat{i}(:,y_index) = NaN*ones(size(srfexc_mat{i},1),1);
           catdef_mat{i}(:,y_index) = NaN*ones(size(srfexc_mat{i},1),1);
        end 

        % loop over months
        for m = 1:12

            mat_index = mat_index + 1;
            mat_index_clm5_mod = mat_index_clm5_mod + 1;

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

            % define file names

            cn51_read_file = [base_path exp_name '/output/SMAP_EASEv2_M36/cat/ens0000/Y' y_str ...
                           '/M' m_str '/' exp_name '.tavg24_1d_lnd_Nt.monthly.' y_str m_str '.nc4'];

            transient_read_file = [transient_path transient_name '/output/SMAP_EASEv2_M36/cat/ens0000/Y' y_str ...
                           '/M' m_str '/' transient_name '.tavg24_1d_lnd_Nt.monthly.' y_str m_str '.nc4'];

            if (i==1)
               cn40_read_file = [cn40_read_path 'Y' y_str '/M' m_str '/GEOSldas_CN40_local.tavg24_1d_lnd_Nt.monthly.' ...
                                y_str m_str '.nc4'];
               cn45_read_file = [cn45_read_path 'Y' y_str '/M' m_str '/GEOSldas45_M36_rst_ldas06_15.tavg24_1d_lnd_Nt.monthly.' y_str m_str '.nc4'];
               fluxcom_gpp_read_file = [fluxcom_read_path 'GPP.RS.FP-ALL.MLM-ALL.METEO-NONE.720_360.monthly.' y_str '.nc'];
               fluxsat_gpp_read_file = [fluxsat_read_path 'SatFlux_GPP_' y_str m_str '.nc'];
               maiac_gpp_read_file = [maiac_read_path 'GPP_MAIAC_' y_str m_str '.nc'];
               vpm_gpp_read_file = [vpm_read_path 'VPM_GPP_' y_str m_str '.nc'];
               clm5_read_file = [clm5_read_path 'clm5_test.clm2.h0.' y_index_str '-' m_str '.nc'];
               cn51_noStress_read_file = [cn51_noStress_read_path 'Y' y_str '/M' m_str '/GEOSldas_CN51_local_280.tavg24_1d_lnd_Nt.monthly.' y_str m_str '.nc4'];
               cn51_noStress_efpot_read_file = [cn51_noStress_efpot_read_path 'Y' y_str '/M' m_str '/GEOSldas_CN51_local_280.tavg24_1d_lnd_Nt.monthly.' y_str m_str '.nc4'];
               cn51_sat_read_file = [cn51_sat_read_path 'Y' y_str '/M' m_str '/GEOSldas_CN51_local_280.tavg24_1d_lnd_Nt.monthly.' y_str m_str '.nc4'];
               cn51_xroot_read_file = [cn51_xroot_read_path 'Y' y_str '/M' m_str '/GEOSldas_CN51_local_280.tavg24_1d_lnd_Nt.monthly.' y_str m_str '.nc4'];
               cn51_wetin_read_file = [cn51_wetin_read_path 'Y' y_str '/M' m_str '/GEOSldas_CN51_local_280.tavg24_1d_lnd_Nt.monthly.' y_str m_str '.nc4'];

            end 
 
            % initialize arrays

            if (mat_index == 1)

                cngpp_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                cnnee_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                cncloss_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                cntlai_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                cntotc_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                et_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                evpsoil_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                evptrns_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                evpintr_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                evpsbln_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);

                cngpp_transient_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                cnnee_transient_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                cncloss_transient_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                cntlai_transient_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                cntotc_transient_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                et_transient_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                evpsoil_transient_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                evptrns_transient_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                evpintr_transient_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                evpsbln_transient_mat{i} = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);

                if (i ==1)
                   cngpp_cn40_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cnnee_cn40_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cncloss_cn40_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cntlai_cn40_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cntotc_cn40_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   et_cn40_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   evpsoil_cn40_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   evptrns_cn40_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   evpintr_cn40_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   evpsbln_cn40_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   fc_gpp_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   fs_gpp_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   maiac_gpp_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   vpm_gpp_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);

                   cngpp_cn45_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cnnee_cn45_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cncloss_cn45_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cntlai_cn45_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cntotc_cn45_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   et_cn45_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   evpsoil_cn45_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   evptrns_cn45_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   evpintr_cn45_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   evpsbln_cn45_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);

                   clm5_gpp_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   clm5_lai_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   clm5_closs_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   clm5_nee_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   clm5_elai_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   clm5_vegc_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   clm5_totc_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);

                   clm5_of_pi_gpp_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   clm5_of_mod_gpp_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   clm5_of_pi_lai_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   clm5_of_mod_lai_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   clm5_of_pi_nee_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   clm5_of_mod_nee_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);

                   cngpp_cn51_noStress_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cnnee_cn51_noStress_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cntlai_cn51_noStress_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);

                   cngpp_cn51_noStress_efpot_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cnnee_cn51_noStress_efpot_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cntlai_cn51_noStress_efpot_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);

                   cngpp_cn51_sat_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cnnee_cn51_sat_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cntlai_cn51_sat_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);

                   cngpp_cn51_xroot_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cnnee_cn51_xroot_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cntlai_cn51_xroot_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);

                   cngpp_cn51_wetin_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cnnee_cn51_wetin_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                   cntlai_cn51_wetin_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
                end
            end % mat_index
  
            % read data

            if exist(cn51_read_file)
               cngpp_mat{i}(:,mat_index) = double(ncread(cn51_read_file,'CNGPP'));
               cnnee_mat{i}(:,mat_index) = double(ncread(cn51_read_file,'CNNEE'));
               cncloss_mat{i}(:,mat_index) = double(ncread(cn51_read_file,'CLOSS'));
               cntlai_mat{i}(:,mat_index) = double(ncread(cn51_read_file,'CNTLAI'));
               cntotc_mat{i}(:,mat_index) = double(ncread(cn51_read_file,'CNTOTC'));
               evpsoil_mat{i}(:,mat_index) = double(ncread(cn51_read_file,'EVPSOIL'));
               evptrns_mat{i}(:,mat_index) = double(ncread(cn51_read_file,'EVPTRNS'));
               evpintr_mat{i}(:,mat_index) = double(ncread(cn51_read_file,'EVPINTR'));
               evpsbln_mat{i}(:,mat_index) = double(ncread(cn51_read_file,'EVPSBLN'));
               et_mat{i}(:,mat_index) = evpsoil_mat{i}(:,mat_index) + evptrns_mat{i}(:,mat_index) ...
                                      + evpintr_mat{i}(:,mat_index) + evpsbln_mat{i}(:,mat_index);
            else 
               cngpp_mat{i}(:,mat_index) = NaN*ones(size(cngpp_mat{i},1),1);
               cnnee_mat{i}(:,mat_index) = NaN*ones(size(cnnee_mat{i},1),1);
               cncloss_mat{i}(:,mat_index) = NaN*ones(size(cncloss_mat{i},1),1);
               cntlai_mat{i}(:,mat_index) = NaN*ones(size(cntlai_mat{i},1),1);
               cntotc_mat{i}(:,mat_index) = NaN*ones(size(cntotc_mat{i},1),1);
               evpsoil_mat{i}(:,mat_index) = NaN*ones(size(evpsoil_mat{i},1),1);
               evptrns_mat{i}(:,mat_index) = NaN*ones(size(evptrns_mat{i},1),1);
               evpintr_mat{i}(:,mat_index) = NaN*ones(size(evpintr_mat{i},1),1);
               evpsbln_mat{i}(:,mat_index) = NaN*ones(size(evpsbln_mat{i},1),1);
               et_mat{i}(:,mat_index) = NaN*ones(size(et_mat{i},1),1);
            end

            if exist(transient_read_file)
               cngpp_transient_mat{i}(:,mat_index) = double(ncread(transient_read_file,'CNGPP'));
               cnnee_transient_mat{i}(:,mat_index) = double(ncread(transient_read_file,'CNNEE'));
               cncloss_transient_mat{i}(:,mat_index) = double(ncread(transient_read_file,'CLOSS'));
               cntlai_transient_mat{i}(:,mat_index) = double(ncread(transient_read_file,'CNTLAI'));
               cntotc_transient_mat{i}(:,mat_index) = double(ncread(transient_read_file,'CNTOTC'));
               evpsoil_transient_mat{i}(:,mat_index) = double(ncread(transient_read_file,'EVPSOIL'));
               evptrns_transient_mat{i}(:,mat_index) = double(ncread(transient_read_file,'EVPTRNS'));
               evpintr_transient_mat{i}(:,mat_index) = double(ncread(transient_read_file,'EVPINTR'));
               evpsbln_transient_mat{i}(:,mat_index) = double(ncread(transient_read_file,'EVPSBLN'));
               et_transient_mat{i}(:,mat_index) = evpsoil_transient_mat{i}(:,mat_index) + evptrns_transient_mat{i}(:,mat_index) ...
                                      + evpintr_transient_mat{i}(:,mat_index) + evpsbln_transient_mat{i}(:,mat_index);
            else 
               cngpp_transient_mat{i}(:,mat_index) = NaN*ones(size(cngpp_transient_mat{i},1),1);
               cnnee_transient_mat{i}(:,mat_index) = NaN*ones(size(cnnee_transient_mat{i},1),1);
               cncloss_transient_mat{i}(:,mat_index) = NaN*ones(size(cncloss_transient_mat{i},1),1);
               cntlai_transient_mat{i}(:,mat_index) = NaN*ones(size(cntlai_transient_mat{i},1),1);
               cntotc_transient_mat{i}(:,mat_index) = NaN*ones(size(cntotc_transient_mat{i},1),1);
               evpsoil_transient_mat{i}(:,mat_index) = NaN*ones(size(evpsoil_transient_mat{i},1),1);
               evptrns_transient_mat{i}(:,mat_index) = NaN*ones(size(evptrns_transient_mat{i},1),1);
               evpintr_transient_mat{i}(:,mat_index) = NaN*ones(size(evpintr_transient_mat{i},1),1);
               evpsbln_transient_mat{i}(:,mat_index) = NaN*ones(size(evpsbln_transient_mat{i},1),1);
               et_transient_mat{i}(:,mat_index) = NaN*ones(size(et_transient_mat{i},1),1);
            end 

            if (i == 1)
               if exist(cn40_read_file)
                  cngpp_cn40_mat(:,mat_index) = double(ncread(cn40_read_file,'CNGPP'));
                  cnnee_cn40_mat(:,mat_index) = double(ncread(cn40_read_file,'CNNEE'));
                  cncloss_cn40_mat(:,mat_index) = double(ncread(cn40_read_file,'CLOSS'));
                  cntlai_cn40_mat(:,mat_index) = double(ncread(cn40_read_file,'CNTLAI'));
                  cntotc_cn40_mat(:,mat_index) = double(ncread(cn40_read_file,'CNTOTC'));
                  evpsoil_cn40_mat(:,mat_index) = double(ncread(cn40_read_file,'EVPSOIL'));
                  evptrns_cn40_mat(:,mat_index) = double(ncread(cn40_read_file,'EVPTRNS'));
                  evpintr_cn40_mat(:,mat_index) = double(ncread(cn40_read_file,'EVPINTR'));
                  evpsbln_cn40_mat(:,mat_index) = double(ncread(cn40_read_file,'EVPSBLN'));
                  et_cn40_mat(:,mat_index) = evpsoil_cn40_mat(:,mat_index) + evptrns_cn40_mat(:,mat_index) ...
                                              + evpintr_cn40_mat(:,mat_index) + evpsbln_cn40_mat(:,mat_index);
                else
                  cngpp_cn40_mat(:,mat_index) = NaN*ones(size(cngpp_cn40_mat,1),1);
                  cnnee_cn40_mat(:,mat_index) = NaN*ones(size(cnnee_cn40_mat,1),1);
                  cncloss_cn40_mat(:,mat_index) = NaN*ones(size(cncloss_cn40_mat,1),1);
                  cntlai_cn40_mat(:,mat_index) = NaN*ones(size(cntlai_cn40_mat,1),1);
                  cntotc_cn40_mat(:,mat_index) = NaN*ones(size(cntotc_cn40_mat,1),1);
                  evpsoil_cn40_mat(:,mat_index) = NaN*ones(size(evpsoil_cn40_mat,1),1);
                  evptrns_cn40_mat(:,mat_index) = NaN*ones(size(evptrns_cn40_mat,1),1);
                  evpintr_cn40_mat(:,mat_index) = NaN*ones(size(evpintr_cn40_mat,1),1);
                  evpsbln_cn40_mat(:,mat_index) = NaN*ones(size(evpsbln_cn40_mat,1),1);
                  et_cn40_mat(:,mat_index) = NaN*ones(size(et_cn40_mat,1),1);
                end

               if exist(cn45_read_file)
                  tmp_mat = double(ncread(cn45_read_file,'CNGPP'));
                  cngpp_cn45_mat(:,mat_index) = tmp_mat(data_index);
                  tmp_mat = double(ncread(cn45_read_file,'CNNEE'));
                  cnnee_cn45_mat(:,mat_index) = tmp_mat(data_index);
                  tmp_mat = double(ncread(cn45_read_file,'CLOSS'));
                  cncloss_cn45_mat(:,mat_index) = tmp_mat(data_index);
                  tmp_mat = double(ncread(cn45_read_file,'CNTLAI'));
                  cntlai_cn45_mat(:,mat_index) = tmp_mat(data_index);
                  tmp_mat = double(ncread(cn45_read_file,'CNTOTC'));
                  cntotc_cn45_mat(:,mat_index) = tmp_mat(data_index);
                  tmp_mat = double(ncread(cn45_read_file,'EVPTRNS'));
                  evptrns_cn45_mat(:,mat_index) = tmp_mat(data_index);
                  tmp_mat = double(ncread(cn45_read_file,'EVPSOIL'));
                  evpsoil_cn45_mat(:,mat_index) = tmp_mat(data_index);
                  tmp_mat = double(ncread(cn45_read_file,'EVPINTR'));
                  evpintr_cn45_mat(:,mat_index) = tmp_mat(data_index);
                  tmp_mat = double(ncread(cn45_read_file,'EVPSBLN'));
                  evpsbln_cn45_mat(:,mat_index) = tmp_mat(data_index);
                  et_cn45_mat(:,mat_index) = evpsoil_cn45_mat(:,mat_index) + evptrns_cn45_mat(:,mat_index) ...
                                              + evpintr_cn45_mat(:,mat_index) + evpsbln_cn45_mat(:,mat_index);
                else
                  cngpp_cn45_mat(:,mat_index) = NaN*ones(size(cngpp_cn45_mat,1),1);
                  cnnee_cn45_mat(:,mat_index) = NaN*ones(size(cnnee_cn45_mat,1),1);
                  cncloss_cn45_mat(:,mat_index) = NaN*ones(size(cncloss_cn45_mat,1),1);
                  cntlai_cn45_mat(:,mat_index) = NaN*ones(size(cntlai_cn45_mat,1),1);
                  cntotc_cn45_mat(:,mat_index) = NaN*ones(size(cntotc_cn45_mat,1),1);
                  evpsoil_cn45_mat(:,mat_index) = NaN*ones(size(evpsoil_cn45_mat,1),1);
                  evptrns_cn45_mat(:,mat_index) = NaN*ones(size(evptrns_cn45_mat,1),1);
                  evpintr_cn45_mat(:,mat_index) = NaN*ones(size(evpintr_cn45_mat,1),1);
                  evpsbln_cn45_mat(:,mat_index) = NaN*ones(size(evpsbln_cn45_mat,1),1);
                  et_cn45_mat(:,mat_index) = NaN*ones(size(et_cn45_mat,1),1);
                end

                if exist(cn51_noStress_read_file)
                  cngpp_cn51_noStress_mat(:,mat_index) = double(ncread(cn51_noStress_read_file,'CNGPP'));
                  cnnee_cn51_noStress_mat(:,mat_index) = double(ncread(cn51_noStress_read_file,'CNNEE'));
                  cntlai_cn51_noStress_mat(:,mat_index) = double(ncread(cn51_noStress_read_file,'CNTLAI'));
                else
                  cngpp_cn51_noStress_mat(:,mat_index) = NaN*ones(size(cngpp_cn51_noStress_mat,1),1);
                  cnnee_cn51_noStress_mat(:,mat_index) = NaN*ones(size(cnnee_cn51_noStress_mat,1),1);
                  cntlai_cn51_noStress_mat(:,mat_index) = NaN*ones(size(cntlai_cn51_noStress_mat,1),1);
                end

                if exist(cn51_noStress_efpot_read_file)
                  cngpp_cn51_noStress_efpot_mat(:,mat_index) = double(ncread(cn51_noStress_efpot_read_file,'CNGPP'));
                  cnnee_cn51_noStress_efpot_mat(:,mat_index) = double(ncread(cn51_noStress_efpot_read_file,'CNNEE'));
                  cntlai_cn51_noStress_efpot_mat(:,mat_index) = double(ncread(cn51_noStress_efpot_read_file,'CNTLAI'));
                else
                  cngpp_cn51_noStress_efpot_mat(:,mat_index) = NaN*ones(size(cngpp_cn51_noStress_efpot_mat,1),1);
                  cnnee_cn51_noStress_efpot_mat(:,mat_index) = NaN*ones(size(cnnee_cn51_noStress_efpot_mat,1),1);
                  cntlai_cn51_noStress_efpot_mat(:,mat_index) = NaN*ones(size(cntlai_cn51_noStress_efpot_mat,1),1);
                end

                if exist(cn51_sat_read_file)
                  cngpp_cn51_sat_mat(:,mat_index) = double(ncread(cn51_sat_read_file,'CNGPP'));
                  cnnee_cn51_sat_mat(:,mat_index) = double(ncread(cn51_sat_read_file,'CNNEE'));
                  cntlai_cn51_sat_mat(:,mat_index) = double(ncread(cn51_sat_read_file,'CNTLAI'));
                else
                  cngpp_cn51_sat_mat(:,mat_index) = NaN*ones(size(cngpp_cn51_sat_mat,1),1);
                  cnnee_cn51_sat_mat(:,mat_index) = NaN*ones(size(cnnee_cn51_sat_mat,1),1);
                  cntlai_cn51_sat_mat(:,mat_index) = NaN*ones(size(cntlai_cn51_sat_mat,1),1);
                end

                if exist(cn51_xroot_read_file)
                  cngpp_cn51_xroot_mat(:,mat_index) = double(ncread(cn51_xroot_read_file,'CNGPP'));
                  cnnee_cn51_xroot_mat(:,mat_index) = double(ncread(cn51_xroot_read_file,'CNNEE'));
                  cntlai_cn51_xroot_mat(:,mat_index) = double(ncread(cn51_xroot_read_file,'CNTLAI'));
                else
                  cngpp_cn51_xroot_mat(:,mat_index) = NaN*ones(size(cngpp_cn51_xroot_mat,1),1);
                  cnnee_cn51_xroot_mat(:,mat_index) = NaN*ones(size(cnnee_cn51_xroot_mat,1),1);
                  cntlai_cn51_xroot_mat(:,mat_index) = NaN*ones(size(cntlai_cn51_xroot_mat,1),1);
                end

                if exist(cn51_wetin_read_file)
                  cngpp_cn51_wetin_mat(:,mat_index) = double(ncread(cn51_wetin_read_file,'CNGPP'));
                  cnnee_cn51_wetin_mat(:,mat_index) = double(ncread(cn51_wetin_read_file,'CNNEE'));
                  cntlai_cn51_wetin_mat(:,mat_index) = double(ncread(cn51_wetin_read_file,'CNTLAI'));
                else 
                  cngpp_cn51_wetin_mat(:,mat_index) = NaN*ones(size(cngpp_cn51_wetin_mat,1),1);
                  cnnee_cn51_wetin_mat(:,mat_index) = NaN*ones(size(cnnee_cn51_wetin_mat,1),1);
                  cntlai_cn51_wetin_mat(:,mat_index) = NaN*ones(size(cntlai_cn51_wetin_mat,1),1);
                end

                if exist(clm5_read_file)
                   tmp_gpp = double(ncread(clm5_read_file,'GPP'));
                   tmp_nee = double(ncread(clm5_read_file,'NEE'));
                   tmp_closs = double(ncread(clm5_read_file,'COL_FIRE_CLOSS'));
                   tmp_lai = double(ncread(clm5_read_file,'TLAI'));
                   tmp_elai = double(ncread(clm5_read_file,'ELAI'));
                   tmp_vegc = double(ncread(clm5_read_file,'TOTVEGC'));
                   tmp_totc = double(ncread(clm5_read_file,'TOTCOLC'));
                   for l = 1:length(lat_vec)
                       clm5_gpp_mat(l,mat_index) = tmp_gpp(clm5_lon_index(l),clm5_lat_index(l));
                       clm5_nee_mat(l,mat_index) = tmp_nee(clm5_lon_index(l),clm5_lat_index(l));
                       clm5_closs_mat(l,mat_index) = tmp_closs(clm5_lon_index(l),clm5_lat_index(l));
                       clm5_elai_mat(l,mat_index) = tmp_elai(clm5_lon_index(l),clm5_lat_index(l));
                       clm5_lai_mat(l,mat_index) = tmp_lai(clm5_lon_index(l),clm5_lat_index(l));
                       clm5_vegc_mat(l,mat_index) = tmp_vegc(clm5_lon_index(l),clm5_lat_index(l));
                       clm5_totc_mat(l,mat_index) = tmp_totc(clm5_lon_index(l),clm5_lat_index(l));
                   end
                 else
                   clm5_gpp_mat(:,mat_index) = NaN*ones(length(lat_vec),1);
                   clm5_nee_mat(:,mat_index) = NaN*ones(length(lat_vec),1);
                   clm5_closs_mat(:,mat_index) = NaN*ones(length(lat_vec),1);
                   clm5_lai_mat(:,mat_index) = NaN*ones(length(lat_vec),1);
                   clm5_elai_mat(:,mat_index) = NaN*ones(length(lat_vec),1);
                   clm5_vegc_mat(:,mat_index) = NaN*ones(length(lat_vec),1);
                   clm5_totc_mat(:,mat_index) = NaN*ones(length(lat_vec),1);
                  end

                if (exist(fluxsat_gpp_read_file)~=0)
                   fs_gpp_tmp_mat = ncread(fluxsat_gpp_read_file,'GPP');
                   for k = 1:length(lat_vec)
                       fs_gpp_mat(k,mat_index) = fs_gpp_tmp_mat(fs_lon_index(k),fs_lat_index(k));
                   end
                end

                if (exist(maiac_gpp_read_file)~=0)
                   maiac_gpp_tmp_mat = ncread(maiac_gpp_read_file,'GPP');
                   maiac_gpp_tmp_mat(maiac_gpp_tmp_mat==-999) = NaN;
                   for k = 1:length(lat_vec)
                       maiac_gpp_mat(k,mat_index) = maiac_gpp_tmp_mat(maiac_lon_index(k),maiac_lat_index(k));
                   end
                end

                if (exist(vpm_gpp_read_file)~=0)
                   vpm_gpp_tmp_mat = ncread(vpm_gpp_read_file,'GPP');
                   for k = 1:length(lat_vec)
                       vpm_gpp_mat(k,mat_index) = vpm_gpp_tmp_mat(vpm_lon_index(k),vpm_lat_index(k));
                   end
                end

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

                for l = 1:length(lat_vec)
                    clm5_of_pi_gpp_mat(l,mat_index) = clm5_of_gpp_vec(clm5_of_lon_index(l),clm5_of_lat_index(l),mat_index);
                    clm5_of_pi_lai_mat(l,mat_index) = clm5_of_lai_vec(clm5_of_lon_index(l),clm5_of_lat_index(l),mat_index);
                    clm5_of_pi_nee_mat(l,mat_index) = clm5_of_nee_vec(clm5_of_lon_index(l),clm5_of_lat_index(l),mat_index);
                    if (mat_index_clm5_mod<=size(clm5_of_gpp_vec,3))
                       clm5_of_mod_gpp_mat(l,mat_index) = clm5_of_gpp_vec(clm5_of_lon_index(l),clm5_of_lat_index(l),mat_index_clm5_mod);
                       clm5_of_mod_lai_mat(l,mat_index) = clm5_of_lai_vec(clm5_of_lon_index(l),clm5_of_lat_index(l),mat_index_clm5_mod);
                       clm5_of_mod_nee_mat(l,mat_index) = clm5_of_nee_vec(clm5_of_lon_index(l),clm5_of_lat_index(l),mat_index_clm5_mod);
                    end
                end

            end % i == 1
            
        end % m

        if (i == 1)
           if (exist(fluxcom_gpp_read_file)~=0)
              gpp_tmp_mat = ncread(fluxcom_gpp_read_file,'GPP');
              for k = 1:length(lat_vec)
                  start_index = (y_index-1)*12 + 1;
                  stop_index = y_index*12;
                  fc_gpp_mat(k,start_index:stop_index) = gpp_tmp_mat(fc_lon_index(k),fc_lat_index(k),:);
               end
           end
        end
      end % y
end % i

% convert CLM5 carbon fluxes from gC/m2/s to kgC/m2/s
% and states from gc/m2 to kgC/m2

clm5_gpp_mat = clm5_gpp_mat*10^(-3);
clm5_nee_mat = clm5_nee_mat*10^(-3);
clm5_closs_mat = clm5_closs_mat*10^(-3);
clm5_vegc_mat = clm5_vegc_mat*10^(-3);
clm5_totc_mat = clm5_totc_mat*10^(-3);

clm5_of_pi_gpp_mat = clm5_of_pi_gpp_mat*10^(-3);
clm5_of_mod_gpp_mat = clm5_of_mod_gpp_mat*10^(-3);
clm5_of_pi_nee_mat = clm5_of_pi_nee_mat*10^(-3);
clm5_of_mod_nee_mat = clm5_of_mod_nee_mat*10^(-3);

% convert observed GPP from gC m^-2 day^-1 to kgC m^-2 s^-1

fc_gpp_mat = fc_gpp_mat*10^(-3);
fc_gpp_mat = fc_gpp_mat./86400;

fs_gpp_mat = fs_gpp_mat*10^(-3);
fs_gpp_mat = fs_gpp_mat./86400;

maiac_gpp_mat = maiac_gpp_mat*10^(-3);
maiac_gpp_mat = maiac_gpp_mat./86400;

vpm_gpp_mat = vpm_gpp_mat*10^(-3);
vpm_gpp_mat = vpm_gpp_mat./86400;

for i = 1:nr_of_cycles

    ann_index = 0;

    mean_srfexc_mat{i} = nanmean(srfexc_mat{i},1);
    mean_rzexc_mat{i} = nanmean(rzexc_mat{i},1);
    mean_catdef_mat{i} = nanmean(catdef_mat{i},1);

    for y = start_year:stop_year

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

        tmp = cngpp_mat{i}(:,year_index);
        mean_global_ann_cngpp{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_cngpp{i}(:,ann_index) = nanmean(tmp,2);

        tmp = cnnee_mat{i}(:,year_index);
        mean_global_ann_cnnee{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_cnnee{i}(:,ann_index) = nanmean(tmp,2);

        tmp = cncloss_mat{i}(:,year_index);
        mean_global_ann_cncloss{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_cncloss{i}(:,ann_index) = nanmean(tmp,2);

        tmp = cntlai_mat{i}(:,year_index);
        mean_global_ann_cntlai{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_cntlai{i}(:,ann_index) = nanmean(tmp,2);

        tmp = cntotc_mat{i}(:,year_index);
        mean_global_ann_cntotc{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_cntotc{i}(:,ann_index) = nanmean(tmp,2);

        tmp = evpsoil_mat{i}(:,year_index);
        mean_global_ann_evpsoil{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_evpsoil{i}(:,ann_index) = nanmean(tmp,2);

        tmp = evptrns_mat{i}(:,year_index);
        mean_global_ann_evptrns{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_evptrns{i}(:,ann_index) = nanmean(tmp,2);

        tmp = evpintr_mat{i}(:,year_index);
        mean_global_ann_evpintr{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_evpintr{i}(:,ann_index) = nanmean(tmp,2);

        tmp = evpsbln_mat{i}(:,year_index);
        mean_global_ann_evpsbln{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_evpsbln{i}(:,ann_index) = nanmean(tmp,2);

        tmp = et_mat{i}(:,year_index);
        mean_global_ann_et{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_et{i}(:,ann_index) = nanmean(tmp,2);

        tmp = cngpp_transient_mat{i}(:,year_index);
        mean_global_ann_transient_cngpp{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_transient_cngpp{i}(:,ann_index) = nanmean(tmp,2);

        tmp = cnnee_transient_mat{i}(:,year_index);
        mean_global_ann_transient_cnnee{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_transient_cnnee{i}(:,ann_index) = nanmean(tmp,2);

        tmp = cncloss_transient_mat{i}(:,year_index);
        mean_global_ann_transient_cncloss{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_transient_cncloss{i}(:,ann_index) = nanmean(tmp,2);

        tmp = cntlai_transient_mat{i}(:,year_index);
        mean_global_ann_transient_cntlai{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_transient_cntlai{i}(:,ann_index) = nanmean(tmp,2);

        tmp = cntotc_transient_mat{i}(:,year_index);
        mean_global_ann_transient_cntotc{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_transient_cntotc{i}(:,ann_index) = nanmean(tmp,2);

        tmp = evpsoil_transient_mat{i}(:,year_index);
        mean_global_ann_transient_evpsoil{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_transient_evpsoil{i}(:,ann_index) = nanmean(tmp,2);

        tmp = evptrns_transient_mat{i}(:,year_index);
        mean_global_ann_transient_evptrns{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_transient_evptrns{i}(:,ann_index) = nanmean(tmp,2);

        tmp = evpintr_transient_mat{i}(:,year_index);
        mean_global_ann_transient_evpintr{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_transient_evpintr{i}(:,ann_index) = nanmean(tmp,2);

        tmp = evpsbln_transient_mat{i}(:,year_index);
        mean_global_ann_transient_evpsbln{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_transient_evpsbln{i}(:,ann_index) = nanmean(tmp,2);

        tmp = et_transient_mat{i}(:,year_index);
        mean_global_ann_transient_et{i}(ann_index) = nanmean(tmp(:));
        mean_local_ann_transient_et{i}(:,ann_index) = nanmean(tmp,2);

        if (i == 1)

           tmp = cngpp_cn40_mat(:,year_index);
           mean_global_ann_cn40_cngpp(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn40_cngpp(:,ann_index) = nanmean(tmp,2);

           tmp = cnnee_cn40_mat(:,year_index);
           mean_global_ann_cn40_cnnee(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn40_cnnee(:,ann_index) = nanmean(tmp,2);

           tmp = cncloss_cn40_mat(:,year_index);
           mean_global_ann_cn40_cncloss(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn40_cncloss(:,ann_index) = nanmean(tmp,2);

           tmp = cntlai_cn40_mat(:,year_index);
           mean_global_ann_cn40_cntlai(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn40_cntlai(:,ann_index) = nanmean(tmp,2);

           tmp = cntotc_cn40_mat(:,year_index);
           mean_global_ann_cn40_cntotc(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn40_cntotc(:,ann_index) = nanmean(tmp,2);

           tmp = evpsoil_cn40_mat(:,year_index);
           mean_global_ann_cn40_evpsoil(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn40_evpsoil(:,ann_index) = nanmean(tmp,2);

           tmp = evptrns_cn40_mat(:,year_index);
           mean_global_ann_cn40_evptrns(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn40_evptrns(:,ann_index) = nanmean(tmp,2);

           tmp = evpintr_cn40_mat(:,year_index);
           mean_global_ann_cn40_evpintr(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn40_evpintr(:,ann_index) = nanmean(tmp,2);

           tmp = evpsbln_cn40_mat(:,year_index);
           mean_global_ann_cn40_evpsbln(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn40_evpsbln(:,ann_index) = nanmean(tmp,2);

           tmp = et_cn40_mat(:,year_index);
           mean_global_ann_cn40_et(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn40_et(:,ann_index) = nanmean(tmp,2);

           tmp = fc_gpp_mat(:,year_index);
           mean_global_ann_fc_gpp(ann_index) = nanmean(tmp(:));
           mean_local_ann_fc_gpp(:,ann_index) = nanmean(tmp,2);

           tmp = fs_gpp_mat(:,year_index);
           mean_global_ann_fs_gpp(ann_index) = nanmean(tmp(:));
           mean_local_ann_fs_gpp(:,ann_index) = nanmean(tmp,2);

           tmp = maiac_gpp_mat(:,year_index);
           mean_global_ann_maiac_gpp(ann_index) = nanmean(tmp(:));
           mean_local_ann_maiac_gpp(:,ann_index) = nanmean(tmp,2);

           tmp = vpm_gpp_mat(:,year_index);
           mean_global_ann_vpm_gpp(ann_index) = nanmean(tmp(:));
           mean_local_ann_vpm_gpp(:,ann_index) = nanmean(tmp,2);

           tmp = cngpp_cn45_mat(:,year_index);
           mean_global_ann_cn45_cngpp(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn45_cngpp(:,ann_index) = nanmean(tmp,2);

           tmp = cnnee_cn45_mat(:,year_index);
           mean_global_ann_cn45_cnnee(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn45_cnnee(:,ann_index) = nanmean(tmp,2);

           tmp = cncloss_cn45_mat(:,year_index);
           mean_global_ann_cn45_cncloss(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn45_cncloss(:,ann_index) = nanmean(tmp,2);

           tmp = cntlai_cn45_mat(:,year_index);
           mean_global_ann_cn45_cntlai(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn45_cntlai(:,ann_index) = nanmean(tmp,2);

           tmp = cntotc_cn45_mat(:,year_index);
           mean_global_ann_cn45_cntotc(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn45_cntotc(:,ann_index) = nanmean(tmp,2);

           tmp = evpsoil_cn45_mat(:,year_index);
           mean_global_ann_cn45_evpsoil(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn45_evpsoil(:,ann_index) = nanmean(tmp,2);

           tmp = evptrns_cn45_mat(:,year_index);
           mean_global_ann_cn45_evptrns(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn45_evptrns(:,ann_index) = nanmean(tmp,2);

           tmp = evpintr_cn45_mat(:,year_index);
           mean_global_ann_cn45_evpintr(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn45_evpintr(:,ann_index) = nanmean(tmp,2);

           tmp = evpsbln_cn45_mat(:,year_index);
           mean_global_ann_cn45_evpsbln(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn45_evpsbln(:,ann_index) = nanmean(tmp,2);

           tmp = et_cn45_mat(:,year_index);
           mean_global_ann_cn45_et(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn45_et(:,ann_index) = nanmean(tmp,2);

           tmp = clm5_gpp_mat(:,year_index);
           mean_global_ann_clm5_cngpp(ann_index) = nanmean(tmp(:));
           mean_local_ann_clm5_cngpp(:,ann_index) = nanmean(tmp,2);

           tmp = clm5_nee_mat(:,year_index);
           mean_global_ann_clm5_cnnee(ann_index) = nanmean(tmp(:));
           mean_local_ann_clm5_cnnee(:,ann_index) = nanmean(tmp,2);

           tmp = clm5_closs_mat(:,year_index);
           mean_global_ann_clm5_cncloss(ann_index) = nanmean(tmp(:));
           mean_local_ann_clm5_cncloss(:,ann_index) = nanmean(tmp,2);

           tmp = clm5_lai_mat(:,year_index);
           mean_global_ann_clm5_cntlai(ann_index) = nanmean(tmp(:));
           mean_local_ann_clm5_cntlai(:,ann_index) = nanmean(tmp,2);

           tmp = clm5_totc_mat(:,year_index);
           mean_global_ann_clm5_cntotc(ann_index) = nanmean(tmp(:));
           mean_local_ann_clm5_cntotc(:,ann_index) = nanmean(tmp,2);

           tmp = clm5_of_pi_gpp_mat(:,year_index);
           mean_global_ann_clm5_of_pi_cngpp(ann_index) = nanmean(tmp(:));
           mean_local_ann_clm5_of_pi_cngpp(:,ann_index) = nanmean(tmp,2);

           tmp = clm5_of_mod_gpp_mat(:,year_index);
           mean_global_ann_clm5_of_mod_cngpp(ann_index) = nanmean(tmp(:));
           mean_local_ann_clm5_of_mod_cngpp(:,ann_index) = nanmean(tmp,2);

           tmp = clm5_of_pi_lai_mat(:,year_index);
           mean_global_ann_clm5_of_pi_cntlai(ann_index) = nanmean(tmp(:));
           mean_local_ann_clm5_of_pi_cntlai(:,ann_index) = nanmean(tmp,2);

           tmp = clm5_of_mod_lai_mat(:,year_index);
           mean_global_ann_clm5_of_mod_cntlai(ann_index) = nanmean(tmp(:));
           mean_local_ann_clm5_of_mod_cntlai(:,ann_index) = nanmean(tmp,2);

           tmp = clm5_of_pi_nee_mat(:,year_index);
           mean_global_ann_clm5_of_pi_cnnee(ann_index) = nanmean(tmp(:));
           mean_local_ann_clm5_of_pi_cnnee(:,ann_index) = nanmean(tmp,2);

           tmp = clm5_of_mod_nee_mat(:,year_index);
           mean_global_ann_clm5_of_mod_cnnee(ann_index) = nanmean(tmp(:));
           mean_local_ann_clm5_of_mod_cnnee(:,ann_index) = nanmean(tmp,2);

           tmp = cngpp_cn51_noStress_mat(:,year_index);
           mean_global_ann_cn51_noStress_cngpp(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn51_noStress_cngpp(:,ann_index) = nanmean(tmp,2);

           tmp = cnnee_cn51_noStress_mat(:,year_index);
           mean_global_ann_cn51_noStress_cnnee(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn51_noStress_cnnee(:,ann_index) = nanmean(tmp,2);

           tmp = cntlai_cn51_noStress_mat(:,year_index);
           mean_global_ann_cn51_noStress_cntlai(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn51_noStress_cntlai(:,ann_index) = nanmean(tmp,2);

           tmp = cngpp_cn51_noStress_efpot_mat(:,year_index);
           mean_global_ann_cn51_noStress_efpot_cngpp(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn51_noStress_efpot_cngpp(:,ann_index) = nanmean(tmp,2);

           tmp = cnnee_cn51_noStress_efpot_mat(:,year_index);
           mean_global_ann_cn51_noStress_efpot_cnnee(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn51_noStress_efpot_cnnee(:,ann_index) = nanmean(tmp,2);

           tmp = cntlai_cn51_noStress_efpot_mat(:,year_index);
           mean_global_ann_cn51_noStress_efpot_cntlai(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn51_noStress_efpot_cntlai(:,ann_index) = nanmean(tmp,2);

           tmp = cngpp_cn51_sat_mat(:,year_index);
           mean_global_ann_cn51_sat_cngpp(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn51_sat_cngpp(:,ann_index) = nanmean(tmp,2);

           tmp = cnnee_cn51_sat_mat(:,year_index);
           mean_global_ann_cn51_sat_cnnee(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn51_sat_cnnee(:,ann_index) = nanmean(tmp,2);

           tmp = cntlai_cn51_sat_mat(:,year_index);
           mean_global_ann_cn51_sat_cntlai(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn51_sat_cntlai(:,ann_index) = nanmean(tmp,2);

           tmp = cngpp_cn51_xroot_mat(:,year_index);
           mean_global_ann_cn51_xroot_cngpp(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn51_xroot_cngpp(:,ann_index) = nanmean(tmp,2);

           tmp = cnnee_cn51_xroot_mat(:,year_index);
           mean_global_ann_cn51_xroot_cnnee(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn51_xroot_cnnee(:,ann_index) = nanmean(tmp,2);

           tmp = cntlai_cn51_xroot_mat(:,year_index);
           mean_global_ann_cn51_xroot_cntlai(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn51_xroot_cntlai(:,ann_index) = nanmean(tmp,2);

           tmp = cngpp_cn51_wetin_mat(:,year_index);
           mean_global_ann_cn51_wetin_cngpp(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn51_wetin_cngpp(:,ann_index) = nanmean(tmp,2);

           tmp = cnnee_cn51_wetin_mat(:,year_index);
           mean_global_ann_cn51_wetin_cnnee(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn51_wetin_cnnee(:,ann_index) = nanmean(tmp,2);

           tmp = cntlai_cn51_wetin_mat(:,year_index);
           mean_global_ann_cn51_wetin_cntlai(ann_index) = nanmean(tmp(:));
           mean_local_ann_cn51_wetin_cntlai(:,ann_index) = nanmean(tmp,2);

        end

     end 
end

% compute deltas 

for i = 1:(length(mean_global_ann_cngpp)-1)

    mean_global_ann_cngpp_delta(i) = nanmean(mean_global_ann_cngpp{i+1}(:)) - nanmean(mean_global_ann_cngpp{i}(:));
    mean_global_ann_cntlai_delta(i) = nanmean(mean_global_ann_cntlai{i+1}(:)) - nanmean(mean_global_ann_cntlai{i}(:));
    mean_global_ann_cnnee_delta(i) = nanmean(mean_global_ann_cnnee{i+1}(:)) - nanmean(mean_global_ann_cnnee{i}(:));
    mean_global_ann_cncloss_delta(i) = nanmean(mean_global_ann_cncloss{i+1}(:)) - nanmean(mean_global_ann_cncloss{i}(:));

end

figure
plot(1:(length(mean_global_ann_cngpp)-1),mean_global_ann_cngpp_delta,'LineWidth',2)
grid on
xlabel('cycles completed [-]')
ylabel('\Delta avg. GPP [kg m^{-2} s^{-1}]')
title('change in cycle-average GPP [kg m^{-2} s^{-1}]')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/Delta_CNGPP'],'png')

figure
plot(1:(length(mean_global_ann_cntlai)-1),mean_global_ann_cntlai_delta,'LineWidth',2)
grid on
xlabel('cycles completed [-]')
ylabel('\Delta avg. LAI [-]')
title('change in cycle-average LAI [-]')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/Delta_CNTLAI'],'png')

figure
plot(1:(length(mean_global_ann_cnnee)-1),mean_global_ann_cnnee_delta,'LineWidth',2)
grid on
xlabel('cycles completed [-]')
ylabel('\Delta avg. NEE [kg m^{-2} s^{-1}]')
title('change in cycle-average NEE [kg m^{-2} s^{-1}]')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/Delta_CNNEE'],'png')

figure
plot(1:(length(mean_global_ann_cncloss)-1),mean_global_ann_cncloss_delta,'LineWidth',2)
grid on
xlabel('cycles completed [-]')
ylabel('\Delta avg. CLOSS [kg m^{-2} s^{-1}]')
title('change in cycle-average CLOSS [kg m^{-2} s^{-1}]')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/Delta_CNCLOSS'],'png')


noStress_years = 11;

%mean_global_ann_clm5_cnnee(26) = NaN;

legend_str = {};

figure
for i = 1:nr_of_cycles
    if (i < 8)
        LS = '-';
    else
        LS = '--';
    end
    plot(start_year:stop_year,mean_global_ann_cngpp{i},'LineWidth',2,'LineStyle',LS);
    hold on
    legend_str{i} =  ['CN51 cycle ' num2str(i) ' (@280ppm)'];
end
plot(start_year:stop_year,mean_global_ann_fc_gpp,'LineWidth',2,'LineStyle','-.');
hold on
plot(start_year:stop_year,mean_global_ann_fs_gpp,'LineWidth',2,'LineStyle','-.');
hold on
plot(start_year:stop_year,mean_global_ann_maiac_gpp,'LineWidth',2,'LineStyle','-.');
hold on
plot(start_year:stop_year,mean_global_ann_vpm_gpp,'LineWidth',2,'LineStyle','-.');
hold on 
plot(start_year:stop_year,mean_global_ann_cn40_cngpp,'LineWidth',2,'LineStyle','-','Color','k');
hold on
plot(start_year:stop_year,mean_global_ann_cn45_cngpp,'LineWidth',2,'LineStyle','--','Color','k');
hold on
%plot(start_year:stop_year,mean_global_ann_clm5_cngpp,'LineWidth',2,'LineStyle','--','Color','m');
%hold on
plot(start_year:stop_year,mean_global_ann_clm5_of_pi_cngpp,'LineWidth',2,'LineStyle','-','Color','b');
hold on
plot(start_year:stop_year,mean_global_ann_clm5_of_mod_cngpp,'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{i+1} = 'FluxCom';
legend_str{i+2} = 'FluxSat';
legend_str{i+3} = 'MAIAC';
legend_str{i+4} = 'VPM';
legend_str{i+5} = 'CN40 (@390ppm)';
legend_str{i+6} = 'CN45 (@390ppm)';
%legend_str{i+7} = 'CLM5 Jana';
legend_str{i+7} = 'CLM5 NCAR (1850 - 1890)';
legend_str{i+8} = 'CLM5 NCAR (1980 - 2014)';
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_cycles_280'],'png')

legend_str = {};

figure
for i = 1:nr_of_cycles
    if (i < 8)
        LS = '-';
    else
        LS = '--';
    end
    plot(start_year:stop_year,mean_global_ann_cngpp{i},'LineWidth',2,'LineStyle',LS);
    hold on
    legend_str{i} =  ['CN51 cycle ' num2str(i) ' (@280ppm)'];
end
plot(start_year:(start_year+noStress_years-1),mean_global_ann_cn51_noStress_cngpp(1:noStress_years),'LineWidth',2,'LineStyle','-.');
hold on
plot(start_year:(start_year+noStress_years-1),mean_global_ann_cn51_noStress_efpot_cngpp(1:noStress_years),'LineWidth',2,'LineStyle','-.');
hold on
plot(start_year:stop_year,mean_global_ann_cn40_cngpp,'LineWidth',2,'LineStyle','-','Color','k');
hold on
plot(start_year:stop_year,mean_global_ann_cn45_cngpp,'LineWidth',2,'LineStyle','--','Color','k');
hold on
%plot(start_year:stop_year,mean_global_ann_clm5_cngpp,'LineWidth',2,'LineStyle','--','Color','m');
%hold on
plot(start_year:stop_year,mean_global_ann_clm5_of_pi_cngpp,'LineWidth',2,'LineStyle','-','Color','b');
hold on
plot(start_year:stop_year,mean_global_ann_clm5_of_mod_cngpp,'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{i+1} = 'CN51 (@280ppm, no stress RH)';
legend_str{i+2} = 'CN51 (@280ppm, no stress pot. EF)';
legend_str{i+3} = 'CN40 (@390ppm)';
legend_str{i+4} = 'CN45 (@390ppm)';
legend_str{i+5} = 'CLM5 NCAR (1850 - 1890)';
legend_str{i+6} = 'CLM5 NCAR (1980 - 2014)';
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_cycles_280_noStress'],'png')

legend_str = {};

figure
for i = 1:nr_of_cycles
    if (i < 8)
        LS = '-';
    else
        LS = '--';
    end
    plot(start_year:stop_year,mean_global_ann_cngpp{i},'LineWidth',2,'LineStyle',LS);
    hold on
    legend_str{i} =  ['CN51 cycle ' num2str(i) ' (@280ppm)'];
end
plot(start_year:stop_year,mean_global_ann_cn51_sat_cngpp,'LineWidth',2,'LineStyle','-.');
hold on
plot(start_year:stop_year,mean_global_ann_cn40_cngpp,'LineWidth',2,'LineStyle','-','Color','k');
hold on
plot(start_year:stop_year,mean_global_ann_cn45_cngpp,'LineWidth',2,'LineStyle','--','Color','k');
hold on
%plot(start_year:stop_year,mean_global_ann_clm5_cngpp,'LineWidth',2,'LineStyle','--','Color','m');
%hold on
plot(start_year:stop_year,mean_global_ann_clm5_of_pi_cngpp,'LineWidth',2,'LineStyle','-','Color','b');
hold on
plot(start_year:stop_year,mean_global_ann_clm5_of_mod_cngpp,'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{i+1} = 'CN51 (@280ppm, saturation)';
legend_str{i+2} = 'CN40 (@390ppm)';
legend_str{i+3} = 'CN45 (@390ppm)';
legend_str{i+4} = 'CLM5 NCAR (1850 - 1890)';
legend_str{i+5} = 'CLM5 NCAR (1980 - 2014)';
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_cycles_280_sat'],'png')

legend_str = {};

figure
for i = 1:nr_of_cycles
    if (i < 8) 
        LS = '-'; 
    else 
        LS = '--';
    end  
    plot(start_year:stop_year,mean_global_ann_cngpp{i},'LineWidth',2,'LineStyle',LS);
    hold on
    legend_str{i} =  ['CN51 cycle ' num2str(i) ' (@280ppm)'];
end
plot(start_year:stop_year,mean_global_ann_cn51_xroot_cngpp,'LineWidth',2,'LineStyle','-.');
hold on
plot(start_year:stop_year,mean_global_ann_cn40_cngpp,'LineWidth',2,'LineStyle','-','Color','k');
hold on
plot(start_year:stop_year,mean_global_ann_cn45_cngpp,'LineWidth',2,'LineStyle','--','Color','k');
hold on
%plot(start_year:stop_year,mean_global_ann_clm5_cngpp,'LineWidth',2,'LineStyle','--','Color','m');
%hold on
plot(start_year:stop_year,mean_global_ann_clm5_of_pi_cngpp,'LineWidth',2,'LineStyle','-','Color','b');
hold on
plot(start_year:stop_year,mean_global_ann_clm5_of_mod_cngpp,'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{i+1} = 'CN51 (@280ppm, xroo=-25000MPa)';
legend_str{i+2} = 'CN40 (@390ppm)';
legend_str{i+3} = 'CN45 (@390ppm)';
legend_str{i+4} = 'CLM5 NCAR (1850 - 1890)';
legend_str{i+5} = 'CLM5 NCAR (1980 - 2014)';
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_cycles_280_xroot'],'png')

legend_str = {};

figure
for i = 1:nr_of_cycles
    if (i < 8) 
        LS = '-'; 
    else 
        LS = '--';
    end  
    plot(start_year:stop_year,mean_global_ann_cngpp{i},'LineWidth',2,'LineStyle',LS);
    hold on
    legend_str{i} =  ['CN51 cycle ' num2str(i) ' (@280ppm)'];
end
plot(start_year:stop_year,mean_global_ann_cn51_wetin_cngpp,'LineWidth',2,'LineStyle','-.');
hold on
plot(start_year:stop_year,mean_global_ann_cn40_cngpp,'LineWidth',2,'LineStyle','-','Color','k');
hold on
plot(start_year:stop_year,mean_global_ann_cn45_cngpp,'LineWidth',2,'LineStyle','--','Color','k');
hold on
%plot(start_year:stop_year,mean_global_ann_clm5_cngpp,'LineWidth',2,'LineStyle','--','Color','m');
%hold on
plot(start_year:stop_year,mean_global_ann_clm5_of_pi_cngpp,'LineWidth',2,'LineStyle','-','Color','b');
hold on
plot(start_year:stop_year,mean_global_ann_clm5_of_mod_cngpp,'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{i+1} = 'CN51 (@280ppm, wetness per zone)';
legend_str{i+2} = 'CN40 (@390ppm)';
legend_str{i+3} = 'CN45 (@390ppm)';
legend_str{i+4} = 'CLM5 NCAR (1850 - 1890)';
legend_str{i+5} = 'CLM5 NCAR (1980 - 2014)';
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_cycles_280_wetin'],'png')

%figure
%for i = 1:nr_of_transient_cycles
%    if (i < 8) 
%        LS = '-'; 
%    else 
%        LS = '--';
%    end  
%    plot(start_year:stop_year,mean_global_ann_transient_cngpp{i},'LineWidth',2,'LineStyle',LS);
%    hold on
%    legend_str{i} =  ['CN51 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_fc_gpp,'LineWidth',2,'LineStyle','-.');
%hold on
%plot(start_year:stop_year,mean_global_ann_fs_gpp,'LineWidth',2,'LineStyle','-.');
%hold on
%plot(start_year:stop_year,mean_global_ann_maiac_gpp,'LineWidth',2,'LineStyle','-.');
%hold on
%plot(start_year:stop_year,mean_global_ann_vpm_gpp,'LineWidth',2,'LineStyle','-.');
%hold on 
%plot(start_year:stop_year,mean_global_ann_cn40_cngpp,'LineWidth',2,'LineStyle','-','Color','k');
%hold on
%plot(start_year:stop_year,mean_global_ann_cn45_cngpp,'LineWidth',2,'LineStyle','--','Color','k');
%%hold on
%%plot(start_year:stop_year,mean_global_ann_clm5_cngpp,'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{i+1} = 'FluxCom';
%legend_str{i+2} = 'FluxSat';
%legend_str{i+3} = 'MAIAC';
%legend_str{i+4} = 'VPM';
%legend_str{i+5} = 'CN40 (@390ppm)';
%legend_str{i+6} = 'CN45 (@390ppm)';
%%legend_str{i+7} = 'CLM5';
%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_cycles_transient'],'png')

legend_str = {};

figure
for i = 1:nr_of_cycles
    if (i < 8)
        LS = '-';
    else
        LS = '--';
    end
    plot(start_year:stop_year,mean_global_ann_cnnee{i},'LineWidth',2,'LineStyle',LS);
    hold on
    legend_str{i} =  ['CN51 cycle ' num2str(i) ' (@280ppm)'];
end
hold on
plot(start_year:stop_year,mean_global_ann_cn40_cnnee,'LineWidth',2,'LineStyle','-','Color','k');
hold on
plot(start_year:stop_year,mean_global_ann_cn45_cnnee,'LineWidth',2,'LineStyle','--','Color','k');
%hold on
%plot(start_year:stop_year,mean_global_ann_clm5_cnnee,'LineWidth',2,'LineStyle','--','Color','m');
hold on
plot(start_year:stop_year,mean_global_ann_clm5_of_pi_cnnee,'LineWidth',2,'LineStyle','-','Color','b');
hold on
plot(start_year:stop_year,mean_global_ann_clm5_of_mod_cnnee,'LineWidth',2,'LineStyle','-','Color','m');
hold off
grid on
xlabel('year')
ylabel('NEE [kg m^{-2} s^{-1}]')
title('annual average NEE [kg m^{-2} s^{-1}]')
legend_str{i+1} = 'CN40 (@390ppm)';
legend_str{i+2} = 'CN45 (@390ppm)';
%legend_str{i+3} = 'CLM5';
legend_str{i+3} = 'CLM5 NCAR (1850 - 1890)';
legend_str{i+4} = 'CLM5 NCAR (1980 - 2014)';
legend(legend_str,'Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/NEE_ann_CN51_cycles_280'],'png')

legend_str = {};

figure
for i = 1:nr_of_cycles
    if (i < 8)
        LS = '-';
    else
        LS = '--';
    end
    plot(start_year:stop_year,mean_global_ann_cnnee{i},'LineWidth',2,'LineStyle',LS);
    hold on
    legend_str{i} =  ['CN51 cycle ' num2str(i) ' (@280ppm)'];
end
plot(start_year:(start_year+noStress_years-1),mean_global_ann_cn51_noStress_cnnee(1:noStress_years),'LineWidth',2,'LineStyle','-.');
hold on
plot(start_year:(start_year+noStress_years-1),mean_global_ann_cn51_noStress_efpot_cnnee(1:noStress_years),'LineWidth',2,'LineStyle','-.');
hold on
plot(start_year:stop_year,mean_global_ann_cn40_cnnee,'LineWidth',2,'LineStyle','-','Color','k');
hold on
plot(start_year:stop_year,mean_global_ann_cn45_cnnee,'LineWidth',2,'LineStyle','--','Color','k');
%hold on
%plot(start_year:stop_year,mean_global_ann_clm5_cnnee,'LineWidth',2,'LineStyle','--','Color','m');
hold on
plot(start_year:stop_year,mean_global_ann_clm5_of_pi_cnnee,'LineWidth',2,'LineStyle','-','Color','b');
hold on
plot(start_year:stop_year,mean_global_ann_clm5_of_mod_cnnee,'LineWidth',2,'LineStyle','-','Color','m');
hold off
grid on
xlabel('year')
ylabel('NEE [kg m^{-2} s^{-1}]')
title('annual average NEE [kg m^{-2} s^{-1}]')
legend_str{i+1} = 'CN51 (@280ppm, no stress RH)';
legend_str{i+2} = 'CN51 (@280ppm, no stress pot. EF)';
legend_str{i+3} = 'CN40 (@390ppm)';
legend_str{i+4} = 'CN45 (@390ppm)';
%legend_str{i+3} = 'CLM5';
legend_str{i+5} = 'CLM5 NCAR (1850 - 1890)';
legend_str{i+6} = 'CLM5 NCAR (1980 - 2014)';
legend(legend_str,'Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/NEE_ann_CN51_cycles_280_noStress'],'png')


legend_str = {};

figure
for i = 1:nr_of_cycles
    if (i < 8)
        LS = '-';
    else
        LS = '--';
    end 
    plot(start_year:stop_year,mean_global_ann_cncloss{i},'LineWidth',2,'LineStyle',LS);
    hold on
    legend_str{i} =  ['CN51 cycle ' num2str(i) ' (@280ppm)'];
end
hold on
plot(start_year:stop_year,mean_global_ann_cn40_cncloss,'LineWidth',2,'LineStyle','-','Color','k');
hold on
plot(start_year:stop_year,mean_global_ann_cn45_cncloss,'LineWidth',2,'LineStyle','--','Color','k');
hold on
plot(start_year:stop_year,mean_global_ann_clm5_cncloss,'LineWidth',2,'LineStyle','--','Color','m');
hold off 
grid on
xlabel('year')
ylabel('CLOSS [kg m^{-2} s^{-1}]')
title('annual average CLOSS [kg m^{-2} s^{-1}]')
legend_str{i+1} = 'CN40 (@390ppm)';
legend_str{i+2} = 'CN45 (@390ppm)';
legend_str{i+3} = 'CLM5';
legend(legend_str,'Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/CLOSS_ann_CN51_cycles_280'],'png')

legend_str = {};

figure
for i = 1:nr_of_cycles
    if (i < 8)
        LS = '-';
    else
        LS = '--';
    end 
    plot(start_year:stop_year,mean_global_ann_cntlai{i},'LineWidth',2,'LineStyle',LS);
    hold on
    legend_str{i} =  ['CN51 cycle ' num2str(i) ' (@280ppm)'];
end
hold on
plot(start_year:stop_year,mean_global_ann_cn40_cntlai,'LineWidth',2,'LineStyle','-','Color','k');
hold on
plot(start_year:stop_year,mean_global_ann_cn45_cntlai,'LineWidth',2,'LineStyle','--','Color','k');
%hold on
%plot(start_year:stop_year,mean_global_ann_clm5_cntlai,'LineWidth',2,'LineStyle','--','Color','m');
hold on 
plot(start_year:stop_year,mean_global_ann_clm5_of_pi_cntlai,'LineWidth',2,'LineStyle','-','Color','b');
hold on
plot(start_year:stop_year,mean_global_ann_clm5_of_mod_cntlai,'LineWidth',2,'LineStyle','-','Color','m');
hold off 
grid on
xlabel('year')
ylabel('LAI [-]')
title('annual average LAI [-]')
legend_str{i+1} = 'CN40 (@390ppm)';
legend_str{i+2} = 'CN45 (@390ppm)';
%legend_str{i+3} = 'CLM5';
legend_str{i+3} = 'CLM5 NCAR (1850 - 1890)';
legend_str{i+4} = 'CLM5 NCAR (1980 - 2014)';
legend(legend_str,'Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/CNTLAI_ann_CN51_cycles_280'],'png')

legend_str = {};

figure
for i = 1:nr_of_cycles
    if (i < 8)
        LS = '-';
    else
        LS = '--';
    end
    plot(start_year:stop_year,mean_global_ann_cntlai{i},'LineWidth',2,'LineStyle',LS);
    hold on
    legend_str{i} =  ['CN51 cycle ' num2str(i) ' (@280ppm)'];
end
plot(start_year:(start_year+noStress_years-1),mean_global_ann_cn51_noStress_cntlai(1:noStress_years),'LineWidth',2,'LineStyle','-.');
hold on
plot(start_year:(start_year+noStress_years-1),mean_global_ann_cn51_noStress_efpot_cntlai(1:noStress_years),'LineWidth',2,'LineStyle','-.');
hold on
plot(start_year:stop_year,mean_global_ann_cn40_cntlai,'LineWidth',2,'LineStyle','-','Color','k');
hold on
plot(start_year:stop_year,mean_global_ann_cn45_cntlai,'LineWidth',2,'LineStyle','--','Color','k');
%hold on
%plot(start_year:stop_year,mean_global_ann_clm5_cntlai,'LineWidth',2,'LineStyle','--','Color','m');
hold on
plot(start_year:stop_year,mean_global_ann_clm5_of_pi_cntlai,'LineWidth',2,'LineStyle','-','Color','b');
hold on
plot(start_year:stop_year,mean_global_ann_clm5_of_mod_cntlai,'LineWidth',2,'LineStyle','-','Color','m');
hold off
grid on
xlabel('year')
ylabel('LAI [-]')
title('annual average LAI [-]')
legend_str{i+1} = 'CN51 (@280ppm, no stress RH)';
legend_str{i+2} = 'CN51 (@280ppm, no stress pot. EF)';
legend_str{i+3} = 'CN40 (@390ppm)';
legend_str{i+4} = 'CN45 (@390ppm)';
%legend_str{i+3} = 'CLM5';
legend_str{i+5} = 'CLM5 NCAR (1850 - 1890)';
legend_str{i+6} = 'CLM5 NCAR (1980 - 2014)';
legend(legend_str,'Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/CNTLAI_ann_CN51_cycles_280_noStress'],'png')


legend_str = {};

figure
for i = 1:nr_of_cycles
    if (i < 8)
        LS = '-';
    else
        LS = '--';
    end
    plot(start_year:stop_year,mean_global_ann_cntotc{i},'LineWidth',2,'LineStyle',LS);
    hold on
    legend_str{i} =  ['CN51 cycle ' num2str(i) ' (@280ppm)'];
end
hold on
plot(start_year:stop_year,mean_global_ann_cn40_cntotc,'LineWidth',2,'LineStyle','-','Color','k');
hold on
plot(start_year:stop_year,mean_global_ann_cn45_cntotc,'LineWidth',2,'LineStyle','--','Color','k');
hold on
plot(start_year:stop_year,mean_global_ann_clm5_cntotc,'LineWidth',2,'LineStyle','--','Color','m');
hold off
grid on
xlabel('year')
ylabel('TOTC [kgC m^{-2}]')
title('annual average TOTC [kgC m^{-2}]')
legend_str{i+1} = 'CN40 (@390ppm)';
legend_str{i+2} = 'CN45 (@390ppm)';
legend_str{i+3} = 'CLM5';
legend(legend_str,'Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/CNTOTC_ann_CN51_cycles_280'],'png')

for i = 1:length(lat_vec)

figure
LS = '-';
plot(start_year:stop_year,mean_local_ann_cntotc{1}(i,:),'LineWidth',2,'LineStyle',LS);
hold on
legend_str{1} =  ['CN51 cycle 1'];
hold on
plot(start_year:stop_year,mean_local_ann_cn40_cntotc(i,:),'LineWidth',2,'LineStyle','-','Color','k');
hold on
plot(start_year:stop_year,mean_local_ann_cn45_cntotc(i,:),'LineWidth',2,'LineStyle','--','Color','k');
hold on
plot(start_year:stop_year,mean_local_ann_clm5_cntotc(i,:),'LineWidth',2,'LineStyle','--','Color','m');
hold off
grid on
xlabel('year')
ylabel('TOTC [kgC m^{-2}]')
title('annual average TOTC [kgC m^{-2}] @280ppm')
legend_str{2} = 'CN40';
legend_str{3} = 'CN45';
legend_str{4} = 'CLM5';
legend(legend_str,'Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/CNTOTC_ann_CN51_cycles_280_loc' num2str(i)],'png')

figure
LS = '-';
plot(start_year:stop_year,mean_local_ann_cntlai{1}(i,:),'LineWidth',2,'LineStyle',LS);
hold on
legend_str{1} =  ['CN51 cycle 1'];
hold on
plot(start_year:stop_year,mean_local_ann_cn40_cntlai(i,:),'LineWidth',2,'LineStyle','-','Color','k');
hold on
plot(start_year:stop_year,mean_local_ann_cn45_cntlai(i,:),'LineWidth',2,'LineStyle','--','Color','k');
hold on
plot(start_year:stop_year,mean_local_ann_clm5_cntlai(i,:),'LineWidth',2,'LineStyle','--','Color','m');
hold off
grid on
xlabel('year')
ylabel('LAI [-]')
title('annual average LAI [-] @280ppm')
legend_str{2} = 'CN40';
legend_str{3} = 'CN45';
legend_str{4} = 'CLM5';
legend(legend_str,'Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/CNTLAI_ann_CN51_cycles_280_loc' num2str(i)],'png')
end


legend_str = {};

figure
for i = 1:nr_of_cycles
    if (i < 8)
        LS = '-';
    else
        LS = '--';
    end
    plot(start_year:stop_year,mean_global_ann_et{i},'LineWidth',2,'LineStyle',LS);
    hold on
    legend_str{i} =  ['CN51 cycle ' num2str(i)];
end
hold on
plot(start_year:stop_year,mean_global_ann_cn40_et,'LineWidth',2,'LineStyle','-','Color','k');
hold on
plot(start_year:stop_year,mean_global_ann_cn45_et,'LineWidth',2,'LineStyle','--','Color','k');
hold off
grid on
xlabel('year')
ylabel('total ET [W m^{-2}]')
title('annual average total ET [W m^{-2}] @280ppm')
legend_str{i+1} = 'CN40';
legend_str{i+2} = 'CN45';
legend(legend_str,'Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/ET_ann_CN51_cycles_280'],'png')

legend_str = {};

figure   
for i = 1:nr_of_cycles
    if (i < 8)
        LS = '-';
    else
        LS = '--';
    end
    plot(start_year:stop_year,mean_global_ann_evpsoil{i},'LineWidth',2,'LineStyle',LS);
    hold on
    legend_str{i} =  ['CN51 cycle ' num2str(i)];
end
hold on
plot(start_year:stop_year,mean_global_ann_cn40_evpsoil,'LineWidth',2,'LineStyle','-','Color','k');
hold on
plot(start_year:stop_year,mean_global_ann_cn45_evpsoil,'LineWidth',2,'LineStyle','--','Color','k');
hold off
grid on
xlabel('year')
ylabel('EVPSOIL [W m^{-2}]')
title('annual average EVPSOIL [W m^{-2}] @280ppm')
legend_str{i+1} = 'CN40';
legend_str{i+2} = 'CN45';
legend(legend_str,'Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/EVPSOIL_ann_CN51_cycles_280'],'png')

legend_str = {};

figure
for i = 1:nr_of_cycles
    if (i < 8)
        LS = '-';
    else
        LS = '--';
    end
    plot(start_year:stop_year,mean_global_ann_evptrns{i},'LineWidth',2,'LineStyle',LS);
    hold on
    legend_str{i} =  ['CN51 cycle ' num2str(i)];
end
hold on
plot(start_year:stop_year,mean_global_ann_cn40_evptrns,'LineWidth',2,'LineStyle','-','Color','k');
hold on
plot(start_year:stop_year,mean_global_ann_cn45_evptrns,'LineWidth',2,'LineStyle','--','Color','k');
hold off
grid on
xlabel('year')
ylabel('EVPTRNS [W m^{-2}]')
title('annual average EVPTRNS [W m^{-2}] @280ppm')
legend_str{i+1} = 'CN40';
legend_str{i+2} = 'CN45';
legend(legend_str,'Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/EVPTRNS_ann_CN51_cycles_280'],'png')

legend_str = {};

figure
for i = 1:nr_of_cycles
    if (i < 8)
        LS = '-';
    else
        LS = '--';
    end
    plot(start_year:stop_year,mean_global_ann_evpintr{i},'LineWidth',2,'LineStyle',LS);
    hold on
    legend_str{i} =  ['CN51 cycle ' num2str(i)];
end
hold on
plot(start_year:stop_year,mean_global_ann_cn40_evpintr,'LineWidth',2,'LineStyle','-','Color','k');
hold on
plot(start_year:stop_year,mean_global_ann_cn45_evpintr,'LineWidth',2,'LineStyle','--','Color','k');
hold off
grid on
xlabel('year')
ylabel('EVPINTR [W m^{-2}]')
title('annual average EVPINTR [W m^{-2}] @280ppm')
legend_str{i+1} = 'CN40';
legend_str{i+2} = 'CN45';
legend(legend_str,'Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/EVPINTR_ann_CN51_cycles_280'],'png')

legend_str = {};

figure
for i = 1:nr_of_cycles
    if (i < 8)
        LS = '-';
    else
        LS = '--';
    end
    plot(start_year:stop_year,mean_global_ann_evpsbln{i},'LineWidth',2,'LineStyle',LS);
    hold on
    legend_str{i} =  ['CN51 cycle ' num2str(i)];
end
hold on
plot(start_year:stop_year,mean_global_ann_cn40_evpsbln,'LineWidth',2,'LineStyle','-','Color','k');
hold on
plot(start_year:stop_year,mean_global_ann_cn45_evpsbln,'LineWidth',2,'LineStyle','--','Color','k');
hold off
grid on
xlabel('year')
ylabel('EVPSBLN [W m^{-2}]')
title('annual average EVPSBLN [W m^{-2}] @280ppm')
legend_str{i+1} = 'CN40';
legend_str{i+2} = 'CN45';
legend(legend_str,'Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/EVPSBLN_ann_CN51_cycles_280'],'png')

legend_str = {};  
                  
figure            
for i = 1:nr_of_cycles
    if (i < 8)
        LS = '-';
    else           
        LS = '--'; 
    end                
    plot(start_year:stop_year,mean_srfexc_mat{i},'LineWidth',2,'LineStyle',LS);
    hold on     
    legend_str{i} =  ['CN51 cycle ' num2str(i)];
end             
hold off           
grid on            
xlabel('year')     
ylabel('SRFEXC [kg m^{-2}]')
title('Jan-01 SRFEXC [kg m^{-2}] @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/SRFEXC_ann_CN51_cycles_280'],'png')

legend_str = {};

figure
for i = 1:nr_of_cycles
    if (i < 8)
        LS = '-';
    else
        LS = '--';
    end
    plot(start_year:stop_year,mean_rzexc_mat{i},'LineWidth',2,'LineStyle',LS);
    hold on
    legend_str{i} =  ['CN51 cycle ' num2str(i)];
end
hold off
grid on
xlabel('year')
ylabel('RZEXC [kg m^{-2}]')
title('Jan-01 RZEXC [kg m^{-2}] @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/RZEXC_ann_CN51_cycles_280'],'png')

legend_str = {};

figure
for i = 1:nr_of_cycles
    if (i < 8)
        LS = '-';
    else
        LS = '--';
    end
    plot(start_year:stop_year,mean_catdef_mat{i},'LineWidth',2,'LineStyle',LS);
    hold on
    legend_str{i} =  ['CN51 cycle ' num2str(i)];
end
hold off
grid on
xlabel('year')
ylabel('CATDEF [kg m^{-2}]')
title('Jan-01 CATDEF [kg m^{-2}] @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/CATDEF_ann_CN51_cycles_280'],'png')

% plot PFT-specific plots

for p = 1:length(cn51_pft_name)

    cn51_pft_index = find(cn51_dom_pft==p);

    tmp = mean_local_ann_fc_gpp(cn51_pft_index,:);
    tmp_fc = nanmean(tmp,1);
    tmp = mean_local_ann_fs_gpp(cn51_pft_index,:);
    tmp_fs = nanmean(tmp,1);
    tmp = mean_local_ann_maiac_gpp(cn51_pft_index,:);
    tmp_maiac = nanmean(tmp,1);
    tmp = mean_local_ann_vpm_gpp(cn51_pft_index,:);
    tmp_vpm = nanmean(tmp,1);
    tmp = mean_local_ann_cn40_cngpp(cn51_pft_index,:);
    tmp_cn40 = nanmean(tmp,1);

    tmp = mean_local_ann_cn40_cntlai(cn51_pft_index,:);
    tmp_cn40_lai = nanmean(tmp,1);


    legend_str = {};

    figure
    for i = 1:nr_of_cycles
        tmp = mean_local_ann_cngpp{i}(cn51_pft_index,:);
        tmp_cn51 = nanmean(tmp,1);
        if (i < 8)
            LS = '-';
        else
            LS = '--';
        end
        plot(start_year:stop_year,tmp_cn51,'LineWidth',2,'LineStyle',LS);
        hold on
        legend_str{i} =  ['CN51 cycle ' num2str(i)];
     end
     plot(start_year:stop_year,tmp_fc,'LineWidth',2,'LineStyle','-.');
     hold on
     plot(start_year:stop_year,tmp_fs,'LineWidth',2,'LineStyle','-.');
     hold on
     plot(start_year:stop_year,tmp_maiac,'LineWidth',2,'LineStyle','-.');
     hold on
     plot(start_year:stop_year,tmp_vpm,'LineWidth',2,'LineStyle','-.');
     hold on
     plot(start_year:stop_year,tmp_cn40,'LineWidth',2,'LineStyle','-','Color','k');
     hold off
     grid on
     xlabel('year')
     ylabel('GPP [kg m^{-2} s^{-1}]')
     title(['annual average GPP [kg m^{-2} s^{-1}] @280ppm for ' cn51_pft_name_short{p}])
     legend_str{i+1} = 'FluxCom';
     legend_str{i+2} = 'FluxSat';
     legend_str{i+3} = 'MAIAC';
     legend_str{i+4} = 'VPM';
     legend_str{i+5} = 'CN40';
     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_cycles_280_PFT' num2str(p)],'png')

    legend_str = {}; 

    figure
    for i = 1:nr_of_cycles
        tmp = mean_local_ann_cntlai{i}(cn51_pft_index,:);
        tmp_cn51_lai = nanmean(tmp,1);
        if (i < 8)
            LS = '-';
        else
            LS = '--';
        end 
        plot(start_year:stop_year,tmp_cn51_lai,'LineWidth',2,'LineStyle',LS);
        hold on
        legend_str{i} =  ['CN51 cycle ' num2str(i)];
     end 
     hold on
     plot(start_year:stop_year,tmp_cn40_lai,'LineWidth',2,'LineStyle','-','Color','k');
     hold off
     grid on
     xlabel('year')
     ylabel('LAI [-]')
     title(['annual average LAI [-] @280ppm for ' cn51_pft_name_short{p}])
     legend_str{i+1} = 'CN40';
     legend(legend_str,'Location','bestoutside')
     set(gcf,'paperunits','centimeters')
     set(gcf,'paperposition',[1 1 25 15])
     saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/CNTLAI_ann_CN51_cycles_280_PFT' num2str(p)],'png')
end

% identify locations with 0 GPP at all times (separately for each cycle)

for i = 1:nr_of_cycles

    tmp = mean_local_ann_cngpp{i};
    zero_mat = tmp==0;
    zero_index = sum(zero_mat,2)==size(zero_mat,2);
    zero_tiles = find(sum(zero_mat,2)==size(zero_mat,2));

    zero_index_mat{i} = zero_index;
    zero_tiles_mat{i} = zero_tiles;
    sum_zero_index{i} = sum(zero_index);
    zero_index_pft{i} = cn51_ity(zero_index,:);
    zero_index_frac{i} = cn51_fvg(zero_index,:);

    load coastlines;
    [latcells, loncells] = polysplit(coastlat, coastlon);
    lat_range = [-60 90];
    lon_range = [-180 180];
    latlim = [-60 90];
    lonlim = [-180 180];

    cmap_mat = colormap(jet(15));

    figure
    ax = axesm('miller','MapLatLimit',lat_range,'MapLonLimit',lon_range,...
               'FLatLimit',lat_range,'FLonLimit',lon_range);
    framem
    gridm on;
    setm(ax,'MLineLocation',60,'PLineLocation',60,'MeridianLabel','on',...
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',60,...
    'MLabelParallel','south','FontSize',12)
    geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
    scatterm(lat_vec(zero_index),lon_vec(zero_index),30,cmap_mat(cn51_dom_pft(zero_index),:),"filled")
    caxis([1 15]);
    set(gcf,'Colormap',cmap_mat);
    plotm(coastlat, coastlon,'k')
    set(gca,'XLim',lonlim);
    set(gca,'YLim',latlim);
    tightmap
    title(['locations with zero GPP'])
    h = colorbar;
    set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
    set(gcf,'paperunits','centimeters')
    set(gcf,'paperposition',[1 1 17 15])
    saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/locs_CN51_280_cycle' num2str(i)],'png')


end

% map of differences between first and last cycle

diff_cngpp_cfirst_clast = mean_local_ann_cngpp{nr_of_cycles} - mean_local_ann_cngpp{1};
mean_diff_cngpp_cfirst_clast = nanmean(diff_cngpp_cfirst_clast,2);

load coastlines;
[latcells, loncells] = polysplit(coastlat, coastlon);
lat_range = [-60 90];
lon_range = [-180 180];
latlim = [-60 90];
lonlim = [-180 180];

caxis_min = min(mean_diff_cngpp_cfirst_clast(:));
caxis_max = max(mean_diff_cngpp_cfirst_clast(:));

figure
ax = axesm('miller','MapLatLimit',lat_range,'MapLonLimit',lon_range,...
               'FLatLimit',lat_range,'FLonLimit',lon_range);
framem
gridm on;
setm(ax,'MLineLocation',60,'PLineLocation',60,'MeridianLabel','on',...
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',60,...
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
scatterm(lat_vec,lon_vec,30,mean_diff_cngpp_cfirst_clast,"filled")
caxis([caxis_min caxis_max]);
colormap(jet(10))
plotm(coastlat, coastlon,'k')
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['\Delta GPP c12-c1 [kg m^{-2} s^{-1}]'])
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/GPP_diff_CN51_280_cfirst_clast'],'png')

% map of differences between last cycle and CatchCN40

diff_cngpp_cn40_clast = mean_local_ann_cngpp{nr_of_cycles} - mean_local_ann_cn40_cngpp;
mean_diff_cngpp_cn40_clast = nanmean(diff_cngpp_cn40_clast,2);

caxis_min = min(mean_diff_cngpp_cn40_clast(:));
caxis_max = max(mean_diff_cngpp_cn40_clast(:));

figure
ax = axesm('miller','MapLatLimit',lat_range,'MapLonLimit',lon_range,...
               'FLatLimit',lat_range,'FLonLimit',lon_range);
framem
gridm on;
setm(ax,'MLineLocation',60,'PLineLocation',60,'MeridianLabel','on',...
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',60,...
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
scatterm(lat_vec,lon_vec,30,mean_diff_cngpp_cn40_clast,"filled")
caxis([caxis_min caxis_max]);
colormap(jet(10))
plotm(coastlat, coastlon,'k')
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['\Delta GPP CN51 c12 - CN40 [kg m^{-2} s^{-1}]'])
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/GPP_diff_CN51_280_CN40_clast'],'png')

