clc
clear 
close all 

%% load in the sgp data files for different stations. 
% stations available
addpath('/Users/yang560/Dropbox/PNNL/Matlab_functions/utilities/')

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];
        
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'];
lh_stn = [];

for istn = 1:length(stn_list)
    lh_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]) 
                    lh_tmp    = ncread(stn_name1, 'be_latent_heat_flux');
                    qc_tmp    = ncread(stn_name1, 'qc_be_latent_heat_flux');
                    lh_hour   = -1*lh_tmp(1:end,1);
                    qc_hour   = qc_tmp(1:end,1);             
                    lh_hour(qc_hour ~= 0) = nan;
                    
                    ratio = length(lh_tmp(isnan(lh_tmp)))/length(lh_tmp);

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

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

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

                lh_save = cat(2, lh_save, lh_daily);
            end 
        end % loop month
    
    end % loop year
    
    lh_stn = cat(3,lh_stn,lh_save);
end

lh_stn = squeeze(lh_stn);

%%
%% load in wrf
addpath('./wrf_lhsh/')
lh_wrf_2011 = ncread('noahmp_lhsh_2011.nc','LH');  % runoption, stations, days
lh_wrf_2012 = ncread('noahmp_lhsh_2012.nc','LH');
lh_wrf_2013 = ncread('noahmp_lhsh_2013.nc','LH');
lh_wrf_2014 = ncread('noahmp_lhsh_2014.nc','LH');
lh_wrf_2015 = ncread('noahmp_lhsh_2015.nc','LH');

lh_wrf = cat(3, lh_wrf_2011,lh_wrf_2012, lh_wrf_2013, lh_wrf_2014, lh_wrf_2015) * 0.1;

%%
    close all

    plot_stn_daily(lh_stn, lh_wrf);
 
    lh_wrf_mean = squeeze(mean(lh_wrf,2));

%% CALCULATE RMSE AND LINEAR CORRELATION COEFFICIENT

for ii = 1:13
    [rmse(ii), r(ii),bias(ii)] = calc_rmse_r(lh_stn, lh_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
            
            lh_wrf_mon(:,:,imon) = mean(lh_wrf(:,:,inds:inde),3);
% find the ratio of invalid data, if more than 50% than put the monthly
% value as missing.
            stn_data = lh_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
                    lh_stn_mon(imon,istn) = nan;
                else 
                    lh_stn_mon(imon,istn) = nanmean(lh_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 = lh_stn_mon(imon, :);
    ratio = length(stn_temp(isnan(stn_temp)))/length(stn_temp) ;
    if ratio > 0.5  % if nan value ratio greater than 50%
        lh_stn_mon_mean(imon) = nan;
    else
        lh_stn_mon_mean(imon) = nanmean(lh_stn_mon(imon,:),2);% months, stations 
    end


end 

lh_wrf_mon_mean = squeeze(nanmean(lh_wrf_mon,2) );   % ropt, station, months

%%

addpath('/Users/yang560/Dropbox/PNNL/Matlab_functions/utilities/')
make_it_tight = true;
subplot = @(m,n,p) subtightplot (m, n, p, [0.055 0.055], [0.07 0.07], [0.07 0.07]);
if ~make_it_tight,  clear subplot;  end

close all
figure(1) 
set(gcf,'PaperPositionMode', 'auto')
set(gcf,'Position', [0 0 800 600])
subplot(3,2,1)
plot(1:60, lh_stn_mon_mean ,'k-o'); hold on
plot(1:60, lh_wrf_mon_mean(1,:),'LineWidth', 2);hold on
plot(1:60, lh_wrf_mon_mean(2,:),'LineWidth', 2);hold on
legend('OBS', 'Default (A1)', 'gSSURGO (A2)', 'Location', 'Northwest'); legend boxoff
xlim([1,60]); 
xticks([ 0:6:60]);grid on;
ylabel('LH (W m^{-2})')
title_left('a) Soil type')

subplot(3,2,2)
plot(1:60, lh_stn_mon_mean ,'k-o'); hold on
plot(1:60, lh_wrf_mon_mean(2,:),'LineWidth', 2);hold on
plot(1:60, lh_wrf_mon_mean(3,:),'LineWidth', 2);hold on
legend('OBS', 'old SM table (A2)', 'new SM table (A3)', 'Location', 'Northwest'); legend boxoff
xlim([1,60]); 
xticks([ 0:6:60]);grid on;
title_left('b) Soil table')


subplot(3,2,3)
plot(1:60, lh_stn_mon_mean ,'k-o'); hold on
plot(1:60, lh_wrf_mon_mean(3,:),'LineWidth', 2);hold on
plot(1:60, lh_wrf_mon_mean(4,:),'LineWidth', 2);hold on
plot(1:60, lh_wrf_mon_mean(5,:),'LineWidth', 2);hold on
legend('OBS', 'Default Vegetation (A3)', 'MODIS climatological (A4)', 'MODIS time varying (A5)', 'Location', 'Northwest'); legend boxoff
xlim([1,60]); 
xticks([ 0:6:60]);grid on;
ylabel('LH (W m^{-2})')
title_left('c) Dynamic vegetation')


subplot(3,2,4)
plot(1:60, lh_stn_mon_mean ,'k-o'); hold on
plot(1:60, lh_wrf_mon_mean(5,:),'LineWidth', 2);hold on
plot(1:60, lh_wrf_mon_mean(6,:),'LineWidth', 2);hold on
legend('OBS','NLDAS-2 prcp (A5)','Stage-IV prcp (A6)', 'Location', 'Northwest'); legend boxoff
xlim([1,60]); 
xticks([ 0:6:60]);grid on;
title_left('d) Precipitation forcing')

subplot(3,2,5)
plot(1:60, lh_stn_mon_mean ,'k-o'); hold on
plot(1:60, lh_wrf_mon_mean(6,:),'LineWidth', 2);hold on
% plot(1:60, lh_wrf_mon_mean(7,:));hold on
% plot(1:60, lh_wrf_mon_mean(8,:));hold on
plot(1:60, lh_wrf_mon_mean(9,:),'LineWidth', 2);hold on
legend('OBS', 'no routing (A6)', 'routing (A9)', 'Location', 'Northwest'); legend boxoff
xlim([1,60]); xlabel('Mon')
xticks([ 0:6:60]);grid on;
ylabel('LH (W m^{-2})')
title_left('e) Routing')

subplot(3,2,6)
plot(1:60, lh_stn_mon_mean ,'k-o'); hold on
plot(1:60, lh_wrf_mon_mean(10,:),'LineWidth', 2);hold on
plot(1:60, lh_wrf_mon_mean(11,:),'LineWidth', 2);hold on
plot(1:60, lh_wrf_mon_mean(12,:),'LineWidth', 2);hold on
plot(1:60, lh_wrf_mon_mean(13,:),'LineWidth', 2);hold on
legend('OBS', '100m (B2)', '250m (B3)', '500m (B4)','1000m (B5)','Location', 'Northwest'); legend boxoff
xlim([1,60]); xlabel('Mon')
xticks([ 0:6:60]);grid on;
title_left('f) Routing resolution')

print('-depsc', 'LH_monthly')


%%
lh_stn_mon_tmp = reshape(lh_stn_mon_mean, [12,5]); 
lh_wrf_mon_tmp = reshape(lh_wrf_mon_mean, [13, 12, 5]);

lh_stn_mon_min = min(lh_stn_mon_tmp, [], 2); 
lh_stn_mon_max = max(lh_stn_mon_tmp, [], 2); 
lh_stn_mon_avg = nanmean(lh_stn_mon_tmp, 2);
lh_stn_mon_std = nanstd(lh_stn_mon_tmp,0, 2);

lh_wrf_mon_min = min(lh_wrf_mon_tmp, [], 3);
lh_wrf_mon_max = max(lh_wrf_mon_tmp, [], 3);
lh_wrf_mon_avg = nanmean(lh_wrf_mon_tmp, 3);
lh_wrf_mon_std = nanstd(lh_wrf_mon_tmp,0, 3);


xx = [1:12];


figure(1) 
set(gcf,'PaperPositionMode', 'auto')
set(gcf,'Position', [0 0 800 600])
subplot(3,2,1)
% plot(1:12, lh_stn_mon_avg ,'k-o'); hold on
% fill([xx fliplr(xx)],[lh_stn_mon_max' fliplr(lh_stn_mon_min')] ,'k'); hold on;alpha(0.15)
% plot(1:12, lh_wrf_mon_avg(1,:));hold on
% fill([xx fliplr(xx)],[lh_wrf_mon_max(1,:) fliplr(lh_wrf_mon_min(1,:))] ,'b'); hold on;alpha(0.15)
% plot(1:12, lh_wrf_mon_avg(2,:));hold on
% fill([xx fliplr(xx)],[lh_wrf_mon_max(2,:) fliplr(lh_wrf_mon_min(2,:))] ,'r'); hold on;alpha(0.15)

e=errorbar(1:12, lh_stn_mon_avg,lh_stn_mon_std, '-o'); hold on
e.Color = 'black'; clear e;
e=errorbar(1:12, lh_wrf_mon_avg(1,:),lh_wrf_mon_std(1,:)); hold on
e.Color = 'blue';clear e;
e=errorbar(1:12, lh_wrf_mon_avg(2,:),lh_wrf_mon_std(2,:)); hold on
e.Color = 'red';clear e;

legend('OBS', 'Default', 'gSSURGO', 'Location', 'Northeast')
xlim([1,12]); xlabel('Mon')
xticks([1:1:12]);grid on;
ylabel('LH (W m^{-2})')
title_left('a) Soil type')

subplot(3,2,2)
e=errorbar(1:12, lh_stn_mon_avg,lh_stn_mon_std, '-o'); hold on
e.Color = 'black';
e=errorbar(1:12, lh_wrf_mon_avg(2,:),lh_wrf_mon_std(2,:)); hold on
e.Color = 'blue';
e=errorbar(1:12, lh_wrf_mon_avg(3,:),lh_wrf_mon_std(3,:)); hold on
e.Color = 'red';

legend('OBS', 'old SM table', 'new SM table', 'Location', 'Northeast')
xlim([1,12]); xlabel('Mon')
xticks([1:1:12]);grid on;
ylabel('LH (W m^{-2})')
title_left('b) Soil table')


subplot(3,2,3)
e=errorbar(1:12, lh_stn_mon_avg,lh_stn_mon_std, '-o'); hold on
e.Color = 'black';
e=errorbar(1:12, lh_wrf_mon_avg(3,:),lh_wrf_mon_std(3,:)); hold on
e.Color = 'blue';
e=errorbar(1:12, lh_wrf_mon_avg(4,:),lh_wrf_mon_std(4,:)); hold on
e.Color = 'red';
e=errorbar(1:12, lh_wrf_mon_avg(5,:),lh_wrf_mon_std(5,:)); hold on
e.Color = 'cyan';

legend('OBS', 'Default Vegetation', 'MODIS climatological', 'MODIS time varying', 'Location', 'Northeast')
xlim([1,12]); xlabel('Mon')
xticks([1:1:12]);grid on;
ylabel('LH (W m^{-2})')
title_left('c) Dynamic vegetation')


subplot(3,2,4)
e=errorbar(1:12, lh_stn_mon_avg,lh_stn_mon_std, '-o'); hold on
e.Color = 'black';
e=errorbar(1:12, lh_wrf_mon_avg(5,:),lh_wrf_mon_std(5,:)); hold on
e.Color = 'blue';
e=errorbar(1:12, lh_wrf_mon_avg(6,:),lh_wrf_mon_std(6,:)); hold on
e.Color = 'red';

legend('OBS','NLDAS-2 prcp','Stage-IV prcp', 'Location', 'Northeast')
xlim([1,12]); xlabel('Mon')
xticks([1:1:12]);grid on;
ylabel('LH (W m^{-2})')
title_left('d) Precipitation forcing')

subplot(3,2,5)
e=errorbar(1:12, lh_stn_mon_avg,lh_stn_mon_std, '-o'); hold on
e.Color = 'black';
e=errorbar(1:12, lh_wrf_mon_avg(6,:),lh_wrf_mon_std(6,:)); hold on
e.Color = 'blue';
e=errorbar(1:12, lh_wrf_mon_avg(9,:),lh_wrf_mon_std(9,:)); hold on
e.Color = 'red';
% plot(1:12, lh_wrf_mon_mean(7,:));hold on
% plot(1:12, lh_wrf_mon_mean(8,:));hold on
% plot(1:12, lh_wrf_mon_mean(9,:));hold on
legend('OBS', 'no routing', 'routing', 'Location', 'Northeast')
xlim([1,12]); xlabel('Mon')
xticks([1:1:12]);grid on;
ylabel('LH (W m^{-2})')
title_left('e) Routing')

subplot(3,2,6)
e=errorbar(1:12, lh_stn_mon_avg,lh_stn_mon_std, '-o'); hold on
e.Color = 'black';
e=errorbar(1:12, lh_wrf_mon_avg(10,:),lh_wrf_mon_std(10,:)); hold on
e.Color = 'blue';
e=errorbar(1:12, lh_wrf_mon_avg(11,:),lh_wrf_mon_std(11,:)); hold on
e.Color = 'red';
e=errorbar(1:12, lh_wrf_mon_avg(12,:),lh_wrf_mon_std(12,:)); hold on
e.Color = 'cyan';
e=errorbar(1:12, lh_wrf_mon_avg(13,:),lh_wrf_mon_std(13,:)); hold on
e.Color = 'green';

% plot(1:12, lh_stn_mon_mean ,'k-o'); hold on
% plot(1:12, lh_wrf_mon_mean(10,:));hold on
% plot(1:12, lh_wrf_mon_mean(11,:));hold on
% plot(1:12, lh_wrf_mon_mean(12,:));hold on
% plot(1:12, lh_wrf_mon_mean(13,:));hold on
legend('OBS', '100m', '250m', '500m','1000m','Location', 'Northeast')
xlim([1,12]); xlabel('Mon')
xticks([1:1:12]);grid on;
ylabel('LH (W m^{-2})')
title_left('f) Routing resolution')

print('-depsc', 'LH_monthly_errorbar')


%%
xx = [1:12];
close all

make_it_tight = true;
subplot = @(m,n,p) subtightplot (m, n, p, [0.065 0.055], [0.1 0.045], [0.08 0.05]);
if ~make_it_tight,  clear subplot;  end

figure(1) 
set(gcf,'PaperPositionMode', 'auto')
set(gcf,'Position', [100 100 1000 1000])

subplot(3,2,1)
% plot(1:12, lh_stn_mon_avg ,'k-o'); hold on
% fill([xx fliplr(xx)],[lh_stn_mon_max' fliplr(lh_stn_mon_min')] ,'k'); hold on;alpha(0.15)
% plot(1:12, lh_wrf_mon_avg(1,:));hold on
% fill([xx fliplr(xx)],[lh_wrf_mon_max(1,:) fliplr(lh_wrf_mon_min(1,:))] ,'b'); hold on;alpha(0.15)
% plot(1:12, lh_wrf_mon_avg(2,:));hold on
% fill([xx fliplr(xx)],[lh_wrf_mon_max(2,:) fliplr(lh_wrf_mon_min(2,:))] ,'r'); hold on;alpha(0.15)
boxplot(lh_stn_mon_tmp')

e=errorbar(1:12, lh_stn_mon_avg,lh_stn_mon_std, '-o'); hold on
e.Color = 'black'; clear e;
e=errorbar(1:12, lh_wrf_mon_avg(1,:),lh_wrf_mon_std(1,:)); hold on
e.Color = 'blue';clear e;
e=errorbar(1:12, lh_wrf_mon_avg(2,:),lh_wrf_mon_std(2,:)); hold on
e.Color = 'red';clear e;

legend('OBS', 'Default', 'gSSURGO', 'Location', 'Northeast')
xlim([1,12]); xlabel('Mon')
xticks([1:1:12]);grid on;
ylabel('LH (W m^{-2})')
title('Soil type')

subplot(3,2,2)
e=errorbar(1:12, lh_stn_mon_avg,lh_stn_mon_std, '-o'); hold on
e.Color = 'black';
e=errorbar(1:12, lh_wrf_mon_avg(2,:),lh_wrf_mon_std(2,:)); hold on
e.Color = 'blue';
e=errorbar(1:12, lh_wrf_mon_avg(3,:),lh_wrf_mon_std(3,:)); hold on
e.Color = 'red';

legend('OBS', 'old SM table', 'new SM table', 'Location', 'Northeast')
xlim([1,12]); xlabel('Mon')
xticks([1:1:12]);grid on;
ylabel('LH (W m^{-2})')
title('Soil table')


subplot(3,2,3)
e=errorbar(1:12, lh_stn_mon_avg,lh_stn_mon_std, '-o'); hold on
e.Color = 'black';
e=errorbar(1:12, lh_wrf_mon_avg(3,:),lh_wrf_mon_std(3,:)); hold on
e.Color = 'blue';
e=errorbar(1:12, lh_wrf_mon_avg(4,:),lh_wrf_mon_std(4,:)); hold on
e.Color = 'red';
e=errorbar(1:12, lh_wrf_mon_avg(5,:),lh_wrf_mon_std(5,:)); hold on
e.Color = 'cyan';

legend('OBS', 'Default Vegetation', 'MODIS climatological', 'MODIS time varying', 'Location', 'Northeast')
xlim([1,12]); xlabel('Mon')
xticks([1:1:12]);grid on;
ylabel('LH (W m^{-2})')
title('Dynamic vegetation')


subplot(3,2,4)
e=errorbar(1:12, lh_stn_mon_avg,lh_stn_mon_std, '-o'); hold on
e.Color = 'black';
e=errorbar(1:12, lh_wrf_mon_avg(5,:),lh_wrf_mon_std(5,:)); hold on
e.Color = 'blue';
e=errorbar(1:12, lh_wrf_mon_avg(6,:),lh_wrf_mon_std(6,:)); hold on
e.Color = 'red';

legend('OBS','NLDAS-2 prcp','Stage-IV prcp', 'Location', 'Northeast')
xlim([1,12]); xlabel('Mon')
xticks([1:1:12]);grid on;
ylabel('LH (W m^{-2})')
title('Precipitation forcing')

subplot(3,2,5)
e=errorbar(1:12, lh_stn_mon_avg,lh_stn_mon_std, '-o'); hold on
e.Color = 'black';
e=errorbar(1:12, lh_wrf_mon_avg(6,:),lh_wrf_mon_std(6,:)); hold on
e.Color = 'blue';
e=errorbar(1:12, lh_wrf_mon_avg(9,:),lh_wrf_mon_std(9,:)); hold on
e.Color = 'red';
% plot(1:12, lh_wrf_mon_mean(7,:));hold on
% plot(1:12, lh_wrf_mon_mean(8,:));hold on
% plot(1:12, lh_wrf_mon_mean(9,:));hold on
legend('OBS', 'no routing', 'routing', 'Location', 'Northeast')
xlim([1,12]); xlabel('Mon')
xticks([1:1:12]);grid on;
ylabel('LH (W m^{-2})')
title('Routing')

subplot(3,2,6)
e=errorbar(1:12, lh_stn_mon_avg,lh_stn_mon_std, '-o'); hold on
e.Color = 'black';
e=errorbar(1:12, lh_wrf_mon_avg(10,:),lh_wrf_mon_std(10,:)); hold on
e.Color = 'blue';
e=errorbar(1:12, lh_wrf_mon_avg(11,:),lh_wrf_mon_std(11,:)); hold on
e.Color = 'red';
e=errorbar(1:12, lh_wrf_mon_avg(12,:),lh_wrf_mon_std(12,:)); hold on
e.Color = 'cyan';
e=errorbar(1:12, lh_wrf_mon_avg(13,:),lh_wrf_mon_std(13,:)); hold on
e.Color = 'green';

% plot(1:12, lh_stn_mon_mean ,'k-o'); hold on
% plot(1:12, lh_wrf_mon_mean(10,:));hold on
% plot(1:12, lh_wrf_mon_mean(11,:));hold on
% plot(1:12, lh_wrf_mon_mean(12,:));hold on
% plot(1:12, lh_wrf_mon_mean(13,:));hold on
legend('OBS', '100m', '250m', '500m','1000m','Location', 'Northeast')
xlim([1,12]); xlabel('Mon')
xticks([1:1:12]);grid on;
ylabel('LH (W m^{-2})')
title('Routing resolution')

print('-depsc', 'LH_monthly_errorbar')

%%

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