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/'; 
cn51_c3_read_path = '/discover/nobackup/jkolassa/CN51_test_local/GEOSldas_CN51_local_c3/output/SMAP_EASEv2_M36/cat/ens0000/';

c1_start_file = [cn51_read_path 'Y1980/M01/GEOSldas_CN51_local.tavg24_1d_lnd_Nt.monthly.198001.nc4'];
c1_stop_file = [cn51_read_path 'Y2019/M12/GEOSldas_CN51_local.tavg24_1d_lnd_Nt.monthly.201912.nc4'];
c2_start_file = [cn51_c2_read_path 'Y1980/M01/GEOSldas_CN51_local_c2.tavg24_1d_lnd_Nt.monthly.198001.nc4'];
c2_stop_file = [cn51_c2_read_path 'Y2019/M12/GEOSldas_CN51_local_c2.tavg24_1d_lnd_Nt.monthly.201912.nc4'];
c3_start_file = [cn51_c3_read_path 'Y1980/M01/GEOSldas_CN51_local_c3.tavg24_1d_lnd_Nt.monthly.198001.nc4'];

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

c1_start_gpp = double(ncread(c1_start_file,'CNGPP'));
mean_c1_start_gpp = nanmean(c1_start_gpp);
c1_stop_gpp = double(ncread(c1_stop_file,'CNGPP'));
mean_c1_stop_gpp = nanmean(c1_stop_gpp);
c2_start_gpp = double(ncread(c2_start_file,'CNGPP'));
mean_c2_start_gpp = nanmean(c2_start_gpp);
c2_stop_gpp = double(ncread(c2_stop_file,'CNGPP'));
mean_c2_stop_gpp = nanmean(c2_stop_gpp);
c3_start_gpp = double(ncread(c3_start_file,'CNGPP'));
mean_c3_start_gpp = nanmean(c3_start_gpp);

plot_labels = {'C1 start','C1 stop','C2 start','C2 stop','C3 start'};
gpp_plot_vec = [mean_c1_start_gpp mean_c1_stop_gpp mean_c2_start_gpp mean_c2_stop_gpp mean_c3_start_gpp];

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

figure
bar(X,gpp_plot_vec)
ylabel('GPP [kg m^{-2} s^{-1}]')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/transition_GPP_monthly'],'png')

for d = 1:31

    d_str = num2str(d,'%02i');

    c2_start_file_daily = [cn51_c2_read_path 'Y1980/M01/GEOSldas_CN51_local_c2.tavg24_1d_lnd_Nt.198001' d_str '_1200z.nc4'];
    c2_start_gpp_daily(:,d) = double(ncread(c2_start_file_daily,'CNGPP'));

    c1_stop_file_daily = [cn51_read_path 'Y2019/M12/GEOSldas_CN51_local.tavg24_1d_lnd_Nt.201912' d_str '_1200z.nc4'];
    c1_stop_gpp_daily(:,d) = double(ncread(c1_stop_file_daily,'CNGPP'));

end

mean_daily_gpp_c2_start = nanmean(c2_start_gpp_daily,1);
mean_daily_gpp_c1_stop = nanmean(c1_stop_gpp_daily,1);

figure
plot(1:31,mean_daily_gpp_c2_start,'-*r')
hold on
plot(1:31,mean_daily_gpp_c1_stop,'-*b')
hold off
xlabel('day')
ylabel('GPP [kg m^{-2} s^{-1}]')
legend('C2 start daily','C1 stop daily')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/transition_GPP_daily'],'png')
