clc
clear all
close all

cn51_read_path = '/discover/nobackup/projects/geoscm/fzeng/Catchment-CN5.1/GEOSldas_CN51/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/';

start_year = 1980;
stop_year = 1980; 

rst_file = '/discover/nobackup/projects/geoscm/fzeng/Catchment-CN5.1/GEOSldas_CN51/output/SMAP_EASEv2_M36/rs/ens0000/Y1980/M01/GEOSldas_CN51.catchcnclm51_internal_rst.19800101_0000';
tile_id = ncread(rst_file,'TILE_ID');
ityp = ncread(rst_file,'ITY');
fvg = ncread (rst_file,'FVG');

clm45_rst_file = '/discover/nobackup/projects/geoscm/fzeng/GEOSldas45_M36_rst_ldas06_15/output/SMAP_EASEv2_M36/rs/ens0000/Y1980/M01/GEOSldas45_M36_rst_ldas06_15.catchcnclm45_internal_rst.19800101_0000';
clm45_tile_id = ncread(clm45_rst_file,'TILE_ID');
clm45_ityp = ncread(clm45_rst_file,'ITY');
clm45_fvg = ncread (clm45_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 

for i = 1:length(clm45_ityp)
    [max_value,max_index] = max(clm45_fvg(i,:));
    if (max_value>=0.7)
       clm45_dom_pft(i) = clm45_ityp(i,max_index);
    else
       clm45_dom_pft(i) = NaN;
    end 
end

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


for i = 1:length(pft_name_cn51)

    if (i==10)
       pft_map_index{i} = [10 11];
    elseif (i==11)
       pft_map_index{i} = 12;
    elseif (i==12)
       pft_map_index{i} = 13;
    elseif (i==13)
       pft_map_index{i} = [14 15];
    elseif (i==14)
       pft_map_index{i} = [16 17];
    elseif (i==15)
       pft_map_index{i} = [18 19];
    else
       pft_map_index{i} = i;
    end
    
end

pft_name_short_cn51 = {...
           '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 'Y2001/M01/GEOSldas_CN51.tavg24_1d_lnd_Nt.monthly.200101.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'];

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


for i = 1:length(lat_vec)

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

end

% read data

mat_index = 0;
y_index = 0;

for y = start_year:stop_year

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

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

    % loop over months
    for m = 1:12

        mat_index = mat_index + 1;

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

        cn51_read_file = [cn51_read_path 'Y' y_str '/M' m_str '/GEOSldas_CN51.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'];

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

           cnnee_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           clm45_cnnee_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);

           cnrs_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           clm45_cnrs_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);

           cncloss_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           clm45_cncloss_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);

        end

        cngpp_mat(:,mat_index) = double(ncread(cn51_read_file,'CNGPP'));
        cnnee_mat(:,mat_index) = double(ncread(cn51_read_file,'CNNEE'));      
        cnrs_mat(:,mat_index) = double(ncread(cn51_read_file,'CNSR'));
        cncloss_mat(:,mat_index) = double(ncread(cn51_read_file,'CLOSS'));

        tmp_mat = double(ncread(clm45_read_file,'CNGPP'));
        clm45_cngpp_mat(:,mat_index) = tmp_mat(data_index);
  
        tmp_mat = double(ncread(clm45_read_file,'CNNEE'));
        clm45_cnnee_mat(:,mat_index) = tmp_mat(data_index);

        tmp_mat = double(ncread(clm45_read_file,'CNSR'));
        clm45_cnrs_mat(:,mat_index) = tmp_mat(data_index);

        tmp_mat = double(ncread(clm45_read_file,'CLOSS'));
        clm45_cncloss_mat(:,mat_index) = tmp_mat(data_index);

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

end

% remove extra entries

cngpp_mat(:,(mat_index+1):end) = [];
clm45_cngpp_mat(:,(mat_index+1):end) = [];

cnnee_mat(:,(mat_index+1):end) = [];
clm45_cnnee_mat(:,(mat_index+1):end) = [];

cnrs_mat(:,(mat_index+1):end) = [];
clm45_cnrs_mat(:,(mat_index+1):end) = [];

cncloss_mat(:,(mat_index+1):end) = [];
clm45_cncloss_mat(:,(mat_index+1):end) = [];

% compute differences

diff_gpp_cn51_cn45 = cngpp_mat - clm45_cngpp_mat;
mean_diff_gpp_cn51_cn45 = nanmean(diff_gpp_cn51_cn45,2);
diff_nee_cn51_cn45 = cnnee_mat - clm45_cnnee_mat;
mean_diff_nee_cn51_cn45 = nanmean(diff_nee_cn51_cn45,2);
diff_rs_cn51_cn45 = cnrs_mat - clm45_cnrs_mat;
mean_diff_rs_cn51_cn45 = nanmean(diff_rs_cn51_cn45,2);
diff_closs_cn51_cn45 = cncloss_mat - clm45_cncloss_mat;
mean_diff_closs_cn51_cn45 = nanmean(diff_closs_cn51_cn45,2);

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 = 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 = cnnee_mat(:,year_index);
    mean_global_ann_cnnee(ann_index) = nanmean(tmp(:));
    mean_local_ann_cnnee(:,ann_index) = nanmean(tmp,2);

    tmp = clm45_cnnee_mat(:,year_index);
    mean_global_ann_clm45_cnnee(ann_index) = nanmean(tmp(:));
    mean_local_ann_clm45_cnnee(:,ann_index) = nanmean(tmp,2);

    tmp = cnrs_mat(:,year_index);
    mean_global_ann_cnrs(ann_index) = nanmean(tmp(:));
    mean_local_ann_cnrs(:,ann_index) = nanmean(tmp,2);

    tmp = clm45_cnrs_mat(:,year_index);
    mean_global_ann_clm45_cnrs(ann_index) = nanmean(tmp(:));
    mean_local_ann_clm45_cnrs(:,ann_index) = nanmean(tmp,2);

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

    tmp = clm45_cncloss_mat(:,year_index);
    mean_global_ann_clm45_cncloss(ann_index) = nanmean(tmp(:));
    mean_local_ann_clm45_cncloss(:,ann_index) = nanmean(tmp,2);

end

figure
plot(start_year:stop_year,mean_global_ann_cnnee,'LineWidth',2);
hold on
plot(start_year:stop_year,mean_global_ann_clm45_cnnee,'LineWidth',2);
hold off
grid on
xlabel('year')
ylabel('NEE [kg m^{-2} s^{-1}]')
title('monthly average NEE [kg m^{-2} s^{-1}] (1980 - 2020)')
legend('CN51','CN45','Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/NEE_annual_CN51_CN45'],'png')

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 off
grid on
xlabel('year')
ylabel('GPP [kg m^{-2} s^{-1}]')
title('monthly average GPP [kg m^{-2} s^{-1}] (1980 - 2020)')
legend('CN51','CN45','Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/GPP_annual_CN51_CN45'],'png')

figure
plot(start_year:stop_year,mean_global_ann_cnrs,'LineWidth',2);
hold on
plot(start_year:stop_year,mean_global_ann_clm45_cnrs,'LineWidth',2);
hold off
grid on
xlabel('year')
ylabel('RS [kg m^{-2} s^{-1}]')
title('monthly average RS [kg m^{-2} s^{-1}] (1980 - 2020)')
legend('CN51','CN45','Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/RS_annual_CN51_CN45'],'png')

figure
plot(start_year:stop_year,mean_global_ann_cncloss,'LineWidth',2);
hold on
plot(start_year:stop_year,mean_global_ann_clm45_cncloss,'LineWidth',2);
hold off
grid on
xlabel('year')
ylabel('CLOSS [kg m^{-2} s^{-1}]')
title('monthly average CLOSS [kg m^{-2} s^{-1}] (1980 - 2020)')
legend('CN51','CN45','Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/CLOSS_annual_CN51_CN45'],'png')


for m = 1:12

    month_index = (month_vec==m);

    tmp = cngpp_mat(:,month_index);
    mean_global_mon_cngpp(m) = nanmean(tmp(:));
    mean_local_mon_cngpp(:,m) = nanmean(tmp,2);

    tmp = clm45_cngpp_mat(:,month_index);
    mean_global_mon_clm45_cngpp(m) = nanmean(tmp(:));
    mean_local_mon_clm45_cngpp(:,m) = nanmean(tmp,2);

    tmp = cnnee_mat(:,month_index);
    mean_global_mon_cnnee(m) = nanmean(tmp(:));
    mean_local_mon_cnnee(:,m) = nanmean(tmp,2);

    tmp = clm45_cnnee_mat(:,month_index);
    mean_global_mon_clm45_cnnee(m) = nanmean(tmp(:));
    mean_local_mon_clm45_cnnee(:,m) = nanmean(tmp,2);

    tmp = cnrs_mat(:,month_index);
    mean_global_mon_cnrs(m) = nanmean(tmp(:));
    mean_local_mon_cnrs(:,m) = nanmean(tmp,2);

    tmp = clm45_cnrs_mat(:,month_index);
    mean_global_mon_clm45_cnrs(m) = nanmean(tmp(:));
    mean_local_mon_clm45_cnrs(:,m) = nanmean(tmp,2);

    tmp = cncloss_mat(:,month_index);
    mean_global_mon_cncloss(m) = nanmean(tmp(:));
    mean_local_mon_cncloss(:,m) = nanmean(tmp,2);

    tmp = clm45_cncloss_mat(:,month_index);
    mean_global_mon_clm45_cncloss(m) = nanmean(tmp(:));
    mean_local_mon_clm45_cncloss(:,m) = nanmean(tmp,2);
end

figure
plot(1:12,mean_global_mon_cnnee,'LineWidth',2);
hold on
plot(1:12,mean_global_mon_clm45_cnnee,'LineWidth',2);
hold off
grid on
xlabel('year')
ylabel('NEE [kg m^{-2} s^{-1}]')
title('monthly average NEE [kg m^{-2} s^{-1}] (1980 - 2020)')
legend('CN51','CN45','Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/NEE_monthly_CN51_CN45'],'png')

figure
plot(1:12,mean_global_mon_cngpp,'LineWidth',2);
hold on
plot(1:12,mean_global_mon_clm45_cngpp,'LineWidth',2);
hold off 
grid on
xlabel('year')
ylabel('GPP [kg m^{-2} s^{-1}]')
title('monthly average GPP [kg m^{-2} s^{-1}] (1980 - 2020)')
legend('CN51','CN45','Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/GPP_monthly_CN51_CN45'],'png')

figure
plot(1:12,mean_global_mon_cnrs,'LineWidth',2);
hold on
plot(1:12,mean_global_mon_clm45_cnrs,'LineWidth',2);
hold off
grid on
xlabel('year')
ylabel('RS [kg m^{-2} s^{-1}]')
title('monthly average RS [kg m^{-2} s^{-1}] (1980 - 2020)')
legend('CN51','CN45','Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/RS_monthly_CN51_CN45'],'png')

figure
plot(1:12,mean_global_mon_cncloss,'LineWidth',2);
hold on
plot(1:12,mean_global_mon_clm45_cncloss,'LineWidth',2);
hold off
grid on
xlabel('year')
ylabel('CLOSS [kg m^{-2} s^{-1}]')
title('monthly average CLOSS [kg m^{-2} s^{-1}] (1980 - 2020)')
legend('CN51','CN45','Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/CLOSS_monthly_CN51_CN45'],'png')

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

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

mean_diff_gpp_cn51_cn45_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
mean_diff_nee_cn51_cn45_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
mean_diff_rs_cn51_cn45_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));
mean_diff_closs_cn51_cn45_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2));


for i = 1:length(lon_vec)

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

    mean_diff_gpp_cn51_cn45_mat(lat_index,lon_index) = mean_diff_gpp_cn51_cn45(i);
    mean_diff_nee_cn51_cn45_mat(lat_index,lon_index) = mean_diff_nee_cn51_cn45(i);
    mean_diff_rs_cn51_cn45_mat(lat_index,lon_index) = mean_diff_rs_cn51_cn45(i);
    mean_diff_closs_cn51_cn45_mat(lat_index,lon_index) = mean_diff_closs_cn51_cn45(i);
end

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

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

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_cn45_mat(:),0.01);
s2 = quantile(mean_diff_gpp_cn51_cn45_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_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 diff. GPP (CN51 - CN45) '])
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_CN45'],'png')

s1 = quantile(mean_diff_nee_cn51_cn45_mat(:),0.01);
s2 = quantile(mean_diff_nee_cn51_cn45_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_nee_cn51_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 diff. NEE (CN51 - CN45) '])
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_NEE_CN51_CN45'],'png')

s1 = quantile(mean_diff_rs_cn51_cn45_mat(:),0.01);
s2 = quantile(mean_diff_rs_cn51_cn45_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_rs_cn51_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 diff. RS (CN51 - CN45) '])
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_RS_CN51_CN45'],'png')

s1 = quantile(mean_diff_gpp_cn51_cn45_mat(:),0.01);
s2 = quantile(mean_diff_gpp_cn51_cn45_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_closs_cn51_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 diff. CLOSS (CN51 - CN45) '])
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_CLOSS_CN51_CN45'],'png')

% RMSE per PFT

for i = 1:length(pft_name_cn51)

    pft_index = find(dom_pft==i);
    if (length(pft_map_index{i})>1)
       pft_index_cn45 = find((clm45_dom_pft==pft_map_index{i}(1)) | (clm45_dom_pft==pft_map_index{i}(2)));
    else
       pft_index_cn45 = find(clm45_dom_pft==pft_map_index{i});
    end

    tmp = cngpp_mat(pft_index,:);
    mean_cngpp_pft(i) = nanmean(tmp(:));
    tmp = clm45_cngpp_mat(pft_index_cn45,:);
    mean_clm45_cngpp_pft(i) = nanmean(tmp(:));

    tmp = cnnee_mat(pft_index,:);
    mean_cnnee_pft(i) = nanmean(tmp(:));
    tmp = clm45_cnnee_mat(pft_index_cn45,:);
    mean_clm45_cnnee_pft(i) = nanmean(tmp(:));

    tmp = cnrs_mat(pft_index,:);
    mean_cnrs_pft(i) = nanmean(tmp(:));
    tmp = clm45_cnrs_mat(pft_index_cn45,:);
    mean_clm45_cnrs_pft(i) = nanmean(tmp(:));

    tmp = cncloss_mat(pft_index,:);
    mean_cncloss_pft(i) = nanmean(tmp(:));
    tmp = clm45_cncloss_mat(pft_index_cn45,:);
    mean_clm45_cncloss_pft(i) = nanmean(tmp(:));

end

gpp_pft_plot_mat = [mean_cngpp_pft; mean_clm45_cngpp_pft];
nee_pft_plot_mat = [mean_cnnee_pft; mean_clm45_cnnee_pft];
rs_pft_plot_mat = [mean_cnrs_pft; mean_clm45_cnrs_pft];
closs_pft_plot_mat = [mean_cncloss_pft; mean_clm45_cncloss_pft];

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

figure
bar(X,gpp_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/GPP_CN51_CN45_PFT'],'png')

figure
bar(X,nee_pft_plot_mat)
ylabel('NEE 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/NEE_CN51_CN45_PFT'],'png')

figure
bar(X,rs_pft_plot_mat)
ylabel('RS 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/RS_CN51_CN45_PFT'],'png')

figure
bar(X,closs_pft_plot_mat)
ylabel('CLOSS 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/CLOSS_CN51_CN45_PFT'],'png')

