clear all
close all

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

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

load('station_data.mat');

time_stn2 = data(:,1);
temp_stn = data(:,2);

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;

terrain = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run,'/wrfout_d03_2011-01-01_00:00:00'),'HGT');
mask = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run,'/wrfout_d03_2011-01-01_00:00:00'),'LANDMASK');
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');
land = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run,'/wrfout_d03_2011-01-01_00:00:00'),'IVGTYP');

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

% Plot for comparison.

x_min = 212;   % Valley domain boundaries (find min height).
x_max = 240;
y_min = 290;
y_max = 322;

% Find out what land use type the grid point is...

land_t = land(:,:,1);

station_land_use = [land_t((x_grid(1)-1):(x_grid(1)+1),(y_grid(1)-1):(y_grid(1)+1)); ...
    land_t((x_grid(2)-1):(x_grid(2)+1),(y_grid(2)-1):(y_grid(2)+1)); ...
    land_t((x_grid(3)-1):(x_grid(3)+1),(y_grid(3)-1):(y_grid(3)+1)); ...
    land_t((x_grid(4)-1):(x_grid(4)+1),(y_grid(4)-1):(y_grid(4)+1)); ...
    land_t((x_grid(5)-1):(x_grid(5)+1),(y_grid(5)-1):(y_grid(5)+1)); ...
    land_t((x_grid(6)-1):(x_grid(6)+1),(y_grid(6)-1):(y_grid(6)+1)); ...
    land_t((x_grid(7)-1):(x_grid(7)+1),(y_grid(7)-1):(y_grid(7)+1))];

station_land_use = reshape(station_land_use',[3,3,7]);

stn_unique = unique(station_land_use);

nlcd_lnd_id = {'Evergreen Needleleaf Forest', ...  
          'Evergreen Broadleaf Forest', ...
          'Deciduous Needleleaf Forest', ... 
          'Deciduous Broadleaf Forest', ... 
          'Mixed Forest', ... 
          'Closed Shrubland', ... 
          'Open Shrubland', ... 
          'Woody Savanna', ... 
          'Savanna', ... 
          'Grassland', ... 
          'Permanent Wetland', ... 
          'Cropland', ... 
          'Urban and Built-Up', ...  
          'Cropland / Natural Veg. Mosaic', ... 
          'Permanent Snow', ... 
          'Barren / Sparsely Vegetated', ... 
          'IGBP Water', ... 
          'Unclassified', ...
          'Fill Value', ...
          'Unclassified', ...
          'Open Water', ...
          'Perennial Ice/Snow', ... 
          'Developed Open Space', ...  
          'Developed Low Intensity', ...  
          'Developed Medium Intensity', ...  
          'Developed High Intensity', ...  
          'Barren Land', ... 
          'Deciduous Forest', ... 
          'Evergreen Forest', ...
          'Mixed Forest', ... 
          'Dwarf Scrub', ... 
          'Shrub/Scrub', ... 
          'Grassland/Herbaceous', ... 
          'Sedge/Herbaceous', ... 
          'Lichens', ... 
          'Moss', ... 
          'Pasture/Hay', ... 
          'Cultivated Crops', ... 
          'Woody Wetland', ... 
          'Emergent Herbaceous Wetland' 
          };
      
luc_unique = ['blank', nlcd_lnd_id(stn_unique)];     
      
for i = 1:length(stn_unique)
    station_land_use(find(station_land_use == stn_unique(i))) = i ;
end

station_land_use(find(station_land_use == 8)) = 9 ;

nlcd_lnd_id_10 = {'Pasture/Hay', ...
    'Water', ... 
    'Wetland', ...
    'Developed High Intensity', ...
    'Developed Medium Intensity', ... 
    'Developed Low Intensity', ...
    'Developed Open Space', ...  
    'Evergreen Forest', ...
    'Deciduous Forest', ... 
    'Cultivated Crops', ...
    'Grassland', ...
    'Shrubland', ... 
    'Barren Land' 
    };

cmap_2011 = colormap([0.8 1.0 0.2; ...
    0.0 0.0 1.0; ...
    0.0 1.0 1.0; ...
    0.3 0.3 0.3; ...
    0.4 0.4 0.4; ...
    0.5 0.5 0.5; ...
    0.6 0.6 0.6; ...
    0.0 0.4 0.2; ...
    0.2 0.6 0.2; ...
    0.3 0.7 0.0; ...
    0.0 1.0 0.0; ...
    0.8 0.4 0.2; ...
    0.6 0.4 0.0; ...
    ]);

figure(1)
ax = usamap([lati(x_min,y_min),lati(x_max,y_max)],[long(x_min,y_max),long(x_max,y_min)])
surfm(lati,long,land(:,:,1))
hold on
mod = plotm(lat_grid, lon_grid, 'w.', 'MarkerSize', 20)
hold on
iss = plotm(lat(8:9), lon(8:9), 'k.', 'MarkerSize', 20)
hold on
time_height = plotm(lati(225,302),long(225,302),'k*')
legend([mod iss time_height], {'ISFS Station','ISS Station','Model Time Height'})
title('Location of ISFS Stations and Closest Model Grid Point')


figure(2)
for i = 1:7
   subplot(2,4,i)
   ax(i) = imagesc(station_land_use(:,:,i),[1 8]);
   hold on
   plot(2,2,'k.','MarkerSize',50)
   title(strcat('ISFS ',int2str(i)))
   set(gca,'Ydir','normal')
   set(gca,'YTick',[1.5,2.5,3.5],'YTickLabel',[1.3,2.7,4])
   set(gca,'XTick',[1.5,2.5,3.5],'XTickLabel',[1.3,2.7,4])
   caxis([1 8])
   xlabel('Distance (km)')
   if i == 1
      ylabel('Distance (km)')
   elseif i == 5
      ylabel('Distance (km)')
   end
end

h=colorbar;
set(h,'Position', [.7425 .106 .028 .34]);
set(h,'YTick',0.5:9.5,'YTickLabel',luc_unique)
contourcmap('jet',1)

%%

x_min = 150;    % Boundaries for modification to snow.
x_max = 425;
y_min = 180;
y_max = 425;

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

lx_min = 160; 
lx_max = 240;
ly_min = 290;
ly_max = 405;

%%

gray_map = colormap(gray);
gray_map = [gray_map(30:64,:)];

%%

ew_min  = 215;
ew_max  = 235;
y_pt    = 300;

ns_min  = 290;
ns_max  = 340;
x_pt    = 225;

x_min = 210;   % Valley domain boundaries (find min height).
x_max = 240;
y_min = 290;
y_max = 340;


figure('Position',[100 100 500 600])
ax = usamap([lati(lx_min,ly_min),lati(lx_max,ly_max)],[long(lx_min,ly_max),long(lx_max,ly_min)]);
set(gca, 'XAxisLocation', 'top')
contourfm(lati,long,double(terrain(:,:,1)),[1200:150:3000],'k')
hold on
colormap(gray_map)
h1 = colorbar;
caxis([1200 3000])
set(h1,'YTick',[1200:300:3600],'Location','southoutside')
xlabel(h1, 'Elevation (m ASL)','FontSize',10)
contourm(lati,long,double(terrain(:,:,1)),[1200:300:300],'k','LineWidth',2)
hold on
contourm(lati,long,double(mask(:,:,1)),[1],'k','LineWidth',2)
hold on
mod = plotm(lat_grid, lon_grid, 'w.', 'MarkerSize', 20)
hold on
iss = plotm(lat(8:9), lon(8:9), 'k.', 'MarkerSize', 20)
hold on
slc = plotm(lat(17), lon(17), 'k.', 'MarkerSize', 20)
hold on
log = plotm(lat(23), lon(23), 'k.', 'MarkerSize', 20)
textm(lat(23)-0.0225, lon(23)-0.2,'Logan','Color','White')
hold on
syr = plotm(lat(30), lon(30), 'k.', 'MarkerSize', 20)
textm(lat(30)-0.02, lon(30)-0.28,'Syracuse','Color','White')
textm(41.175,-112.725,{'Great','Salt Lake'},'Color','White')
textm(40.375,-113.23,'(b)','Color','White','FontSize',14)
hold on
v_sub = plotm(lati(vx_min:vx_max,vy_min),long(vx_min:vx_max,vy_min),'k-','LineWidth',2)
hold on
v_sub = plotm(lati(vx_min:vx_max,vy_max),long(vx_min:vx_max,vy_max),'k-','LineWidth',2)
hold on
v_sub = plotm(lati(vx_min,vy_min:vy_max),long(vx_min,vy_min:vy_max),'k-','LineWidth',2)
hold on
v_sub = plotm(lati(vx_max,vy_min:vy_max),long(vx_max,vy_min:vy_max),'k-','LineWidth',2)

hold on
v_sub = plotm(lati(lx_min:lx_max,ly_min),long(lx_min:lx_max,ly_min),'k-','LineWidth',2)
hold on
v_sub = plotm(lati(lx_min:lx_max,ly_max),long(lx_min:lx_max,ly_max),'k-','LineWidth',2)
hold on
v_sub = plotm(lati(lx_min,ly_min:ly_max),long(lx_min,ly_min:ly_max),'k-','LineWidth',2)
hold on
v_sub = plotm(lati(lx_max,ly_min:ly_max),long(lx_max,ly_min:ly_max),'k-','LineWidth',2)

hold on
time_height_x = plotm(lati(ew_min:ew_max,y_pt),long(ew_min:ew_max,y_pt),'g-','LineWidth',1.5)
hold on
time_height_y = plotm(lati(x_pt,ns_min:ns_max),long(x_pt,ns_min:ns_max),'b-','LineWidth',1.5)

textm(lati(ew_min,y_pt),long(ew_min,y_pt),'A1','Color','Green','FontWeight','Bold')
textm(lati(ew_max,y_pt),long(ew_max,y_pt),'B1','Color','Green','FontWeight','Bold')

textm(lati(x_pt,ns_min),long(x_pt,ns_min),'A2','Color','Blue','FontWeight','Bold')
textm(lati(x_pt,ns_max),long(x_pt,ns_max),'B2','Color','Blue','FontWeight','Bold')
%% Pull albedo data from grid points of stations.

clear alb

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

alb = zeros(size(lati));
alb2 = zeros(size(lati));
temp = zeros(size(lati));
temp2 = zeros(size(lati));
temp3 = zeros(size(lati));
snowh = 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))
      
      alb_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'),'ALBEDO');
      
      alb = cat(3, alb,alb_temp);    
      
      alb2_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');
      
      alb2 = cat(3, alb2,alb2_temp);    
      
      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);
      
      temp3_temp = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run3,'/wrfout_d03_',int2str(year(d)),'-',mon,'-',day,'_',hour,':00:00'),'T2');
      
      temp3 = cat(3, temp3,temp3_temp);
      
      snowh_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'),'SNOWH');
      
      snowh = cat(3, snowh,snowh_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
   
alb(:,:,1) = []; % Remove first time step of zeros.
alb2(:,:,1) = [];
temp(:,:,1) = [];
temp2(:,:,1) = [];
temp3(:,:,1) = [];
snowh(:,:,1) = [];
time_model(1) = [];   % Remove first time (repeat).


%%

avg_alb = mean(alb,3);
flg = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run,'/wrfout_d03_2011-01-01_00:00:00'),'LANDMASK');
flg = double(flg(:,:,1));
terrain = double(terrain(:,:,1));

figure(1)
subplot(1,2,1)
title('USGS Albedo Hour 0')
ax1 = usamap([lati(lx_min,ly_min),lati(lx_max,ly_max)],[long(lx_min,ly_max),long(lx_max,ly_min)]);
surfm(lati,long,alb(:,:,1));
colorbar('YTick',[0.3:0.025:0.7])
contourcmap('jet',[0.3:0.005:0.7])
caxis([0.3 0.7])
contourm(lati,long,terrain,'k','LevelStep',125,'LineWidth',0.5)

subplot(1,2,2)
title('USGS Average Albedo')
ax1 = usamap([lati(lx_min,ly_min),lati(lx_max,ly_max)],[long(lx_min,ly_max),long(lx_max,ly_min)]);
surfm(lati,long,avg_alb);
colorbar('YTick',[0.3:0.025:0.7])
contourcmap('jet',[0.3:0.005:0.7])
caxis([0.3 0.7])
contourm(lati,long,terrain,'k','LevelStep',125,'LineWidth',0.5)


%%

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

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

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

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

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

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

s1 = subplot(3,1,1)

p = get(s1,'position');
p(3) = p(3)*0.85; 
set(s1, 'position', p);

plot(time_model,snowh_model,'linewidth',2)
set(gca,'XTickLabel','')
ylabel('Snow Depth (m)')

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

s2 = subplot(3,1,2)

p = get(s2,'position');
p(3) = p(3)*0.85; 
set(s2, 'position', p);

plot(time_model,alb_model,'linewidth',2)
hold on
plot(iop05(:,1),iop05(:,3:9),'.','markersize',10)
axis([min(time_model) max(time_model) 0.4 1])
set(gca,'XTickLabel','')
ylabel('Albedo')

fig_nudge(gca,0,0.06)

l = legend('Model ISFS1','Model ISFS2','Model ISFS3','Model ISFS4','Model ISFS5','Model ISFS6','Model ISFS7','ISFS1','ISFS2','ISFS3','ISFS4','ISFS5','ISFS6','ISFS7','location','northeastoutside')
set(l,'Position',[0.825 0.735 0.1 0.1])

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

s3 = subplot(3,1,3)

p = get(s3,'position');
p(3) = p(3)*0.85; 
set(s3, 'position', p);

plot(time_model,alb2_model,'linewidth',2)
hold on
plot(iop05(:,1),iop05(:,3:9),'.','markersize',10)
axis([min(time_model) max(time_model) 0.4 1])
datetick('x','mm/dd')
xlabel('Date (Ticks at 00 UTC)')
ylabel('Albedo')
set(gca,'YTick',[0.4:0.1:1],'YTickLabel',[0.4:0.1:1])

fig_nudge(gca,0,0.12)

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

%% Single ISFS Plots

close all

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

% % SYRACUSE % % 

time_temp = [time_stn2{21}];

for i = 1:length(time_temp)
    time_syr(i) = strcat(time_temp{i}(1:10),{' '},time_temp{i}(12:19));
end

s1 = subplot(4,1,1)

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

plot(datenum(time_syr),cell2mat(temp_stn(21)),'k','LineWidth',1.5)
hold on
plot(time_model,temp_model(:,11),'LineWidth',1.5)
hold on
plot(time_model,temp2_model(:,11),'LineWidth',1.5)
axis([min(time_model) max(time_model) -22.5 2.5])
set(gca,'XTickLabel','')
legend('Syracuse','BASE','DFLT','location','southeast')
ylabel('Temperature (^oC)')

text(min(time_model)+0.01,1.25,'(a)','FontSize',12)

s2 = subplot(4,1,2)

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

plot(data(:,1),data(:,2),'k-','linewidth',1.5)
hold on
plot(time_model,temp_model(:,1),'LineWidth',1.5)
hold on
plot(time_model,temp2_model(:,1),'LineWidth',1.5)
axis([min(time_model) max(time_model) -22.5 2.5])
set(gca,'XTickLabel','')
legend('ISFS 1','BASE','DFLT','location','southeast')
ylabel('Temperature (^oC)')

text(min(time_model)+0.01,1.25,'(b)','FontSize',12)

s3 = subplot(4,1,3)

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

plot(data(:,1),data(:,3),'k-','linewidth',1.5)
hold on
plot(time_model,temp_model(:,2),'LineWidth',1.5)
hold on
plot(time_model,temp2_model(:,2),'LineWidth',1.5)

set(gca,'xlim',[time_model(1) time_model(end)])
set(gca,'ylim',[ -22.5 2.5])
datetick('x','mm/dd','keeplimits')
set(gca,'XTickLabel','')

legend('ISFS 2','BASE','DFLT','location','southeast')
ylabel('Temperature (^oC)')

text(min(time_model)+0.01,1.25,'(c)','FontSize',12)

s4 = subplot(4,1,4)

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

plot(data(:,1),data(:,6),'k-','linewidth',1.5)
hold on
plot(time_model,temp_model(:,6),'LineWidth',1.5)
hold on
plot(time_model,temp2_model(:,6),'LineWidth',1.5)

set(gca,'xlim',[time_model(1) time_model(end)])
set(gca,'ylim',[ -22.5 2.5])
datetick('x','mm/dd','keeplimits')

legend('ISFS 6','BASE','DFLT','location','southeast')
xlabel('Date (Ticks at 00 UTC)')
ylabel('Temperature (^oC)')

text(min(time_model)+0.01,1.25,'(d)','FontSize',12)

%%

n = [1,2,5,6,3,4,7]

figure(5)
for i = 1:7
   subaxis(2,4,n(i))
   plot(data(:,1),data(:,1+i),'k-','linewidth',0.5)
   hold on
   plot(time_model,temp_model(:,i))
   hold on
   plot(time_model,temp2_model(:,i))
   %hold on
   %plot(time_model,temp3_model(:,i),'linewidth',1.25)
   axis([min(time_model) max(time_model) -22.5 2.5])
   set(gca,'XTick',linspace(min(time_model),max(time_model),4))
   datetick('x','dd-hh UTC','keeplimits','keepticks')
   title(strcat('ISFS Station ',int2str(i)))
   if i == 1
       ylabel('Temperature (^oC)')
   end
   if i == 5
       ylabel('Temperature (^oC)')
   end
   if i == 7
       ylabel('Temperature (^oC)')
   end
   if i == 3
       xlabel('Time (UTC)')
       ylabel('Temperature (^oC)')
   end
   if i == 4 
       xlabel('Time (UTC)')
   end
   if i == 6 
       xlabel('Time (UTC)')
   end
   if i == 7 
      xlabel('Time (UTC)')
   end
end
l = legend('Station Data','BASE','USGS','location','eastoutside');
%l = legend('Station Data','Dec 30','Dec 31','Jan 01','location','southeast')
set(l,'Position', [.76 .334 .02 .2]);

% % LOGAN % % 

time_temp = [time_stn2{7}];

for i = 1:length(time_temp)
    time_logan(i) = strcat(time_temp{i}(1:10),{' '},time_temp{i}(12:19));
end

figure(6)
plot(datenum(time_logan),cell2mat(temp_stn(7)),'k')
hold on
plot(time_model,temp_model(:,8))
hold on
plot(time_model,temp2_model(:,8))
hold on
plot(time_model,temp3_model(:,8))
axis([min(time_model) max(time_model) -25 0])
datetick('x','mm/dd','keeplimits','keepticks')
legend('Logan','DEC30','BASE','JAN01','location','southeast')
xlabel('Time (UTC)')
ylabel('Temperature (^oC)')

% % VERNAL % % 

time_temp = [time_stn2{8}];

for i = 1:length(time_temp)
    time_vern(i) = strcat(time_temp{i}(1:10),{' '},time_temp{i}(12:19));
end

figure(7)
plot(datenum(time_vern),cell2mat(temp_stn(8)),'k')
hold on
plot(time_model,temp_model(:,9))
hold on
plot(time_model,temp2_model(:,9))
%hold on
%plot(time_model,temp3_model(:,9))
axis([min(time_model) max(time_model) -27.5 -5])
datetick('x','dd','keeplimits','keepticks')
title('Station and Model Data Time Series, Vernal, UT')
legend('Station Data','NLCD 2011','USGS','location','southeast')

% % ROOSEVELT % % 

time_temp = [time_stn2{9}];

for i = 1:length(time_temp)
    time_roose(i) = strcat(time_temp{i}(1:10),{' '},time_temp{i}(12:19));
end
figure(8)
plot(datenum(time_roose),cell2mat(temp_stn(9)),'k')
hold on
plot(time_model,temp_model(:,10))
hold on
plot(time_model,temp2_model(:,10))
%hold on
%plot(time_model,temp3_model(:,10))
axis([min(time_model) max(time_model) -32.5 -5])
datetick('x','dd','keeplimits','keepticks')
title('Station and Model Data Time Series, Roosevelt, UT')
legend('Station Data','NLCD 2011','USGS','location','southeast')

% % SYRACUSE % % 

time_temp = [time_stn2{21}];

for i = 1:length(time_temp)
    time_syr(i) = strcat(time_temp{i}(1:10),{' '},time_temp{i}(12:19));
end

figure(9)
plot(datenum(time_syr),cell2mat(temp_stn(21)),'k','LineWidth',1.5)
hold on
plot(time_model,temp_model(:,11),'LineWidth',1.5)
hold on
plot(time_model,temp2_model(:,11),'LineWidth',1.5)
%hold on
%plot(time_model,temp3_model(:,11))
axis([min(time_model) max(time_model) -22.5 2.5])
set(gca,'XTick',linspace(min(time_model),max(time_model),9))
datetick('x','dd-hh UTC','keeplimits','keepticks')
legend('Syracuse','BASE','USGS','location','southeast')
xlabel('Time (UTC)')
ylabel('Temperature (^oC)')

%% Find the average model albedo at each station and the actual to compare.

mean_stn = mean(iop05(:,2:8))
mean_mod = mean(alb_model)
