clc
clear 
close all 

%% load in the sgp data files for different stations. 
% stations available
addpath('./sgpswats/')
stn_list = {'E11','E21','E12', 'E13','E15','E31', 'E32', 'E33'...
            'E34','E35','E36','E37','E38','E39','E40', 'E9'};

year = 2011;

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
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'];
soil_stn = [];
for istn = 1:16
    soil_save = [];
    nvalid = 0;
    for imon = 1:12
        for iday = 1:days_in_mon(imon)
            mon  = char(mon_list(imon));
            day  = sprintf('%02i',iday);
            stn_name = ['./sgpswats/','sgpswats', char(stn_list(istn)), '.b1.' , num2str(year) , mon, day, '.000700.cdf'];
             if isfile(stn_name)
                 
                 try 
                    soilm_tmp = ncread(stn_name, 'watcont_W');
                    qc_tmp    = ncread(stn_name, 'qc_watcont_W');
                 catch 
                 	soilm_tmp = ncread(stn_name, 'watcont_w');
                    qc_tmp    = ncread(stn_name, 'qc_watcont_w');  
                 end
                 
                soil_hour = soilm_tmp(1,1:end);
                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);
                nvalid = nvalid + 1;
             else
                stn_name = ['./sgpswats/','sgpswats', char(stn_list(istn)), '.b1.' , num2str(year) , mon, day, '.000700.nc'];
                if isfile(stn_name)
 
                     try 
                        soilm_tmp = ncread(stn_name, 'watcont_W');
                        qc_tmp    = ncread(stn_name, 'qc_watcont_W');
                     catch 
                        soilm_tmp = ncread(stn_name, 'watcont_w');
                        qc_tmp    = ncread(stn_name, 'qc_watcont_w');  
                     end

                    soil_hour = soilm_tmp(1,1:end);
                    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);
                    nvalid = nvalid + 1;
                else
    %                 disp([stn_name, ' not exist'])
                    soil_save = cat(2, soil_save, nan);
                end
                
             end
             
             
             
             
        end 
    end
    
    if (nvalid/365 > 0.8)
        disp([char(stn_list(istn)),' is valid'])
    end
    
    soil_stn = cat(3,soil_stn,soil_save);
end

%% load in wrf
wrf_title = ['noahmp_soilm_',num2str(year),'.nc'];
soil_wrf = ncread(wrf_title,'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)
%% 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);