clear all
close all

set(0,'defaultfigureposition',[500 500 500 420]);

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

run  = 'usgs-dec31';

load('station_data.mat');

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

%% Get radiation data from wrf output.

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));
swd = zeros(size(lati));
swu = 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);      
      
      swd_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'),'SWDNB');
      
      swd = cat(3,swd,swd_temp);
      
      swu_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'),'SWUPB');
      
      swu = cat(3,swu,swu_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.
swd(:,:,1) = [];
swu(:,:,1) = [];
time_model(1) = [];   % Remove first time (repeat).

%% Calculate albedo.

alb_rad = swu ./ swd;

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

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

%% Plot stuff.

figure('Position', [100, 100, 1000, 800])
plot(time_model,alb_rad_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')
legend('ISFS1','ISFS2','ISFS3','ISFS4','ISFS5','ISFS6','ISFS7','location','southeast')

latlim = [min(lati(:)) max(lati(:))];
lonlim = [min(long(:)) max(long(:))];

figure(2)
title('ALB')
ax1 = usamap(latlim,lonlim);
surfm(lati,long,alb_rad(:,:,20));
colorbar('YTick',[0:1])
contourcmap('jet',[0:0.1:1])
caxis([0 1])
contourm(lati,long,double(terrain(:,:,1)),'k','LevelStep',250)


