clc
clear all
close all

cn51_read_path = '/discover/nobackup/jkolassa/CN51_test_local/GEOSldas_CN51_local/output/SMAP_EASEv2_M36/cat/ens0000/';
cn51_c2_read_path = '/discover/nobackup/jkolassa/CN51_test_local/GEOSldas_CN51_local_c2/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/';
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/';

start_year = 2001;
stop_year = 2016; 

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_local.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'];
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'];

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

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

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

    dist_vec = sqrt((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);
    fs_lon_index(i) = find(lon_dist==min(lon_dist));

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

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

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

    lon_dist = sqrt((vpm_lon_vec-lon_vec(i)).^2);
    vpm_lon_index(i) = find(lon_dist==min(lon_dist));
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_local.tavg24_1d_lnd_Nt.monthly.' y_str m_str '.nc4'];
        cn51_c2_read_file = [cn51_c2_read_path 'Y' y_str '/M' m_str '/GEOSldas_CN51_local_c2.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'];
        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'];

        if (mat_index == 1)

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

        end

        cngpp_mat(:,mat_index) = double(ncread(cn51_read_file,'CNGPP'));
        cngpp_c2_mat(:,mat_index) = double(ncread(cn51_c2_read_file,'CNGPP'));      

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

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

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 = cngpp_c2_mat(:,year_index);
    mean_global_ann_cngpp_c2(ann_index) = nanmean(tmp(:));
    mean_local_ann_cngpp_c2(:,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 = 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);

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_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 off
grid on
xlabel('year')
ylabel('GPP [kg m^{-2} s^{-1}]')
title('annual average GPP [kg m^{-2} s^{-1}]')
legend('CN51','CN45','FluxCom','FluxSat','MAIAC','VPM','Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/CN51_eval_Fluxcom_GPP'],'png')

figure
plot(start_year:stop_year,mean_global_ann_cngpp,'LineWidth',2);
hold on
plot(start_year:stop_year,mean_global_ann_cngpp_c2,'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 off
grid on
xlabel('year')
ylabel('GPP [kg m^{-2} s^{-1}]')
title('annual average GPP [kg m^{-2} s^{-1}]')
legend('CN51 cycle 1','CN51 cycle 2','FluxCom','FluxSat','MAIAC','VPM','Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/GPP_ann_CN51_cycles'],'png')

%amazon_index = find(tile_id==21274);
amazon_index = 97;

figure
plot(start_year:stop_year,mean_local_ann_cngpp(amazon_index,:),'LineWidth',2);
hold on
plot(start_year:stop_year,mean_local_ann_cngpp_c2(amazon_index,:),'LineWidth',2);
hold on
%plot(start_year:stop_year,mean_global_ann_cn40_gpp,'LineWidth',2);
%hold on
plot(start_year:stop_year,mean_local_ann_fc_gpp(amazon_index,:),'LineWidth',2,'LineStyle','--');
hold on
plot(start_year:stop_year,mean_local_ann_fs_gpp(amazon_index,:),'LineWidth',2,'LineStyle','--');
hold on
plot(start_year:stop_year,mean_local_ann_maiac_gpp(amazon_index,:),'LineWidth',2,'LineStyle','--');
hold on
plot(start_year:stop_year,mean_local_ann_vpm_gpp(amazon_index,:),'LineWidth',2,'LineStyle','--');
hold off
grid on
xlabel('year')
ylabel('GPP [kg m^{-2} s^{-1}]')
title('annual average GPP [kg m^{-2} s^{-1}]')
legend('CN51 cycle 1','CN51 cycle 2','FluxCom','FluxSat','MAIAC','VPM','Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/GPP_ann_CN51_cycles_amazon'],'png')

tmp = start_year:stop_year;
x_vec = [tmp tmp];
mean_local_ann_cngpp_all = [mean_local_ann_cngpp mean_local_ann_cngpp_c2];
mean_global_ann_cngpp_all = [mean_global_ann_cngpp mean_global_ann_cngpp_c2];

figure
plot(x_vec,mean_local_ann_cngpp_all(amazon_index,:),'LineWidth',2);
hold off
grid on
xlabel('year')
ylabel('GPP [kg m^{-2} s^{-1}]')
title('annual average GPP [kg m^{-2} s^{-1}]')
legend('CN51 cycle 1','Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/GPP_ann_CN51_cycles_amazon_seq'],'png')

figure
plot(x_vec,mean_global_ann_cngpp_all,'LineWidth',2);
grid on
xlabel('year')
ylabel('GPP [kg m^{-2} s^{-1}]')
title('annual average GPP [kg m^{-2} s^{-1}]')
legend('CN51 cycle 1','Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/GPP_ann_CN51_cycles_seq'],'png')
