clc
clear 
close all 

%% load in the sgp data files for different stations. 
% stations available
data_dir = './ebbr_ecor_data/215758/';
addpath('./ebbr_ecor_data/215758/')
stn_list = {'E32','E34','E35', 'E36','E21','E31', 'E33', 'E37','E38'};
% stn_list = {'E32','E34','E35', 'E36'};    
days_in_mon = [31, 28, 31, 30, 31, 30, 31, 31,30, 31, 30, 31];
        
year = '2012';
mon_list  = {'01', '02', '03','04', '05', '06', '07', '08', '09', '10', '11', '12'};

%     stn_name = ['sgpstamp', char(stn_list(1)), '.b1.' , year , mon, day, '.00000.nc'];
sh_stn = [];

for istn = 1:length(stn_list)
    sh_save = [];
    
    for year = 2011:2015
        if leapyear(year)
            days_in_mon = [31, 29, 31, 30, 31, 30, 31, 31,30, 31, 30, 31];
        else
            days_in_mon = [31, 28, 31, 30, 31, 30, 31, 31,30, 31, 30, 31];
        end
        
        for imon = 1:12
            for iday = 1:days_in_mon(imon)
                mon  = char(mon_list(imon));
                day  = sprintf('%02i',iday);

                stn_name1 = ['sgp30baebbr', char(stn_list(istn)), '.s1.' , num2str(year) , mon, day, '.000000.cdf'];
                stn_name2 = ['sgp30qcecor', char(stn_list(istn)), '.s1.' , num2str(year) , mon, day, '.000000.cdf'];

%                 stn_name1 = 'sgp30baebbrE32.s1.20120101.000000.cdf';
                if isfile([data_dir,stn_name1]) 
                    sh_tmp    = ncread(stn_name1, 'be_sensible_heat_flux');
                    qc_tmp    = ncread(stn_name1, 'qc_be_sensible_heat_flux');
                    sh_hour   = -1*sh_tmp(1:end,1);
                    qc_hour   = qc_tmp(1:end,1);             
                    sh_hour(qc_hour ~= 0) = nan;
                    
                    ratio = length(sh_tmp(isnan(sh_tmp)))/length(sh_tmp);

                    if ratio > 0.2    % if ratio of nan values is greater than 0.2 than set daily mean to zero.
                        sh_daily = nan;
                    % CONVERT HALF HOURLY TO DAILY MEAN. 
                    else
                        sh_daily = nanmean(sh_hour);
                    end
                    
                elseif isfile([data_dir,stn_name2]) 
                    sh_tmp    = ncread(stn_name2, 'sensible_heat_flux');
                    qc_tmp    = ncread(stn_name2, 'qc_sensible_heat_flux');
                    sh_hour   = sh_tmp(1:end,1);
                    qc_hour   = qc_tmp(1:end,1);             
                    sh_hour(qc_hour ~= 0) = nan;

                    ratio = length(sh_tmp(isnan(sh_tmp)))/length(sh_tmp);

                    if ratio > 0.2
                        sh_daily = nan;
                    % CONVERT HALF HOURLY TO DAILY MEAN. 
                    else                    
                    % CONVERT HALF HOURLY TO DAILY MEAN. 
                        sh_daily = nanmean(sh_hour);
                    end
                else 
                    sh_daily = nan;
                end

                sh_save = cat(2, sh_save, sh_daily);
            end 
        end % loop month
    
    end % loop year
    
    sh_stn = cat(3,sh_stn,sh_save);
end

sh_stn = squeeze(sh_stn);

%%
%% load in wrf
addpath('./wrf_lhsh/')
sh_wrf_2011 = ncread('noahmp_lhsh_2011.nc','SH');  % runoption, stations, days
sh_wrf_2012 = ncread('noahmp_lhsh_2012.nc','SH');
sh_wrf_2013 = ncread('noahmp_lhsh_2013.nc','SH');
sh_wrf_2014 = ncread('noahmp_lhsh_2014.nc','SH');
sh_wrf_2015 = ncread('noahmp_lhsh_2015.nc','SH');

sh_wrf = cat(3, sh_wrf_2011,sh_wrf_2012, sh_wrf_2013, sh_wrf_2014, sh_wrf_2015) * 0.1;

%%
    close all

    plot_stn_daily(sh_stn, sh_wrf);
 
    sh_wrf_mean = squeeze(mean(sh_wrf,2));

%% CALCULATE RMSE AND LINEAR CORRELATION COEFFICIENT

for ii = 1:13
    [rmse(ii), r(ii),bias(ii)] = calc_rmse_r(sh_stn, sh_wrf, ii);
end


%% Convert daily to monthly 
imon = 1;  day_tot = 0;
for year = 2011:2015
        if leapyear(year)
            days_in_mon = [31, 29, 31, 30, 31, 30, 31, 31,30, 31, 30, 31];
        else
            days_in_mon = [31, 28, 31, 30, 31, 30, 31, 31,30, 31, 30, 31];
        end

        for mon = 1:12
            if mon == 1
                inds = 1 + day_tot;
                inde = days_in_mon(mon)+ inds-1;              
            else 
               inds = sum(days_in_mon(1:mon-1))+1+ day_tot;
               inde = days_in_mon(mon) + inds-1;
            end
            
            sh_wrf_mon(:,:,imon) = mean(sh_wrf(:,:,inds:inde),3);
% find the ratio of invalid data, if more than 50% than put the monthly
% value as missing.
            stn_data = sh_stn(inds:inde,:);
            for istn = 1:size(stn_data,2)
                stn_data_istn = stn_data(:,istn);
                nan_num = length(stn_data_istn(isnan(stn_data_istn)));
                ratio = nan_num/(inde - inds + 1)
                if ratio > .5
                    sh_stn_mon(imon,istn) = nan;
                else 
                    sh_stn_mon(imon,istn) = nanmean(sh_stn(inds:inde,istn),1);
                end 
            end
        imon = imon + 1;
        end
        
        day_tot = inde;
        
end
    

%% plot monthly values
% ---- for each station1 ----
% ---- for the mean stations -----
for imon = 1:60
    stn_temp = sh_stn_mon(imon, :);
    ratio = length(stn_temp(isnan(stn_temp)))/length(stn_temp) ;
    if ratio > 0.5  % if nan value ratio greater than 50%
        sh_stn_mon_mean(imon) = nan;
    else
        sh_stn_mon_mean(imon) = nanmean(sh_stn_mon(imon,:),2);% months, stations 
    end


end 

sh_wrf_mon_mean = squeeze(nanmean(sh_wrf_mon,2) );   % ropt, station, months

xx = 1:60;
figure
plot(xx, sh_stn_mon_mean, 'k-o'); hold on;
plot(xx, sh_wrf_mon_mean(1,:));hold on;
plot(xx, sh_wrf_mon_mean(2,:));hold on;
plot(xx, sh_wrf_mon_mean(3,:));hold on;


%%
close all
figure(1) 
set(gcf,'PaperPositionMode', 'auto')
set(gcf,'Position', [100 100 1000 1000])
subplot(321)
plot(1:60, sh_stn_mon_mean ,'k-o'); hold on
plot(1:60, sh_wrf_mon_mean(1,:));hold on
plot(1:60, sh_wrf_mon_mean(2,:));hold on
legend('OBS', 'Default', 'gSSURGO', 'Location', 'Northwest')
xlim([13,60]); xlabel('Mon')
xticks([ 0:6:60]);grid on;
ylabel('SH (W m^{-2})')
title('Soil type')

subplot(322)
plot(1:60, sh_stn_mon_mean ,'k-o'); hold on
plot(1:60, sh_wrf_mon_mean(2,:));hold on
plot(1:60, sh_wrf_mon_mean(3,:));hold on
legend('OBS', 'old SM table', 'new SM table', 'Location', 'Northwest')
xlim([13,60]); xlabel('Mon')
xticks([ 0:6:60]);grid on;
ylabel('SH (W m^{-2})')
title('Soil table')


subplot(323)
plot(1:60, sh_stn_mon_mean ,'k-o'); hold on
plot(1:60, sh_wrf_mon_mean(3,:));hold on
plot(1:60, sh_wrf_mon_mean(4,:));hold on
plot(1:60, sh_wrf_mon_mean(5,:));hold on
legend('OBS', 'Default Vegetation', 'MODIS climatological', 'MODIS time varying', 'Location', 'Northwest')
xlim([13,60]); xlabel('Mon')
xticks([ 0:6:60]);grid on;
ylabel('SH (W m^{-2})')
title('Dynamic vegetation')


subplot(324)
plot(1:60, sh_stn_mon_mean ,'k-o'); hold on
plot(1:60, sh_wrf_mon_mean(5,:));hold on
plot(1:60, sh_wrf_mon_mean(6,:));hold on
legend('OBS','NLDAS-2 prcp','Stage-IV prcp', 'Location', 'Northwest')
xlim([13,60]); xlabel('Mon')
xticks([ 0:6:60]);grid on;
ylabel('SH (W m^{-2})')
title('Precipitation forcing')

subplot(325)
plot(1:60, sh_stn_mon_mean ,'k-o'); hold on
plot(1:60, sh_wrf_mon_mean(6,:));hold on
% plot(1:60, sh_wrf_mon_mean(7,:));hold on
% plot(1:60, sh_wrf_mon_mean(8,:));hold on
plot(1:60, sh_wrf_mon_mean(9,:));hold on
legend('OBS', 'no routing', 'routing', 'Location', 'Northwest')
xlim([13,60]); xlabel('Mon')
xticks([ 0:6:60]);grid on;
ylabel('SH (W m^{-2})')
title('Routing')

subplot(326)
plot(1:60, sh_stn_mon_mean ,'k-o'); hold on
plot(1:60, sh_wrf_mon_mean(10,:));hold on
plot(1:60, sh_wrf_mon_mean(11,:));hold on
plot(1:60, sh_wrf_mon_mean(12,:));hold on
plot(1:60, sh_wrf_mon_mean(13,:));hold on
legend('OBS', '100m', '250m', '500m','1000m','Location', 'Northwest')
xlim([13,60]); xlabel('Mon')
xticks([ 0:6:60]);grid on;
ylabel('SH (W m^{-2})')
title('Routing resolution')

print('-depsc', 'SH_monthly')

%%

for ii = 1:13
    [rmse_mon(ii), r_mon(ii),bias_mon(ii)] = calc_rmse_r(sh_stn_mon(13:60,:), sh_wrf_mon(:,:,13:60), ii);
end