clc
clear all
close all

cn_read_path = '/discover/nobackup/jkolassa/test_dir/GEOSldas_CN51/output/SMAP_EASEv2_M36/cat/ens0000/';
clm_read_path = '/discover/nobackup/projects/geoscm/fzeng/CLM/archive/clm5_test/lnd/hist/';

start_year = 1980;
stop_year = 1980;


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

    % loop over months
    for m = 1:12

        mat_index = mat_index + 1;

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

        cn_read_file = [cn_read_path 'Y' y_str '/M' m_str '/GEOSldas_CN51.tavg24_1d_lnd_Nt.monthly.' y_str m_str '.nc4'];
        clm_read_file = [clm_read_path 'clm5_test.clm2.h0.0001-' m_str '.nc'];

        if (mat_index == 1)
           lat_vec = double(ncread(cn_read_file,'lat'));
           lon_vec = double(ncread(cn_read_file,'lon'));

           clm_lat_data = ncread(clm_read_file,'lat');
           clm_lon_data = ncread(clm_read_file,'lon');
           clm_lon_data(clm_lon_data>180) = clm_lon_data(clm_lon_data>180) - 360;

           cntlai_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           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);
           cnrootc_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           cnvegc_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           cntotc_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           sfmc_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           rzmc_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);

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

           for i = 1:length(lat_vec)
               lat_dist = sqrt((clm_lat_data-lat_vec(1)).^2);
               lat_index = find(lat_dist==min(lat_dist));
               lon_dist = sqrt((clm_lon_data-lon_vec(1)).^2);
               lon_index = find(lon_dist==min(lon_dist));
           end
        end

        cntlai_mat(:,mat_index) = double(ncread(cn_read_file,'CNTLAI'));
        cngpp_mat(:,mat_index) = double(ncread(cn_read_file,'CNGPP'));
        cnnee_mat(:,mat_index) = double(ncread(cn_read_file,'CNNEE'));
        cnroot_mat(:,mat_index) = double(ncread(cn_read_file,'CNROOT'));
        cnvegc_mat(:,mat_index) = double(ncread(cn_read_file,'CNVEGC'));
        cntotc_mat(:,mat_index) = double(ncread(cn_read_file,'CNTOTC'));
        sfmc_mat(:,mat_index) = double(ncread(cn_read_file,'SFMC'));
        rzmc_mat(:,mat_index) = double(ncread(cn_read_file,'RZMC'));

        tmp_mat = ncread(clm_read_file,'GPP');
        clmgpp_mat(:,m) = tmp_mat(lon_index,lat_index);
        tmp_mat_1 = ncread(clm_read_file,'LAISHA');
        tmp_mat_2 = ncread(clm_read_file,'LAISUN');
        clmtlai_mat(:,m) = tmp_mat_1(lon_index,lat_index) + tmp_mat_2(lon_index,lat_index);


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

cntlai_mat(:,mat_index+1:end) = [];
cngpp_mat(:,mat_index+1:end) = [];
cnnee_mat(:,mat_index+1:end) = [];
cnroot_mat(:,mat_index+1:end) = [];
cnvegc_mat(:,mat_index+1:end) = [];
cntotc_mat(:,mat_index+1:end) = [];
sfmc_mat(:,mat_index+1:end) = [];
rzmc_mat(:,mat_index+1:end) = [];

% convert CLM GPP from gc/m^2/s to kgC/m^2/s

clmgpp_mat = clmgpp_mat*10^(-3);
