clc
clear all 
close all 

%exp1_path = '/discover/nobackup/projects/geoscm/fzeng/Catchment-CN5.1/local_runs`/';
%exp1_name = 'GEOSldas_CN51_local_280';

exp1_path = '/discover/nobackup/jkolassa/CN51_test_local/';
exp1_name = 'GEOSldas_CN51_local_280';

exp2_path = '/discover/nobackup/jkolassa/new/';
exp2_name = 'GEOSldas_CN51_local_280';

exp3_path = '/discover/nobackup/jkolassa/CN51_test_local/';
exp3_name = 'GEOSldas_CN51_local_280_fcur';

exp4_path = '/discover/nobackup/jkolassa/new/';
exp4_name = 'GEOSldas_CN51_280_global';

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/CN51_test_local/GEOSldas_CN51_local/output/SMAP_EASEv2_M36/cat/ens0000/Y1980/M01/GEOSldas_CN51_local.tavg24_1d_lnd_Nt.monthly.198001.nc4'];
lat_vec = double(ncread(cn51_read_file,'lat'));
lon_vec = double(ncread(cn51_read_file,'lon'));

exp2_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'];
exp2_lat_vec = double(ncread(exp2_read_file,'lat'));
exp2_lon_vec = double(ncread(exp2_read_file,'lon'));

exp2_rst_file = ['/discover/nobackup/jkolassa/new/GEOSldas_CN51_local_280/output/SMAP_EASEv2_M36/rs/ens0000/Y1980/M01/GEOSldas_CN51_local_280.catchcnclm51_internal_rst.19800101_0000'];
exp2_ityp = double(ncread(exp2_rst_file,'ITY'));
exp2_fveg = double(ncread(exp2_rst_file,'FVG'));

%exp3_lat_vec = exp2_lat_vec;
%exp3_lon_vec = exp2_lon_vec;

exp3_read_file = ['/discover/nobackup/jkolassa/CN51_test_local/GEOSldas_CN51_local_280_fcur/output/SMAP_EASEv2_M36/cat/ens0000/Y1980/M01/GEOSldas_CN51_local_280_fcur.tavg24_1d_lnd_Nt.monthly.198001.nc4'];
exp3_lat_vec = double(ncread(exp3_read_file,'lat'));
exp3_lon_vec = double(ncread(exp3_read_file,'lon'));

exp4_read_file = ['/discover/nobackup/projects/geoscm/fzeng/Catchment-CN5.1/GEOSldas_CN51/output/SMAP_EASEv2_M36/cat/ens0000/Y1980/M01/GEOSldas_CN51.tavg24_1d_lnd_Nt.monthly.198001.nc4'];
exp4_lat_vec = double(ncread(exp4_read_file,'lat'));
exp4_lon_vec = double(ncread(exp4_read_file,'lon'));

for i = 1:length(exp2_lat_vec)

    dist_vec = sqrt((exp4_lat_vec-exp2_lat_vec(i)).^2 + (exp4_lon_vec-exp2_lon_vec(i)).^2);
    exp4_loc_index(i) = find(dist_vec==min(dist_vec));
end

    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

    for m = 1:12

        mat_index = mat_index + 1;

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

        % define file names

        exp1_read_file = [exp1_path exp1_name '/output/SMAP_EASEv2_M36/cat/ens0000/Y' y_str ...
                         '/M' m_str '/' exp1_name '.tavg24_1d_lnd_Nt.monthly.' y_str m_str '.nc4'];

        exp2_read_file = [exp2_path exp2_name '/output/SMAP_EASEv2_M36/cat/ens0000/Y' y_str ...
                         '/M' m_str '/' exp2_name '.tavg24_1d_lnd_Nt.monthly.' y_str m_str '.nc4'];

        exp3_read_file = [exp3_path exp3_name '/output/SMAP_EASEv2_M36/cat/ens0000/Y' y_str ...
                         '/M' m_str '/' exp3_name '.tavg24_1d_lnd_Nt.monthly.' y_str m_str '.nc4'];

        exp4_read_file = [exp4_path exp4_name '/output/SMAP_EASEv2_M36/cat/ens0000/Y' y_str ...
                         '/M' m_str '/' exp4_name '.tavg24_1d_lnd_Nt.monthly.' y_str m_str '.nc4'];


        if (mat_index==1)
           exp1_cngpp_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           exp1_cnnee_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           exp1_cncloss_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           exp1_cntlai_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           exp1_cntotc_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           exp1_et_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           exp1_evpsoil_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           exp1_evptrns_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           exp1_evpintr_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           exp1_evpsbln_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);

           exp2_cngpp_mat = NaN*ones(length(exp2_lat_vec),(stop_year-start_year+1)*12);
           exp2_cnnee_mat = NaN*ones(length(exp2_lat_vec),(stop_year-start_year+1)*12);
           exp2_cncloss_mat = NaN*ones(length(exp2_lat_vec),(stop_year-start_year+1)*12);
           exp2_cntlai_mat = NaN*ones(length(exp2_lat_vec),(stop_year-start_year+1)*12);
           exp2_cntotc_mat = NaN*ones(length(exp2_lat_vec),(stop_year-start_year+1)*12);
           exp2_et_mat = NaN*ones(length(exp2_lat_vec),(stop_year-start_year+1)*12);
           exp2_evpsoil_mat = NaN*ones(length(exp2_lat_vec),(stop_year-start_year+1)*12);
           exp2_evptrns_mat = NaN*ones(length(exp2_lat_vec),(stop_year-start_year+1)*12);
           exp2_evpintr_mat = NaN*ones(length(exp2_lat_vec),(stop_year-start_year+1)*12);
           exp2_evpsbln_mat = NaN*ones(length(exp2_lat_vec),(stop_year-start_year+1)*12);

           exp3_cngpp_mat = NaN*ones(length(exp3_lat_vec),(stop_year-start_year+1)*12);
           exp3_cnnee_mat = NaN*ones(length(exp3_lat_vec),(stop_year-start_year+1)*12);
           exp3_cncloss_mat = NaN*ones(length(exp3_lat_vec),(stop_year-start_year+1)*12);
           exp3_cntlai_mat = NaN*ones(length(exp3_lat_vec),(stop_year-start_year+1)*12);
           exp3_cntotc_mat = NaN*ones(length(exp3_lat_vec),(stop_year-start_year+1)*12);
           exp3_et_mat = NaN*ones(length(exp3_lat_vec),(stop_year-start_year+1)*12);
           exp3_evpsoil_mat = NaN*ones(length(exp3_lat_vec),(stop_year-start_year+1)*12);
           exp3_evptrns_mat = NaN*ones(length(exp3_lat_vec),(stop_year-start_year+1)*12);
           exp3_evpintr_mat = NaN*ones(length(exp3_lat_vec),(stop_year-start_year+1)*12);
           exp3_evpsbln_mat = NaN*ones(length(exp3_lat_vec),(stop_year-start_year+1)*12);

           exp4_cngpp_mat = NaN*ones(length(exp2_lat_vec),(stop_year-start_year+1)*12);
           exp4_cnnee_mat = NaN*ones(length(exp2_lat_vec),(stop_year-start_year+1)*12);
           exp4_cncloss_mat = NaN*ones(length(exp2_lat_vec),(stop_year-start_year+1)*12);
           exp4_cntlai_mat = NaN*ones(length(exp2_lat_vec),(stop_year-start_year+1)*12);
           exp4_cntotc_mat = NaN*ones(length(exp2_lat_vec),(stop_year-start_year+1)*12);
           exp4_et_mat = NaN*ones(length(exp2_lat_vec),(stop_year-start_year+1)*12);
           exp4_evpsoil_mat = NaN*ones(length(exp2_lat_vec),(stop_year-start_year+1)*12);
           exp4_evptrns_mat = NaN*ones(length(exp2_lat_vec),(stop_year-start_year+1)*12);
           exp4_evpintr_mat = NaN*ones(length(exp2_lat_vec),(stop_year-start_year+1)*12);
           exp4_evpsbln_mat = NaN*ones(length(exp2_lat_vec),(stop_year-start_year+1)*12);

        end

        if exist(exp1_read_file)
           
           exp1_cngpp_mat(:,mat_index) = double(ncread(exp1_read_file,'CNGPP'));
           exp1_cnnee_mat(:,mat_index) = double(ncread(exp1_read_file,'CNNEE'));
           exp1_cncloss_mat(:,mat_index) = double(ncread(exp1_read_file,'CLOSS'));
           exp1_cntlai_mat(:,mat_index) = double(ncread(exp1_read_file,'CNTLAI'));
           exp1_cntotc_mat(:,mat_index) = double(ncread(exp1_read_file,'CNTOTC'));
           exp1_evpsoil_mat(:,mat_index) = double(ncread(exp1_read_file,'EVPSOIL'));
           exp1_evptrns_mat(:,mat_index) = double(ncread(exp1_read_file,'EVPTRNS'));
           exp1_evpintr_mat(:,mat_index) = double(ncread(exp1_read_file,'EVPINTR'));
           exp1_evpsbln_mat(:,mat_index) = double(ncread(exp1_read_file,'EVPSBLN'));
           exp1_et_mat(:,mat_index) = exp1_evpsoil_mat(:,mat_index) + exp1_evptrns_mat(:,mat_index) ...
                                      + exp1_evpintr_mat(:,mat_index) + exp1_evpsbln_mat(:,mat_index);
         else
           exp1_cngpp_mat(:,mat_index) = NaN*ones(size(exp1_cngpp_mat,1),1);
           exp1_cnnee_mat(:,mat_index) = NaN*ones(size(exp1_cnnee_mat,1),1);
           exp1_cncloss_mat(:,mat_index) = NaN*ones(size(exp1_cncloss_mat,1),1);
           exp1_cntlai_mat(:,mat_index) = NaN*ones(size(exp1_cntlai_mat,1),1);
           exp1_cntotc_mat(:,mat_index) = NaN*ones(size(exp1_cntotc_mat,1),1);
           exp1_evpsoil_mat(:,mat_index) = NaN*ones(size(exp1_evpsoil_mat,1),1);
           exp1_evptrns_mat(:,mat_index) = NaN*ones(size(exp1_evptrns_mat,1),1);
           exp1_evpintr_mat(:,mat_index) = NaN*ones(size(exp1_evpintr_mat,1),1);
           exp1_evpsbln_mat(:,mat_index) = NaN*ones(size(exp1_evpsbln_mat,1),1);
           exp1_et_mat(:,mat_index) = NaN*ones(size(exp1_et_mat,1),1);
        end

        if exist(exp2_read_file)
           disp('read exp2')
           exp2_cngpp_mat(:,mat_index) = double(ncread(exp2_read_file,'CNGPP'));
           exp2_cnnee_mat(:,mat_index) = double(ncread(exp2_read_file,'CNNEE'));
           exp2_cncloss_mat(:,mat_index) = double(ncread(exp2_read_file,'CLOSS'));
           exp2_cntlai_mat(:,mat_index) = double(ncread(exp2_read_file,'CNTLAI'));
           exp2_cntotc_mat(:,mat_index) = double(ncread(exp2_read_file,'CNTOTC'));
           exp2_evpsoil_mat(:,mat_index) = double(ncread(exp2_read_file,'EVPSOIL'));
           exp2_evptrns_mat(:,mat_index) = double(ncread(exp2_read_file,'EVPTRNS'));
           exp2_evpintr_mat(:,mat_index) = double(ncread(exp2_read_file,'EVPINTR'));
           exp2_evpsbln_mat(:,mat_index) = double(ncread(exp2_read_file,'EVPSBLN'));
           exp2_et_mat(:,mat_index) = exp2_evpsoil_mat(:,mat_index) + exp2_evptrns_mat(:,mat_index) ...
                                      + exp2_evpintr_mat(:,mat_index) + exp2_evpsbln_mat(:,mat_index);
         else
           exp2_cngpp_mat(:,mat_index) = NaN*ones(size(exp2_cngpp_mat,1),1);
           exp2_cnnee_mat(:,mat_index) = NaN*ones(size(exp2_cnnee_mat,1),1);
           exp2_cncloss_mat(:,mat_index) = NaN*ones(size(exp2_cncloss_mat,1),1);
           exp2_cntlai_mat(:,mat_index) = NaN*ones(size(exp2_cntlai_mat,1),1);
           exp2_cntotc_mat(:,mat_index) = NaN*ones(size(exp2_cntotc_mat,1),1);
           exp2_evpsoil_mat(:,mat_index) = NaN*ones(size(exp2_evpsoil_mat,1),1);
           exp2_evptrns_mat(:,mat_index) = NaN*ones(size(exp2_evptrns_mat,1),1);
           exp2_evpintr_mat(:,mat_index) = NaN*ones(size(exp2_evpintr_mat,1),1);
           exp2_evpsbln_mat(:,mat_index) = NaN*ones(size(exp2_evpsbln_mat,1),1);
           exp2_et_mat(:,mat_index) = NaN*ones(size(exp2_et_mat,1),1);
        end

        if exist(exp3_read_file)
       
           exp3_cngpp_mat(:,mat_index) = double(ncread(exp3_read_file,'CNGPP'));
           exp3_cnnee_mat(:,mat_index) = double(ncread(exp3_read_file,'CNNEE'));
           exp3_cncloss_mat(:,mat_index) = double(ncread(exp3_read_file,'CLOSS'));
           exp3_cntlai_mat(:,mat_index) = double(ncread(exp3_read_file,'CNTLAI'));
           exp3_cntotc_mat(:,mat_index) = double(ncread(exp3_read_file,'CNTOTC'));
           exp3_evpsoil_mat(:,mat_index) = double(ncread(exp3_read_file,'EVPSOIL'));
           exp3_evptrns_mat(:,mat_index) = double(ncread(exp3_read_file,'EVPTRNS'));
           exp3_evpintr_mat(:,mat_index) = double(ncread(exp3_read_file,'EVPINTR'));
           exp3_evpsbln_mat(:,mat_index) = double(ncread(exp3_read_file,'EVPSBLN'));
           exp3_et_mat(:,mat_index) = exp3_evpsoil_mat(:,mat_index) + exp3_evptrns_mat(:,mat_index) ...
                                      + exp3_evpintr_mat(:,mat_index) + exp3_evpsbln_mat(:,mat_index);
         else
           exp3_cngpp_mat(:,mat_index) = NaN*ones(size(exp3_cngpp_mat,1),1);
           exp3_cnnee_mat(:,mat_index) = NaN*ones(size(exp3_cnnee_mat,1),1);
           exp3_cncloss_mat(:,mat_index) = NaN*ones(size(exp3_cncloss_mat,1),1);
           exp3_cntlai_mat(:,mat_index) = NaN*ones(size(exp3_cntlai_mat,1),1);
           exp3_cntotc_mat(:,mat_index) = NaN*ones(size(exp3_cntotc_mat,1),1);
           exp3_evpsoil_mat(:,mat_index) = NaN*ones(size(exp3_evpsoil_mat,1),1);
           exp3_evptrns_mat(:,mat_index) = NaN*ones(size(exp3_evptrns_mat,1),1);
           exp3_evpintr_mat(:,mat_index) = NaN*ones(size(exp3_evpintr_mat,1),1);
           exp3_evpsbln_mat(:,mat_index) = NaN*ones(size(exp3_evpsbln_mat,1),1);
           exp3_et_mat(:,mat_index) = NaN*ones(size(exp3_et_mat,1),1);
        end

        if exist(exp4_read_file)

           tmp = double(ncread(exp4_read_file,'CNGPP'));
           exp4_cngpp_mat(:,mat_index) = tmp(exp4_loc_index);
           tmp = double(ncread(exp4_read_file,'CNNEE'));
           exp4_cnnee_mat(:,mat_index) = tmp(exp4_loc_index);
           tmp = double(ncread(exp4_read_file,'CLOSS'));
           exp4_cncloss_mat(:,mat_index) = tmp(exp4_loc_index);
           tmp = double(ncread(exp4_read_file,'CNTLAI'));
           exp4_cntlai_mat(:,mat_index) = tmp(exp4_loc_index);
           tmp = double(ncread(exp4_read_file,'CNTOTC'));
           exp4_cntotc_mat(:,mat_index) = tmp(exp4_loc_index);
           tmp = double(ncread(exp4_read_file,'EVPSOIL'));
           exp4_evpsoil_mat(:,mat_index) = tmp(exp4_loc_index);
           tmp = double(ncread(exp4_read_file,'EVPTRNS'));
           exp4_evptrns_mat(:,mat_index) = tmp(exp4_loc_index);
           tmp = double(ncread(exp4_read_file,'EVPINTR'));
           exp4_evpintr_mat(:,mat_index) = tmp(exp4_loc_index);
           tmp = double(ncread(exp4_read_file,'EVPSBLN'));
           exp4_evpsbln_mat(:,mat_index) = tmp(exp4_loc_index);
           exp4_et_mat(:,mat_index) = exp4_evpsoil_mat(:,mat_index) + exp4_evptrns_mat(:,mat_index) ...
                                      + exp4_evpintr_mat(:,mat_index) + exp4_evpsbln_mat(:,mat_index);
         else
           exp4_cngpp_mat(:,mat_index) = NaN*ones(size(exp4_cngpp_mat,1),1);
           exp4_cnnee_mat(:,mat_index) = NaN*ones(size(exp4_cnnee_mat,1),1);
           exp4_cncloss_mat(:,mat_index) = NaN*ones(size(exp4_cncloss_mat,1),1);
           exp4_cntlai_mat(:,mat_index) = NaN*ones(size(exp4_cntlai_mat,1),1);
           exp4_cntotc_mat(:,mat_index) = NaN*ones(size(exp4_cntotc_mat,1),1);
           exp4_evpsoil_mat(:,mat_index) = NaN*ones(size(exp4_evpsoil_mat,1),1);
           exp4_evptrns_mat(:,mat_index) = NaN*ones(size(exp4_evptrns_mat,1),1);
           exp4_evpintr_mat(:,mat_index) = NaN*ones(size(exp4_evpintr_mat,1),1);
           exp4_evpsbln_mat(:,mat_index) = NaN*ones(size(exp4_evpsbln_mat,1),1);
           exp4_et_mat(:,mat_index) = NaN*ones(size(exp4_et_mat,1),1);
        end

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

    end 

    datenum_vec_y(y_index) = datenum(y,12,31);
end

ann_index = 0;

for y = start_year:stop_year

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

    tmp = exp1_cngpp_mat(:,year_index);
    mean_global_ann_cngpp_exp1(ann_index) = nanmean(tmp(:));
    mean_local_ann_cngpp_exp1(:,ann_index) = nanmean(tmp,2);

    tmp = exp1_cnnee_mat(:,year_index);
    mean_global_ann_cnnee_exp1(ann_index) = nanmean(tmp(:));
    mean_local_ann_cnnee_exp1(:,ann_index) = nanmean(tmp,2);

    tmp = exp1_cncloss_mat(:,year_index);
    mean_global_ann_cncloss_exp1(ann_index) = nanmean(tmp(:));
    mean_local_ann_cncloss_exp1(:,ann_index) = nanmean(tmp,2);

    tmp = exp1_cntlai_mat(:,year_index);
    mean_global_ann_cntlai_exp1(ann_index) = nanmean(tmp(:));
    mean_local_ann_cntlai_exp1(:,ann_index) = nanmean(tmp,2);

    tmp = exp1_cntotc_mat(:,year_index);
    mean_global_ann_cntotc_exp1(ann_index) = nanmean(tmp(:));
    mean_local_ann_cntotc_exp1(:,ann_index) = nanmean(tmp,2);

    tmp = exp1_evpsoil_mat(:,year_index);
    mean_global_ann_evpsoil_exp1(ann_index) = nanmean(tmp(:));
    mean_local_ann_evpsoil_exp1(:,ann_index) = nanmean(tmp,2);

    tmp = exp1_evptrns_mat(:,year_index);
    mean_global_ann_evptrns_exp1(ann_index) = nanmean(tmp(:));
    mean_local_ann_evptrns_exp1(:,ann_index) = nanmean(tmp,2);

    tmp = exp1_evpintr_mat(:,year_index);
    mean_global_ann_evpintr_exp1(ann_index) = nanmean(tmp(:));
    mean_local_ann_evpintr_exp1(:,ann_index) = nanmean(tmp,2);

    tmp = exp1_evpsbln_mat(:,year_index);
    mean_global_ann_evpsbln_exp1(ann_index) = nanmean(tmp(:));
    mean_local_ann_evpsbln_exp1(:,ann_index) = nanmean(tmp,2);

    tmp = exp1_et_mat(:,year_index);
    mean_global_ann_et_exp1(ann_index) = nanmean(tmp(:));
    mean_local_ann_et_exp1(:,ann_index) = nanmean(tmp,2);

    tmp = exp2_cngpp_mat(:,year_index);
    mean_global_ann_cngpp_exp2(ann_index) = nanmean(tmp(:));
    mean_local_ann_cngpp_exp2(:,ann_index) = nanmean(tmp,2);

    tmp = exp2_cnnee_mat(:,year_index);
    mean_global_ann_cnnee_exp2(ann_index) = nanmean(tmp(:));
    mean_local_ann_cnnee_exp2(:,ann_index) = nanmean(tmp,2);

    tmp = exp2_cncloss_mat(:,year_index);
    mean_global_ann_cncloss_exp2(ann_index) = nanmean(tmp(:));
    mean_local_ann_cncloss_exp2(:,ann_index) = nanmean(tmp,2);

    tmp = exp2_cntlai_mat(:,year_index);
    mean_global_ann_cntlai_exp2(ann_index) = nanmean(tmp(:));
    mean_local_ann_cntlai_exp2(:,ann_index) = nanmean(tmp,2);

    tmp = exp2_cntotc_mat(:,year_index);
    mean_global_ann_cntotc_exp2(ann_index) = nanmean(tmp(:));
    mean_local_ann_cntotc_exp2(:,ann_index) = nanmean(tmp,2);

    tmp = exp2_evpsoil_mat(:,year_index);
    mean_global_ann_evpsoil_exp2(ann_index) = nanmean(tmp(:));
    mean_local_ann_evpsoil_exp2(:,ann_index) = nanmean(tmp,2);

    tmp = exp2_evptrns_mat(:,year_index);
    mean_global_ann_evptrns_exp2(ann_index) = nanmean(tmp(:));
    mean_local_ann_evptrns_exp2(:,ann_index) = nanmean(tmp,2);

    tmp = exp2_evpintr_mat(:,year_index);
    mean_global_ann_evpintr_exp2(ann_index) = nanmean(tmp(:));
    mean_local_ann_evpintr_exp2(:,ann_index) = nanmean(tmp,2);

    tmp = exp2_evpsbln_mat(:,year_index);
    mean_global_ann_evpsbln_exp2(ann_index) = nanmean(tmp(:));
    mean_local_ann_evpsbln_exp2(:,ann_index) = nanmean(tmp,2);

    tmp = exp2_et_mat(:,year_index);
    mean_global_ann_et_exp2(ann_index) = nanmean(tmp(:));
    mean_local_ann_et_exp2(:,ann_index) = nanmean(tmp,2);

    tmp = exp3_cngpp_mat(:,year_index);
    mean_global_ann_cngpp_exp3(ann_index) = nanmean(tmp(:));
    mean_local_ann_cngpp_exp3(:,ann_index) = nanmean(tmp,2);

    tmp = exp3_cnnee_mat(:,year_index);
    mean_global_ann_cnnee_exp3(ann_index) = nanmean(tmp(:));
    mean_local_ann_cnnee_exp3(:,ann_index) = nanmean(tmp,2);

    tmp = exp3_cncloss_mat(:,year_index);
    mean_global_ann_cncloss_exp3(ann_index) = nanmean(tmp(:));
    mean_local_ann_cncloss_exp3(:,ann_index) = nanmean(tmp,2);

    tmp = exp3_cntlai_mat(:,year_index);
    mean_global_ann_cntlai_exp3(ann_index) = nanmean(tmp(:));
    mean_local_ann_cntlai_exp3(:,ann_index) = nanmean(tmp,2);

    tmp = exp3_cntotc_mat(:,year_index);
    mean_global_ann_cntotc_exp3(ann_index) = nanmean(tmp(:));
    mean_local_ann_cntotc_exp3(:,ann_index) = nanmean(tmp,2);

    tmp = exp3_evpsoil_mat(:,year_index);
    mean_global_ann_evpsoil_exp3(ann_index) = nanmean(tmp(:));
    mean_local_ann_evpsoil_exp3(:,ann_index) = nanmean(tmp,2);

    tmp = exp3_evptrns_mat(:,year_index);
    mean_global_ann_evptrns_exp3(ann_index) = nanmean(tmp(:));
    mean_local_ann_evptrns_exp3(:,ann_index) = nanmean(tmp,2);

    tmp = exp3_evpintr_mat(:,year_index);
    mean_global_ann_evpintr_exp3(ann_index) = nanmean(tmp(:));
    mean_local_ann_evpintr_exp3(:,ann_index) = nanmean(tmp,2);

    tmp = exp3_evpsbln_mat(:,year_index);
    mean_global_ann_evpsbln_exp3(ann_index) = nanmean(tmp(:));
    mean_local_ann_evpsbln_exp3(:,ann_index) = nanmean(tmp,2);

    tmp = exp3_et_mat(:,year_index);
    mean_global_ann_et_exp3(ann_index) = nanmean(tmp(:));
    mean_local_ann_et_exp3(:,ann_index) = nanmean(tmp,2);

    tmp = exp4_cngpp_mat(:,year_index);
    mean_global_ann_cngpp_exp4(ann_index) = nanmean(tmp(:));
    mean_local_ann_cngpp_exp4(:,ann_index) = nanmean(tmp,2);

    tmp = exp4_cnnee_mat(:,year_index);
    mean_global_ann_cnnee_exp4(ann_index) = nanmean(tmp(:));
    mean_local_ann_cnnee_exp4(:,ann_index) = nanmean(tmp,2);

    tmp = exp4_cncloss_mat(:,year_index);
    mean_global_ann_cncloss_exp4(ann_index) = nanmean(tmp(:));
    mean_local_ann_cncloss_exp4(:,ann_index) = nanmean(tmp,2);

    tmp = exp4_cntlai_mat(:,year_index);
    mean_global_ann_cntlai_exp4(ann_index) = nanmean(tmp(:));
    mean_local_ann_cntlai_exp4(:,ann_index) = nanmean(tmp,2);

    tmp = exp4_cntotc_mat(:,year_index);
    mean_global_ann_cntotc_exp4(ann_index) = nanmean(tmp(:));
    mean_local_ann_cntotc_exp4(:,ann_index) = nanmean(tmp,2);

    tmp = exp4_evpsoil_mat(:,year_index);
    mean_global_ann_evpsoil_exp4(ann_index) = nanmean(tmp(:));
    mean_local_ann_evpsoil_exp4(:,ann_index) = nanmean(tmp,2);

    tmp = exp4_evptrns_mat(:,year_index);
    mean_global_ann_evptrns_exp4(ann_index) = nanmean(tmp(:));
    mean_local_ann_evptrns_exp4(:,ann_index) = nanmean(tmp,2);

    tmp = exp4_evpintr_mat(:,year_index);
    mean_global_ann_evpintr_exp4(ann_index) = nanmean(tmp(:));
    mean_local_ann_evpintr_exp4(:,ann_index) = nanmean(tmp,2);

    tmp = exp4_evpsbln_mat(:,year_index);
    mean_global_ann_evpsbln_exp4(ann_index) = nanmean(tmp(:));
    mean_local_ann_evpsbln_exp4(:,ann_index) = nanmean(tmp,2);

    tmp = exp4_et_mat(:,year_index);
    mean_global_ann_et_exp4(ann_index) = nanmean(tmp(:));
    mean_local_ann_et_exp4(:,ann_index) = nanmean(tmp,2);

end

full_index = [4 114 118];
index_ity = exp2_ityp(full_index,:);
index_fveg = exp2_fveg(full_index,:);


for i = 1:length(full_index)
    figure
    LS = '-';
    plot(start_year:stop_year,mean_local_ann_cntlai_exp1(full_index(i),:),'LineWidth',2);
    hold on
    plot(start_year:stop_year,mean_local_ann_cntlai_exp2(full_index(i),:),'LineWidth',2);
    hold on
    plot(start_year:stop_year,mean_local_ann_cntlai_exp3(full_index(i),:),'LineWidth',2);
    hold off
    grid on
    xlabel('year')
    ylabel('LAI [-]')
    title('annual average LAI [-] @280ppm')
    legend_str{1} = 'baseline';
    legend_str{2} = 'GC fix';
    legend_str{3} = 'new fcur';
    legend(legend_str,'Location','bestoutside')
    set(gcf,'paperunits','centimeters')
    set(gcf,'paperposition',[1 1 25 15])
    saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/tests/CNTLAI_ann_CN51_cycles_280_loc' num2str(i)],'png')

    color_tmp = get(gca,'ColorOrder');
    nyears = 28;

%    figure
%    LS = '-';
%    plot(datenum_vec(1:nyears*12),exp1_cntlai_mat(full_index(i),1:nyears*12),'LineWidth',2,'Color',color_tmp(1,:));
%    hold on 
%    plot(datenum_vec_y(1:nyears),mean_local_ann_cntlai_exp1(full_index(i),1:nyears),'LineWidth',2,'LineStyle','--','Color',color_tmp(1,:));
%    hold on
%    plot(datenum_vec(1:nyears*12),exp2_cntlai_mat(full_index(i),1:nyears*12),'LineWidth',2,'Color',color_tmp(2,:));
%    hold on
%    plot(datenum_vec_y(1:nyears),mean_local_ann_cntlai_exp2(full_index(i),1:nyears),'LineWidth',2,'LineStyle','--','Color',color_tmp(2,:));
%   % hold on
%    %plot(datenum_vec,exp3_cntlai_mat(full_index(i),:),'LineWidth',2);
%   % hold on
%   % plot(datenum_vec,exp4_cntlai_mat(full_index(i),:),'LineWidth',2);
%    hold off
%    grid on
%    datetick('x','mm-yyyy')
%    xlabel('date [mm-yyyy]')
%    ylabel('LAI [-]')
%    title('monthly average LAI [-] @280ppm')
%    legend_str{1} = 'baseline monthly avg.';
%    legend_str{2} = 'baseline yearly avg.';
%    legend_str{3} = 'GC fix monthly avg.';
%    legend_str{4} = 'GC fix yearly avg.';
%   % legend_str{3} = 'new fcur';
%    %legend_str{4} = 'global';
%    legend(legend_str,'Location','bestoutside')
%    set(gcf,'paperunits','centimeters')
%    set(gcf,'paperposition',[1 1 25 15])
%    saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/tests/CNTLAI_mon_CN51_cycles_280_loc' num2str(i)],'png')

    figure
    LS = '-';
    plot(datenum_vec,exp1_cntlai_mat(full_index(i),:),'LineWidth',2,'Color',color_tmp(1,:));
    hold on
    plot(datenum_vec_y,mean_local_ann_cntlai_exp1(full_index(i),:),'LineWidth',2,'LineStyle','--','Color',color_tmp(1,:));
    hold on
    plot(datenum_vec,exp2_cntlai_mat(full_index(i),:),'LineWidth',2,'Color',color_tmp(2,:));
    hold on
    plot(datenum_vec_y,mean_local_ann_cntlai_exp2(full_index(i),:),'LineWidth',2,'LineStyle','--','Color',color_tmp(2,:));
   % hold on
    %plot(datenum_vec,exp3_cntlai_mat(full_index(i),:),'LineWidth',2);
   % hold on
   % plot(datenum_vec,exp4_cntlai_mat(full_index(i),:),'LineWidth',2);
    hold off
    grid on
    datetick('x','mm-yyyy')
    xlabel('date [mm-yyyy]')
    ylabel('LAI [-]')
    title('monthly average LAI [-] @280ppm')
    legend_str{1} = 'baseline monthly avg.';
    legend_str{2} = 'baseline yearly avg.';
    legend_str{3} = 'GC fix monthly avg.';
    legend_str{4} = 'GC fix yearly avg.';
   % legend_str{3} = 'new fcur';
    %legend_str{4} = 'global';
    legend(legend_str,'Location','bestoutside')
    set(gcf,'paperunits','centimeters')
    set(gcf,'paperposition',[1 1 25 15])
    saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/tests/CNTLAI_mon_CN51_cycles_280_loc' num2str(i)],'png')


    figure
    LS = '-';
    plot(datenum_vec,exp1_cntotc_mat(full_index(i),:),'LineWidth',2);
    hold on
    plot(datenum_vec,exp2_cntotc_mat(full_index(i),:),'LineWidth',2);
    hold on
    plot(datenum_vec,exp3_cntotc_mat(full_index(i),:),'LineWidth',2);
    hold off
    grid on
    datetick('x','mm-yyyy')
    xlabel('date [mm-yyyy]')
    ylabel('TOTC [kgC m^{-2}]')
    title('monthly average total carbon [kgC m^{-2}] @280ppm')
    legend_str{1} = 'baseline';
    legend_str{2} = 'high-precip';
    legend_str{3} = 'new fcur';
    legend(legend_str,'Location','bestoutside')
    set(gcf,'paperunits','centimeters')
    set(gcf,'paperposition',[1 1 25 15])
    saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/tests/CNTOTC_mon_CN51_cycles_280_loc' num2str(i)],'png')


    for y = 2:size(mean_local_ann_cntlai_exp2,2)
        cntlai_diff(i,y-1) = mean_local_ann_cntlai_exp2(full_index(i),y) - mean_local_ann_cntlai_exp2(full_index(i),y-1);
    end
end

figure
LS = '-';
plot(start_year:stop_year,mean_global_ann_cntlai_exp1,'LineWidth',2);
hold on
plot(start_year:stop_year,mean_global_ann_cntlai_exp2,'LineWidth',2);
%hold on
%plot(start_year:stop_year,mean_global_ann_cntlai_exp3,'LineWidth',2);
hold off
grid on
xlabel('year')
ylabel('LAI [-]')
title('annual average LAI [-] @280ppm')
legend_str{1} = 'baseline';
legend_str{2} = 'GC fix';
%legend_str{3} = 'new fcur';
legend(legend_str,'Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/tests/CNTLAI_ann_CN51_cycles_280'],'png')

