close all
clear all

load('station_data.mat')

run = 'usgs-dec31';

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.

x_stns = ones(1,length(loc));
y_stns = ones(1,length(loc));

%% Find points in domain that correspond to each station location.

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

%% Get data from model.

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));
temp = 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))
            
      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);
          
      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) = [];

temp = temp - 273.14;

%% Find model average temperature at each station.

temp_avg = mean(temp,3);

%% Find observed average temperature at each station.

for i = 1:length(data)
   T = cell2mat(data(i,2))
   stns_avg(i) = mean(T)
end

%% Plot cool stuff.

x_min = 160; 
x_max = 240;
y_min = 290;
y_max = 405;

figure(1)
title('USGS Dec31 Model vs Observed Temperature (1 Jan - 10 Jan 2011)')
ax1 = usamap([lati(x_min,y_min),lati(x_max,y_max)],[long(x_min,y_max),long(x_max,y_min)]);
temp_ax = surfm(lati(x_min:x_max,y_min:y_max),long(x_min:x_max,y_min:y_max),temp_avg(x_min:x_max,y_min:y_max))
hold on
stns_ax = scatterm(lat_stns,lon_stns,50,double(stns_avg),'d','filled','MarkerEdgeColor','k')
colorbar
contourcmap('jet',[-12:0.1:-4])
contourm(lati,long,double(terrain(:,:,1)),'k','LevelStep',300,'LineWidth',0.5)