clc
clear all
close all

cn51_read_path = '/discover/nobackup/jkolassa/CN51_global/GEOSldas_CN51_280_global/output/SMAP_EASEv2_M36/cat/ens0000/';
clm45_read_path = '/discover/nobackup/projects/geoscm/fzeng/GEOSldas45_M36_rst_ldas06_15/output/SMAP_EASEv2_M36/cat/ens0000/';
cn40_read_path = '/discover/nobackup/projects/geoscm/fzeng/Catchment-CN40_9km/GEOSldas_1981_present/GEOSldas_CNCLM40_9km/output/SMAP_EASEv2_M09_GLOBAL/cat/ens0000/';
cn45_280_read_path = '/discover/nobackup/jkolassa/test_new/GEOSldas_CN45_global_280/output/SMAP_EASEv2_M36/cat/ens0000/';
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_path = '/discover/nobackup/projects/geoscm/fzeng/CLM/archive/clm5_test/lnd/hist/';
clm5_of_path = '/discover/nobackup/jkolassa/CLM5_sim/';

start_year = 1981;
stop_year = 2016; 
nyears_280 = 9;

rst_file = '/discover/nobackup/jkolassa/CN51_global/GEOSldas_CN51_280_global/output/SMAP_EASEv2_M36/rs/ens0000/Y1980/M01/GEOSldas_CN51_280_global.catchcnclm51_internal_rst.19800101_0000';
tile_id = ncread(rst_file,'TILE_ID');
ityp = ncread(rst_file,'ITY');
fvg = ncread (rst_file,'FVG');

for i = 1:length(ityp)
    [max_value,max_index] = max(fvg(i,:));
    if (max_value>=0.7)
       dom_pft(i) = ityp(i,max_index);
    else
       dom_pft(i) = NaN;
    end
end 

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'...
           };

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

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

% compute mapping indices

cn51_read_file = [cn51_read_path 'Y1980/M01/GEOSldas_CN51_280_global.tavg24_1d_lnd_Nt.monthly.198001.nc4'];
clm45_read_file = [clm45_read_path 'Y2001/M01/GEOSldas45_M36_rst_ldas06_15.tavg24_1d_lnd_Nt.monthly.200101.nc4'];
cn40_read_file = [cn40_read_path 'Y2001/M01/GEOSldas_CNCLM40_9km.tavg24_1d_lnd_Nt.monthly.200101.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'];

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

clm45_lat_vec = double(ncread(clm45_read_file,'lat'));
clm45_lon_vec = double(ncread(clm45_read_file,'lon'));

cn40_lat_vec = double(ncread(clm45_read_file,'lat'));
cn40_lon_vec = double(ncread(clm45_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_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;

% 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');
%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);
    tmp = find(lon_dist==min(lon_dist));
    fc_lon_index(i) = tmp(1);

    dist_vec = sqrt((clm45_lat_vec-lat_vec(i)).^2 + (clm45_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);
    tmp = find(lon_dist==min(lon_dist));
    fs_lon_index(i) = tmp(1);

    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);
    tmp = find(lon_dist==min(lon_dist));
    maiac_lon_index(i) = tmp(1);

    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);
    tmp = find(lon_dist==min(lon_dist));
    vpm_lon_index(i) = tmp(1);

    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

% read data

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

    % 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');

        cn51_read_file = [cn51_read_path 'Y' y_str '/M' m_str '/GEOSldas_CN51_280_global.tavg24_1d_lnd_Nt.monthly.' y_str m_str '.nc4'];
        clm45_read_file = [clm45_read_path 'Y' y_str '/M' m_str '/GEOSldas45_M36_rst_ldas06_15.tavg24_1d_lnd_Nt.monthly.' y_str m_str '.nc4'];
        cn40_read_file = [cn40_read_path 'Y' y_str '/M' m_str '/GEOSldas_CNCLM40_9km.tavg24_1d_lnd_Nt.monthly.' y_str m_str '.nc4'];

        cn45_280_read_file = [cn45_280_read_path 'Y' y_str '/M' m_str '/GEOSldas_CN45_global_280.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_path 'clm5_test.clm2.h0.' y_index_str '-' m_str '.nc'];


        if (mat_index == 1)

           cngpp_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           clm45_cngpp_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           cn40_gpp_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           cn45_280_gpp_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);
           clm5_gpp_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);
        end

        cngpp_mat(:,mat_index) = double(ncread(cn51_read_file,'CNGPP'));
      
        if (exist(cn45_280_read_file)~=0)
           cn45_280_gpp_mat(:,mat_index) = double(ncread(cn45_280_read_file,'CNGPP'));
        end

        tmp_mat = double(ncread(clm45_read_file,'CNGPP'));
        clm45_cngpp_mat(:,mat_index) = tmp_mat(data_index);
      
        tmp_mat = double(ncread(cn40_read_file,'CNGPP'));
        cn40_gpp_mat(:,mat_index) = tmp_mat(cn40_data_index);

        if exist(clm5_read_file)
           tmp_gpp = double(ncread(clm5_read_file,'GPP'));
           for l = 1:length(lat_vec)
               clm5_gpp_mat(l,mat_index) = tmp_gpp(clm5_lon_index(l),clm5_lat_index(l));
           end
         else
           clm5_gpp_mat(:,mat_index) = NaN*ones(length(lat_vec),1);
        end
  
        if (y>=2003)

           if (exist(fluxsat_gpp_read_file)~=0)
              fs_gpp_tmp_mat = ncread(fluxsat_gpp_read_file,'GPP');
              for i = 1:length(lat_vec)
                  fs_gpp_mat(i,mat_index) = fs_gpp_tmp_mat(fs_lon_index(i),fs_lat_index(i));
              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 i = 1:length(lat_vec)
                  maiac_gpp_mat(i,mat_index) = maiac_gpp_tmp_mat(maiac_lon_index(i),maiac_lat_index(i));
              end
           end

           if (exist(vpm_gpp_read_file)~=0)
              vpm_gpp_tmp_mat = ncread(vpm_gpp_read_file,'GPP');
              for i = 1:length(lat_vec)
                  vpm_gpp_mat(i,mat_index) = vpm_gpp_tmp_mat(vpm_lon_index(i),vpm_lat_index(i));
              end
           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);
            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);
            end
        end


    end

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

% remove extra entries

cngpp_mat(:,(mat_index+1):end) = [];
clm45_cngpp_mat(:,(mat_index+1):end) = [];
cn40_gpp_mat(:,(mat_index+1):end) = [];
cn45_280_gpp_mat(:,(mat_index+1):end) = [];
fc_gpp_mat(:,(mat_index+1):end) = [];
fs_gpp_mat(:,(mat_index+1):end) = [];
maiac_gpp_mat(:,(mat_index+1):end) = [];
vpm_gpp_mat(:,(mat_index+1):end) = [];
clm5_of_pi_gpp_mat(:,(mat_index+1):end) = [];
clm5_of_mod_gpp_mat(:,(mat_index+1):end) = [];

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

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

clm5_gpp_mat = clm5_gpp_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);

ann_index = 0;

for y = start_year:stop_year

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

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

    tmp = cn40_gpp_mat(:,year_index);
    mean_global_ann_cn40_gpp(ann_index) = nanmean(tmp(:));
    mean_local_ann_cn40_gpp(:,ann_index) = nanmean(tmp,2);

    tmp = cn45_280_gpp_mat(:,year_index);
    mean_global_ann_cn45_280_gpp(ann_index) = nanmean(tmp(:));
    mean_local_ann_cn45_280_gpp(:,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 = clm45_cngpp_mat(:,year_index);
    mean_global_ann_clm45_cngpp(ann_index) = nanmean(tmp(:));
    mean_local_ann_clm45_cngpp(:,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_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);

end

figure
plot(start_year:stop_year,mean_global_ann_cngpp,'LineWidth',2);
hold on
plot(start_year:stop_year,mean_global_ann_clm45_cngpp,'LineWidth',2);
hold on
plot(start_year:stop_year,mean_global_ann_cn45_280_gpp,'LineWidth',2);
hold on
%plot(start_year:stop_year,mean_global_ann_cn40_gpp,'LineWidth',2);
%hold on
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_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('CN51 (@280ppm)','CN45 (@390ppm)','CN45 (@280ppm)','FluxCom','FluxSat','MAIAC','VPM','CLM5','Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/CN51_eval_Fluxcom_GPP'],'png')

% compute and plot global difference maps

mean_gpp_cn51 = nanmean(cngpp_mat,2);
mean_gpp_cn45 = nanmean(clm45_cngpp_mat,2);
mean_gpp_fc = nanmean(fc_gpp_mat,2);
mean_gpp_clm5 = nanmean(clm5_gpp_mat,2);
mean_gpp_clm5_of_pi = nanmean(clm5_of_pi_gpp_mat,2);
mean_gpp_clm5_of_mod = nanmean(clm5_of_mod_gpp_mat,2);

cngpp_mat(cngpp_mat==0) = NaN;

tmp = cngpp_mat(:,1:120)*(36000^2)*(31*24*60*60)./10;
tmp_year = nansum(tmp,2);
tot_ann_gpp_cn51_2 = nansum(tmp_year)*10^(-12);
tmp = clm45_cngpp_mat(1:120)*(36000^2)*(31*24*60*60)./10;
tmp_year = nansum(tmp,2);
tot_ann_gpp_cn45_2 = nansum(tmp_year)*10^(-12);
tmp = cn45_280_gpp_mat(1:120)*(36000^2)*(31*24*60*60)./10;
tmp_year = nansum(tmp,2);
tot_ann_gpp_cn45_280_2 = nansum(tmp_year)*10^(-12);


tmp = mean_gpp_cn51*(36000^2)*(365*24*60*60);  % convert from kgC m^-2 s^-1 to kgC y^-1
tot_ann_gpp_cn51 = nansum(tmp)*(10^(-12));     % compute global sum and convert to PgC y^-1
tmp = mean_gpp_cn45*(36000^2)*(365*24*60*60);  % convert from kgC m^-2 s^-1 to kgC y^-1
tot_ann_gpp_cn45 = nansum(tmp)*(10^(-12));     % compute global sum and convert to PgC y^-1
tmp = mean_gpp_fc*(36000^2)*(365*24*60*60);  % convert from kgC m^-2 s^-1 to kgC y^-1
tot_ann_gpp_fc = nansum(tmp)*(10^(-12));     % compute global sum and convert to PgC y^-1
tmp = mean_gpp_clm5*(36000^2)*(365*24*60*60);  % convert from kgC m^-2 s^-1 to kgC y^-1
tot_ann_gpp_clm5 = nansum(tmp)*(10^(-12));     % compute global sum and convert to PgC y^-1
tmp = mean_gpp_clm5_of_pi*(36000^2)*(365*24*60*60);  % convert from kgC m^-2 s^-1 to kgC y^-1
tot_ann_gpp_clm5_of_pi = nansum(tmp)*(10^(-12));     % compute global sum and convert to PgC y^-1
tmp = mean_gpp_clm5_of_mod*(36000^2)*(365*24*60*60);  % convert from kgC m^-2 s^-1 to kgC y^-1
tot_ann_gpp_clm5_of_mod = nansum(tmp)*(10^(-12));     % compute global sum and convert to PgC y^-1


mean_gpp_cn51_short = nanmean(mean_local_ann_cngpp(:,1:nyears_280),2);
mean_gpp_cn45_short = nanmean(mean_local_ann_clm45_cngpp(:,1:nyears_280),2);
mean_gpp_cn45_280_short = nanmean(mean_local_ann_cn45_280_gpp(:,1:nyears_280),2);

tmp = mean_gpp_cn51_short*(36000^2)*(365*24*60*60);  % convert from kgC m^-2 s^-1 to kgC y^-1
tot_ann_gpp_cn51_short = nansum(tmp)*(10^(-12));     % compute global sum and convert to PgC y^-1
tmp = mean_gpp_cn45_short*(36000^2)*(365*24*60*60);  % convert from kgC m^-2 s^-1 to kgC y^-1
tot_ann_gpp_cn45_short = nansum(tmp)*(10^(-12));     % compute global sum and convert to PgC y^-1
tmp = mean_gpp_cn45_280_short*(36000^2)*(365*24*60*60);  % convert from kgC m^-2 s^-1 to kgC y^-1
tot_ann_gpp_cn45_280_short = nansum(tmp)*(10^(-12));     % compute global sum and convert to PgC y^-1

mean_gpp_cn51_short_gcy = mean_gpp_cn51_short*1000*(365*24*60*60); % gC m^-2 y^-1
mean_gpp_cn45_short_gcy = mean_gpp_cn45_short*1000*(365*24*60*60); % gC m^-2 y^-1
mean_gpp_cn45_280_short_gcy = mean_gpp_cn45_280_short*1000*(365*24*60*60); % gC m^-2 y^-1

diff_gpp_cn51_fc = cngpp_mat - fc_gpp_mat;
mean_diff_gpp_cn51_fc = nanmean(diff_gpp_cn51_fc,2);
tmp = (fc_gpp_mat-cngpp_mat).^2;
mse_gpp_cn51_fc = nanmean(tmp,2);
rmse_gpp_cn51_fc = sqrt(mse_gpp_cn51_fc);

diff_gpp_cn45_fc = clm45_cngpp_mat - fc_gpp_mat;
mean_diff_gpp_cn45_fc = nanmean(diff_gpp_cn45_fc,2);
tmp = (fc_gpp_mat-clm45_cngpp_mat).^2;
mse_gpp_cn45_fc = nanmean(tmp,2);
rmse_gpp_cn45_fc = sqrt(mse_gpp_cn45_fc);

closeness_gpp_fc = rmse_gpp_cn51_fc - rmse_gpp_cn45_fc;

test_lat = unique(lat_vec);
test_lon = unique(lon_vec);

[lon_mat,lat_mat] = meshgrid(test_lon,test_lat);

mean_diff_gpp_cn51_fc_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
mean_diff_gpp_cn45_fc_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
rmse_gpp_cn51_fc_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
rmse_gpp_cn45_fc_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
closeness_gpp_fc_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
R_gpp_cn51_fc_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
R_gpp_cn45_fc_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));


mean_gpp_cn51_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
mean_gpp_cn45_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
mean_gpp_fc_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
mean_gpp_clm5_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
mean_gpp_clm5_of_pi_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
mean_gpp_clm5_of_mod_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));

mean_gpp_cn51_year1_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
mean_gpp_cn45_280_year1_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
mean_gpp_cn45_year1_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
mean_gpp_fc_year1_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));

mean_gpp_cn51_short_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
mean_gpp_cn45_280_short_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
mean_gpp_cn45_short_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));

mean_gpp_cn51_short_gcy_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
mean_gpp_cn45_280_short_gcy_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
mean_gpp_cn45_short_gcy_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));

for i = 1:length(lon_vec)

    cngpp_vec = cngpp_mat(i,:);
    fc_gpp_vec = fc_gpp_mat(i,:);
    clm45_cngpp_vec = clm45_cngpp_mat(i,:);

    data_index_cn51 = isnan(cngpp_vec) + isnan(fc_gpp_vec);
    data_index_cn51 = data_index_cn51 == 0;

    R = corrcoef(cngpp_vec(data_index_cn51),fc_gpp_vec(data_index_cn51));
    R_gpp_cn51_fc(i) = R(1,2);

    data_index_cn45 = isnan(clm45_cngpp_vec) + isnan(fc_gpp_vec);
    data_index_cn45 = data_index_cn45 == 0;

    R = corrcoef(clm45_cngpp_vec(data_index_cn45),fc_gpp_vec(data_index_cn45));
    R_gpp_cn45_fc(i) = R(1,2);

    % map tile data to plottable grid

    lat_index = test_lat == lat_vec(i);
    lon_index = test_lon == lon_vec(i);

    mean_gpp_cn51_mat(lat_index,lon_index) = mean_gpp_cn51(i); 
    mean_gpp_cn45_mat(lat_index,lon_index) = mean_gpp_cn45(i);
    mean_gpp_fc_mat(lat_index,lon_index) = mean_gpp_fc(i);
    mean_gpp_clm5_mat(lat_index,lon_index) = mean_gpp_clm5(i);
    mean_gpp_clm5_of_pi_mat(lat_index,lon_index) = mean_gpp_clm5_of_pi(i);
    mean_gpp_clm5_of_mod_mat(lat_index,lon_index) = mean_gpp_clm5_of_mod(i);

    mean_diff_gpp_cn51_fc_mat(lat_index,lon_index) = mean_diff_gpp_cn51_fc(i);
    mean_diff_gpp_cn45_fc_mat(lat_index,lon_index) = mean_diff_gpp_cn45_fc(i);
    rmse_gpp_cn51_fc_mat(lat_index,lon_index) = rmse_gpp_cn51_fc(i);
    rmse_gpp_cn45_fc_mat(lat_index,lon_index) = rmse_gpp_cn45_fc(i);
    R_gpp_cn51_fc_mat(lat_index,lon_index) = R_gpp_cn51_fc(i);
    R_gpp_cn45_fc_mat(lat_index,lon_index) = R_gpp_cn45_fc(i);

    mean_gpp_cn51_year1_mat(lat_index,lon_index) = mean_local_ann_cngpp(i,1);
    mean_gpp_cn45_280_year1_mat(lat_index,lon_index) = mean_local_ann_cn45_280_gpp(i,1);
    mean_gpp_cn45_year1_mat(lat_index,lon_index) = mean_local_ann_clm45_cngpp(i,1);
    mean_gpp_fc_year1_mat(lat_index,lon_index) = mean_local_ann_fc_gpp(i,1);

    mean_gpp_cn51_short_mat(lat_index,lon_index) = mean_gpp_cn51_short(i);
    mean_gpp_cn45_280_short_mat(lat_index,lon_index) = mean_gpp_cn45_280_short(i);
    mean_gpp_cn45_short_mat(lat_index,lon_index) = mean_gpp_cn45_short(i);

    mean_gpp_cn51_short_gcy_mat(lat_index,lon_index) = mean_gpp_cn51_short_gcy(i);
    mean_gpp_cn45_280_short_gcy_mat(lat_index,lon_index) = mean_gpp_cn45_280_short_gcy(i);
    mean_gpp_cn45_short_gcy_mat(lat_index,lon_index) = mean_gpp_cn45_short_gcy(i);
end

load coastlines;
[latcells, loncells] = polysplit(coastlat, coastlon);

latlim = [-60 90];
lonlim = [-180 180];

cmap_mat = cmocean('speed',12);
%cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
cmap = cmap_mat;
s1 = quantile(mean_gpp_cn51_mat(:),0.01);
s2 = quantile(mean_gpp_cn51_mat(:),0.99);
caxis_min = s1;
caxis_max = s2;

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,mean_gpp_cn51_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['mean GPP CN51 [kgC 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/global/GPP_CN51'],'png')

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,mean_gpp_cn45_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['mean GPP CN45 [kgC 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/global/GPP_CN45'],'png')

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,mean_gpp_fc_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['mean GPP FluxCom [kgC 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/global/GPP_FC'],'png')

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,mean_gpp_clm5_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['mean GPP CLM5 [kgC 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/global/GPP_CLM5'],'png')

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,mean_gpp_clm5_of_pi_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['mean GPP CLM5 NCAR [kgC m^{-2} s^{-1}] (1850 - 1890)'])
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/global/GPP_CLM5_NCAR_1850_1890'],'png')

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,mean_gpp_clm5_of_mod_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['mean GPP CLM5 NCAR [kgC m^{-2} s^{-1}] (1980 - 2014))'])
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/global/GPP_CLM5_NCAR_1980_2014'],'png')

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,mean_gpp_cn51_year1_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['mean GPP CN51 [kgC 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/global/GPP_CN51_year1'],'png')

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)     
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,mean_gpp_cn45_280_year1_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['mean GPP CN45 @280ppm [kgC 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/global/GPP_CN45_280_year1'],'png')

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,mean_gpp_cn45_year1_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['mean GPP CN45 @390ppm [kgC 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/global/GPP_CN45_year1'],'png')

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)     
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,mean_gpp_fc_year1_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['mean GPP FluxCom [kgC 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/global/GPP_FC_year1'],'png')

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,mean_gpp_cn51_short_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['mean GPP CN51 [kgC 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/global/GPP_CN51_short'],'png')

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,mean_gpp_cn45_short_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['mean GPP CN45 [kgC 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/global/GPP_CN45_short'],'png')

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,mean_gpp_cn45_280_short_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['mean GPP CN45 (@280ppm) [kgC 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/global/GPP_CN45_280_short'],'png')

cmap_mat = cmocean('speed',10);
%cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
cmap = cmap_mat;
s1 = quantile(mean_gpp_cn51_mat(:),0.01);
s2 = quantile(mean_gpp_cn51_mat(:),0.99);
caxis_min = 0;
caxis_max = 4000;

newColormap = cmap(1,:);
newColormap = [newColormap; repmat(cmap(2,:),1,1)];
newColormap = [newColormap; repmat(cmap(3,:),2,1)];
newColormap = [newColormap; repmat(cmap(4,:),4,1)];
newColormap = [newColormap; repmat(cmap(5,:),4,1)];
newColormap = [newColormap; repmat(cmap(6,:),4,1)];
newColormap = [newColormap; repmat(cmap(7,:),4,1)];
newColormap = [newColormap; repmat(cmap(8,:),5,1)];
newColormap = [newColormap; repmat(cmap(9,:),5,1)];
newColormap = [newColormap; repmat(cmap(10,:),10,1)];


figure
ax = axesm('pcarree','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,mean_gpp_cn51_short_gcy_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',newColormap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['mean GPP CN51 (@280ppm) [kgC m^{-2} s^{-1}]'])
h = colorbar;
h.Ticks = [0 100 200 400 800 1200 1600 2000 2500 3500 4000];
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/global/GPP_CN51_short_gcy'],'png')

figure
ax = axesm('pcarree','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,mean_gpp_cn45_short_gcy_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',newColormap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['mean GPP CN45 (@390ppm) [kgC m^{-2} s^{-1}]'])
h = colorbar;
h.Ticks = [0 100 200 400 800 1200 1600 2000 2500 3500 4000];
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/global/GPP_CN45_short_gcy'],'png')

figure
ax = axesm('pcarree','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,mean_gpp_cn45_280_short_gcy_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',newColormap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['mean GPP CN45 (@280ppm) [kgC m^{-2} s^{-1}]'])
h = colorbar;
h.Ticks = [0 100 200 400 800 1200 1600 2000 2500 3500 4000];
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/global/GPP_CN45_280_short_gcy'],'png')

cmap_mat = cmocean('balance',12);
%cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
cmap = cmap_mat; 
s1 = quantile(mean_diff_gpp_cn51_fc_mat(:),0.01); 
s2 = quantile(mean_diff_gpp_cn51_fc_mat(:),0.99);
caxis_min = -1*max(abs(s1),abs(s2));
caxis_max = max(abs(s1),abs(s2));

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,mean_diff_gpp_cn51_fc_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['mean diff. GPP (CN51 - FluxCom) [kgC m^{-2} s^{-1}] avg.: ' num2str(nanmean(mean_diff_gpp_cn51_fc_mat(:)),2)])
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/global/diff_GPP_CN51_FC'],'png')

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on; 
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,mean_diff_gpp_cn45_fc_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['mean diff. GPP (CN45 - FluxCom) [kgC m^{-2} s^{-1}] avg.: ' num2str(nanmean(mean_diff_gpp_cn45_fc_mat(:)),2)])
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/global/diff_GPP_CN45_FC'],'png')

cmap_mat = cmocean('matter',12);
%cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
cmap = cmap_mat;
s1 = quantile(rmse_gpp_cn51_fc_mat(:),0.01);
s2 = quantile(rmse_gpp_cn51_fc_mat(:),0.99);
caxis_min = s1;
caxis_max = s2;

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,rmse_gpp_cn51_fc_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['RMSE GPP CN51 vs. FluxCom [kgC m^{-2} s^{-1}] avg.: ' num2str(nanmean(rmse_gpp_cn51_fc_mat(:)),2)])
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/global/RMSE_GPP_CN51_FC'],'png')

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,rmse_gpp_cn45_fc_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['RMSE GPP CN45 vs. FluxCom [kgC m^{-2} s^{-1}] avg.: ' num2str(nanmean(rmse_gpp_cn45_fc_mat(:)),2)])
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/global/RMSE_GPP_CN45_FC'],'png')

cmap_mat = cmocean('matter',12);
%cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
cmap = cmap_mat;
s1 = quantile(R_gpp_cn51_fc_mat(:),0.01);
s2 = quantile(R_gpp_cn51_fc_mat(:),0.99);
caxis_min = 0;
caxis_max = 1;

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,R_gpp_cn51_fc_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['R GPP CN51 vs. FluxCom [-] avg.: ' num2str(nanmean(R_gpp_cn51_fc_mat(:)),2)])
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/global/R_GPP_CN51_FC'],'png')

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                       
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off;
gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                       
    'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                            
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,R_gpp_cn45_fc_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                             
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title(['R GPP CN45 vs. FluxCom [-] avg.: ' num2str(nanmean(R_gpp_cn45_fc_mat(:)),2)])
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/global/R_GPP_CN45_FC'],'png')

% RMSE per PFT

for i = 1:length(pft_name)

    pft_index = find(dom_pft==i);

    rmse_gpp_cn51_fc_pft(i) = nanmean(rmse_gpp_cn51_fc(pft_index));
    rmse_gpp_cn45_fc_pft(i) = nanmean(rmse_gpp_cn45_fc(pft_index));

end

rmse_gpp_fc_pft_plot_mat = [rmse_gpp_cn51_fc_pft; rmse_gpp_cn45_fc_pft];

X = categorical(pft_name_short);
X = reordercats(X,pft_name_short);

figure
bar(X,rmse_gpp_fc_pft_plot_mat)
ylabel('GPP RMSE [kg m^{-2} s^{-1}]')
legend('CN51','CN45')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/RMSE_GPP_CN51_CN45_FC_PFT'],'png')
