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];
mont =  [1,1,1,1,1,1,1,1];

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

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

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


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.
      
      temp_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'),'T');
      
      run1_temp = cat(4,run1_temp,temp_temp);      
      
      temp_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'),'T');
      
      run2_temp = cat(4,run2_temp,temp_temp);
      
      % Q Cloud from WRF.
      
%       cloud_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'),'QCLOUD');
%       
%       run1_cloud = cat(4,run1_cloud,cloud_temp);      
%       
%       cloud_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'),'QCLOUD');
%       
%       run2_cloud = cat(4,run2_cloud,cloud_temp);
      
      % Albedo from WRF.
      
      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);
      
      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_temp(:,:,:,1:6) = []; % Remove first time (repeat).
run2_temp(:,:,:,1:6) = [];
% run1_cloud(:,:,:,1:6) = [];
% run2_cloud(:,:,:,1:6) = [];
run1_alb(:,:,1:6) = [];
run2_alb(:,:,1:6) = [];

time_model(1) = [];   % Remove first time (repeat).

%% Calculate spacial averages of variables.

run1_temp_lev1 = reshape(run1_temp(:,:,1,:),[474,498,192]);
run2_temp_lev1 = reshape(run2_temp(:,:,1,:),[474,498,192]);

run1_avg_temp = mean(run1_temp_lev1,3);
run2_avg_temp = mean(run2_temp_lev1,3);

temp_diff = run1_avg_temp - run2_avg_temp;

% run1_cloud_lev10 = reshape(run1_cloud(:,:,1:10,:),[474,498,10,192]);
% run2_cloud_lev10 = reshape(run2_cloud(:,:,1:10,:),[474,498,10,192]);
% 
% run1_avg_cloud = mean(run1_cloud_lev10,3);
% run2_avg_cloud = mean(run2_cloud_lev10,3);
% 
% cloud_diff = run1_avg_cloud - run2_avg_cloud;

run1_avg_alb = mean(run1_alb,3);
run2_avg_alb = mean(run2_alb,3);

alb_diff = run1_avg_alb - run2_avg_alb;

%% 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_temp_day = reshape(run1_temp(:,:,1,day),[474,498,48]);
run2_temp_day = reshape(run2_temp(:,:,1,day),[474,498,48]);

run1_temp_night = reshape(run1_temp(:,:,1,night),[474,498,48]);
run2_temp_night = reshape(run2_temp(:,:,1,night),[474,498,48]);

% run1_cloud_day = reshape(run1_cloud(:,:,1:10,day),[474,498,10,48]);
% run2_cloud_day = reshape(run2_cloud(:,:,1:10,day),[474,498,10,48]);
% 
% run1_cloud_night = reshape(run1_cloud(:,:,1:10,night),[474,498,10,48]);
% run2_cloud_night = reshape(run2_cloud(:,:,1:10,night),[474,498,10,48]);

% Average

run1_avg_temp_day = mean(run1_temp_day,3);
run2_avg_temp_day = mean(run2_temp_day,3);

run1_avg_temp_night = mean(run1_temp_night,3);
run2_avg_temp_night = mean(run2_temp_night,3);

% run1_avg_cloud_day = mean(mean(run1_cloud_day,4),3);
% run2_avg_cloud_day = mean(mean(run2_cloud_day,4),3);
% 
% run1_avg_cloud_night = mean(mean(run1_cloud_night,4),3);
% run2_avg_cloud_night = mean(mean(run2_cloud_night,4),3);

% Difference

temp_diff_day   = run1_avg_temp_day - run2_avg_temp_day;
temp_diff_night = run1_avg_temp_night - run2_avg_temp_night;

% cloud_diff_day   = run1_avg_cloud_day - run2_avg_cloud_day;
% cloud_diff_night = run1_avg_cloud_night - run2_avg_cloud_night;

%% Plot comparison.

close all

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

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

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

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

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

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

basin_temp_diff = temp_diff(x_min:x_max,y_min:y_max);
basin_temp_diff_day = temp_diff_day(x_min:x_max,y_min:y_max);
basin_temp_diff_night = temp_diff_night(x_min:x_max,y_min:y_max);
basin_cloud_diff = cloud_diff(x_min:x_max,y_min:y_max);
basin_terrain = terrain(x_min:x_max,y_min:y_max);
basin_flg = flg(x_min:x_max,y_min:y_max);

lake_temp_diff = temp_diff(lx_min:lx_max,ly_min:ly_max);
lake_temp_diff_day = temp_diff_day(lx_min:lx_max,ly_min:ly_max);
lake_temp_diff_night = temp_diff_night(lx_min:lx_max,ly_min:ly_max);
lake_cloud_diff = cloud_diff(lx_min:lx_max,ly_min:ly_max);
lake_cloud_diff_day = cloud_diff_day(lx_min:lx_max,ly_min:ly_max);
lake_cloud_diff_night = cloud_diff_night(lx_min:lx_max,ly_min:ly_max);
lake_terrain = terrain(lx_min:lx_max,ly_min:ly_max);
lake_flg = flg(lx_min:lx_max,ly_min:ly_max);


%%

title2 = {'Average Temperature Difference 01/01/11-01/10/11,','NLCD 31, Dec 2010 - NLCD 01, Jan 2011'}

figure(5)
title(title2)
ax2 = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
surfm(lat(lx_min:lx_max,ly_min:ly_max),lon(lx_min:lx_max,ly_min:ly_max),lake_temp_diff);
colorbar('YTick',[-2:0.5:2])
contourcmap('jet',[-2:0.25:2])
caxis([-2 2])
contourm(lat(lx_min:lx_max,ly_min:ly_max),lon(lx_min:lx_max,ly_min:ly_max),lake_flg,'k','LevelStep',1,'LineWidth',2)
contourm(lat(lx_min:lx_max,ly_min:ly_max),lon(lx_min:lx_max,ly_min:ly_max),lake_terrain,'k','LevelStep',125,'LineWidth',0.5)

title3= {'Average Temperature Difference 01/01/11-01/10/11,','NLCD 31, Dec 2010 - NLCD 01, Jan 2011'}

figure(6)
title(title3)
ax2 = usamap([lat(x_min,y_min),lat(x_max,y_max)],[lon(x_min,y_max),lon(x_max,y_min)]);
surfm(lat(x_min:x_max,y_min:y_max),lon(x_min:x_max,y_min:y_max),basin_temp_diff);
colorbar('YTick',[-2:0.5:2])
contourcmap('jet',[-2:0.25:2])
caxis([-2 2])
contourm(lat(x_min:x_max,y_min:y_max),lon(x_min:x_max,y_min:y_max),basin_flg,'k','LevelStep',1,'LineWidth',2)
contourm(lat(x_min:x_max,y_min:y_max),lon(x_min:x_max,y_min:y_max),basin_terrain,'k','LevelStep',125,'LineWidth',0.5)

figure(7)
ax3 = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
surfm(lat(lx_min:lx_max,ly_min:ly_max),lon(lx_min:lx_max,ly_min:ly_max),lake_cloud_diff);
colorbar('YTick',[-2.5e-5:0.5e-5:2.5e-5])
contourcmap('jet',[-2.5e-5:0.5e-5:2.5e-5])
contourm(lat(lx_min:lx_max,ly_min:ly_max),lon(lx_min:lx_max,ly_min:ly_max),lake_flg,'k','LevelStep',1,'LineWidth',2)
contourm(lat(lx_min:lx_max,ly_min:ly_max),lon(lx_min:lx_max,ly_min:ly_max),lake_terrain,'k','LevelStep',125,'LineWidth',0.5)
str = {'Cloud Water Mixing Ratio Difference (qcloud),',...
    '01/01/11-01/10/11, NLCD 31, Dec 2010 - USGS 31, Dec 2010'}
title(str)

figure(8)
title('Average Cloud Difference , NLCD 2011 - USGS')
ax1 = usamap(latlim,lonlim);
temp_ax = surfm(lat,lon,cloud_diff);
geoshow(ax1, states, 'FaceColor', 'none');
colorbar('YTick',[-4:4])
contourcmap('jet',[-4:0.5:4])
caxis([-4 4])
contourm(lat,lon,terrain,'k','LevelStep',250)
contourm(lat,lon,flg,'k','LevelStep',1,'LineWidth',2)

%%

close all

title4 = {'Day'}

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

subplot(1,2,1)
title(title4,'FontSize',12)
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(temp_diff_day),[-4:0.25:4],'LineStyle','none');
%colorbar('YTick',[-4:0.5: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)

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

title5 = {'Night'}

subplot(1,2,2)
title(title5,'FontSize',12)
ax2 = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
plabel('off')
contourfm(double(lat),double(lon),double(temp_diff_night),[-4:0.25:4],'LineStyle','none');
%colorbar('YTick',[-4:0.5:4])
fig_nudge(ax2,-0.08,0)
caxis([-4 4])
contourm(lat,lon,flg,'k','LevelStep',1,'LineWidth',2)
contourm(lat,lon,terrain,'k','LevelStep',125,'LineWidth',0.5)

c = colorbar;
contourcmap('jet',[-4:0.25:4])
hL = title(c,'^oC','FontSize',10); 
set(c,'position',[0.84 0.245 0.012 0.55]) 

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

%%

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

ax3 = 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(alb_diff),[-0.5:0.05:0.5],'LineStyle','none');
colorbar('YTick',[-0.5:0.1:0.5],'FontSize',12)
contourcmap('jet')
caxis([-0.5 0.5])
contourm(lat,lon,flg,'k','LevelStep',1,'LineWidth',2)
contourm(lat,lon,terrain,'k','LevelStep',125,'LineWidth',0.5)
c = colorbar;
set(c,'fontsize',10);
cpos = c.Position;
cpos(3) = 0.5*cpos(3);
c.Position = cpos;
set(c,'position',[0.84 0.195 0.02 0.655]) 