clc
clear all
close all

cn51_read_path = '/discover/nobackup/jkolassa/CN51_test_local/GEOSldas_CN51_local/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/projects/geoscm/fzeng/Fluxcom/';
alexi_read_path = '/discover/nobackup/projects/geoscm/fzeng/ALEXI/monthly/';

start_year = 2001;
stop_year = 2019; 

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_et_read_file = [fluxcom_read_path 'LE.RS.EBC-ALL.MLM-ALL.METEO-NONE.720_360.monthly.2001.nc'];
alexi_et_read_file = [alexi_read_path 'ALEXI_ET_201501.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_et_read_file,'lat');
fc_lon_vec = ncread(fluxcom_et_read_file,'lon');

alexi_lat_vec = ncread(alexi_et_read_file,'latitude');
alexi_lon_vec = ncread(alexi_et_read_file,'longitude');

for i = 1:length(lat_vec)
    lat_dist = sqrt((fc_lat_vec-lat_vec(i)).^2);
    lon_dist = sqrt((fc_lon_vec-lon_vec(i)).^2);
    if (min(lat_dist)<=0.25) & (min(lon_dist)<=0.25)
       fc_lat_index(i) = find(lat_dist==min(lat_dist));
       fc_lon_index(i) = find(lon_dist==min(lon_dist));
    end

    lat_dist = sqrt((alexi_lat_vec-lat_vec(i)).^2);
    lon_dist = sqrt((alexi_lon_vec-lon_vec(i)).^2);
    if (min(lat_dist)<=0.25) & (min(lon_dist)<=0.25)
       alexi_lat_index(i) = find(lat_dist==min(lat_dist));
       alexi_lon_index(i) = find(lon_dist==min(lon_dist));
    end


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

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'];
        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_et_read_file = [fluxcom_read_path 'LE.RS.EBC-ALL.MLM-ALL.METEO-NONE.720_360.monthly.' y_str '.nc'];
        fluxcom_rn_read_file = [fluxcom_read_path 'Rn.RS.EBC-NONE.MLM-ALL.METEO-NONE.720_360.monthly.' y_str '.nc'];

        alexi_et_read_file = [alexi_read_path 'ALEXI_ET_' y_str m_str '.nc'];

        cn51_forc_file = [cn51_read_path 'Y' y_str '/M' m_str '/GEOSldas_CN51_local.tavg24_1d_lfs_Nt.monthly.' y_str m_str '.nc4'];

        if (mat_index == 1)

           cn51_et_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           cn45_et_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           cn40_et_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);

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

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

           cn51_precip_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           
           cn51_lh_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           cn51_ln_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);           
           cn51_sn_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           cn51_rn_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);

        end

        evpsoil_tmp_mat = double(ncread(cn51_read_file,'EVPSOIL'));
        evptrns_tmp_mat = double(ncread(cn51_read_file,'EVPTRNS'));
        evpintr_tmp_mat = double(ncread(cn51_read_file,'EVPINTR'));      
        evpsbln_tmp_mat = double(ncread(cn51_read_file,'EVPSBLN'));
        cn51_et_mat(:,mat_index) = evpsoil_tmp_mat + evptrns_tmp_mat + ...
                                   evpintr_tmp_mat + evpsbln_tmp_mat;

        cn51_precip_mat(:,mat_index) = double(ncread(cn51_forc_file,'RainfSnowf'));
        cn51_lh_mat(:,mat_index) = double(ncread(cn51_read_file,'LHLAND'));
        cn51_ln_mat(:,mat_index) = double(ncread(cn51_read_file,'LWLAND'));
        cn51_sn_mat(:,mat_index) = double(ncread(cn51_read_file,'SWLAND'));
        cn51_rn_mat(:,mat_index) = cn51_ln_mat(:,mat_index) + cn51_sn_mat(:,mat_index);

        if (y<=2016)
           evpsoil_tmp_mat = double(ncread(clm45_read_file,'EVPSOIL'));
           evptrns_tmp_mat = double(ncread(clm45_read_file,'EVPTRNS'));
           evpintr_tmp_mat = double(ncread(clm45_read_file,'EVPINTR'));    
           evpsbln_tmp_mat = double(ncread(clm45_read_file,'EVPSBLN'));
           cn45_et_mat(:,mat_index) = evpsoil_tmp_mat(data_index) + evptrns_tmp_mat(data_index) + ... 
                                   evpintr_tmp_mat(data_index) + evpsbln_tmp_mat(data_index);
        end        

        evpsoil_tmp_mat = double(ncread(cn40_read_file,'EVPSOIL'));
        evptrns_tmp_mat = double(ncread(cn40_read_file,'EVPTRNS'));
        evpintr_tmp_mat = double(ncread(cn40_read_file,'EVPINTR'));
        evpsbln_tmp_mat = double(ncread(cn40_read_file,'EVPSBLN'));
        cn40_et_mat(:,mat_index) = evpsoil_tmp_mat(cn40_data_index) + evptrns_tmp_mat(cn40_data_index) + ... 
                                   evpintr_tmp_mat(cn40_data_index) + evpsbln_tmp_mat(cn40_data_index);
  
        if (y>=2015)
           alexi_et_tmp_mat = ncread(alexi_et_read_file,'ET');

           for i = 1:length(lat_vec)
               if (alexi_lat_index(i)~=0)
                  alexi_et_mat(i,mat_index) = alexi_et_tmp_mat(alexi_lat_index(i),alexi_lon_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 (y<=2015)

       et_tmp_mat = ncread(fluxcom_et_read_file,'LE');
       et_tmp_mat(et_tmp_mat==-9999) = NaN;

       rn_tmp_mat = ncread(fluxcom_rn_read_file,'Rn');
       rn_tmp_mat(rn_tmp_mat==-9999) = NaN;

       for i = 1:length(lat_vec)

           start_index = (y_index-1)*12 + 1;
           stop_index = y_index*12;

           fc_et_mat(i,start_index:stop_index) = et_tmp_mat(fc_lon_index(i),fc_lat_index(i),:);
           fc_rn_mat(i,start_index:stop_index) = rn_tmp_mat(fc_lon_index(i),fc_lat_index(i),:);

       end
    end
end

% convert Fluxcom variables from MJ m^-2 day^-1 to W m^-2

fc_et_mat = fc_et_mat./86400;
fc_et_mat = fc_et_mat*1000000;
fc_rn_mat = fc_rn_mat./86400;
fc_rn_mat = fc_rn_mat*1000000;

% convert ALEXI ET data from mm/day to W m^-2

alexi_et_mm_mat = alexi_et_mat;
alexi_et_mat = alexi_et_mat./(0.0864*0.408);

% convert CN51 precipitation data to mm/day and create ET matrix in mm/day

cn51_precip_mat = cn51_precip_mat * 86400;
cn51_et_mm_mat = cn51_et_mat*0.0864*0.408;
cn45_et_mm_mat = cn45_et_mat*0.0864*0.408;
fc_et_mm_mat = fc_et_mat*0.0864*0.408;

ann_index = 0;

for y = start_year:stop_year

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

    tmp = cn51_et_mat(:,year_index);
    mean_global_ann_cn51_et(ann_index) = nanmean(tmp(:));
    mean_local_ann_cn51_et(:,ann_index) = nanmean(tmp,2);

    tmp = cn51_lh_mat(:,year_index);
    mean_global_ann_cn51_lh(ann_index) = nanmean(tmp(:));
    mean_local_ann_cn51_lh(:,ann_index) = nanmean(tmp,2);

    tmp = cn51_rn_mat(:,year_index);
    mean_global_ann_cn51_rn(ann_index) = nanmean(tmp(:));
    mean_local_ann_cn51_rn(:,ann_index) = nanmean(tmp,2);

    tmp = cn40_et_mat(:,year_index);
    mean_global_ann_cn40_et(ann_index) = nanmean(tmp(:));
    mean_local_ann_cn40_et(:,ann_index) = nanmean(tmp,2);

    tmp = fc_et_mat(:,year_index);
    mean_global_ann_fc_et(ann_index) = nanmean(tmp(:));
    mean_local_ann_fc_et(:,ann_index) = nanmean(tmp,2);

    tmp = fc_rn_mat(:,year_index);
    mean_global_ann_fc_rn(ann_index) = nanmean(tmp(:));
    mean_local_ann_fc_rn(:,ann_index) = nanmean(tmp,2);

    tmp = alexi_et_mat(:,year_index);
    mean_global_ann_alexi_et(ann_index) = nanmean(tmp(:));
    mean_local_ann_alexi_et(:,ann_index) = nanmean(tmp,2);

    tmp = cn45_et_mat(:,year_index);
    mean_global_ann_cn45_et(ann_index) = nanmean(tmp(:));
    mean_local_ann_cn45_et(:,ann_index) = nanmean(tmp,2);

    tmp = cn51_precip_mat(:,year_index);
    mean_global_ann_cn51_precip(ann_index) = nanmean(tmp(:));
    mean_local_ann_cn51_precip(:,ann_index) = nanmean(tmp,2);

    tmp = cn51_et_mm_mat(:,year_index);
    mean_global_ann_cn51_et_mm(ann_index) = nanmean(tmp(:));
    mean_local_ann_cn51_et_mm(:,ann_index) = nanmean(tmp,2);

    tmp = cn45_et_mm_mat(:,year_index);
    mean_global_ann_cn45_et_mm(ann_index) = nanmean(tmp(:));
    mean_local_ann_cn45_et_mm(:,ann_index) = nanmean(tmp,2);

    tmp = fc_et_mm_mat(:,year_index);
    mean_global_ann_fc_et_mm(ann_index) = nanmean(tmp(:));
    mean_local_ann_fc_et_mm(:,ann_index) = nanmean(tmp,2);

    tmp = alexi_et_mm_mat(:,year_index);
    mean_global_ann_alexi_et_mm(ann_index) = nanmean(tmp(:));
    mean_local_ann_alexi_et_mm(:,ann_index) = nanmean(tmp,2);

end

% ET plots

figure
plot(start_year:stop_year,mean_global_ann_cn51_et,'LineWidth',2);
hold on
plot(start_year:stop_year,mean_global_ann_cn45_et,'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_et,'LineWidth',2,'LineStyle','--');
hold on 
plot(start_year:stop_year,mean_global_ann_alexi_et,'LineWidth',2,'LineStyle','--');
hold off
grid on
xlabel('year')
ylabel('ET [W m^{-2}]')
title('annual average ET [W m^{-2}]')
legend('CN51','CN45','FluxCom','ALEXI','Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 30 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/CN51_eval_Fluxcom_ET'],'png')

figure
plot(start_year:stop_year,mean_global_ann_cn51_et_mm,'LineWidth',2);
hold on
plot(start_year:stop_year,mean_global_ann_cn45_et_mm,'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_et_mm,'LineWidth',2,'LineStyle','--');
hold on
plot(start_year:stop_year,mean_global_ann_alexi_et_mm,'LineWidth',2,'LineStyle','--');
hold off
grid on
xlabel('year')
ylabel('ET [mm day^{-1}]')
title('annual average ET [mm day^{-1}]')
legend('CN51','CN45','FluxCom','ALEXI','Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 30 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/CN51_eval_Fluxcom_ETmm'],'png')

cn51_diff_et_precip = mean_global_ann_cn51_et_mm - mean_global_ann_cn51_precip;
cn45_diff_et_precip = mean_global_ann_cn45_et_mm - mean_global_ann_cn51_precip;
fc_diff_et_precip = mean_global_ann_fc_et_mm - mean_global_ann_cn51_precip;

figure
plot(start_year:stop_year,cn51_diff_et_precip,'LineWidth',2);
hold on
plot(start_year:stop_year,cn45_diff_et_precip,'LineWidth',2);
hold on
%plot(start_year:stop_year,mean_global_ann_cn40_gpp,'LineWidth',2);
%hold on
plot(start_year:stop_year,fc_diff_et_precip,'LineWidth',2,'LineStyle','--');
hold off
grid on
xlabel('year')
ylabel('ET - P  [mm day^{-1}]')
title('annual average ET - P [mm day^{-1}]')
legend('CN51','CN45','FluxCom','Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 30 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/CN51_eval_ETmP'],'png')

% latent heat

figure
plot(start_year:stop_year,mean_global_ann_cn51_lh,'LineWidth',2);
hold on
plot(start_year:stop_year,mean_global_ann_cn51_et,'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_et,'LineWidth',2,'LineStyle','--');
hold off
grid on
xlabel('year')
ylabel('LH [W m^{-2}]')
title('annual average LH [W m^{-2}]')
legend('CN51','CN51 ET','FluxCom','Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 30 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/CN51_eval_Fluxcom_LH'],'png')

% evaporative fraction

mean_global_ann_cn51_ef = mean_global_ann_cn51_lh./mean_global_ann_cn51_rn;
mean_global_ann_fc_ef = mean_global_ann_fc_et./mean_global_ann_fc_rn;

figure
plot(start_year:stop_year,mean_global_ann_cn51_ef,'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_ef,'LineWidth',2,'LineStyle','--');
hold off
grid on
xlabel('year')
ylabel('EF [W m^{-2}]')
title('annual average EF [W m^{-2}]')
legend('CN51','FluxCom','Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 30 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/local/CN51_eval_Fluxcom_EF'],'png')
