clear all
close all

% Grab all variables of interest & store them in appropriately named arrays.

run2 = 'nlcd2011-dec31';
run = 'usgs-dec31';

load('station_data.mat');

time_stn2 = data(:,1);
temp_stn = data(:,2);
spd_stn = data(:,3);
dir_stn = data(:,4);

for i = 1:length(loc)
    lat_stns(i) = str2num(char(loc(i,1)));
    lon_stns(i) = str2num(char(loc(i,2)));
end

load('data/iop5.mat');
load('albedo_isfs.mat');

data(find(data==0)) = nan;

long = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run,'/wrfout_d03_2011-01-01_00:00:00'),'XLONG');
lati = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run,'/wrfout_d03_2011-01-01_00:00:00'),'XLAT');

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

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

% Station location variables.

id  = {'ISFS1',   'ISFS2',   'ISFS3',   'ISFS4',   'ISFS5',   'ISFS6',   'ISFS7',   'ISS-N',  'ISS-S'};
lon = [-112.1031, -111.9136, -111.8448, -111.9946, -111.8045, -112.0606, -111.9729, -111.925, -111.9252, lon_stns];
lat = [40.7785,   40.7374,   40.6980,    40.6810,   40.6127,   40.5847,  40.5140,   40.61,    40.6006,   lat_stns];

% Find closest grid point to station locations

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

% Find the lon/lat associated with those grid points.

lon_grid = [lons(y_grid(1),(x_grid(1))), lons(y_grid(2),(x_grid(2))), lons(y_grid(3),(x_grid(3))), lons(y_grid(4),(x_grid(4))), lons(y_grid(5),(x_grid(5))), lons(y_grid(6),(x_grid(6))), lons(y_grid(7),(x_grid(7)))];
lat_grid = [lats(y_grid(1),(x_grid(1))), lats(y_grid(2),(x_grid(2))), lats(y_grid(3),(x_grid(3))), lats(y_grid(4),(x_grid(4))), lats(y_grid(5),(x_grid(5))), lats(y_grid(6),(x_grid(6))), lats(y_grid(7),(x_grid(7)))];

% Grab the model data.

days = [1, 2, 3, 4, 5, 6, 7, 8];
year = [2011,2011,2011,2011,2011,2011,2011,2011];
mont = [1,1,1,1,1,1,1,1];

temp = zeros(size(lati));
temp2 = zeros(size(lati));
u = zeros(size(lati));
u2 = zeros(size(lati));
v = zeros(size(lati));
v2 = 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-group4/model/cfoster/output/',run,'/wrfout_d03_',int2str(year(d)),'-',mon,'-',day,'_',hour,':00:00'),'T2');
      
      temp = cat(3, temp,temp_temp);
      
      temp2_temp = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run2,'/wrfout_d03_',int2str(year(d)),'-',mon,'-',day,'_',hour,':00:00'),'T2');
      
      temp2 = cat(3, temp2,temp2_temp);
      
      u_temp = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run,'/wrfout_d03_',int2str(year(d)),'-',mon,'-',day,'_',hour,':00:00'),'U');
      
      u = cat(3, u,squeeze(u_temp(2:end,:,1,:)));
      
      u2_temp = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run2,'/wrfout_d03_',int2str(year(d)),'-',mon,'-',day,'_',hour,':00:00'),'U');
      
      u2 = cat(3, u2,squeeze(u2_temp(2:end,:,1,:)));
      
      v_temp = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run,'/wrfout_d03_',int2str(year(d)),'-',mon,'-',day,'_',hour,':00:00'),'V');
      
      v = cat(3, v,squeeze(v_temp(:,2:end,1,:)));
      
      v2_temp = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run2,'/wrfout_d03_',int2str(year(d)),'-',mon,'-',day,'_',hour,':00:00'),'V');
      
      v2 = cat(3, v2,squeeze(v2_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) = [];
temp2(:,:,1) = [];
u(:,:,1) = [];
u2(:,:,1) = [];
v(:,:,1) = [];
v2(:,:,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]), ...
    reshape([temp(x_grid(6),y_grid(6),:)],[length(temp(1,1,:)),1]), ...
    reshape([temp(x_grid(7),y_grid(7),:)],[length(temp(1,1,:)),1]), ...
    reshape([temp(x_grid(16),y_grid(16),:)],[length(temp(1,1,:)),1]), ... % Logan 
    reshape([temp(x_grid(17),y_grid(17),:)],[length(temp(1,1,:)),1]), ... % Vernal 
    reshape([temp(x_grid(18),y_grid(18),:)],[length(temp(1,1,:)),1]), ... % Roosevelt
    reshape([temp(x_grid(30),y_grid(30),:)],[length(temp(1,1,:)),1])] - 273.14; % Syracuse

temp2_model = [reshape([temp2(x_grid(1),y_grid(1),:)],[length(temp2(1,1,:)),1]), ...
    reshape([temp2(x_grid(2),y_grid(2),:)],[length(temp2(1,1,:)),1]), ...
    reshape([temp2(x_grid(3),y_grid(3),:)],[length(temp2(1,1,:)),1]), ...
    reshape([temp2(x_grid(4),y_grid(4),:)],[length(temp2(1,1,:)),1]), ...
    reshape([temp2(x_grid(5),y_grid(5),:)],[length(temp2(1,1,:)),1]), ...
    reshape([temp2(x_grid(6),y_grid(6),:)],[length(temp2(1,1,:)),1]), ...
    reshape([temp2(x_grid(7),y_grid(7),:)],[length(temp2(1,1,:)),1]), ...
    reshape([temp2(x_grid(16),y_grid(16),:)],[length(temp2(1,1,:)),1]), ... % Logan
    reshape([temp2(x_grid(17),y_grid(17),:)],[length(temp2(1,1,:)),1]), ... % Vernal     
    reshape([temp2(x_grid(18),y_grid(18),:)],[length(temp2(1,1,:)),1]), ... % Roosevelt
    reshape([temp2(x_grid(30),y_grid(30),:)],[length(temp2(1,1,:)),1])] - 273.14; % Syracuse

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]), ...
    reshape([u(x_grid(6),y_grid(6),:)],[length(u(1,1,:)),1]), ...
    reshape([u(x_grid(7),y_grid(7),:)],[length(u(1,1,:)),1]), ...
    reshape([u(x_grid(16),y_grid(16),:)],[length(u(1,1,:)),1]), ... % Logan 
    reshape([u(x_grid(17),y_grid(17),:)],[length(u(1,1,:)),1]), ... % Vernal 
    reshape([u(x_grid(18),y_grid(18),:)],[length(u(1,1,:)),1]), ... % Roosevelt
    reshape([u(x_grid(30),y_grid(30),:)],[length(u(1,1,:)),1])]; % Syracuse

u2_model = [reshape([u2(x_grid(1),y_grid(1),:)],[length(u2(1,1,:)),1]), ...
    reshape([u2(x_grid(2),y_grid(2),:)],[length(u2(1,1,:)),1]), ...
    reshape([u2(x_grid(3),y_grid(3),:)],[length(u2(1,1,:)),1]), ...
    reshape([u2(x_grid(4),y_grid(4),:)],[length(u2(1,1,:)),1]), ...
    reshape([u2(x_grid(5),y_grid(5),:)],[length(u2(1,1,:)),1]), ...
    reshape([u2(x_grid(6),y_grid(6),:)],[length(u2(1,1,:)),1]), ...
    reshape([u2(x_grid(7),y_grid(7),:)],[length(u2(1,1,:)),1]), ...
    reshape([u2(x_grid(16),y_grid(16),:)],[length(u2(1,1,:)),1]), ... % Logan
    reshape([u2(x_grid(17),y_grid(17),:)],[length(u2(1,1,:)),1]), ... % Vernal     
    reshape([u2(x_grid(18),y_grid(18),:)],[length(u2(1,1,:)),1]), ... % Roosevelt
    reshape([u2(x_grid(30),y_grid(30),:)],[length(u2(1,1,:)),1])]; % Syracuse

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]), ...
    reshape([v(x_grid(6),y_grid(6),:)],[length(v(1,1,:)),1]), ...
    reshape([v(x_grid(7),y_grid(7),:)],[length(v(1,1,:)),1]), ...
    reshape([v(x_grid(16),y_grid(16),:)],[length(v(1,1,:)),1]), ... % Logan 
    reshape([v(x_grid(17),y_grid(17),:)],[length(v(1,1,:)),1]), ... % Vernal 
    reshape([v(x_grid(18),y_grid(18),:)],[length(v(1,1,:)),1]), ... % Roosevelt
    reshape([v(x_grid(30),y_grid(30),:)],[length(v(1,1,:)),1])]; % Syracuse

v2_model = [reshape([v2(x_grid(1),y_grid(1),:)],[length(v2(1,1,:)),1]), ...
    reshape([v2(x_grid(2),y_grid(2),:)],[length(v2(1,1,:)),1]), ...
    reshape([v2(x_grid(3),y_grid(3),:)],[length(v2(1,1,:)),1]), ...
    reshape([v2(x_grid(4),y_grid(4),:)],[length(v2(1,1,:)),1]), ...
    reshape([v2(x_grid(5),y_grid(5),:)],[length(v2(1,1,:)),1]), ...
    reshape([v2(x_grid(6),y_grid(6),:)],[length(v2(1,1,:)),1]), ...
    reshape([v2(x_grid(7),y_grid(7),:)],[length(v2(1,1,:)),1]), ...
    reshape([v2(x_grid(16),y_grid(16),:)],[length(v2(1,1,:)),1]), ... % Logan
    reshape([v2(x_grid(17),y_grid(17),:)],[length(v2(1,1,:)),1]), ... % Vernal     
    reshape([v2(x_grid(18),y_grid(18),:)],[length(v2(1,1,:)),1]), ... % Roosevelt
    reshape([v2(x_grid(30),y_grid(30),:)],[length(v2(1,1,:)),1])]; % Syracuse

%% Get station data.

% ISFS times.

time_isfs = datevec(data(:,1));

time_temp1 = [time_stn2{1}]; % ISFS 1
time_temp = time_temp1(1:length(spd_stn{1}));
for t = 1:length(time_temp)
   time_isfs1(t,1:6) = datevec(strcat(time_temp{t}(1:10),{' '},time_temp{t}(12:19)));
end

time_temp = [time_stn2{2}]; % ISFS 2
for t = 1:length(time_temp)
   time_isfs2(t,1:6) = datevec(strcat(time_temp{t}(1:10),{' '},time_temp{t}(12:19)));
end

time_temp = [time_stn2{6}]; % ISFS 6
for t = 1:length(time_temp)
   time_isfs6(t,1:6) = datevec(strcat(time_temp{t}(1:10),{' '},time_temp{t}(12:19)));
end

time_temp = [time_stn2{21}]; % Syracuse
for t = 1:length(time_temp)
   time_syr(t,1:6) = datevec(strcat(time_temp{t}(1:10),{' '},time_temp{t}(12:19)));
end


% Parse data to hourly.

hour = find(time_isfs(:,5) == 3);
isfs_hour = time_isfs(hour,:);

hour_isfs1 = find(time_isfs1(:,5) == 3);
isfs1_hour = time_isfs1(hour_isfs1,:);

hour_isfs2 = find(time_isfs2(:,5) == 3);
isfs2_hour = time_isfs2(hour_isfs2,:);

hour_isfs6 = find(time_isfs6(:,5) == 0 | time_isfs6(:,5) == 3);
isfs6_hour = time_isfs6(hour_isfs6,:);

hour_syr = find(time_syr(:,5) == 0);
syr_hour = time_syr(hour_syr,:);

% Match observation with hourly wrf output.

model_vec = datevec(time_model);

%%

% ISFS 

n = 1; 

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

% ISFS 1

n = 1; 

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

% ISFS 2

n = 1; 

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

% ISFS 6

n = 1; 

for i = 1:length(isfs6_hour)
    if isfs6_hour(i,2:5) == [model_vec(n,2:4) model_vec(n,5)+3] | isfs6_hour(i,2:5) == [model_vec(n,2:5)]
        final_index6(n) = i;
        n = n + 1;
    end
end

% Syracuse

n = 1; 

for i = 1:length(syr_hour)
    if syr_hour(i,2:5) == model_vec(n,2:5)  
        final_syr(n) = i;
        n = n + 1;
    end
end

%%

temp_int1 = data(:,2);
temp_stn1 = temp_int1(hour);
spd_int1 = spd_stn{1};
spd_stn1 = spd_int1(hour_isfs1);
dir_int1 = dir_stn{1};
dir_stn1 = dir_int1(hour_isfs1);
isfs1_final = [temp_model(:,1) temp2_model(:,1) temp_stn1(final_index) ...
    u_model(:,1) u2_model(:,1) v_model(:,1) v2_model(:,1) ...
    spd_stn1(final_index1) dir_stn1(final_index1)];

temp_int2 = data(:,3);
temp_stn2 = temp_int2(hour);
spd_int2 = spd_stn{2};
spd_stn2 = spd_int2(hour_isfs2);
dir_int2 = dir_stn{2};
dir_stn2 = dir_int2(hour_isfs2);
isfs2_final = [temp_model(:,2) temp2_model(:,2) temp_stn2(final_index) ...
    u_model(:,2) u2_model(:,2) v_model(:,2) v2_model(:,2)...
    spd_stn2(final_index2) dir_stn2(final_index2)];

temp_int3 = data(:,4);
temp_stn3 = temp_int3(hour);
isfs3_final = [temp_model(:,3) temp2_model(:,3) temp_stn3(final_index) ...
    u_model(:,1) u2_model(:,1) v_model(:,1) v2_model(:,1)];

temp_int4 = data(:,5);
temp_stn4 = temp_int4(hour);
isfs4_final = [temp_model(:,4) temp2_model(:,4) temp_stn4(final_index) ...
    u_model(:,1) u2_model(:,1) v_model(:,1) v2_model(:,1)];

temp_int5 = data(:,6);
temp_stn5 = temp_int5(hour);
isfs5_final = [temp_model(:,5) temp2_model(:,5) temp_stn5(final_index) ...
    u_model(:,1) u2_model(:,1) v_model(:,1) v2_model(:,1)];

temp_int6 = data(:,7);
temp_stn6 = temp_int6(hour);
isfs6_final = [temp_model(:,6) temp2_model(:,6) temp_stn6(final_index) ...
    u_model(:,6) u2_model(:,6) v_model(:,6) v2_model(:,6)];

temp_int7 = data(:,8);
temp_stn7 = temp_int7(hour);
isfs7_final = [temp_model(:,7) temp2_model(:,7) temp_stn7(final_index) ...
    u_model(:,1) u2_model(:,1) v_model(:,1) v2_model(:,1)];

temp_syr = temp_stn{21};
spd_syr = spd_stn{21};
spd_syr_fin = spd_syr(hour_syr);
dir_syr = dir_stn{21};
dir_syr_fin = dir_syr(hour_syr);
syr_final = [temp_model(final_syr,11) temp2_model(final_syr,11) temp_syr(hour_syr(1:length(final_syr))) ...   
    u_model(:,11) u2_model(:,11) v_model(:,11) v2_model(:,11) ...
    spd_stn1(final_syr) dir_stn1(final_syr)];

%% Calculate RMSE.

base_rmse1 = sqrt(mean((isfs1_final(24:end,2) - isfs1_final(24:end,3)).^2));   
base_rmse2 = sqrt(mean((isfs2_final(:,2) - isfs2_final(:,3)).^2));   
base_rmse3 = sqrt(mean((isfs3_final(24:end,2) - isfs3_final(24:end,3)).^2));   
base_rmse4 = sqrt(mean((isfs4_final(24:end,2) - isfs4_final(24:end,3)).^2));   
base_rmse5 = sqrt(mean((isfs5_final(:,2) - isfs5_final(:,3)).^2));   
base_rmse6 = sqrt(mean((isfs6_final(:,2) - isfs6_final(:,3)).^2));   
base_rmse7 = sqrt(mean((isfs7_final(:,2) - isfs7_final(:,3)).^2));   
base_rmse_syr = sqrt(mean((syr_final(:,2) - syr_final(:,3)).^2)); 

usgs_rmse1 = sqrt(mean((isfs1_final(24:end,1) - isfs1_final(24:end,3)).^2));   
usgs_rmse2 = sqrt(mean((isfs2_final(:,1) - isfs2_final(:,3)).^2));   
usgs_rmse3 = sqrt(mean((isfs3_final(24:end,1) - isfs3_final(24:end,3)).^2));   
usgs_rmse4 = sqrt(mean((isfs4_final(24:end,1) - isfs4_final(24:end,3)).^2));   
usgs_rmse5 = sqrt(mean((isfs5_final(:,1) - isfs5_final(:,3)).^2));   
usgs_rmse6 = sqrt(mean((isfs6_final(:,1) - isfs6_final(:,3)).^2));
usgs_rmse7 = sqrt(mean((isfs7_final(:,1) - isfs7_final(:,3)).^2));
usgs_rmse_syr = sqrt(mean((syr_final(:,1) - syr_final(:,3)).^2)); 

%% Calculate model wind speed. 

% [ang1_usgs,spd1_usgs] = cart2compass(isfs1_final(:,4),isfs1_final(:,6));
% [ang1_base,spd1_base] = cart2compass(isfs1_final(:,5),isfs1_final(:,7));
% 
% [ang2_usgs,spd2_usgs] = cart2compass(isfs2_final(:,4),isfs2_final(:,6));
% [ang2_base,spd2_base] = cart2compass(isfs2_final(:,5),isfs2_final(:,7));
% 
% [ang3_usgs,spd3_usgs] = cart2compass(isfs3_final(:,4),isfs3_final(:,6));
% [ang3_base,spd3_base] = cart2compass(isfs3_final(:,5),isfs3_final(:,7));
% 
% [ang4_usgs,spd4_usgs] = cart2compass(isfs4_final(:,4),isfs4_final(:,6));
% [ang4_base,spd4_base] = cart2compass(isfs4_final(:,5),isfs4_final(:,7));
% 
% [ang5_usgs,spd5_usgs] = cart2compass(isfs5_final(:,4),isfs5_final(:,6));
% [ang5_base,spd5_base] = cart2compass(isfs5_final(:,5),isfs5_final(:,7));
% 
% [ang6_usgs,spd6_usgs] = cart2compass(isfs6_final(:,4),isfs6_final(:,6));
% [ang6_base,spd6_base] = cart2compass(isfs6_final(:,5),isfs6_final(:,7));
% 
% [ang_syr_usgs,spd_syr_usgs] = cart2compass(syr_final(:,4),syr_final(:,6));
% [ang_syr_base,spd_syr_base] = cart2compass(syr_final(:,5),syr_final(:,7));
% 

spd_int6 = spd_stn{6};
spd_stn6 = spd_int6(hour_isfs6);
dir_int6 = dir_stn{6};
dir_stn6 = dir_int6(hour_isfs6);

rad = 4.0*atan(1.0)/180 
u_stn_1 = -1.*isfs1_final(:,8).*sin(rad.*(isfs1_final(:,9))); 
v_stn_1 = -1.*isfs1_final(:,8).*cos(rad.*(isfs1_final(:,9))); 

u_stn_2 = -1.*isfs2_final(:,8).*sin(rad.*(isfs2_final(:,9))); 
v_stn_2 = -1.*isfs2_final(:,8).*cos(rad.*(isfs2_final(:,9))); 

u_stn_6 = -1.*spd_stn6(:).*sin(rad.*(dir_stn6(:)));
v_stn_6 = -1.*spd_stn6(:).*cos(rad.*(dir_stn6(:)));

u_stn_syr = -1.*syr_final(:,8).*sin(rad.*(syr_final(:,9))); 
v_stn_syr = -1.*syr_final(:,8).*cos(rad.*(syr_final(:,9))); 
%% Plots

figure(1)
subplot(2,4,1)
plot(datetime(isfs_hour(final_index,:)),isfs1_final(:,1),'r-', ...
     datetime(isfs_hour(final_index,:)),isfs1_final(:,2),'b-', ...
     datetime(isfs_hour(final_index,:)),isfs1_final(:,3),'k-')
title('ISFS 1') 
ylim([-25 5])
text(max(datenum(isfs_hour(final_index,:)))-3.5,-21,'BASE RMSE = 3.49')
text(max(datenum(isfs_hour(final_index,:)))-3.5,-23,'DFLT RMSE = 4.49')

subplot(2,4,2)
plot(datetime(isfs_hour(final_index,:)),isfs2_final(:,1),'r-', ...
     datetime(isfs_hour(final_index,:)),isfs2_final(:,2),'b-', ...
     datetime(isfs_hour(final_index,:)),isfs2_final(:,3),'k-')
title('ISFS 2') 
ylim([-25 5])
text(max(datenum(isfs_hour(final_index,:)))-3.5,-21,'BASE RMSE = 1.73')
text(max(datenum(isfs_hour(final_index,:)))-3.5,-23,'DFLT RMSE = 3.61') 

subplot(2,4,3)
plot(datetime(isfs_hour(final_index,:)),isfs3_final(:,1),'r-', ...
     datetime(isfs_hour(final_index,:)),isfs3_final(:,2),'b-', ...
     datetime(isfs_hour(final_index,:)),isfs3_final(:,3),'k-')
title('ISFS 3') 
ylim([-25 5])
text(max(datenum(isfs_hour(final_index,:)))-3.5,-21,'BASE RMSE = 1.78')
text(max(datenum(isfs_hour(final_index,:)))-3.5,-23,'DFLT RMSE = 2.21')
 
subplot(2,4,4)
plot(datetime(isfs_hour(final_index,:)),isfs4_final(:,1),'r-', ...
     datetime(isfs_hour(final_index,:)),isfs4_final(:,2),'b-', ...
     datetime(isfs_hour(final_index,:)),isfs4_final(:,3),'k-')
title('ISFS 4') 
ylim([-25 5])
text(max(datenum(isfs_hour(final_index,:)))-3.5,-21,'BASE RMSE = 2.12')
text(max(datenum(isfs_hour(final_index,:)))-3.5,-23,'DFLT RMSE = 2.89')

subplot(2,4,5)
plot(datetime(isfs_hour(final_index,:)),isfs5_final(:,1),'r-', ...
     datetime(isfs_hour(final_index,:)),isfs5_final(:,2),'b-', ...
     datetime(isfs_hour(final_index,:)),isfs5_final(:,3),'k-')
title('ISFS 5') 
ylim([-25 5])
text(max(datenum(isfs_hour(final_index,:)))-3.5,-21,'BASE RMSE = 2.11')
text(max(datenum(isfs_hour(final_index,:)))-3.5,-23,'DFLT RMSE = 2.63')
 
subplot(2,4,6)
plot(datetime(isfs_hour(final_index,:)),isfs6_final(:,1),'r-', ...
     datetime(isfs_hour(final_index,:)),isfs6_final(:,2),'b-', ...
     datetime(isfs_hour(final_index,:)),isfs6_final(:,3),'k-')
title('ISFS 6') 
ylim([-25 5])
text(max(datenum(isfs_hour(final_index,:)))-3.5,-21,'BASE RMSE = 3.32')
text(max(datenum(isfs_hour(final_index,:)))-3.5,-23,'DFLT RMSE = 4.61')

subplot(2,4,7)
plot(datetime(isfs_hour(final_index,:)),isfs7_final(:,1),'r-', ...
     datetime(isfs_hour(final_index,:)),isfs7_final(:,2),'b-', ...
     datetime(isfs_hour(final_index,:)),isfs7_final(:,3),'k-')
title('ISFS 7') 
ylim([-25 5])
text(max(datenum(isfs_hour(final_index,:)))-3.5,-21,'BASE RMSE = 1.88')
text(max(datenum(isfs_hour(final_index,:)))-3.5,-23,'DFLT RMSE = 2.71')

subplot(2,4,8)
plot(datetime(syr_hour(1:length(final_syr),:)),syr_final(:,1),'r-', ...
     datetime(syr_hour(1:length(final_syr),:)),syr_final(:,2),'b-', ...
     datetime(syr_hour(1:length(final_syr),:)),syr_final(:,3),'k-')
title('Syracuse') 
ylim([-25 5])
text(max(datenum(syr_hour(1:length(final_syr),:)))-3.5,-21,'BASE RMSE = 2.59')
text(max(datenum(syr_hour(1:length(final_syr),:)))-3.5,-23,'DFLT RMSE = 4.85')




%%

figure('Position',[100, 100, 900, 800])

s1 = subplot(3,1,1)

p = get(s1,'position');
p(4) = p(4)*1.25; 
set(s1, 'position', p);

quiver(datenum(isfs_hour(final_index(24:1:96),:)),1.5+zeros(length(datenum(isfs1_hour(24:1:96,:))),1), ...
    u_stn_1(24:1:96),v_stn_1(24:1:96),'k','LineWidth',1.25)
hold on
quiver(datenum(isfs_hour(final_index(24:1:96),:)),1+zeros(length(datenum(isfs_hour(24:1:96,:))),1), ...
    isfs1_final(24:1:96,4),isfs1_final(24:1:96,6),'b','LineWidth',1.25)
hold on
quiver(datenum(isfs_hour(final_index(24:1:96),:)),0.5+zeros(length(datenum(isfs_hour(24:1:96,:))),1), ...
    isfs1_final(24:1:96,5),isfs1_final(24:1:96,7),'r','LineWidth',1.25)

datetick('x','mmm dd-hhZ')
axis([datenum(isfs_hour(final_index(20),:)) datenum(isfs_hour(final_index(100),:)) 0 2])
set(gca,'YTickLabel','')

text(datenum(isfs_hour(final_index(20),:))+0.01,1.9,'(a)','FontSize',12)

s2 = subplot(3,1,2)

p = get(s2,'position');
p(4) = p(4)*1.25; 
set(s2, 'position', p);

quiver(datenum(isfs_hour(final_index(24:1:96),:)),1.5+zeros(length(datenum(isfs1_hour(24:1:96,:))),1), ...
    u_stn_2(24:1:96),v_stn_2(24:1:96),'k','LineWidth',1.25)
hold on
quiver(datenum(isfs_hour(final_index(24:1:96),:)),1+zeros(length(datenum(isfs_hour(24:1:96,:))),1), ...
    isfs2_final(24:1:96,4),isfs2_final(24:1:96,6),'b','LineWidth',1.25)
hold on
quiver(datenum(isfs_hour(final_index(24:1:96),:)),0.5+zeros(length(datenum(isfs_hour(24:1:96,:))),1), ...
    isfs2_final(24:1:96,5),isfs2_final(24:1:96,7),'r','LineWidth',1.25)

datetick('x','mmm dd-hhZ')
axis([datenum(isfs_hour(final_index(20),:)) datenum(isfs_hour(final_index(100),:)) 0 2])
set(gca,'YTickLabel','')

text(datenum(isfs_hour(final_index(20),:))+0.01,1.9,'(b)','FontSize',12)

s3 = subplot(3,1,3)

p = get(s3,'position');
p(4) = p(4)*1.25; 
set(s3, 'position', p);

quiver(datenum(isfs6_hour((24:1:96),:)),1.5+zeros(length(datenum(isfs6_hour((24:1:96),:))),1), ...
    u_stn_6(24:1:96),v_stn_6(24:1:96),'k','LineWidth',1.25)
hold on
quiver(datenum(isfs_hour(final_index(24:1:96),:)),1+zeros(length(datenum(isfs_hour(24:1:96,:))),1), ...
    isfs6_final(24:1:96,4),isfs6_final(24:1:96,6),'b','LineWidth',1.25)
hold on
quiver(datenum(isfs_hour(final_index(24:1:96),:)),0.5+zeros(length(datenum(isfs_hour(24:1:96,:))),1), ...
    isfs6_final(24:1:96,5),isfs6_final(24:1:96,7),'r','LineWidth',1.25)


% quiver(datenum(isfs_hour(final_index(24:1:96),:)),1.5+zeros(length(datenum(isfs1_hour(24:1:96,:))),1), ...
%     u_stn_syr(24:1:96),v_stn_syr(24:1:96),'k','LineWidth',1.25)
% hold on
% quiver(datenum(isfs_hour(final_index(24:1:96),:)),1+zeros(length(datenum(isfs_hour(24:1:96,:))),1), ...
%     syr_final(24:1:96,4),syr_final(24:1:96,6),'b','LineWidth',1.25)
% hold on
% quiver(datenum(isfs_hour(final_index(24:1:96),:)),0.5+zeros(length(datenum(isfs_hour(24:1:96,:))),1), ...
%     syr_final(24:1:96,5),syr_final(24:1:96,7),'r','LineWidth',1.25)

datetick('x','mmm dd-hhZ')
axis([datenum(isfs_hour(final_index(20),:)) datenum(isfs_hour(final_index(100),:)) 0 2])
set(gca,'YTickLabel','')

legend('Station Data','BASE','DFLT','location','southwest')
xlabel('Date (Ticks at 00 UTC)')

text(datenum(isfs_hour(final_index(20),:))+0.01,1.9,'(c)','FontSize',12)
%%

figure(2)
[hax,ha1,ha2] = plotyy([datetime(isfs1_hour(1:length(final_index1),:)), ...
    datetime(isfs1_hour(1:length(final_index1),:)), ...
    datetime(isfs1_hour(1:length(final_index1),:))], ...
    [isfs1_final(:,8),spd1_usgs,spd1_base], ...
    [datetime(isfs1_hour(1:length(final_index1),:)), ...
    datetime(isfs1_hour(1:length(final_index1),:)), ...
    datetime(isfs1_hour(1:length(final_index1),:))], ...
    [isfs1_final(:,9),ang1_usgs,ang1_base])
set(ha1(1),'LineStyle','-','Color','k')
set(ha1(2),'LineStyle','-','Color','r')
set(ha1(3),'LineStyle','-','Color','b')
set(ha2(1),'LineStyle','none','Marker','*','Color','k')
set(ha2(2),'LineStyle','none','Marker','.','Color','r')
set(ha2(3),'LineStyle','none','Marker','.','Color','b')

figure(2)
[hax,ha1,ha2] = plotyy([datetime(isfs1_hour(1:length(final_index1),:))], [isfs1_final(:,8)], ...
    [datetime(isfs1_hour(1:length(final_index1),:))],[isfs1_final(:,9)])
set(ha1(1),'LineStyle','-','Color','k')
set(ha2(1),'LineStyle','none','Marker','*','Color','k')
