clear all
close all

%%

load('../stats_project/data/double_blend.mat');
load('../stats_project/data/pcaps_timeheight.mat');

master = master';

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


days = [1    2    3    4    5    6    7    8 ];%[10   11   12   13   14   15   16   17];
year = [2011 2011 2011 2011 2011 2011 2011 2011];%[2016 2016 2016 2016 2016 2016 2016 2016];
mont = [1    1    1    1    1    1    1    1 ];%[2    2    2    2    2    2    2    2];

temp  = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run,'/wrfout_d03_2011-01-01_00:00:00'),'T');
geo   = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run,'/wrfout_d03_2011-01-01_00:00:00'),'PH');
geo_p = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run,'/wrfout_d03_2011-01-01_00:00:00'),'PHB');
cloud = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run,'/wrfout_d03_2011-01-01_00:00:00'),'QCLOUD');

temp2  = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run2,'/wrfout_d03_2011-01-01_00:00:00'),'T');
geo2   = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run2,'/wrfout_d03_2011-01-01_00:00:00'),'PH');
geo_p2 = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run2,'/wrfout_d03_2011-01-01_00:00:00'),'PHB');
cloud2 = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run2,'/wrfout_d03_2011-01-01_00:00:00'),'QCLOUD');


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/',run,'/wrfout_d03_',int2str(year(d)),'-',mon,'-',day,'_',hour,':00:00'),'T');
      
      temp = cat(4,temp,temp_temp);  
         
      temp_geo = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run,'/wrfout_d03_',int2str(year(d)),'-',mon,'-',day,'_',hour,':00:00'),'PH');
      
      geo = cat(4,geo,temp_geo); 
      
      temp_geo_p = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run,'/wrfout_d03_',int2str(year(d)),'-',mon,'-',day,'_',hour,':00:00'),'PHB');
      
      geo_p = cat(4,geo_p,temp_geo_p);   
      
      temp_cloud = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run,'/wrfout_d03_',int2str(year(d)),'-',mon,'-',day,'_',hour,':00:00'),'QCLOUD');
      
      cloud = cat(4,cloud,temp_cloud);   
       
      temp_temp2 = 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');
      
      temp2 = cat(4,temp2,temp_temp2);  
         
      temp_geo2 = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run2,'/wrfout_d03_',int2str(year(d)),'-',mon,'-',day,'_',hour,':00:00'),'PH');
      
      geo2 = cat(4,geo2,temp_geo2); 
      
      temp_geo_p2 = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run2,'/wrfout_d03_',int2str(year(d)),'-',mon,'-',day,'_',hour,':00:00'),'PHB');
      
      geo_p2 = cat(4,geo_p2,temp_geo_p2);   
      
      temp_cloud2 = 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');
      
      cloud2 = cat(4,cloud2,temp_cloud2); 
      
      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:6) = []; % Remove first time (repeat).
cloud(:,:,:,1:6) = []; % Remove first time (repeat).

temp = temp + 300;
cloud = cloud .* 1000;

geo(:,:,:,1:6) = [];
geo_p(:,:,:,1:6) = [];

hgt = (geo + geo_p) ./ 9.8;

temp2(:,:,:,1:6) = []; % Remove first time (repeat).
cloud2(:,:,:,1:6) = []; % Remove first time (repeat).

temp2 = temp2 + 300;
cloud2 = cloud2 .* 1000;

geo2(:,:,:,1:6) = [];
geo_p2(:,:,:,1:6) = [];

hgt2 = (geo2 + geo_p2) ./ 9.8;

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

%% Set up time height.
 
% x_pt = 225;
% y_pt = 302;

%x_pt = 219;
%y_pt = 343;

x_pt = 225;
y_pt = 301;

th_temp = temp(x_pt-1:x_pt+1,y_pt-1:y_pt+1,1:40,:);
th_cloud = cloud(x_pt-1:x_pt+1,y_pt-1:y_pt+1,1:40,:);
th_hgt = hgt(x_pt-1:x_pt+1,y_pt-1:y_pt+1,1:40,:);

th_avg_temp = mean(mean(th_temp,1),2);
th_avg_cloud = mean(mean(th_cloud,1),2);
th_avg_hgt = mean(mean(th_hgt,1),2);

temp_fin = reshape(th_avg_temp,[40,length(time_model)]);
cloud_fin = reshape(th_avg_cloud,[40,length(time_model)]);
hgt_fin = reshape(th_avg_hgt,[40,length(time_model)]);

th_temp2 = temp2(x_pt-1:x_pt+1,y_pt-1:y_pt+1,1:40,:);
th_cloud2 = cloud2(x_pt-1:x_pt+1,y_pt-1:y_pt+1,1:40,:);
th_hgt2 = hgt2(x_pt-1:x_pt+1,y_pt-1:y_pt+1,1:40,:);

th_avg_temp2 = mean(mean(th_temp2,1),2);
th_avg_cloud2 = mean(mean(th_cloud2,1),2);
th_avg_hgt2 = mean(mean(th_hgt2,1),2);

temp_fin2 = reshape(th_avg_temp2,[40,length(time_model)]);
cloud_fin2 = reshape(th_avg_cloud2,[40,length(time_model)]);
hgt_fin2 = reshape(th_avg_hgt2,[40,length(time_model)]);

%%

start_index = find(master == time_model(1));
end_index = find(master == time_model(end));

temp_obs = theta(:,start_index:end_index);
rh_obs = RHth(:,start_index:end_index);

hgt_rnd = round(hgt_fin(1:26,:),-1);
temp_rnd = zeros(size(hgt(hgt_rnd)));

for i = 1:length(hgt_rnd(:,1))
    for j = 1:length(hgt_rnd(1,:))
        temp_rnd(i,j) = temp_obs(find(hgts == hgt_rnd(i,j)),j);
        rh_rnd(i,j) = rh_obs(find(hgts == hgt_rnd(i,j)),j);
    end 
end

hgt_rnd2 = round(hgt_fin2(1:26,:),-1);
temp_rnd2 = zeros(size(hgt(hgt_rnd2)));

for i = 1:length(hgt_rnd2(:,1))
    for j = 1:length(hgt_rnd2(1,:))
        temp_rnd2(i,j) = temp_obs(find(hgts == hgt_rnd2(i,j)),j);
    end 
end

x = zeros(size(hgt_rnd));

for i = 1:26
   x(i,:) = time_model';
end

temp_diff = temp_fin(1:26,:) - temp_rnd;

temp_diff2 = temp_fin2(1:26,:) - temp_rnd2;

%%

cntrs = [251,252,253,254,256,257,258,259, ...
    261,262,263,264,266,267,268,269, ...
    271,272,273,274,276,277,278,279, ...
    281,282,283,284,286,287,288,289, ...
    291,292,293,294,296,297,298,299, ...
    301,302,303,304];

%%



%%

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

s1 = subplot(5,1,1)

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

contourf(x,hgt_rnd(1:26,:),temp_rnd(1:26,:),cntrs)
colormap(s1,'jet')
hold on;
[C1,h1] = contour(x,hgt_rnd(1:26,:),temp_rnd(1:26,:),[250:5:305],'k-','linewidth',2);
clabel(C1,h1);
caxis([250 305])
hold on;
[C2,h2] = contour(x,hgt_rnd(1:26,:),rh_rnd(1:26,:),[90 90],'b-.','linewidth',2);
set(gca,'xlim',[time_model(1) time_model(end)])
set(gca,'ylim',[min(min(hgt_fin)) min(hgt_fin(25,:))])
set(gca,'XTickLabel','')
ylabel('Height ASL (m)')

fl = fig_label('a','ul',18,0.0025,0.0525,0.205,230);
set(fl,'color','k');

s2 = subplot(5,1,2)

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

contourf(x,hgt_fin(1:26,:),temp_fin(1:26,:),cntrs)
hold on
[C1,h1] = contour(x,hgt_fin(1:26,:),temp_fin(1:26,:),[250:5:305],'k-','linewidth',2);
clabel(C1,h1);
colormap(s2,'jet')
caxis([250 305])
hold on;
[C2, h2] = contour(x,hgt_fin(1:26,:),cloud_fin(1:26,:),[0.1 0.1],'b-.','linewidth',2);
set(gca,'xlim',[time_model(1) time_model(end)])
set(gca,'ylim',[min(min(hgt_fin)) min(hgt_fin(25,:))])
set(gca,'XTickLabel','')
ylabel('Height ASL (m)')

fl = fig_label('b','ul',14,0.004,0.08,0.205,230);
set(fl,'color','k');

s3 = subplot(5,1,3)

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

contourf(x,hgt_rnd,temp_diff,[-8:1:8],'LineStyle','none')
hold on;
[C1, h1] = contour(x,hgt_rnd,temp_diff,[-8:1:8],'k-');
clabel(C1,h1,[-8:2:8]);
colormap(s3,redblue)
caxis([-8 8])
axis([min(time_model) max(time_model) min(hgt_rnd(1,:)) min(hgt_rnd(26,:))])
set(gca,'XTickLabel','')
ylabel('Height ASL (m)')

fl = fig_label('c','ul',18,0.003935,0.0525,0.205,230);
set(fl,'color','k');

s4 = subplot(5,1,4)

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

contourf(x,hgt_fin2(1:26,:),temp_fin2(1:26,:),cntrs)
hold on
[C1,h1] = contour(x,hgt_fin2(1:26,:),temp_fin2(1:26,:),[250:5:305],'k-','linewidth',2);
clabel(C1,h1);
colormap(s4,'jet')
caxis([250 305])
hold on;
[C2, h2] = contour(x,hgt_fin2(1:26,:),cloud_fin2(1:26,:),[0.1 0.1],'b-.','linewidth',2);
set(gca,'xlim',[time_model(1) time_model(end)])
set(gca,'ylim',[min(min(hgt_fin)) min(hgt_fin(25,:))])
set(gca,'XTickLabel','')
ylabel('Height ASL (m)')

fl = fig_label('d','ul',14,0.004,0.08,0.205,230);
set(fl,'color','k');

s5 = subplot(5,1,5)

p = get(s5,'position');
p(4) = p(4)*1.3;
set(s5, 'position', p);

contourf(x,hgt_rnd2,temp_diff2,[-8:1:8],'LineStyle','none')
hold on
[C1, h1] = contour(x,hgt_rnd2,temp_diff2,[-8:1:8],'k-');
clabel(C1,h1,[-8:2:8]);
colormap(s5,redblue)
caxis([-8 8])
datetick('x','mm/DD','keepticks')
axis([min(time_model) max(time_model) min(hgt_rnd2(1,:)) min(hgt_rnd2(26,:))])
ylabel('Height ASL (m)')

%fig_nudge(gca,0,0.055)

fl = fig_label('e','ul',14,0.004,0.07,0.190,230);
set(fl,'color','k');

