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,2011,2011,2011];
mont =  [1,1,1,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'),'U');
run2_temp = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run2,'/wrfout_d03_2011-01-01_00:00:00'),'U');

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

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'),'U');
      
      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'),'U');
      
      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'),'V');
      
      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'),'V');
      
      run2_cloud = cat(4,run2_cloud,cloud_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) = [];
time_model(1) = []; 

%run1_cloud_lev10 = mean(reshape(run1_cloud(:,:,1:10,:),[474,498,10,192]),3); % Average first 10 levels.
%run2_cloud_lev10 = mean(reshape(run2_cloud(:,:,1:10,:),[474,498,10,192]),3);

run1_u = reshape(run1_temp(:,:,1,:),[474,499,192]);
run2_u = reshape(run2_temp(:,:,1,:),[474,499,192]);

run1_v = reshape(run1_cloud(:,:,1,:),[474,499,192]);
run2_v = reshape(run2_cloud(:,:,1,:),[474,499,192]);

%%

% lx_min = 155; 
% lx_max = 235;
% ly_min = 285;
% ly_max = 400;
% 
% lake_temp1 = run1_temp(lx_min:lx_max,ly_min:ly_max,:);% - 273.14;
% lake_temp2 = run2_temp(lx_min:lx_max,ly_min:ly_max,:);
% lake_cloud1 = run1_cloud(lx_min:lx_max,ly_min:ly_max,:);
% lake_cloud2 = run2_cloud(lx_min:lx_max,ly_min:ly_max,:);

wind = (run1_u.^2 + run1_v.^2).^(1/2)

% lake_terrain = terrain(lx_min:lx_max,ly_min:ly_max);
% lake_flg = flg(lx_min:lx_max,ly_min:ly_max);

%%

run = 'NLCD2011'

for i = 1:length(time_model)
   
   figure(1);
   ax = 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_temp1(:,:,i));
   %colorbar('YTick',[0:0.25e-5:2.5e-5])
   %contourcmap('pink',[0:0.25e-5:2.5e-5])
   colorbar('YTick',[-20:1:10])
   contourcmap('jet',[-20:0.25:10])
   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 = {'Model Temperature', datestr(time_model(i))}
   title(str)
   
   filename = strcat('/uufs/chpc.utah.edu/common/home/u0875930/public_html/wrf_output/iop5/temperature/hourly/',run,datestr(time_model(i)));
   saveas(1,filename,'png')
   
%    drawnow
%    frame = getframe(1);
%    im = frame2im(frame);
%    [imind,cm] = rgb2ind(im,256);
% 
%       if i == 1;
%          imwrite(imind,cm,filename,'gif','Loopcount',inf);
%       else
%          imwrite(imind,cm,filename,'gif','WriteMode','append')
%       end
%    
end
