clc 
clear all
close all

exp_path = '/discover/nobackup/jkolassa/SMAP_TC_surf_test/';
exp_name_louis = 'Idai_Louis';
exp_name_helf = 'Idai_Helfand';
                                                                                        
month_length = [31 28 14 30 31 30 31 31 30 31 30 31];                                                         

louis_et_vec = [];
helf_et_vec = [];
louis_wind_vec = [];
helf_wind_vec = [];
louis_sfmc_vec = [];
helf_sfmc_vec = [];
louis_lai_vec = [];
helf_lai_vec = [];
louis_tsurf_vec = [];
helf_tsurf_vec = [];

t_vec = 30:100:2330;

for m = 1:3
    m
    m_str = num2str(m,'%02i');
    for d = 1:month_length(m)

        d_str = num2str(d,'%02i');
        d
        for t = 1:length(t_vec)

        t_str = num2str(t_vec(t),'%04i');
        louis_lnd_file = [exp_path exp_name_louis '/output/SMAP_EASEv2_M09_GLOBAL/cat/ens0000/Y2019/M' ...
                          m_str '/' exp_name_louis '.tavg24_1d_lnd_Nt.2019' m_str d_str '_' t_str 'z.nc4'];

        helf_lnd_file = [exp_path exp_name_helf '/output/SMAP_EASEv2_M09_GLOBAL/cat/ens0000/Y2019/M' ...
                          m_str '/' exp_name_helf '.tavg24_1d_lnd_Nt.2019' m_str d_str '_' t_str 'z.nc4'];

        louis_lfs_file = [exp_path exp_name_louis '/output/SMAP_EASEv2_M09_GLOBAL/cat/ens0000/Y2019/M' ...
                          m_str '/' exp_name_louis '.tavg24_1d_lfs_Nt.2019' m_str d_str '_' t_str 'z.nc4'];

         helf_lfs_file = [exp_path exp_name_helf '/output/SMAP_EASEv2_M09_GLOBAL/cat/ens0000/Y2019/M' ...
                          m_str '/' exp_name_helf '.tavg24_1d_lfs_Nt.2019' m_str d_str '_' t_str 'z.nc4'];

         lat_data = ncread(louis_lnd_file,'lat');
         lon_data = ncread(louis_lnd_file,'lon');

         louis_wind_data = ncread(louis_lfs_file,'Wind');
         louis_wind_vec = [louis_wind_vec; louis_wind_data(:)];
         helf_wind_data = ncread(helf_lfs_file,'Wind');
         helf_wind_vec = [helf_wind_vec; helf_wind_data(:)];

         louis_et_data = ncread(louis_lnd_file,'EVLAND');
         louis_et_data = louis_et_data*86400;
         louis_et_vec = [louis_et_vec; louis_et_data(:)];

         helf_et_data = ncread(helf_lnd_file,'EVLAND');
         helf_et_data = helf_et_data*86400;
         helf_et_vec = [helf_et_vec; helf_et_data(:)];

         louis_sfmc_data = ncread(louis_lnd_file,'SFMC');
         louis_sfmc_vec = [louis_sfmc_vec; louis_sfmc_data(:)];
         helf_sfmc_data = ncread(helf_lnd_file,'SFMC');
         helf_sfmc_vec = [helf_sfmc_vec; helf_sfmc_data(:)];

         louis_lai_data = ncread(louis_lnd_file,'CNTLAI');
         louis_lai_vec = [louis_lai_vec; louis_lai_data(:)];
         helf_lai_data = ncread(helf_lnd_file,'CNTLAI');
         helf_lai_vec = [helf_lai_vec; helf_lai_data(:)];

         louis_tsurf_data = ncread(louis_lnd_file,'TSURF');                                         
         louis_tsurf_vec = [louis_tsurf_vec; louis_tsurf_data(:)];                                       
         helf_tsurf_data = ncread(helf_lnd_file,'TSURF');                                           
         helf_tsurf_vec = [helf_tsurf_vec; helf_tsurf_data(:)];  

       end
    end
end

lai_bin_low = [0 2.5 5 7.5];
lai_bin_high = [2.5 5 7.5 10];

sfmc_bin_low = 0:0.2:0.6;
sfmc_bin_high = 0.2:0.2:0.8;

wind_bin_low = 0:5:(max(louis_wind_vec(:))-5);
wind_bin_high = 5:5:max(louis_wind_vec(:));

tsurf_bin_low = 290:10:310;
tsurf_bin_high = 300:10:320;

CM = jet(length(sfmc_bin_low));
surf_scheme = 'Louis';

for l = 1:length(lai_bin_low)

    louis_lai_index = (louis_lai_vec>lai_bin_low(l)) + (louis_lai_vec<=lai_bin_high(l));
    louis_lai_index = louis_lai_index == 2;

    helf_lai_index = (helf_lai_vec>lai_bin_low(l)) + (helf_lai_vec<=lai_bin_high(l));
    helf_lai_index = helf_lai_index == 2;
    
    for t = 1:length(tsurf_bin_low)

        louis_tsurf_index = (louis_tsurf_vec>tsurf_bin_low(t)) + (louis_tsurf_vec<=tsurf_bin_high(t));           
        louis_tsurf_index = louis_tsurf_index == 2;                                                        

        helf_tsurf_index = (helf_tsurf_vec>tsurf_bin_low(t)) + (helf_tsurf_vec<=tsurf_bin_high(t));              
        helf_tsurf_index = helf_tsurf_index == 2;  
    figure

    for s = 1:length(sfmc_bin_low)
        
        louis_sfmc_index = (louis_sfmc_vec>sfmc_bin_low(s)) + (louis_sfmc_vec<=sfmc_bin_high(s));
        louis_sfmc_index = louis_sfmc_index == 2;
    
        helf_sfmc_index = (helf_sfmc_vec>sfmc_bin_low(s)) + (helf_sfmc_vec<=sfmc_bin_high(s));
        helf_sfmc_index = helf_sfmc_index == 2;
        
        louis_data_index = louis_lai_index + louis_sfmc_index + louis_tsurf_index;
        louis_data_index = louis_data_index == 3;

        helf_data_index = helf_lai_index + helf_sfmc_index + helf_tsurf_index;
        helf_data_index = helf_data_index == 3;

        if strcmp(surf_scheme,'Louis')
        for i = 1:length(wind_bin_low)
            wind_index = (louis_wind_vec>wind_bin_low(i)) + (louis_wind_vec<=wind_bin_high(i));
            wind_index = wind_index == 2;
            new_data_index = louis_data_index + wind_index;
            new_data_index = new_data_index == 2;
            et_plot_vec(i) = nanmean(louis_et_vec(new_data_index));
            error_vec(i) = nanstd(louis_et_vec(new_data_index),0);
            wind_plot_vec(i) = nanmean([wind_bin_low(i) wind_bin_high(i)]);   
        end
        errorbar(wind_plot_vec,et_plot_vec,error_vec,'Color',CM(s,:),'LineWidth',2)
        elseif strcmp(surf_scheme,'Helfand')
        for i = 1:length(wind_bin_low)
            wind_index = (helf_wind_vec>wind_bin_low(i)) + (helf_wind_vec<=wind_bin_high(i));    
            wind_index = wind_index == 2;                                                          
            new_data_index = helf_data_index + wind_index;                                        
            new_data_index = new_data_index == 2;
            et_plot_vec(i) = nanmean(helf_et_vec(new_data_index));
            error_vec(i) = nanstd(helf_et_vec(new_data_index),0);                                 
            wind_plot_vec(i) = nanmean([wind_bin_low(i) wind_bin_high(i)]);                        
        end                                                                                        
        errorbar(wind_plot_vec,et_plot_vec,error_vec,'Color',CM(s,:),'LineWidth',2) 
        end
        hold on
    end
    hold off
    grid on
    axis([0 max(louis_wind_vec) 0 15])
    xlabel('Wind speed at reference height [m s^-1]')
    ylabel('Evapotranspiration [mm day^-1]')
    title(['ET vs Wind speed for LAI range [' num2str(lai_bin_low(l)) ' ' num2str(lai_bin_high(l)) '] and Tsurf range [' num2str(tsurf_bin_low(t)) ' ' num2str(tsurf_bin_high(t)) ']'])
    legend('SM: 0 - 0.2','SM: 0.2 - 0.4','SM: 0.4 - 0.6','SM: 0.6 - 0.8')
    print('-dpng','-r300',['ET_vs_Wind_LAI_' num2str(lai_bin_low(l)) '_' num2str(lai_bin_high(l)) '_' ...
    'Tsurf_' num2str(tsurf_bin_low(t)) '_' num2str(tsurf_bin_high(t)) '_' surf_scheme '.png']); 
end
end
