clc
clear all 
close all

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 = 1981;
stop_year = 2020; 
    
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];

% read data

cn40_read_file = [cn40_read_path 'Y1981/M01/GEOSldas_CNCLM40_9km.tavg24_1d_lnd_Nt.monthly.198101.nc4'];
lat_vec = double(ncread(cn40_read_file,'lat'));
lon_vec = double(ncread(cn40_read_file,'lon'));

mat_index = 0;
mat_index_d = 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');

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


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

        if exist(cn40_read_file)
           cnnee_mat(:,mat_index) = double(ncread(cn40_read_file,'CNNEE'));
        else
           cnnee_mat(:,mat_index) = NaN*ones(size(cnnee_mat,1),1);
        end

    end % m
end % y

ann_index = 0;

for y = start_year:stop_year

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

    tmp = cnnee_mat(:,year_index);
    mean_global_ann_cnnee(ann_index) = nanmean(tmp(:));
    mean_local_ann_cnnee(:,ann_index) = nanmean(tmp,2);

end 

figure
plot(start_year:stop_year,mean_global_ann_cnnee,'LineWidth',2)
grid on
xlabel('year')
ylabel('NEE [kg m^{-2} s^{-1}]')
title('annual average NEE [kg m^{-2} s^{-1}]')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN40_ann_glob_NEE'],'png')


