clc
clear 
close all 

%% load in the sgp data files for different stations. 
% stations available
addpath('./sgpstamp/')
stn_list = {'E11','E21','E12', 'E13','E15','E31', 'E32', 'E33'...
            'E34','E35','E36','E37','E38','E39','E40', 'E9'};
        
days_in_mon = [30, 31, 31];
        
year = '2016';
mon_list  = {'06', '07', '08'};

%     stn_name = ['sgpstamp', char(stn_list(1)), '.b1.' , year , mon, day, '.00000.nc'];
soil_stn = [];
for istn = 1:16
    soil_save = [];
    for imon = 1:3
        for iday = 1:days_in_mon(imon)
            mon  = char(mon_list(imon));
            day  = sprintf('%02i',iday);
            stn_name = ['sgpstamp', char(stn_list(istn)), '.b1.' , year , mon, day, '.000000.nc'];

            soilm_tmp = ncread(stn_name, 'soil_specific_water_content_west');
            qc_tmp    = ncread(stn_name, 'qc_soil_specific_water_content_west');

            soil_hour = soilm_tmp(1,1:end)*0.01;
            qc_hour   = qc_tmp(1,1:end);              % CALCULATE THE FIRST LEVEL

            soil_hour(qc_hour ~= 0) = nan;
                
            % CONVERT HALF HOURLY TO DAILY MEAN. 
            soil_daily = nanmean(soil_hour);
            
            soil_save = cat(2, soil_save, soil_daily);
        end 
    end
    
    soil_stn = cat(3,soil_stn,soil_save);
end

%% load in wrf
soil_wrf = ncread('noahmp_soilm_2016.nc','SOIL_M');  % lev, day, stn, runoption

%%
% close all
% for ii = 1:16
%     plot_stn_daily(soil_stn, soil_wrf, ii)
% end 
%%
plot_stn_daily_composite(soil_stn, soil_wrf)
print('-depsc', 'soil_moisture_comparison')
%% CALCULATE RMSE AND LINEAR CORRELATION COEFFICIENT
for ii = 1:16
    [rmse(ii,:), r(ii,:),bias(ii,:)] = calc_rmse_r(soil_stn, soil_wrf, ii);
end 

[rmse, r,bias] = calc_rmse_r_composite(soil_stn, soil_wrf);