clear all
close all

% Get static data of interest.

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

terrain = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run1,'/wrfout_d03_2011-01-01_00:00:00'),'HGT');
mask = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run1,'/wrfout_d03_2011-01-01_00:00:00'),'LANDMASK');
lon = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run1,'/wrfout_d03_2011-01-01_00:00:00'),'XLONG');
lat = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run1,'/wrfout_d03_2011-01-01_00:00:00'),'XLAT');
land = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run1,'/wrfout_d03_2011-01-01_00:00:00'),'IVGTYP');
flg = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run1,'/wrfout_d03_2011-01-01_00:00:00'),'LANDMASK');

lat = double(lat(:,:,1));
lon = double(lon(:,:,1));
terrain = double(terrain(:,:,1));
flg = double(flg(:,:,1));

latlim = [min(lat(:)) max(lat(:))];
lonlim = [min(lon(:)) max(lon(:))];


% Loop through WRF output files to get variables of interest.

% Vectors of days and then the year and month associated with each day.

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

run1_u = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run1,'/wrfout_d03_2010-12-31_00:00:00'),'U');
run2_u = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run2,'/wrfout_d03_2011-01-01_00:00:00'),'U');

run1_v = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run1,'/wrfout_d03_2010-12-31_00:00:00'),'V');
run2_v = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run2,'/wrfout_d03_2011-01-01_00:00:00'),'V');

geo   = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run1,'/wrfout_d03_2011-01-01_00:00:00'),'PH');
geo_p = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run1,'/wrfout_d03_2011-01-01_00:00:00'),'PHB');

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 from WRF.
      
      u_temp = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run1,'/wrfout_d03_',int2str(year(d)),'-',mon,'-',day,'_',hour,':00:00'),'U');
      
      run1_u = cat(4,run1_u,u_temp);      
      
      u_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');
      
      run2_u = cat(4,run2_u,u_temp);
      
      % Q Cloud from WRF.
      
      v_temp = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run1,'/wrfout_d03_',int2str(year(d)),'-',mon,'-',day,'_',hour,':00:00'),'V');
      
      run1_v = cat(4,run1_v,v_temp);      
      
      v_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');
      
      run2_v = cat(4,run2_v,v_temp);
      
      temp_geo = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run1,'/wrfout_d03_',int2str(year(d)),'-',mon,'-',day,'_',hour,':00:00'),'PH');
      
      geo = cat(4,geo,temp_geo); 
      
      temp_geo_p = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run1,'/wrfout_d03_',int2str(year(d)),'-',mon,'-',day,'_',hour,':00:00'),'PHB');
      
      geo_p = cat(4,geo_p,temp_geo_p);   
      
      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
   
run1_u(:,:,:,1:6) = []; % Remove first time (repeat).
run2_u(:,:,:,1:6) = [];
run1_v(:,:,:,1:6) = [];
run2_v(:,:,:,1:6) = [];
time_model(1) = [];   % Remove first time (repeat).

geo(:,:,:,1:6) = [];
geo_p(:,:,:,1:6) = [];

%% Time heights

hgt = (geo + geo_p) ./ 9.8;

x_pt = 225;
y_pt = 302;
top = 25;
  
% USGS

th_u = run2_u(x_pt-1:x_pt+1,y_pt-1:y_pt+1,1:25,:);
th_v = run2_v(x_pt-1:x_pt+1,y_pt-1:y_pt+1,1:25,:);
th_hgt = hgt(x_pt-1:x_pt+1,y_pt-1:y_pt+1,1:25,:);

th_avg_u = mean(mean(th_u,1),2);
th_avg_v = mean(mean(th_v,1),2);
th_avg_hgt = mean(mean(th_hgt,1),2);

u_fin = reshape(th_avg_u,[top,length(time_model)]);
v_fin = reshape(th_avg_v,[top,length(time_model)]);
hgt_fin = reshape(th_avg_hgt,[top,length(time_model)]);

% NLCD

th_u1 = run1_u(x_pt-1:x_pt+1,y_pt-1:y_pt+1,1:25,:);
th_v1 = run1_v(x_pt-1:x_pt+1,y_pt-1:y_pt+1,1:25,:);

th_avg_u1 = mean(mean(th_u1,1),2);
th_avg_v1 = mean(mean(th_v1,1),2);

u_fin1 = reshape(th_avg_u1,[top,length(time_model)]);
v_fin1 = reshape(th_avg_v1,[top,length(time_model)]);

x = zeros(size(hgt_fin));

for i = 1:25
   x(i,:) = time_model';
end

%% Plot time heights.

close all

load('../stats_project/data/double_blend.mat');
load('../stats_project/data/pcaps_timeheight.mat');

master = master';

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

subplot(4,1,1)

contourf(master,height,v,[-12:1:12],'LineStyle','none')
contourcmap('jet',[-12:12])
hold on
contour(master,height,v,[0:12],'k-');
hold on
contour(master,height,v,[-12:-1],'k--');
hold on
plot(master,zeros(length(master))+1800,'k--','LineWidth',2)
set(gca,'xlim',[time_model(1) time_model(end)])
set(gca,'ylim',[min(min(hgt_fin)) min(hgt_fin(25,:))])
set(gca,'XTickLabel','')
ylabel('Height ASL (m)')

fl = fig_label('a','ul',18,0.0025,0.0525,0.205,230);
set(fl,'color','k');

subplot(4,1,2)

fig_nudge(gca,0,0.05)

contourf(x,hgt_fin,v_fin,[-12:1:12],'LineStyle','none')
contourcmap('jet',[-12:12])
hold on
contour(x,hgt_fin,v_fin,[0:12],'k-');
hold on
contour(x,hgt_fin,v_fin,[-12:-1],'k--');
hold on
plot(time_model,zeros(length(time_model))+1800,'k--','LineWidth',2)
set(gca,'xlim',[time_model(1) time_model(end)])
set(gca,'ylim',[min(min(hgt_fin)) min(hgt_fin(25,:))])
set(gca,'XTickLabel','')
ylabel('Height ASL (m)')

fl = fig_label('b','ul',14,0.004,0.08,0.205,230);
set(fl,'color','k');

subplot(4,1,3)

fig_nudge(gca,0,0.1)

contourf(x,hgt_fin,v_fin1,[-12:1:12],'LineStyle','none')
hold on
contour(x,hgt_fin,v_fin1,[0:12],'k-');
hold on
contour(x,hgt_fin,v_fin1,[-12:-1],'k--');
hold on
plot(time_model,zeros(length(time_model))+1800,'k--','LineWidth',2)
set(gca,'xlim',[time_model(1) time_model(end)])
set(gca,'ylim',[min(min(hgt_fin)) min(hgt_fin(25,:))])
set(gca,'XTickLabel','')
ylabel('Height ASL (m)')

fl = fig_label('c','ul',18,0.0025,0.0525,0.205,230);
set(fl,'color','k');

g = subplot(4,1,4)

fig_nudge(g,0,0.245)

p = get(g,'position');
p(4) = p(4)*0.375; % Add 10 percent to height
set(g, 'position', p);

contourf(x,hgt_fin,v_fin1-v_fin,[-12:1:12],'LineStyle','none')
caxis([-12 12])

c = colorbar;
contourcmap('jet',[-12:12])
hL = title(c, '    m s^-^1','FontSize',12)
set(c,'position',[0.915 0.365 0.01 0.55]) 

hold on
contour(x,hgt_fin,v_fin1-v_fin,[0:12],'k-');
hold on
contour(x,hgt_fin,v_fin1-v_fin,[-12:-1],'k--');
set(gca,'xlim',[time_model(1) time_model(end)])
set(gca,'ylim',[min(min(hgt_fin)) 1800])
datetick('x','mm/dd')
xlabel('Date (Ticks at 00 Z)','FontSize',12)

fl = fig_label('d','ul',14,0.004,0.21,0.205,180);
set(fl,'color','k');


%% Calculate spacial averages of variables.

run1_u_lev1 = reshape(run1_u(1:474,:,1,:),[474,498,192]);
run2_u_lev1 = reshape(run2_u(1:474,:,1,:),[474,498,192]);

run1_avg_u = mean(run1_u_lev1,3);
run2_avg_u = mean(run2_u_lev1,3);

u_diff = run1_avg_u - run2_avg_u;

run1_v_lev1 = reshape(run1_v(:,1:498,1,:),[474,498,192]);
run2_v_lev1 = reshape(run2_v(:,1:498,1,:),[474,498,192]);

run1_avg_v = mean(mean(run1_v_lev1,4),3);
run2_avg_v = mean(mean(run2_v_lev1,4),3);

v_diff = run1_avg_v - run2_avg_v;

run1_mag = ((run1_u_lev1.^2) + (run1_v_lev1.^2)) .^ (1/2);
run1_avg_mag = mean(run1_mag,3);

run2_mag = ((run2_u_lev1.^2) + (run2_v_lev1.^2)) .^ (1/2);
run2_avg_mag = mean(run2_mag,3);

mag_diff = run1_avg_mag - run2_avg_mag;

%% Plot vector differences.

vx_min = 212; 
vx_max = 240;
vy_min = 290;
vy_max = 322;

figure(1)
subplot(1,2,1)
ax = usamap([lat(vx_min,vy_min),lat(vx_max,vy_max)],[lon(vx_min,vy_max),lon(vx_max,vy_min)]);
contourfm(lat(vx_min:vx_max,vy_min:vy_max),lon(vx_min:vx_max,vy_min:vy_max),double(run1_mag(vx_min:vx_max,vy_min:vy_max,119)),[0:0.2:5],'LineStyle','none')
contourcmap('jet',[0:0.2:5],'colorbar','on')
caxis([0 5])
hold on
contourm(lat(vx_min:vx_max,vy_min:vy_max),lon(vx_min:vx_max,vy_min:vy_max),terrain(vx_min:vx_max,vy_min:vy_max),'k','LevelStep',150)
hold on
quivermc(lat(vx_min:2:vx_max,vy_min:2:vy_max),lon(vx_min:2:vx_max,vy_min:2:vy_max),run1_u(vx_min:2:vx_max,vy_min:2:vy_max,1,119),run1_v(vx_min:2:vx_max,vy_min:2:vy_max,1,119))
title('BASE 6 January 0000 UTC')

subplot(1,2,2)
ax = usamap([lat(vx_min,vy_min),lat(vx_max,vy_max)],[lon(vx_min,vy_max),lon(vx_max,vy_min)]);
mag1=contourfm(lat(vx_min:vx_max,vy_min:vy_max),lon(vx_min:vx_max,vy_min:vy_max),double(run2_mag(vx_min:vx_max,vy_min:vy_max,119)),[0:0.2:5],'LineStyle','none')
contourcmap('jet',[0:0.2:5],'colorbar','on')
hold on
contourm(lat(vx_min:vx_max,vy_min:vy_max),lon(vx_min:vx_max,vy_min:vy_max),terrain(vx_min:vx_max,vy_min:vy_max),'k','LevelStep',150)
hold on
quivermc(lat(vx_min:2:vx_max,vy_min:2:vy_max),lon(vx_min:2:vx_max,vy_min:2:vy_max),run2_u(vx_min:2:vx_max,vy_min:2:vy_max,1,119),run2_v(vx_min:2:vx_max,vy_min:2:vy_max,1,119))
caxis([0 5])
title('USGS 6 January 0000 UTC')

subplot(1,2,1)
ax = usamap([lat(vx_min,vy_min),lat(vx_max,vy_max)],[lon(vx_min,vy_max),lon(vx_max,vy_min)]);
contourfm(lat(vx_min:vx_max,vy_min:vy_max),lon(vx_min:vx_max,vy_min:vy_max),double(run1_mag(vx_min:vx_max,vy_min:vy_max,125)),[0:0.2:5],'LineStyle','none')
contourcmap('jet',[0:0.2:5],'colorbar','on')
caxis([0 5])
hold on
contourm(lat(vx_min:vx_max,vy_min:vy_max),lon(vx_min:vx_max,vy_min:vy_max),terrain(vx_min:vx_max,vy_min:vy_max),'k','LevelStep',150)
hold on
quivermc(lat(vx_min:2:vx_max,vy_min:2:vy_max),lon(vx_min:2:vx_max,vy_min:2:vy_max),run1_u(vx_min:2:vx_max,vy_min:2:vy_max,1,125),run1_v(vx_min:2:vx_max,vy_min:2:vy_max,1,125))
title('BASE 6 January 0600 UTC')

subplot(1,2,2)
ax = usamap([lat(vx_min,vy_min),lat(vx_max,vy_max)],[lon(vx_min,vy_max),lon(vx_max,vy_min)]);
mag1=contourfm(lat(vx_min:vx_max,vy_min:vy_max),lon(vx_min:vx_max,vy_min:vy_max),double(run2_mag(vx_min:vx_max,vy_min:vy_max,125)),[0:0.2:5],'LineStyle','none')
contourcmap('jet',[0:0.2:5],'colorbar','on')
hold on
contourm(lat(vx_min:vx_max,vy_min:vy_max),lon(vx_min:vx_max,vy_min:vy_max),terrain(vx_min:vx_max,vy_min:vy_max),'k','LevelStep',150)
hold on
quivermc(lat(vx_min:2:vx_max,vy_min:2:vy_max),lon(vx_min:2:vx_max,vy_min:2:vy_max),run2_u(vx_min:2:vx_max,vy_min:2:vy_max,1,125),run2_v(vx_min:2:vx_max,vy_min:2:vy_max,1,125))
caxis([0 5])
title('USGS 6 January 0600 UTC')

figure(3)
subplot(1,2,1)
ax = usamap([lat(vx_min,vy_min),lat(vx_max,vy_max)],[lon(vx_min,vy_max),lon(vx_max,vy_min)]);
contourfm(lat(vx_min:vx_max,vy_min:vy_max),lon(vx_min:vx_max,vy_min:vy_max),double(run1_mag(vx_min:vx_max,vy_min:vy_max,131)),[0:0.2:5],'LineStyle','none')
contourcmap('jet',[0:0.2:5],'colorbar','on')
caxis([0 5])
hold on
contourm(lat(vx_min:vx_max,vy_min:vy_max),lon(vx_min:vx_max,vy_min:vy_max),terrain(vx_min:vx_max,vy_min:vy_max),'k','LevelStep',150)
hold on
quivermc(lat(vx_min:2:vx_max,vy_min:2:vy_max),lon(vx_min:2:vx_max,vy_min:2:vy_max),run1_u(vx_min:2:vx_max,vy_min:2:vy_max,1,131),run1_v(vx_min:2:vx_max,vy_min:2:vy_max,1,131))
title('BASE 6 January 1200 UTC')

subplot(1,2,2)
ax = usamap([lat(vx_min,vy_min),lat(vx_max,vy_max)],[lon(vx_min,vy_max),lon(vx_max,vy_min)]);
mag1=contourfm(lat(vx_min:vx_max,vy_min:vy_max),lon(vx_min:vx_max,vy_min:vy_max),double(run2_mag(vx_min:vx_max,vy_min:vy_max,131)),[0:0.2:5],'LineStyle','none')
contourcmap('jet',[0:0.2:5],'colorbar','on')
hold on
contourm(lat(vx_min:vx_max,vy_min:vy_max),lon(vx_min:vx_max,vy_min:vy_max),terrain(vx_min:vx_max,vy_min:vy_max),'k','LevelStep',150)
hold on
quivermc(lat(vx_min:2:vx_max,vy_min:2:vy_max),lon(vx_min:2:vx_max,vy_min:2:vy_max),run2_u(vx_min:2:vx_max,vy_min:2:vy_max,1,131),run2_v(vx_min:2:vx_max,vy_min:2:vy_max,1,131))
caxis([0 5])
title('USGS 6 January 1200 UTC')

figure(4)
subplot(1,2,1)
ax = usamap([lat(vx_min,vy_min),lat(vx_max,vy_max)],[lon(vx_min,vy_max),lon(vx_max,vy_min)]);
contourfm(lat(vx_min:vx_max,vy_min:vy_max),lon(vx_min:vx_max,vy_min:vy_max),double(run1_mag(vx_min:vx_max,vy_min:vy_max,137)),[0:0.2:5],'LineStyle','none')
contourcmap('jet',[0:0.2:5],'colorbar','on')
caxis([0 5])
hold on
contourm(lat(vx_min:vx_max,vy_min:vy_max),lon(vx_min:vx_max,vy_min:vy_max),terrain(vx_min:vx_max,vy_min:vy_max),'k','LevelStep',150)
hold on
quivermc(lat(vx_min:2:vx_max,vy_min:2:vy_max),lon(vx_min:2:vx_max,vy_min:2:vy_max),run1_u(vx_min:2:vx_max,vy_min:2:vy_max,1,137),run1_v(vx_min:2:vx_max,vy_min:2:vy_max,1,137))
title('BASE 6 January 1800 UTC')

subplot(1,2,2)
ax = usamap([lat(vx_min,vy_min),lat(vx_max,vy_max)],[lon(vx_min,vy_max),lon(vx_max,vy_min)]);
mag1=contourfm(lat(vx_min:vx_max,vy_min:vy_max),lon(vx_min:vx_max,vy_min:vy_max),double(run2_mag(vx_min:vx_max,vy_min:vy_max,137)),[0:0.2:5],'LineStyle','none')
contourcmap('jet',[0:0.2:5],'colorbar','on')
hold on
contourm(lat(vx_min:vx_max,vy_min:vy_max),lon(vx_min:vx_max,vy_min:vy_max),terrain(vx_min:vx_max,vy_min:vy_max),'k','LevelStep',150)
hold on
quivermc(lat(vx_min:2:vx_max,vy_min:2:vy_max),lon(vx_min:2:vx_max,vy_min:2:vy_max),run2_u(vx_min:2:vx_max,vy_min:2:vy_max,1,137),run2_v(vx_min:2:vx_max,vy_min:2:vy_max,1,137))
caxis([0 5])
title('USGS 6 January 1800 UTC')

figure(5)
ax = usamap([lat(vx_min,vy_min),lat(vx_max,vy_max)],[lon(vx_min,vy_max),lon(vx_max,vy_min)]);
contourm(lat(vx_min:vx_max,vy_min:vy_max),lon(vx_min:vx_max,vy_min:vy_max),terrain(vx_min:vx_max,vy_min:vy_max),'k','LevelStep',150)
hold on
quivermc(lat(vx_min:2:vx_max,vy_min:2:vy_max),lon(vx_min:2:vx_max,vy_min:2:vy_max),u_diff(vx_min:2:vx_max,vy_min:2:vy_max),v_diff(vx_min:2:vx_max,vy_min:2:vy_max))
title('BASE - USGS Vector Difference')

lx_min = 155+13; 
lx_max = 235;
ly_min = 285;
ly_max = 400;


figure(6)
ax = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
contourm(lat(lx_min:lx_max,ly_min:ly_max),lon(lx_min:lx_max,ly_min:ly_max),terrain(lx_min:lx_max,ly_min:ly_max),'k','LevelStep',150)
hold on
contourm(lat(lx_min:lx_max,ly_min:ly_max),lon(lx_min:lx_max,ly_min:ly_max),flg(lx_min:lx_max,ly_min:ly_max),[1],'k')
hold on
quivermc(lat(lx_min:5:lx_max,ly_min:5:ly_max),lon(lx_min:5:lx_max,ly_min:5:ly_max),u_diff(lx_min:5:lx_max,ly_min:5:ly_max),v_diff(lx_min:5:lx_max,ly_min:5:ly_max))
title('BASE - USGS Vector Difference')

%% Day and night averages.

night = [ 7:12,31:36,55:60,79:84,103:108,127:132,151:156,175:180];

day   = [19:24,43:48,67:72,91:96,115:120,139:144,163:168,187:192];

% Data

run1_u_day = reshape(run1_u_lev1(:,:,day),[474,498,192/4]);
run2_u_day = reshape(run2_u_lev1(:,:,day),[474,498,192/4]);

run1_u_night = reshape(run1_u_lev1(:,:,night),[474,498,192/4]);
run2_u_night = reshape(run2_u_lev1(:,:,night),[474,498,192/4]);

run1_v_day = reshape(run1_v_lev1(:,:,day),[474,498,192/4]);
run2_v_day = reshape(run2_v_lev1(:,:,day),[474,498,192/4]);

run1_v_night = reshape(run1_v_lev1(:,:,night),[474,498,192/4]);
run2_v_night = reshape(run2_v_lev1(:,:,night),[474,498,192/4]);

% Average

run1_avg_u_day = mean(run1_u_day,3);
run2_avg_u_day = mean(run2_u_day,3);

run1_avg_u_night = mean(run1_u_night,3);
run2_avg_u_night = mean(run2_u_night,3);

run1_avg_v_day = mean(run1_v_day,3);
run2_avg_v_day = mean(run2_v_day,3);

run1_avg_v_night = mean(run1_v_night,3);
run2_avg_v_night = mean(run2_v_night,3);

run1_avg_mag_day = ((run1_avg_u_day.^2) + (run1_avg_v_day.^2)).^ 0.5;
run2_avg_mag_day = ((run2_avg_u_day.^2) + (run2_avg_v_day.^2)) .^ 0.5;

run1_avg_mag_night = ((run1_avg_u_night.^2) + (run1_avg_v_night.^2)) .^ 0.5;
run2_avg_mag_night = ((run2_avg_u_night.^2) + (run2_avg_v_night.^2)) .^ 0.5;

% Difference

u_diff_day   = run1_avg_u_day - run2_avg_u_day;
u_diff_night = run1_avg_u_night - run2_avg_u_night;

v_diff_day   = run1_avg_v_day - run2_avg_v_day;
v_diff_night = run1_avg_v_night - run2_avg_v_night;

mag_diff_day = run1_avg_mag_day - run2_avg_mag_day;
mag_diff_night = run1_avg_mag_night - run2_avg_mag_night;

%%

jet = colormap('jet')
jet_r = jet(64-48:64,:)

figure(7)
subplot(1,2,1)
ax = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
mag1=contourfm(lat(lx_min:lx_max,ly_min:ly_max),lon(lx_min:lx_max,ly_min:ly_max),double(run1_avg_mag_day(lx_min:lx_max,ly_min:ly_max)),[0:0.1:2],'LineStyle','none')
contourcmap('jet',[0:0.1:2],'colorbar','on')
caxis([0 2])
hold on
contourm(lat(lx_min:lx_max,ly_min:ly_max),lon(lx_min:lx_max,ly_min:ly_max),terrain(lx_min:lx_max,ly_min:ly_max),'k','LevelStep',150)
hold on
contourm(lat(lx_min:lx_max,ly_min:ly_max),lon(lx_min:lx_max,ly_min:ly_max),flg(lx_min:lx_max,ly_min:ly_max),[1],'k')
hold on
quivermc(lat(lx_min:5:lx_max,ly_min:5:ly_max),lon(lx_min:5:lx_max,ly_min:5:ly_max),run1_avg_u_day(lx_min:5:lx_max,ly_min:5:ly_max),run1_avg_v_day(lx_min:5:lx_max,ly_min:5:ly_max))
title('BASE Day')

subplot(1,2,2)
ax = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
mag1=contourfm(lat(lx_min:lx_max,ly_min:ly_max),lon(lx_min:lx_max,ly_min:ly_max),double(run2_avg_mag_day(lx_min:lx_max,ly_min:ly_max)),[0:0.1:2],'LineStyle','none')
contourcmap('jet',[0:0.1:2],'colorbar','on')
caxis([0 2])
hold on
contourm(lat(lx_min:lx_max,ly_min:ly_max),lon(lx_min:lx_max,ly_min:ly_max),terrain(lx_min:lx_max,ly_min:ly_max),'k','LevelStep',150)
hold on
contourm(lat(lx_min:lx_max,ly_min:ly_max),lon(lx_min:lx_max,ly_min:ly_max),flg(lx_min:lx_max,ly_min:ly_max),[1],'k')
hold on
quivermc(lat(lx_min:5:lx_max,ly_min:5:ly_max),lon(lx_min:5:lx_max,ly_min:5:ly_max),run2_avg_u_day(lx_min:5:lx_max,ly_min:5:ly_max),run2_avg_v_day(lx_min:5:lx_max,ly_min:5:ly_max))
title('USGS Day')

%%

cmap = colormap(jet)
cmap_short = colormap(cmap(34:60,:))

figure('Position', [100, 100, 600, 450])

subplot(1,2,1)
ax = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
contourfm(lat,lon,double(run1_avg_mag_night),[0:0.2:4],'LineStyle','none')
colormap(cmap_short)
caxis([0 4])
hold on
contourm(lat,lon,terrain,'k','LevelStep',150)
hold on
contourm(lat,lon,flg,[1],'k')
hold on
quivermc(lat(lx_min-20:5:lx_max+20,ly_min-20:5:ly_max+20),lon(lx_min-20:5:lx_max+20,ly_min-20:5:ly_max+20),run1_avg_u_night(lx_min-20:5:lx_max+20,ly_min-20:5:ly_max+20),run1_avg_v_night(lx_min-20:5:lx_max+20,ly_min-20:5:ly_max+20),'scale',3,'Color','k');
title('BASE')

fl = figm_label('a','ul',14,0.035,0.04,13500,13500);
set(fl,'color','k');

subplot(1,2,2)
ax2 = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
plabel('off')
fig_nudge(ax2,-0.08,0)
contourfm(lat,lon,double(run2_avg_mag_night),[0:0.2:4],'LineStyle','none')
colormap(cmap_short)
caxis([0 4])
hold on
contourm(lat,lon,terrain,'k','LevelStep',150)
hold on
contourm(lat,lon,flg,[1],'k')
hold on
quivermc(lat(lx_min-20:5:lx_max+20,ly_min-20:5:ly_max+20),lon(lx_min-20:5:lx_max+20,ly_min-20:5:ly_max+20),run2_avg_u_night(lx_min-20:5:lx_max+20,ly_min-20:5:ly_max+20),run2_avg_v_night(lx_min-20:5:lx_max+20,ly_min-20:5:ly_max+20),'scale',3,'Color','k');
title('DFLT')
c = colorbar;
%contourcmap('winter',[-4:0.25:4])
hL = title(c, '    m s^-^1','FontSize',12)
set(c,'position',[0.84 0.245 0.012 0.55]) 

fl = figm_label('b','ul',14,0.03,0.045,13500,13500);
set(fl,'color','k');
%% Plot comparison.

close all

states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);

figure(1)
title('Average Wind Speed Difference 01/01/1-01/10/11, NLCD 31, Dec 2010 - USGS 31, Dec 2010')
ax1 = usamap(latlim,lonlim);
temp_ax = surfm(lat,lon,mag_diff);
geoshow(ax1, states, 'FaceColor', 'none');
colorbar('YTick',[-1:0.25:1])
contourcmap('jet',[-1:0.05:1])
caxis([-1 1])
contourm(lat,lon,terrain,'k','LevelStep',250)
contourm(lat,lon,flg,'k','LevelStep',1,'LineWidth',2)

figure(2)
title('Average Daytime Wind Speed 01/01/1-01/09/11, NLCD 31, Dec 2010 - USGS 31, Dec 2010')
ax1 = usamap(latlim,lonlim);
temp_ax = surfm(lat,lon,mag_diff_day);
geoshow(ax1, states, 'FaceColor', 'none');
colorbar('YTick',[-1:0.25:1])
contourcmap('jet',[-1:0.05:1])
caxis([-1 1])
contourm(lat,lon,terrain,'k','LevelStep',250)
contourm(lat,lon,flg,'k','LevelStep',1,'LineWidth',2)

figure(3)
title('Average Nighttime Wind Speed 01/01/1-01/09/11, NLCD 31, Dec 2010 - USGS 31, Dec 2010')
ax1 = usamap(latlim,lonlim);
temp_ax = surfm(lat,lon,mag_diff_night);
geoshow(ax1, states, 'FaceColor', 'none');
colorbar('YTick',[-1:0.25:1],'FontSize',12)
contourcmap('jet',[-1:0.05:1])
caxis([-1 1])
contourm(lat,lon,terrain,'k','LevelStep',250)
contourm(lat,lon,flg,'k','LevelStep',1,'LineWidth',2)

%%

lx_min = 155; 
lx_max = 235;
ly_min = 285;
ly_max = 400;

close all

title4 = {'Average Daytime Wind Speed 01/01/11-01/09/11,','NLCD 31, Dec 2010 - USGS 31, Dec 2010'};

figure(4)

subplot(1,2,1)
title(title4,'FontSize',10)
ax2 = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
surfm(lat,lon,mag_diff_day);
colorbar('YTick',[-1:0.25:1])
contourcmap('jet',[-1:0.05:1])
caxis([-1 1])
contourm(lat,lon,flg,'k','LevelStep',1,'LineWidth',2)
contourm(lat,lon,terrain,'k','LevelStep',125,'LineWidth',0.5)

title5 = {'Average Nighttime Wind Speed 01/01/11-01/09/11,','NLCD 31, Dec 2010 - USGS 31, Dec 2010'}

subplot(1,2,2)
title(title5,'FontSize',10)
ax2 = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
surfm(lat,lon,mag_diff_night);
colorbar('YTick',[-1:0.25:1],'FontSize',12)
contourcmap('jet',[-1:0.05:1])
caxis([-1 1])
contourm(lat,lon,flg,'k','LevelStep',1,'LineWidth',2)
contourm(lat,lon,terrain,'k','LevelStep',125,'LineWidth',0.5)

figure(5)

subplot(1,2,1)
title(title4,'FontSize',10)
ax2 = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
contourfm(double(lat),double(lon),double(mag_diff_day),[-1:0.2:1],'LineStyle','none');
colorbar('YTick',[-1:0.25:1])
contourcmap('jet',[-1:0.05:1])
caxis([-1 1])
contourm(lat,lon,flg,'k','LevelStep',1,'LineWidth',2)
contourm(lat,lon,terrain,'k','LevelStep',125,'LineWidth',0.5)

title5 = {'Average Nighttime Wind Speed 01/01/11-01/09/11,','NLCD 31, Dec 2010 - USGS 31, Dec 2010'}

subplot(1,2,2)
title(title5,'FontSize',10)
ax2 = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
contourfm(double(lat),double(lon),double(mag_diff_night),[-1:0.2:1],'LineStyle','none');
colorbar('YTick',[-1:0.25:1],'FontSize',12)
contourcmap('jet',[-1:0.05:1])
caxis([-1 1])
contourm(lat,lon,flg,'k','LevelStep',1,'LineWidth',2)
contourm(lat,lon,terrain,'k','LevelStep',125,'LineWidth',0.5)

%%

close all

title4 = {'NLCD Average Daytime u Wind 01/01/11-01/09/11'}%,'NLCD 31, Dec 2010 - USGS 31, Dec 2010'};

figure(6)

subplot(1,2,1)
title(title4,'FontSize',10)
ax2 = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
contourfm(double(lat),double(lon),double(run2_avg_v_day),[-4:0.25:4],'LineStyle','none');
colorbar('YTick',[-4:1:4])
contourcmap('jet',[-4:0.25:4])
caxis([-4 4])
contourm(lat,lon,flg,'k','LevelStep',1,'LineWidth',2)
contourm(lat,lon,terrain,'k','LevelStep',125,'LineWidth',0.5)

title5 = {'NLCD Average Nighttime u Wind 01/01/11-01/09/11'}%,'NLCD 31, Dec 2010 - USGS 31, Dec 2010'}

subplot(1,2,2)
title(title5,'FontSize',10)
ax2 = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
contourfm(double(lat),double(lon),double(run2_avg_v_night),[-4:0.25:45],'LineStyle','none');
colorbar('YTick',[-4:1:4],'FontSize',12)
contourcmap('jet',[-4:0.25:4])
caxis([-4 4])
contourm(lat,lon,flg,'k','LevelStep',1,'LineWidth',2)
contourm(lat,lon,terrain,'k','LevelStep',125,'LineWidth',0.5)

%%

close all

run = 'NLCD2011'

for i = 1:length(time_model)
   
  

   figure(7)
   
   title4 = {'NLCD u Wind',datestr(time_model(i))};

   subplot(1,2,1)
   title(title4,'FontSize',10)
   ax2 = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
   surfm(lat,lon,run1_u(:,:,1,i));
   colorbar('YTick',[-5:0.5:5])
   contourcmap('jet',[-5:0.25:5])
   caxis([-5 5])
   contourm(lat,lon,flg,'k','LevelStep',1,'LineWidth',2)
   contourm(lat,lon,terrain,'k','LevelStep',125,'LineWidth',0.5)

   title4 = {'NLCD v Wind',datestr(time_model(i))};

   subplot(1,2,2)
   title(title4,'FontSize',10)
   ax2 = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
   surfm(lat,lon,run1_v(:,:,1,i));
   colorbar('YTick',[-5:0.5:5])
   contourcmap('jet',[-5:0.25:5])
   caxis([-5 5])
   contourm(lat,lon,flg,'k','LevelStep',1,'LineWidth',2)
   contourm(lat,lon,terrain,'k','LevelStep',125,'LineWidth',0.5)
   
   filename = strcat('/uufs/chpc.utah.edu/common/home/u0875930/public_html/wrf_output/iop5/wind/hourly/',run,'_',datestr(time_model(i)))
   saveas(7,filename,'png')
   
%    drawnow
%    frame = getframe(1);
%    im = frame2im(frame);
%    [imind,cm] = rgb2ind(im,256);
% 
%       if i == 1;
%          imwrite(imind,cm,filename,'gif','Loopcount',inf);
%       else
%          imwrite(imind,cm,filename,'gif','WriteMode','append')
%       end
%    
end

