clc
clear all
close all

ghg_path = '/discover/nobackup/projects/geoscm/fzeng/Catchment-CN40_9km/GHG_center_data_milan/';
ghg_exp_name = 'GEOSldas_CN40_9km';

ref_path = '/discover/nobackup/projects/geoscm/fzeng/Catchment-CN40_9km/GEOSldas_1981_present/';
ref_exp_name = 'GEOSldas_CNCLM40_9km';

start_year = 1981;
stop_year = 1982;

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

mat_index = 0;

for y = start_year:stop_year

    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

    for m = 1:12

        mat_index = mat_index + 1;

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

        ghg_read_file = [ghg_path ghg_exp_name '/output/SMAP_EASEv2_M09_GLOBAL/cat/ens0000/Y' y_str ...
                         '/M' m_str '/' ghg_exp_name '.tavg24_1d_lnd_Nt.monthly.' y_str m_str '.nc4'];

        ref_read_file = [ref_path ref_exp_name '/output/SMAP_EASEv2_M09_GLOBAL/cat/ens0000/Y' y_str ...
                         '/M' m_str '/' ref_exp_name '.tavg24_1d_lnd_Nt.monthly.' y_str m_str '.nc4'];

        if (mat_index==1)

           lat_vec = ncread(ghg_read_file,'lat');
           lon_vec = ncread(ghg_read_file,'lon');  

           ghg_cngpp_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           ghg_cnnee_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           ghg_cncloss_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           ghg_cntlai_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           ghg_cntotc_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           ghg_et_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           ghg_evpsoil_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           ghg_evptrns_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           ghg_evpintr_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           ghg_evpsbln_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);

           ref_cngpp_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           ref_cnnee_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           ref_cncloss_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           ref_cntlai_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           ref_cntotc_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           ref_et_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           ref_evpsoil_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           ref_evptrns_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           ref_evpintr_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           ref_evpsbln_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);

           diff_cngpp_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           diff_cnnee_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           diff_cncloss_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           diff_cntlai_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           diff_cntotc_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           diff_et_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           diff_evpsoil_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           diff_evptrns_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           diff_evpintr_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           diff_evpsbln_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);

       end

       ghg_cngpp_mat(:,mat_index) = double(ncread(ghg_read_file,'CNGPP'));
       ghg_cnnee_mat(:,mat_index) = double(ncread(ghg_read_file,'CNNEE'));
       ghg_cncloss_mat(:,mat_index) = double(ncread(ghg_read_file,'CLOSS'));
       ghg_cntlai_mat(:,mat_index) = double(ncread(ghg_read_file,'CNTLAI'));
       ghg_cntotc_mat(:,mat_index) = double(ncread(ghg_read_file,'CNTOTC'));
       ghg_evpsoil_mat(:,mat_index) = double(ncread(ghg_read_file,'EVPSOIL'));
       ghg_evptrns_mat(:,mat_index) = double(ncread(ghg_read_file,'EVPTRNS'));
       ghg_evpintr_mat(:,mat_index) = double(ncread(ghg_read_file,'EVPINTR'));
       ghg_evpsbln_mat(:,mat_index) = double(ncread(ghg_read_file,'EVPSBLN'));
       ghg_et_mat(:,mat_index) = ghg_evpsoil_mat(:,mat_index) + ghg_evptrns_mat(:,mat_index) ...
                                 + ghg_evpintr_mat(:,mat_index) + ghg_evpsbln_mat(:,mat_index);

       ref_cngpp_mat(:,mat_index) = double(ncread(ref_read_file,'CNGPP'));
       ref_cnnee_mat(:,mat_index) = double(ncread(ref_read_file,'CNNEE'));
       ref_cncloss_mat(:,mat_index) = double(ncread(ref_read_file,'CLOSS'));
       ref_cntlai_mat(:,mat_index) = double(ncread(ref_read_file,'CNTLAI'));
       ref_cntotc_mat(:,mat_index) = double(ncread(ref_read_file,'CNTOTC'));
       ref_evpsoil_mat(:,mat_index) = double(ncread(ref_read_file,'EVPSOIL'));
       ref_evptrns_mat(:,mat_index) = double(ncread(ref_read_file,'EVPTRNS'));
       ref_evpintr_mat(:,mat_index) = double(ncread(ref_read_file,'EVPINTR'));
       ref_evpsbln_mat(:,mat_index) = double(ncread(ref_read_file,'EVPSBLN'));
       ref_et_mat(:,mat_index) = ref_evpsoil_mat(:,mat_index) + ref_evptrns_mat(:,mat_index) ... 
                                 + ref_evpintr_mat(:,mat_index) + ref_evpsbln_mat(:,mat_index);

       diff_cngpp_mat(:,mat_index) = ghg_cngpp_mat(:,mat_index) - ref_cngpp_mat(:,mat_index);
       diff_cnnee_mat(:,mat_index) = ghg_cnnee_mat(:,mat_index) - ref_cnnee_mat(:,mat_index);
       diff_cncloss_mat(:,mat_index) = ghg_cncloss_mat(:,mat_index) - ref_cncloss_mat(:,mat_index);
       diff_cntlai_mat(:,mat_index) = ghg_cntlai_mat(:,mat_index) - ref_cntlai_mat(:,mat_index);
       diff_cntotc_mat(:,mat_index) = ghg_cntotc_mat(:,mat_index) - ref_cntotc_mat(:,mat_index);
       diff_evpsoil_mat(:,mat_index) = ghg_evpsoil_mat(:,mat_index) - ref_evpsoil_mat(:,mat_index);
       diff_evptrns_mat(:,mat_index) = ghg_evptrns_mat(:,mat_index) - ref_evptrns_mat(:,mat_index);
       diff_evpintr_mat(:,mat_index) = ghg_evpintr_mat(:,mat_index) - ref_evpintr_mat(:,mat_index);
       diff_evpsbln_mat(:,mat_index) = ghg_evpsbln_mat(:,mat_index) - ref_evpsbln_mat(:,mat_index);
       diff_et_mat(:,mat_index) = ghg_et_mat(:,mat_index) - ref_et_mat(:,mat_index);


    end
end

mean_diff_cngpp = nanmean(diff_cngpp_mat,2);
mean_diff_cnnee = nanmean(diff_cnnee_mat,2);
mean_diff_cncloss = nanmean(diff_cncloss_mat,2);
mean_diff_cntlai = nanmean(diff_cntlai_mat,2);
mean_diff_cntotc = nanmean(diff_cntotc_mat,2);
mean_diff_evpsoil = nanmean(diff_evpsoil_mat,2);
mean_diff_evptrns = nanmean(diff_evptrns_mat,2);
mean_diff_evpintr = nanmean(diff_evpintr_mat,2);
mean_diff_evpsbln = nanmean(diff_evpsbln_mat,2);
mean_diff_et = nanmean(diff_et_mat,2);

max_diff_cngpp = max(diff_cngpp_mat,[],2);
max_diff_cnnee = max(diff_cnnee_mat,[],2);
max_diff_cncloss = max(diff_cncloss_mat,[],2);
max_diff_cntlai = max(diff_cntlai_mat,[],2);
max_diff_cntotc = max(diff_cntotc_mat,[],2);
max_diff_evpsoil = max(diff_evpsoil_mat,[],2);
max_diff_evptrns = max(diff_evptrns_mat,[],2);
max_diff_evpintr = max(diff_evpintr_mat,[],2);
max_diff_evpsbln = max(diff_evpsbln_mat,[],2);
max_diff_et = max(diff_et_mat,[],2);

