clear all
close all

data_QLN = readtable('jan2011_QLN.csv');
data_QHW = readtable('jan2011_QHW.csv');
data_VEL = readtable('jan2011_KVEL.csv');
data_SLC = readtable('jan2011_KSLC.csv');
data_LGU = readtable('jan2011_KLGU.csv');

date_QLN = table2array(data_QLN(2:end,'Date_Time'));
date_QHW = table2array(data_QHW(2:end,'Date_Time'));
date_VEL = table2array(data_VEL(2:end,'Date_Time'));
date_SLC = table2array(data_SLC(2:end,'Date_Time'));
date_LGU = table2array(data_LGU(2:end,'Date_Time'));

ts_datenum_QLN = datenum(date_QLN, 'yyyy-mm-ddThh:MM:SSZ');
ts_datenum_QHW = datenum(date_QHW, 'yyyy-mm-ddThh:MM:SSZ');
ts_datenum_VEL = datenum(date_VEL, 'yyyy-mm-ddThh:MM:SSZ');
ts_datenum_SLC = datenum(date_SLC, 'yyyy-mm-ddThh:MM:SSZ');
ts_datenum_LGU = datenum(date_LGU, 'yyyy-mm-ddThh:MM:SSZ');

temp_QLN = table2array(data_QLN(2:end,'air_temp_set_1'));
temp_QHW = table2array(data_QHW(2:end,'air_temp_set_1'));
temp_VEL = table2array(data_VEL(2:end,'air_temp_set_1'));
temp_SLC = table2array(data_SLC(2:end,'air_temp_set_1'));
temp_LGU = table2array(data_LGU(2:end,'air_temp_set_1'));

wspd_QLN = table2array(data_QLN(2:end,'wind_speed_set_1'));
wspd_QHW = table2array(data_QHW(2:end,'wind_speed_set_1'));
wspd_VEL = table2array(data_VEL(2:end,'wind_speed_set_1'));
wspd_SLC = table2array(data_SLC(2:end,'wind_speed_set_1'));
wspd_LGU = table2array(data_LGU(2:end,'wind_speed_set_1'));

wdir_QLN = table2array(data_QLN(2:end,'wind_direction_set_1'));
wdir_QHW = table2array(data_QHW(2:end,'wind_direction_set_1'));
wdir_VEL = table2array(data_VEL(2:end,'wind_direction_set_1'));
wdir_SLC = table2array(data_SLC(2:end,'wind_direction_set_1'));
wdir_LGU = table2array(data_LGU(2:end,'wind_direction_set_1'));

index_QLN = find(strcmp(temp_QLN,''));
temp_QLN(index_QLN) = {'NaN'};
index_QHW = find(strcmp(temp_QHW,''));
temp_QHW(index_QHW) = {'NaN'};

ts_temp_QLN = cellfun(@str2num,temp_QLN);
ts_temp_QHW = cellfun(@str2num,temp_QHW);
ts_temp_VEL = cellfun(@str2num,temp_VEL);
ts_temp_SLC = cellfun(@str2num,temp_SLC);
ts_temp_LGU = cellfun(@str2num,temp_LGU);

index_QLN = find(strcmp(wspd_QLN,''));
wspd_QLN(index_QLN) = {'NaN'};
index_QHW = find(strcmp(wspd_QHW,''));
wspd_QHW(index_QHW) = {'NaN'};
index_LGU = find(strcmp(wspd_LGU,''));
wspd_LGU(index_LGU) = {'NaN'};
index_VEL = find(strcmp(wspd_VEL,''));
wspd_VEL(index_VEL) = {'NaN'};
wspd_LGU(index_LGU) = {'NaN'};
index_SLC = find(strcmp(wspd_SLC,''));
wspd_SLC(index_SLC) = {'NaN'};

ts_wspd_QLN = cellfun(@str2num,wspd_QLN);
ts_wspd_QHW = cellfun(@str2num,wspd_QHW);
ts_wspd_VEL = cellfun(@str2num,wspd_VEL);
ts_wspd_SLC = cellfun(@str2num,wspd_SLC);
ts_wspd_LGU = cellfun(@str2num,wspd_LGU);

index_QLN = find(strcmp(wdir_QLN,''));
wdir_QLN(index_QLN) = {'NaN'};
index_QHW = find(strcmp(wdir_QHW,''));
wdir_QHW(index_QHW) = {'NaN'};
index_VEL = find(strcmp(wdir_VEL,''));
wdir_VEL(index_VEL) = {'NaN'};
index_SLC = find(strcmp(wdir_SLC,''));
wdir_SLC(index_SLC) = {'NaN'};
index_LGU = find(strcmp(wdir_LGU,''));
wdir_LGU(index_LGU) = {'NaN'};

ts_wdir_QLN = cellfun(@str2num,wdir_QLN);
ts_wdir_QHW = cellfun(@str2num,wdir_QHW);
ts_wdir_VEL = cellfun(@str2num,wdir_VEL);
ts_wdir_SLC = cellfun(@str2num,wdir_SLC);
ts_wdir_LGU = cellfun(@str2num,wdir_LGU);

%%

ts_datevec_QLN = datevec(ts_datenum_QLN);
ts_datevec_QHW = datevec(ts_datenum_QHW);
ts_datevec_VEL = datevec(ts_datenum_VEL);
ts_datevec_SLC = datevec(ts_datenum_SLC);
ts_datevec_LGU = datevec(ts_datenum_LGU);

hour_QLN = find(ts_datevec_QLN(:,5)==0);
hour_QHW = find(ts_datevec_QHW(:,5)==45);
hour_VEL = find(ts_datevec_VEL(:,5)==53);
hour_SLC = find(ts_datevec_SLC(:,5)==0);
hour_LGU = find(ts_datevec_LGU(:,5)==51);

data_hour_QLN = [ts_datenum_QLN(hour_QLN) ts_datevec_QLN(hour_QLN,3:5) ts_temp_QLN(hour_QLN) ts_wspd_QLN(hour_QLN) ts_wdir_QLN(hour_QLN)];
data_hour_QHW = [ts_datenum_QHW(hour_QHW) ts_datevec_QHW(hour_QHW,3:5) ts_temp_QHW(hour_QHW) ts_wspd_QHW(hour_QHW) ts_wdir_QHW(hour_QHW)];
data_hour_VEL = [ts_datenum_VEL(hour_VEL) ts_datevec_VEL(hour_VEL,3:5) ts_temp_VEL(hour_VEL) ts_wspd_VEL(hour_VEL) ts_wdir_VEL(hour_VEL)];
data_hour_SLC = [ts_datenum_SLC(hour_SLC) ts_datevec_SLC(hour_SLC,3:5) ts_temp_SLC(hour_SLC) ts_wspd_SLC(hour_SLC) ts_wdir_SLC(hour_SLC)];
data_hour_LGU = [ts_datenum_LGU(hour_LGU) ts_datevec_LGU(hour_LGU,3:5) ts_temp_LGU(hour_LGU) ts_wspd_LGU(hour_LGU) ts_wdir_LGU(hour_LGU)];


%%

lati = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group3/horel_data/cfoster/wrf_daq_runs/jan2011_alldoms/wrfout_d03_2011-01-01_00:00:00'),'XLAT');
long = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group3/horel_data/cfoster/wrf_daq_runs/jan2011_alldoms/wrfout_d03_2011-01-01_00:00:00'),'XLONG');

lati = double(lati(:,:,1));
long = double(long(:,:,1));

lats = lati(:,:,1)';
lons = long(:,:,1)';

id  = {'QLN',     'QHW',     'VEL',     'SLC',     'LGU'};
lon = [-111.7133, -111.8717, -109.5127, -111.9650, -111.8519];
lat = [40.3388,   40.7335,   40.4430,    40.7707,   40.7865];

for y = 1:length(lons(:,1))-1
   for x = 1:length(lons(1,:))-1
      for i = 1:length(lon)
         if (lons(y,x) < lon(i) & lon(i) < lons(y,x+1)) & (lats(y,x) < lat(i) & lat(i) < lats(y+1,x))
            x_grid(i) = x;
            y_grid(i) = y;
         end
      end
   end
end

%%

run  = 'jan2011_alldoms';
days = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11];
year = [2011,2011,2011,2011,2011,2011,2011,2011,2011,2011,2011];
mont = [1,1,1,1,1,1,1,1,1,1,1,1];

% run  = 'dec2013_alldoms';
% days = [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20];
% year = [2013,2013,2013,2013,2013,2013,2013,2013,2013,2013,2013,2013,2013,2013];
% mont = [12,12,12,12,12,12,12,12,12,12,12,12,12,12];

% run  = 'feb2016_alldoms';
% days = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17];
% year = [2016,2016,2016,2016,2016,2016,2016,2016,2016,2016,2016,2016,2016,2016,2016,2016,2016];
% mont = [2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2];

temp = zeros(size(lati));
u = zeros(size(lati));
v = zeros(size(lati));

time_model = datenum(year(1),mont(1),days(1),0,0,0);

for d = 1:length(days)
 
% Deal with time to read in the correct file.

   if days(d) < 10
      day = strcat('0',int2str(days(d)));
   elseif days(d) >= 10
      day = int2str(days(d));
   end
   
   if mont(d) < 10
      mon = strcat('0',int2str(mont(d)));
   elseif mont(d) >= 10
      mon = int2str(mont(d));
   end
   
   for h = 0:3
  
      if h <= 1
         hour = strcat('0',int2str(h*6));
      elseif h >= 1
         hour = int2str(h*6);
      end
      
      disp(strcat(day,',',hour))
      
      temp_temp = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group3/horel_data/cfoster/wrf_daq_runs/',run,'/wrfout_d03_',int2str(year(d)),'-',mon,'-',day,'_',hour,':00:00'),'T2');
      
      temp = cat(3, temp,temp_temp);
      
      u_temp = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group3/horel_data/cfoster/wrf_daq_runs/',run,'/wrfout_d03_',int2str(year(d)),'-',mon,'-',day,'_',hour,':00:00'),'U');
      
      u = cat(3, u,squeeze(u_temp(2:end,:,1,:)));
      
      v_temp = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group3/horel_data/cfoster/wrf_daq_runs/',run,'/wrfout_d03_',int2str(year(d)),'-',mon,'-',day,'_',hour,':00:00'),'V');
      
      v = cat(3, v,squeeze(v_temp(:,2:end,1,:)));
      
      for hr = 0:5
          
         time_model = cat(1,time_model,datenum(year(d),mont(d),days(d),((h*6)+hr),0,0));
         
      end
   end   
end
   
temp(:,:,1) = [];
u(:,:,1) = [];
v(:,:,1) = [];
time_model(1) = [];   % Remove first time (repeat).

%%

temp_model = [reshape([temp(x_grid(1),y_grid(1),:)],[length(temp(1,1,:)),1]), ...
    reshape([temp(x_grid(2),y_grid(2),:)],[length(temp(1,1,:)),1]), ...
    reshape([temp(x_grid(3),y_grid(3),:)],[length(temp(1,1,:)),1]), ...
    reshape([temp(x_grid(4),y_grid(4),:)],[length(temp(1,1,:)),1]), ...
    reshape([temp(x_grid(5),y_grid(5),:)],[length(temp(1,1,:)),1])];

u_model = [reshape([u(x_grid(1),y_grid(1),:)],[length(u(1,1,:)),1]), ...
    reshape([u(x_grid(2),y_grid(2),:)],[length(u(1,1,:)),1]), ...
    reshape([u(x_grid(3),y_grid(3),:)],[length(u(1,1,:)),1]), ...
    reshape([u(x_grid(4),y_grid(4),:)],[length(u(1,1,:)),1]), ...
    reshape([u(x_grid(5),y_grid(5),:)],[length(u(1,1,:)),1])];
    
v_model = [reshape([v(x_grid(1),y_grid(1),:)],[length(v(1,1,:)),1]), ...
    reshape([v(x_grid(2),y_grid(2),:)],[length(v(1,1,:)),1]), ...
    reshape([v(x_grid(3),y_grid(3),:)],[length(v(1,1,:)),1]), ...
    reshape([v(x_grid(4),y_grid(4),:)],[length(v(1,1,:)),1]), ...
    reshape([v(x_grid(5),y_grid(5),:)],[length(v(1,1,:)),1])];

wspd_model = (u_model.^2 + v_model.^2).^(0.5);
wdir_model = cart2compass(u_model,v_model);

%% Bias

temp_model_mean = mean(temp_model-273.14)

ts_temp_QLN(find(isnan(ts_temp_QLN)))=0;
ts_temp_QHW(find(isnan(ts_temp_QHW)))=0;
ts_temp_VEL(find(isnan(ts_temp_VEL)))=0;
ts_temp_SLC(find(isnan(ts_temp_SLC)))=0;
ts_temp_LGU(find(isnan(ts_temp_LGU)))=0;

temp_stn_mean = [mean(ts_temp_QLN) mean(ts_temp_QHW) mean(ts_temp_VEL) ...
    mean(ts_temp_SLC) mean(ts_temp_LGU)];

temp_bias = temp_model_mean - temp_stn_mean

wspd_model_mean = mean(wspd_model)

ts_wspd_QLN(find(isnan(ts_wspd_QLN)))=0;
ts_wspd_QHW(find(isnan(ts_wspd_QHW)))=0;
ts_wspd_VEL(find(isnan(ts_wspd_VEL)))=0;
ts_wspd_SLC(find(isnan(ts_wspd_SLC)))=0;
ts_wspd_LGU(find(isnan(ts_wspd_LGU)))=0;

wspd_stn_mean = [mean(ts_wspd_QLN) mean(ts_wspd_QHW) mean(ts_wspd_VEL) ...
    mean(ts_wspd_SLC) mean(ts_wspd_LGU)];

wspd_bias = wspd_model_mean - wspd_stn_mean

%%

model_vec = datevec(time_model);

n = 1; 

for i = 1:length(time_model)
    if data_hour_QLN(n,2:4) == [model_vec(i,3:5)] 
        final_index_QLN(n) = i;
        n = n + 1;
    end
end

n = 1; 

for i = 1:length(time_model)
    if data_hour_QHW(n,2:4) == [model_vec(i,3:4) model_vec(i,5)+45] 
        final_index_QHW(n) = i;
        n = n + 1;
    end
end

n = 1; 

for i = 1:length(time_model)
    
    if data_hour_VEL(n,2:4) == [model_vec(i,3:4) model_vec(i,5)+53] 
        final_index_VEL(n) = i;
        n = n + 1;
    end
end

n = 1; 

for i = 1:length(time_model)
    if data_hour_SLC(n,2:4) == [model_vec(i,3:5)] 
        final_index_SLC(n) = i;
        n = n + 1;
    end
end

n = 1; 

for i = 1:length(time_model)
    if data_hour_LGU(n,2:4) == [model_vec(i,3:4) model_vec(i,5)+51] 
        final_index_LGU(n) = i;
        n = n + 1;
    end
end

%%
data_hour_QLN(find(isnan(data_hour_QLN(:,5))),5)=0;
data_hour_QHW(find(isnan(data_hour_QHW(:,5))),5)=0;

temp_hour_QLN = [data_hour_QLN(1:end-1,5) temp_model(final_index_QLN,1)-273.14];
temp_hour_QHW = [data_hour_QHW(:,5) temp_model(final_index_QHW,2)-273.14];
temp_hour_VEL = [data_hour_VEL(:,5) temp_model(final_index_VEL,3)-273.14];
temp_hour_SLC = [data_hour_SLC(:,5) temp_model(final_index_SLC,4)-273.14];
temp_hour_LGU = [data_hour_LGU(:,5) temp_model(final_index_LGU,5)-273.14];

u_hour_QLN = [u_model(final_index_QLN,1)];
u_hour_QHW = [u_model(final_index_QHW,2)];
u_hour_VEL = [u_model(final_index_VEL,3)];
u_hour_SLC = [u_model(final_index_SLC,4)];
u_hour_LGU = [u_model(final_index_LGU,5)];

v_hour_QLN = [v_model(final_index_QLN,1)];
v_hour_QHW = [v_model(final_index_QHW,2)];
v_hour_VEL = [v_model(final_index_VEL,3)];
v_hour_SLC = [v_model(final_index_SLC,4)];
v_hour_LGU = [v_model(final_index_LGU,5)];

dir_hour_QLN = cart2compass(u_hour_QLN,v_hour_QLN)
dir_hour_QHW = cart2compass(u_hour_QHW,v_hour_QHW)
dir_hour_VEL = cart2compass(u_hour_VEL,v_hour_VEL)
dir_hour_SLC = cart2compass(u_hour_SLC,v_hour_SLC)
dir_hour_LGU = cart2compass(u_hour_LGU,v_hour_LGU)

%%

temp_rmse = [(sqrt(mean((temp_hour_QLN(find(temp_hour_QLN(:,1)~=0),2) - temp_hour_QLN(find(temp_hour_QLN(:,1)~=0),1)).^2))) ...
    (sqrt(mean((temp_hour_QHW(find(temp_hour_QHW(:,1)~=0),2) - temp_hour_QHW(find(temp_hour_QHW(:,1)~=0),1)).^2))) ...
    (sqrt(mean((temp_hour_VEL(find(temp_hour_VEL(:,1)~=0),2) - temp_hour_VEL(find(temp_hour_VEL(:,1)~=0),1)).^2))) ...
    (sqrt(mean((temp_hour_SLC(find(temp_hour_SLC(:,1)~=0),2) - temp_hour_SLC(find(temp_hour_SLC(:,1)~=0),1)).^2))) ...
    (sqrt(mean((temp_hour_LGU(find(temp_hour_LGU(:,1)~=0),2) - temp_hour_LGU(find(temp_hour_LGU(:,1)~=0),1)).^2)))]

wspd_hour_QLN = (u_hour_QLN.^2 + v_hour_QLN.^2).^(0.5);
wspd_hour_QHW = (u_hour_QHW.^2 + v_hour_QHW.^2).^(0.5);
wspd_hour_VEL = (u_hour_VEL.^2 + v_hour_VEL.^2).^(0.5);
wspd_hour_SLC = (u_hour_SLC.^2 + v_hour_SLC.^2).^(0.5);
wspd_hour_LGU = (u_hour_LGU.^2 + v_hour_LGU.^2).^(0.5);

data_hour_QLN(find(isnan(data_hour_QLN(:,6))),6)=0;
data_hour_QHW(find(isnan(data_hour_QHW(:,6))),6)=0;
data_hour_VEL(find(isnan(data_hour_VEL(:,6))),6)=0;

wspd_rmse = [(sqrt(mean((wspd_hour_QLN - data_hour_QLN(1:end-1,6)).^2))) ...
    (sqrt(mean((wspd_hour_QHW - data_hour_QHW(:,6)).^2))) ...
    (sqrt(mean((wspd_hour_VEL - data_hour_VEL(:,6)).^2))) ...
    (sqrt(mean((wspd_hour_SLC - data_hour_SLC(:,6)).^2))) ...
    (sqrt(mean((wspd_hour_LGU - data_hour_LGU(:,6)).^2)))]

%%

figure(1)

subplot(2,3,1)
plot(time_model,temp_model(:,1)-273.15,'r-',data_hour_QLN(:,1),data_hour_QLN(:,5),'k-')
text(min(time_model)+7,-14,strcat('RMSE: ',num2str(temp_rmse(1))))
title('Lindon, Utah Valley')
xlim([min(time_model) max(time_model)])
ylim([-15 5])
datetick('x','keeplimits','keepticks')
ylabel('Temperature ^oC')

subplot(2,3,2)
plot(time_model,temp_model(:,2)-273.15,'r-',data_hour_QHW(:,1),data_hour_QHW(:,5),'k-')
text(min(time_model)+7,-14,strcat('RMSE: ',num2str(temp_rmse(2))))
title('Hawthorne, Salt Lake Valley')
xlim([min(time_model) max(time_model)])
ylim([-15 5])
datetick('x','keeplimits','keepticks')

subplot(2,3,3)
plot(time_model,temp_model(:,3)-273.15,'r-',data_hour_VEL(:,1),data_hour_VEL(:,5),'k-')
text(min(time_model)+7,-29,strcat('RMSE: ',num2str(temp_rmse(3))))
title('Vernal, Uintah Basin')
xlim([min(time_model) max(time_model)])
ylim([-30 0])
datetick('x','keeplimits','keepticks')
xlabel('Date')

subplot(2,3,4)
plot(time_model,temp_model(:,4)-273.15,'r-',data_hour_SLC(:,1),data_hour_SLC(:,5),'k-')
text(min(time_model)+7,-14,strcat('RMSE: ',num2str(temp_rmse(4))))
title('Salt Lake City, Salt Lake Valley')
xlim([min(time_model) max(time_model)])
ylim([-15 5])
datetick('x','keeplimits','keepticks')
ylabel('Temperature ^oC')
xlabel('Date')

subplot(2,3,5)
plot(time_model,temp_model(:,5)-273.15,'r-',data_hour_LGU(:,1),data_hour_LGU(:,5),'k-')
text(min(time_model)+7,-24,strcat('RMSE: ',num2str(temp_rmse(5))))
title('Logan, Cache Valley')
xlim([min(time_model) max(time_model)])
ylim([-25 5])
datetick('x','keeplimits','keepticks')
legend('Model','Station','Location','NorthWest')
xlabel('Date')


%%

figure(2)

subplot(2,3,1)
plot(time_model,temp_model(:,1)-273.15,'r-',data_hour_QLN(:,1),data_hour_QLN(:,5),'k-')
text(min(time_model)+15,-9,strcat('RMSE: ',num2str(temp_rmse(1))))
title('Lindon, Utah Valley')
xlim([min(time_model)+9 max(time_model)])
ylim([-10 15])
datetick('x','keeplimits','keepticks')
ylabel('Temperature ^oC')

subplot(2,3,2)
plot(time_model,temp_model(:,2)-273.15,'r-',data_hour_QHW(:,1),data_hour_QHW(:,5),'k-')
text(min(time_model)+15,-9,strcat('RMSE: ',num2str(temp_rmse(2))))
title('Hawthorne, Salt Lake Valley')
xlim([min(time_model)+9 max(time_model)])
ylim([-10 15])
datetick('x','keeplimits','keepticks')

subplot(2,3,3)
plot(time_model,temp_model(:,3)-273.15,'r-',data_hour_VEL(:,1),data_hour_VEL(:,5),'k-')
text(min(time_model)+15,-14,strcat('RMSE: ',num2str(temp_rmse(3))))
title('Vernal, Uintah Basin')
xlim([min(time_model)+9 max(time_model)])
ylim([-15 10])
datetick('x','keeplimits','keepticks')
xlabel('Date')

subplot(2,3,4)
plot(time_model,temp_model(:,4)-273.15,'r-',data_hour_SLC(:,1),data_hour_SLC(:,5),'k-')
text(min(time_model)+15,-9,strcat('RMSE: ',num2str(temp_rmse(4))))
title('Salt Lake City, Salt Lake Valley')
xlim([min(time_model)+9 max(time_model)])
ylim([-10 15])
datetick('x','keeplimits','keepticks')
ylabel('Temperature ^oC')
xlabel('Date')

subplot(2,3,5)
plot(time_model,temp_model(:,5)-273.15,'r-',data_hour_LGU(:,1),data_hour_LGU(:,5),'k-')
text(min(time_model)+15,-19,strcat('RMSE: ',num2str(temp_rmse(5))))
title('Logan, Cache Valley')
xlim([min(time_model)+9 max(time_model)])
ylim([-20 15])
datetick('x','keeplimits','keepticks')
legend('Model','Station','Location','NorthWest')
xlabel('Date')


%%

figure(3)

subplot(2,3,1)
plot(time_model,wspd_model(:,1),'r-',data_hour_QLN(:,1),data_hour_QLN(:,6),'k-')
text(min(time_model)+7,4.75,strcat('RMSE: ',num2str(wspd_rmse(1))))
title('Lindon, Utah Valley')
xlim([min(time_model) max(time_model)])
ylim([0 5])
datetick('x','keeplimits','keepticks')
ylabel('wind Speed m/s')

subplot(2,3,2)
plot(time_model,wspd_model(:,2),'r-',data_hour_QHW(:,1),data_hour_QHW(:,6),'k-')
text(min(time_model)+7,4.75,strcat('RMSE: ',num2str(wspd_rmse(2))))
title('Hawthorne, Salt Lake Valley')
xlim([min(time_model) max(time_model)])
ylim([0 5])
datetick('x','keeplimits','keepticks')

subplot(2,3,3)
plot(time_model,wspd_model(:,3),'r-',data_hour_VEL(:,1),data_hour_VEL(:,6),'k-')
text(min(time_model)+7,4.75,strcat('RMSE: ',num2str(wspd_rmse(3))))
title('Vernal, Uintah Basin')
xlim([min(time_model) max(time_model)])
ylim([0 5])
datetick('x','keeplimits','keepticks')
xlabel('Date')

subplot(2,3,4)
plot(time_model,wspd_model(:,4),'r-',data_hour_SLC(:,1),data_hour_SLC(:,6),'k-')
text(min(time_model)+7,4.75,strcat('RMSE: ',num2str(wspd_rmse(4))))
title('Salt Lake City, Salt Lake Valley')
xlim([min(time_model) max(time_model)])
ylim([0 5])
datetick('x','keeplimits','keepticks')
ylabel('Wind Speed m/s')
xlabel('Date')

subplot(2,3,5)
plot(time_model,wspd_model(:,5),'r-',data_hour_LGU(:,1),data_hour_LGU(:,6),'k-')
text(min(time_model)+7,4.75,strcat('RMSE: ',num2str(wspd_rmse(5))))
title('Logan, Cache Valley')
xlim([min(time_model) max(time_model)])
ylim([0 5])
datetick('x','keeplimits','keepticks')
legend('Model','Station','Location','NorthWest')
xlabel('Date')

%%

figure(4)

subplot(2,3,1)
plot(time_model,wspd_model(:,1),'r-',data_hour_QLN(:,1),data_hour_QLN(:,6),'k-')
text(min(time_model)+15,4.75,strcat('RMSE: ',num2str(wspd_rmse(1))))
title('Lindon, Utah Valley')
xlim([min(time_model)+9 max(time_model)])
ylim([0 5])
datetick('x','keeplimits','keepticks')
ylabel('wind Speed m/s')

subplot(2,3,2)
plot(time_model,wspd_model(:,2),'r-',data_hour_QHW(:,1),data_hour_QHW(:,6),'k-')
text(min(time_model)+15,4.75,strcat('RMSE: ',num2str(wspd_rmse(2))))
title('Hawthorne, Salt Lake Valley')
xlim([min(time_model)+9 max(time_model)])
ylim([0 5])
datetick('x','keeplimits','keepticks')

subplot(2,3,3)
plot(time_model,wspd_model(:,3),'r-',data_hour_VEL(:,1),data_hour_VEL(:,6),'k-')
text(min(time_model)+15,4.75,strcat('RMSE: ',num2str(wspd_rmse(3))))
title('Vernal, Uintah Basin')
xlim([min(time_model)+9 max(time_model)])
ylim([0 5])
datetick('x','keeplimits','keepticks')
xlabel('Date')

subplot(2,3,4)
plot(time_model,wspd_model(:,4),'r-',data_hour_SLC(:,1),data_hour_SLC(:,6),'k-')
text(min(time_model)+15,4.75,strcat('RMSE: ',num2str(wspd_rmse(4))))
title('Salt Lake City, Salt Lake Valley')
xlim([min(time_model)+9 max(time_model)])
ylim([0 5])
datetick('x','keeplimits','keepticks')
ylabel('Wind Speed m/s')
xlabel('Date')

subplot(2,3,5)
plot(time_model,wspd_model(:,5),'r-',data_hour_LGU(:,1),data_hour_LGU(:,6),'k-')
text(min(time_model)+15,4.75,strcat('RMSE: ',num2str(wspd_rmse(5))))
title('Logan, Cache Valley')
xlim([min(time_model)+9 max(time_model)])
ylim([0 5])
datetick('x','keeplimits','keepticks')
legend('Model','Station','Location','NorthWest')
xlabel('Date')

%%

figure(5)

subplot(2,3,1)
plot(time_model,wdir_model(:,1),'r.',data_hour_QLN(:,1),data_hour_QLN(:,7),'k.')
title('Lindon, Utah Valley')
axis([min(time_model) max(time_model) 0 360])
datetick('x','keepticks')
ylabel('wind Speed m/s')

subplot(2,3,2)
plot(time_model,wdir_model(:,2),'r.',data_hour_QHW(:,1),data_hour_QHW(:,7),'k.')
title('Hawthorne, Salt Lake Valley')
axis([min(time_model) max(time_model) 0 360])
datetick('x','keepticks')

subplot(2,3,3)
plot(time_model,wdir_model(:,3),'r.',data_hour_VEL(:,1),data_hour_VEL(:,7),'k.')
title('Vernal, Uintah Basin')
axis([min(time_model) max(time_model) 0 360])
datetick('x','keepticks')
xlabel('Date')

subplot(2,3,4)
plot(time_model,wdir_model(:,4),'r.',data_hour_SLC(:,1),data_hour_SLC(:,7),'k.')
title('Salt Lake City, Salt Lake Valley')
axis([min(time_model) max(time_model) 0 360])
datetick('x','keepticks')
ylabel('Wind Speed m/s')
xlabel('Date')

subplot(2,3,5)
plot(time_model,wdir_model(:,5),'r.',data_hour_LGU(:,1),data_hour_LGU(:,7),'k.')
title('Logan, Cache Valley')
axis([min(time_model) max(time_model) 0 360])
datetick('x','keepticks')
legend('Model','Station','Location','NorthWest')
xlabel('Date')
