

% 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_grdflx = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run1,'/wrfout_d03_2010-12-31_00:00:00'),'GRDFLX');
run2_grdflx = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run2,'/wrfout_d03_2011-01-01_00:00:00'),'GRDFLX');

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

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

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

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

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

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.
      
      grdflx_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'),'GRDFLX');
      
      run1_grdflx = cat(3,run1_grdflx,grdflx_temp);      
      
      grdflx_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'),'GRDFLX');
      
      run2_grdflx = cat(3,run2_grdflx,grdflx_temp);
      
      hfx_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'),'HFX');
      
      run1_hfx = cat(3,run1_hfx,hfx_temp);      
      
      hfx_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'),'HFX');
      
      run2_hfx = cat(3,run2_hfx,hfx_temp);
      
      glw_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'),'GLW');
      
      run1_glw = cat(3,run1_glw,glw_temp);      
      
      glw_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'),'GLW');
      
      run2_glw = cat(3,run2_glw,glw_temp);
      
      swd_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'),'SWDOWN');
      
      run1_swd = cat(3,run1_swd,swd_temp);      
      
      swd_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'),'SWDOWN');
      
      run2_swd = cat(3,run2_swd,swd_temp);
      
      alb_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'),'ALBEDO');
      
      run1_alb = cat(3,run1_alb,alb_temp);      
      
      alb_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'),'ALBEDO');
      
      run2_alb = cat(3,run2_alb,alb_temp);
      
      lhf_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'),'LH');
      
      run1_lhf = cat(3,run1_lhf,lhf_temp);      
      
      lhf_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'),'LH');
      
      run2_lhf = cat(3,run2_lhf,lhf_temp);
      
      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_grdflx(:,:,1:6) = []; % Remove first time (repeat).
run2_grdflx(:,:,1:6) = [];
run1_hfx(:,:,1:6) = [];
run2_hfx(:,:,1:6) = [];
run1_glw(:,:,1:6) = [];
run2_glw(:,:,1:6) = [];
run1_swd(:,:,1:6) = [];
run2_swd(:,:,1:6) = [];
run1_alb(:,:,1:6) = [];
run2_alb(:,:,1:6) = [];
run1_lhf(:,:,1:6) = [];
run2_lhf(:,:,1:6) = [];
time_model(1) = [];   % Remove first time (repeat).


%% Calculate spacial averages of variables.

run1_grdflx_lev1 = reshape(run1_grdflx(1:474,:,:),[474,498,192]);
run2_grdflx_lev1 = reshape(run2_grdflx(1:474,:,:),[474,498,192]);

run1_avg_grdflx = mean(run1_grdflx_lev1,3);
run2_avg_grdflx = mean(run2_grdflx_lev1,3);

grdflx_diff = run1_avg_grdflx - run2_avg_grdflx;

run1_hfx_lev1 = reshape(run1_hfx(:,1:498,:),[474,498,192]);
run2_hfx_lev1 = reshape(run2_hfx(:,1:498,:),[474,498,192]);

run1_avg_hfx = mean(mean(run1_hfx_lev1,4),3);
run2_avg_hfx = mean(mean(run2_hfx_lev1,4),3);

hfx_diff = run1_avg_hfx - run2_avg_hfx;

%% 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_grdflx_day = reshape(run1_grdflx_lev1(:,:,day),[474,498,192/4]);
run2_grdflx_day = reshape(run2_grdflx_lev1(:,:,day),[474,498,192/4]);

run1_grdflx_night = reshape(run1_grdflx_lev1(:,:,night),[474,498,192/4]);
run2_grdflx_night = reshape(run2_grdflx_lev1(:,:,night),[474,498,192/4]);

run1_hfx_day = reshape(run1_hfx_lev1(:,:,day),[474,498,192/4]);
run2_hfx_day = reshape(run2_hfx_lev1(:,:,day),[474,498,192/4]);

run1_hfx_night = reshape(run1_hfx_lev1(:,:,night),[474,498,192/4]);
run2_hfx_night = reshape(run2_hfx_lev1(:,:,night),[474,498,192/4]);

% Average

run1_avg_grdflx_day = mean(run1_grdflx_day,3);
run2_avg_grdflx_day = mean(run2_grdflx_day,3);

run1_avg_grdflx_night = mean(run1_grdflx_night,3);
run2_avg_grdflx_night = mean(run2_grdflx_night,3);

run1_avg_hfx_day = mean(run1_hfx_day,3);
run2_avg_hfx_day = mean(run2_hfx_day,3);

run1_avg_hfx_night = mean(run1_hfx_night,3);
run2_avg_hfx_night = mean(run2_hfx_night,3);

% Difference

grdflx_diff_day   = run1_avg_grdflx_day - run2_avg_grdflx_day;
grdflx_diff_night = run1_avg_grdflx_night - run2_avg_grdflx_night;

hfx_diff_day   = run1_avg_hfx_day - run2_avg_hfx_day;
hfx_diff_night = run1_avg_hfx_night - run2_avg_hfx_night;

%% Plot comparison.

close all

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

figure(1)
title('Average Differenced Ground Heat Flux 01/01/1-01/08/11, NLCD 31, Dec 2010 - USGS 31, Dec 2010')
ax1 = usamap(latlim,lonlim);
temp_ax = surfm(lat,lon,grdflx_diff);
geoshow(ax1, states, 'FaceColor', 'none');
colorbar('YTick',[-20:2:20])
contourcmap('jet',[-20:1:20])
caxis([-20 20])
contourm(lat,lon,terrain,'k','LevelStep',250)
contourm(lat,lon,flg,'k','LevelStep',1,'LineWidth',2)

figure(2)
title('Average Differenced Upward Sensible Heat Flux 01/01/1-01/08/11, NLCD 31, Dec 2010 - USGS 31, Dec 2010')
ax1 = usamap(latlim,lonlim);
temp_ax = surfm(lat,lon,hfx_diff_day);
geoshow(ax1, states, 'FaceColor', 'none');
colorbar('YTick',[-20:2:20])
contourcmap('jet',[-20:1:20])
caxis([-20 20])
contourm(lat,lon,terrain,'k','LevelStep',250)
contourm(lat,lon,flg,'k','LevelStep',1,'LineWidth',2)

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

title3 = {'Average Ground Heat Flux 01/01/11-01/08/11,','NLCD 31, Dec 2010 - USGS 31, Dec 2010'};

figure(3)

title(title3,'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,grdflx_diff);
geoshow(ax2,states,'FaceColor','none');
colorbar('YTick',[-20:2:20])
contourcmap('jet',[-20:1:20])
caxis([-20 20])
contourm(lat,lon,terrain,'k','LevelStep',250)
contourm(lat,lon,flg,'k','LevelStep',1,'LineWidth',2)

%%

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

subplot(1,2,1)
title('Day','FontSize',12)
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,hfx_diff_day);
contourcmap('jet',[-20:1:20])
contourm(lat,lon,terrain,'k','LevelStep',250)
contourm(lat,lon,flg,'k','LevelStep',1,'LineWidth',2)
textm(41.4,-112,'4','Color','white','FontSize',18,'FontWeight','bold')

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

subplot(1,2,2)
fig_nudge(gca,-0.08,0)
title('Night','FontSize',12)
ax3 = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
surfm(lat,lon,hfx_diff_night);
caxis([-20 20])
contourcmap('jet',[-20:1:20])
contourm(lat,lon,terrain,'k','LevelStep',250)
contourm(lat,lon,flg,'k','LevelStep',1,'LineWidth',2)
textm(41.575,-112.9,'3','Color','white','FontSize',18,'FontWeight','bold')
textm(41.1,-112.4,'3','Color','white','FontSize',18,'FontWeight','bold')
textm(40.75,-112.4,'3','Color','white','FontSize',18,'FontWeight','bold')
textm(40.94,-112.125,'1','Color','white','FontSize',18,'FontWeight','bold')
textm(40.65,-111.95,'2','Color','black','FontSize',18,'FontWeight','bold')
plabel('off')

c = colorbar;
hL = title(c,'W m^-^2','FontSize',10);

cpos = c.Position;
cpos(3) = 0.5*cpos(3);
c.Position = cpos;
set(c,'position',[0.84 0.245 0.012 0.55]) 

fl = figm_label('b','ul',14,0.04,0.05,13500,13500);
set(fl,'color','k');


%% Time Series

load('station_data_flux.mat');

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

time_stn = data(:,1);
sw_stn = data(:,2);
lw_stn = data(:,3);

long = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run1,'/wrfout_d03_2011-01-01_00:00:00'),'XLONG');
lati = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run1,'/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];
lat = [40.7785,   40.7374,   40.6980,    40.6810,   40.6127,   40.5847,  40.5140,   40.61,    40.6006];

% 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)))];

run1_grdflx_model = [reshape([run1_grdflx(x_grid(1),y_grid(1),:)],[length(run1_grdflx(1,1,:)),1]), ...
    reshape([run1_grdflx(x_grid(2),y_grid(2),:)],[length(run1_grdflx(1,1,:)),1]), ...
    reshape([run1_grdflx(x_grid(3),y_grid(3),:)],[length(run1_grdflx(1,1,:)),1]), ...
    reshape([run1_grdflx(x_grid(4),y_grid(4),:)],[length(run1_grdflx(1,1,:)),1]), ...
    reshape([run1_grdflx(x_grid(5),y_grid(5),:)],[length(run1_grdflx(1,1,:)),1]), ...
    reshape([run1_grdflx(x_grid(6),y_grid(6),:)],[length(run1_grdflx(1,1,:)),1]), ...
    reshape([run1_grdflx(x_grid(7),y_grid(7),:)],[length(run1_grdflx(1,1,:)),1])];

run1_hfx_model = [reshape([run1_hfx(x_grid(1),y_grid(1),:)],[length(run1_hfx(1,1,:)),1]), ...
    reshape([run1_hfx(x_grid(2),y_grid(2),:)],[length(run1_hfx(1,1,:)),1]), ...
    reshape([run1_hfx(x_grid(3),y_grid(3),:)],[length(run1_hfx(1,1,:)),1]), ...
    reshape([run1_hfx(x_grid(4),y_grid(4),:)],[length(run1_hfx(1,1,:)),1]), ...
    reshape([run1_hfx(x_grid(5),y_grid(5),:)],[length(run1_hfx(1,1,:)),1]), ...
    reshape([run1_hfx(x_grid(6),y_grid(6),:)],[length(run1_hfx(1,1,:)),1]), ...
    reshape([run1_hfx(x_grid(7),y_grid(7),:)],[length(run1_hfx(1,1,:)),1])];

run1_glw_model = [reshape([run1_glw(x_grid(1),y_grid(1),:)],[length(run1_glw(1,1,:)),1]), ...
    reshape([run1_glw(x_grid(2),y_grid(2),:)],[length(run1_glw(1,1,:)),1]), ...
    reshape([run1_glw(x_grid(3),y_grid(3),:)],[length(run1_glw(1,1,:)),1]), ...
    reshape([run1_glw(x_grid(4),y_grid(4),:)],[length(run1_glw(1,1,:)),1]), ...
    reshape([run1_glw(x_grid(5),y_grid(5),:)],[length(run1_glw(1,1,:)),1]), ...
    reshape([run1_glw(x_grid(6),y_grid(6),:)],[length(run1_glw(1,1,:)),1]), ...
    reshape([run1_glw(x_grid(7),y_grid(7),:)],[length(run1_glw(1,1,:)),1])];

run1_swd_model = [reshape([run1_swd(x_grid(1),y_grid(1),:)],[length(run1_swd(1,1,:)),1]), ...
    reshape([run1_swd(x_grid(2),y_grid(2),:)],[length(run1_swd(1,1,:)),1]), ...
    reshape([run1_swd(x_grid(3),y_grid(3),:)],[length(run1_swd(1,1,:)),1]), ...
    reshape([run1_swd(x_grid(4),y_grid(4),:)],[length(run1_swd(1,1,:)),1]), ...
    reshape([run1_swd(x_grid(5),y_grid(5),:)],[length(run1_swd(1,1,:)),1]), ...
    reshape([run1_swd(x_grid(6),y_grid(6),:)],[length(run1_swd(1,1,:)),1]), ...
    reshape([run1_swd(x_grid(7),y_grid(7),:)],[length(run1_swd(1,1,:)),1])];

run1_alb_model = [reshape([run1_alb(x_grid(1),y_grid(1),:)],[length(run1_alb(1,1,:)),1]), ...
    reshape([run1_alb(x_grid(2),y_grid(2),:)],[length(run1_alb(1,1,:)),1]), ...
    reshape([run1_alb(x_grid(3),y_grid(3),:)],[length(run1_alb(1,1,:)),1]), ...
    reshape([run1_alb(x_grid(4),y_grid(4),:)],[length(run1_alb(1,1,:)),1]), ...
    reshape([run1_alb(x_grid(5),y_grid(5),:)],[length(run1_alb(1,1,:)),1]), ...
    reshape([run1_alb(x_grid(6),y_grid(6),:)],[length(run1_alb(1,1,:)),1]), ...
    reshape([run1_alb(x_grid(7),y_grid(7),:)],[length(run1_alb(1,1,:)),1])];

run1_lhf_model = [reshape([run1_lhf(x_grid(1),y_grid(1),:)],[length(run1_lhf(1,1,:)),1]), ...
    reshape([run1_lhf(x_grid(2),y_grid(2),:)],[length(run1_lhf(1,1,:)),1]), ...
    reshape([run1_lhf(x_grid(3),y_grid(3),:)],[length(run1_lhf(1,1,:)),1]), ...
    reshape([run1_lhf(x_grid(4),y_grid(4),:)],[length(run1_lhf(1,1,:)),1]), ...
    reshape([run1_lhf(x_grid(5),y_grid(5),:)],[length(run1_lhf(1,1,:)),1]), ...
    reshape([run1_lhf(x_grid(6),y_grid(6),:)],[length(run1_lhf(1,1,:)),1]), ...
    reshape([run1_lhf(x_grid(7),y_grid(7),:)],[length(run1_lhf(1,1,:)),1])];

run2_grdflx_model = [reshape([run2_grdflx(x_grid(1),y_grid(1),:)],[length(run2_grdflx(1,1,:)),1]), ...
    reshape([run2_grdflx(x_grid(2),y_grid(2),:)],[length(run2_grdflx(1,1,:)),1]), ...
    reshape([run2_grdflx(x_grid(3),y_grid(3),:)],[length(run2_grdflx(1,1,:)),1]), ...
    reshape([run2_grdflx(x_grid(4),y_grid(4),:)],[length(run2_grdflx(1,1,:)),1]), ...
    reshape([run2_grdflx(x_grid(5),y_grid(5),:)],[length(run2_grdflx(1,1,:)),1]), ...
    reshape([run2_grdflx(x_grid(6),y_grid(6),:)],[length(run2_grdflx(1,1,:)),1]), ...
    reshape([run2_grdflx(x_grid(7),y_grid(7),:)],[length(run2_grdflx(1,1,:)),1])];

run2_hfx_model = [reshape([run2_hfx(x_grid(1),y_grid(1),:)],[length(run2_hfx(1,1,:)),1]), ...
    reshape([run2_hfx(x_grid(2),y_grid(2),:)],[length(run2_hfx(1,1,:)),1]), ...
    reshape([run2_hfx(x_grid(3),y_grid(3),:)],[length(run2_hfx(1,1,:)),1]), ...
    reshape([run2_hfx(x_grid(4),y_grid(4),:)],[length(run2_hfx(1,1,:)),1]), ...
    reshape([run2_hfx(x_grid(5),y_grid(5),:)],[length(run2_hfx(1,1,:)),1]), ...
    reshape([run2_hfx(x_grid(6),y_grid(6),:)],[length(run2_hfx(1,1,:)),1]), ...
    reshape([run2_hfx(x_grid(7),y_grid(7),:)],[length(run2_hfx(1,1,:)),1])];

run2_glw_model = [reshape([run2_glw(x_grid(1),y_grid(1),:)],[length(run2_glw(1,1,:)),1]), ...
    reshape([run2_glw(x_grid(2),y_grid(2),:)],[length(run2_glw(1,1,:)),1]), ...
    reshape([run2_glw(x_grid(3),y_grid(3),:)],[length(run2_glw(1,1,:)),1]), ...
    reshape([run2_glw(x_grid(4),y_grid(4),:)],[length(run2_glw(1,1,:)),1]), ...
    reshape([run2_glw(x_grid(5),y_grid(5),:)],[length(run2_glw(1,1,:)),1]), ...
    reshape([run2_glw(x_grid(6),y_grid(6),:)],[length(run2_glw(1,1,:)),1]), ...
    reshape([run2_glw(x_grid(7),y_grid(7),:)],[length(run2_glw(1,1,:)),1])];

run2_swd_model = [reshape([run2_swd(x_grid(1),y_grid(1),:)],[length(run2_swd(1,1,:)),1]), ...
    reshape([run2_swd(x_grid(2),y_grid(2),:)],[length(run2_swd(1,1,:)),1]), ...
    reshape([run2_swd(x_grid(3),y_grid(3),:)],[length(run2_swd(1,1,:)),1]), ...
    reshape([run2_swd(x_grid(4),y_grid(4),:)],[length(run2_swd(1,1,:)),1]), ...
    reshape([run2_swd(x_grid(5),y_grid(5),:)],[length(run2_swd(1,1,:)),1]), ...
    reshape([run2_swd(x_grid(6),y_grid(6),:)],[length(run2_swd(1,1,:)),1]), ...
    reshape([run2_swd(x_grid(7),y_grid(7),:)],[length(run2_swd(1,1,:)),1])];

run2_alb_model = [reshape([run2_alb(x_grid(1),y_grid(1),:)],[length(run2_alb(1,1,:)),1]), ...
    reshape([run2_alb(x_grid(2),y_grid(2),:)],[length(run2_alb(1,1,:)),1]), ...
    reshape([run2_alb(x_grid(3),y_grid(3),:)],[length(run2_alb(1,1,:)),1]), ...
    reshape([run2_alb(x_grid(4),y_grid(4),:)],[length(run2_alb(1,1,:)),1]), ...
    reshape([run2_alb(x_grid(5),y_grid(5),:)],[length(run2_alb(1,1,:)),1]), ...
    reshape([run2_alb(x_grid(6),y_grid(6),:)],[length(run2_alb(1,1,:)),1]), ...
    reshape([run2_alb(x_grid(7),y_grid(7),:)],[length(run2_alb(1,1,:)),1])];

run2_lhf_model = [reshape([run2_lhf(x_grid(1),y_grid(1),:)],[length(run2_lhf(1,1,:)),1]), ...
    reshape([run2_lhf(x_grid(2),y_grid(2),:)],[length(run2_lhf(1,1,:)),1]), ...
    reshape([run2_lhf(x_grid(3),y_grid(3),:)],[length(run2_lhf(1,1,:)),1]), ...
    reshape([run2_lhf(x_grid(4),y_grid(4),:)],[length(run2_lhf(1,1,:)),1]), ...
    reshape([run2_lhf(x_grid(5),y_grid(5),:)],[length(run2_lhf(1,1,:)),1]), ...
    reshape([run2_lhf(x_grid(6),y_grid(6),:)],[length(run2_lhf(1,1,:)),1]), ...
    reshape([run2_lhf(x_grid(7),y_grid(7),:)],[length(run2_lhf(1,1,:)),1])];

%%

run1_swa_model = run1_swd_model .* (1 - run1_alb_model); % absorbed short wave

run1_lwu_model = run1_swa_model + run1_glw_model - run1_hfx_model - run1_lhf_model + run1_grdflx_model;

run1_netlw_model = run1_glw_model - run1_lwu_model;


run2_swa_model = run2_swd_model .* (1 - run2_alb_model); % absorbed short wave

run2_lwu_model = run2_swa_model + run2_glw_model - run2_hfx_model - run2_lhf_model + run2_grdflx_model;

run2_netlw_model = run2_glw_model - run2_lwu_model;

%%

time_stn1 = time_stn{1};
sw_stn1 = sw_stn{1};
lw_stn1 = lw_stn{1};

time_temp = time_stn1;
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_stn2 = time_stn{2};
sw_stn2 = sw_stn{2};
lw_stn2 = lw_stn{2};

time_temp = time_stn2;
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_stn3 = time_stn{3};
sw_stn3 = sw_stn{3};
lw_stn3 = lw_stn{3};

time_temp = time_stn3;
for t = 1:length(time_temp)
   time_isfs3(t,1:6) = datevec(strcat(time_temp{t}(1:10),{' '},time_temp{t}(12:19)));
end

time_stn4 = time_stn{4};
sw_stn4 = sw_stn{4};
lw_stn4 = lw_stn{4};

time_temp = time_stn4;
for t = 1:length(time_temp)
   time_isfs4(t,1:6) = datevec(strcat(time_temp{t}(1:10),{' '},time_temp{t}(12:19)));
end

time_stn5 = time_stn{5};
sw_stn5 = sw_stn{5};
lw_stn5 = lw_stn{5};

time_temp = time_stn5;
for t = 1:length(time_temp)
   time_isfs5(t,1:6) = datevec(strcat(time_temp{t}(1:10),{' '},time_temp{t}(12:19)));
end

time_stn6 = time_stn{6};
sw_stn6 = sw_stn{6};
lw_stn6 = lw_stn{6};

time_temp = time_stn6;
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_stn7 = time_stn{7};
sw_stn7 = sw_stn{7};
lw_stn7 = lw_stn{7};

time_temp = time_stn7;
for t = 1:length(time_temp)
   time_isfs7(t,1:6) = datevec(strcat(time_temp{t}(1:10),{' '},time_temp{t}(12:19)));
end

%%

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);

plot(datenum(time_isfs1),lw_stn1,'k-', ...
    time_model,run1_netlw_model(:,1),'b-', ...
    time_model,run2_netlw_model(:,1), 'r-', ...
    'LineWidth',2)
hold on
plot(datenum(time_isfs1),sw_stn1, 'k-', ...
    time_model,run1_swa_model(:,1),'b-', ...
    time_model,run2_swa_model(:,1), 'r-', ...
    'LineWidth',1)
set(gca, 'XTickLabel', '')
axis([min(time_model) max(time_model) -160 260])
ylabel('W m^-^2')

fig_label('(a)','ul',12,0.005,0.05)

s2 = subplot(3,1,2)

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

plot(datenum(time_isfs2),lw_stn2,'k-', ...
    time_model,run1_netlw_model(:,2),'b-', ...
    time_model,run2_netlw_model(:,2), 'r-', ...
    'LineWidth',2)
hold on
plot(datenum(time_isfs2),sw_stn2, 'k-', ...
    time_model,run1_swa_model(:,2),'b-', ...
    time_model,run2_swa_model(:,2), 'r', ...
    'LineWidth',1)
set(gca, 'XTickLabel', '')
axis([min(time_model) max(time_model) -160 260])
ylabel('W m^-^2')

fig_label('(b)','ul',12,0.005,0.05)

s3 = subplot(3,1,3)

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

plot(datenum(time_isfs6),lw_stn6,'k-', ...
    time_model,run1_netlw_model(:,6),'b-', ...
    time_model,run2_netlw_model(:,6), 'r-', ...
    'LineWidth',2)
hold on
plot(datenum(time_isfs6),sw_stn6, 'k-', ...
    time_model,run1_swa_model(:,6),'b-', ...
    time_model,run2_swa_model(:,6), 'r-', ...
    'LineWidth',1)
datetick('x','mm/dd','keepticks')
axis([min(time_model) max(time_model) -160 260])

legend('Station Data','BASE','DFLT','location','southwest')
xlabel('Date (Ticks at 00 UTC)')
ylabel('W m^-^2')

fig_label('(c)','ul',12,0.005,0.05)

%%
figure(1)

subplot(2,3,1)
plot(time_model,run1_netlw_model(:,1),'b-', ...
    time_model,run2_netlw_model(:,1), 'r-', ...
    datenum(time_isfs1),lw_stn1, 'k-')
title('ISFS 1')
datetick('x','mm/dd','keepticks')
axis([min(time_model) max(time_model) -160 20])

subplot(2,3,2)
plot(time_model,run1_netlw_model(:,2),'b-', ...
    time_model,run2_netlw_model(:,2),'r-', ...
    datenum(time_isfs2),lw_stn2,'k-')
title('ISFS 2')
datetick('x','mm/dd','keepticks')
axis([min(time_model) max(time_model) -160 20])

subplot(2,3,3)
plot(time_model,run1_netlw_model(:,3),'b-', ...
    time_model,run2_netlw_model(:,3),'r-', ...
    datenum(time_isfs3),lw_stn3,'k-')
title('ISFS 3')
datetick('x','mm/dd','keepticks')
axis([min(time_model) max(time_model) -160 20])

subplot(2,3,4)
plot(time_model,run1_netlw_model(:,4),'b-', ...
    time_model,run2_netlw_model(:,4),'r-', ...
    datenum(time_isfs4),lw_stn4,'k-')
title('ISFS 4')
datetick('x','mm/dd','keepticks')
axis([min(time_model) max(time_model) -160 20])

subplot(2,3,5)
plot(time_model,run1_netlw_model(:,5),'b-', ...
    time_model,run2_netlw_model(:,5),'r-', ...
    datenum(time_isfs5),lw_stn5,'k-')
title('ISFS 5')
datetick('x','mm/dd','keepticks')
axis([min(time_model) max(time_model) -160 20])

subplot(2,3,6)
plot(time_model,run1_netlw_model(:,6),'b-', ...
    time_model,run2_netlw_model(:,6),'r-', ...
    datenum(time_isfs6),lw_stn6,'k-')
title('ISFS 6')
datetick('x','mm/dd','keepticks')
axis([min(time_model) max(time_model) -160 20])


figure(2)

subplot(2,3,1)
plot(time_model,run1_swa_model(:,1),'b-', ...
    time_model,run2_swa_model(:,1), 'r-', ...
    datenum(time_isfs1),sw_stn1, 'k-')
title('ISFS 1')
datetick('x','mm/dd','keepticks')
axis([min(time_model) max(time_model) -5 250])

subplot(2,3,2)
plot(time_model,run1_swa_model(:,2),'b-', ...
    time_model,run2_swa_model(:,2),'r-', ...
    datenum(time_isfs2),sw_stn2,'k-')
title('ISFS 2')
datetick('x','mm/dd','keepticks')
axis([min(time_model) max(time_model) -5 250])

subplot(2,3,3)
plot(time_model,run1_swa_model(:,3),'b-', ...
    time_model,run2_swa_model(:,3),'r-', ...
    datenum(time_isfs3),sw_stn3,'k-')
title('ISFS 3')
datetick('x','mm/dd','keepticks')
axis([min(time_model) max(time_model) -5 250])

subplot(2,3,4)
plot(time_model,run1_swa_model(:,4),'b-', ...
    time_model,run2_swa_model(:,4),'r-', ...
    datenum(time_isfs4),sw_stn4,'k-')
title('ISFS 4')
datetick('x','mm/dd','keepticks')
axis([min(time_model) max(time_model) -5 250])

subplot(2,3,5)
plot(time_model,run1_swa_model(:,5),'b-', ...
    time_model,run2_swa_model(:,5),'r-', ...
    datenum(time_isfs5),sw_stn5,'k-')
title('ISFS 5')
datetick('x','mm/dd','keepticks')
axis([min(time_model) max(time_model) -5 250])

subplot(2,3,6)
plot(time_model,run1_swa_model(:,6),'b-', ...
    time_model,run2_swa_model(:,6),'r-', ...
    datenum(time_isfs6),sw_stn6,'k-')
title('ISFS 6')
datetick('x','mm/dd','keepticks')
axis([min(time_model) max(time_model) -5 250])

% % ISFS % %

for i = 1:7
figure(i)
subplot(3,1,1)
% plot(time_model,run1_hfx_model(:,i),'r-','linewidth',2)
% hold on
% plot(time_model,run1_lhf_model(:,i),'b-','linewidth',2)
% hold on
% plot(time_model,run1_grdflx_model(:,i),'k-','linewidth',2)
% title('Flux Model Data Time Series, Logan, UT')
% datetick('x','mm/dd','keeplimits','keepticks')
% ylabel('Energy Flux (W/m^2)')
% ylim([-100 100])
% legend('HFX BASE','LHF BASE','GRDFLX BASE','location','southeast')
plot(time_model,run1_swd_model(:,i),'r-','linewidth',2)
hold on
plot(time_model,run1_swa_model(:,i),'b-','linewidth',2)
hold on
plot(time_model,run1_netlw_model(:,i),'k-','linewidth',2)
datetick('x','mm/dd','keeplimits','keepticks')
xlabel('Date')
ylabel('Energy Flux (W/m^2)')
ylim([-200 700])
legend('SWD BASE','SWA BASE','NETLW BASE','location','northeast')
subplot(3,1,2)
plot(time_model,run2_swd_model(:,i),'r-','linewidth',2)
hold on
plot(time_model,run2_swa_model(:,i),'b-','linewidth',2)
hold on
plot(time_model,run2_netlw_model(:,i),'k-','linewidth',2)
datetick('x','mm/dd','keeplimits','keepticks')
xlabel('Date')
ylabel('Energy Flux (W/m^2)')
ylim([-200 700])
legend('SWD USGS','SWA USGS','NETLW USGS','location','northeast')
subplot(3,1,3)
plot(time_model,run1_swd_model(:,i)-run2_swd_model(:,i),'r-','linewidth',2)
hold on
plot(time_model,run1_swa_model(:,i)-run2_swa_model(:,i),'b-','linewidth',2)
hold on
plot(time_model,run1_netlw_model(:,i)-run2_netlw_model(:,i),'k-','linewidth',2)
datetick('x','mm/dd','keeplimits','keepticks')
xlabel('Date')
ylabel('Energy Flux (W/m^2)')
ylim([-200 200])
legend('SWD BASE-USGS','SWA BASE-USGS','NETLW BASE-USGS','location','northeast')
saveas(gcf,['/uufs/chpc.utah.edu/common/home/u0875930/public_html/paper_figs/energy_budget/ISFS_3_Panel' int2str(i) '.png'])
end

