clear all
close all

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

master = master';

%%

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

days = [1, 2, 3, 4, 5, 6, 7, 8];
year = [2011,2011,2011,2011,2011,2011,2011,2011,2011];
mont =  [1,1,1,1,1,1,1,1,1];

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

temp1  = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run1,'/wrfout_d03_2011-01-01_00:00:00'),'T');
geo1   = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run1,'/wrfout_d03_2011-01-01_00:00:00'),'PH');
geo_p1 = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run1,'/wrfout_d03_2011-01-01_00:00:00'),'PHB');
cloud1 = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run1,'/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_temp1 = 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');
      
      temp1 = cat(4,temp1,temp_temp1); 
      
      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_cloud1 = 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');
      
      cloud1 = cat(4,cloud1,temp_cloud1); 
      
      
      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).
temp1(:,:,:,1:6) = []; % Remove first time (repeat).

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

temp = temp + 300;
temp1 = temp1 + 300;

cloud = cloud .* 1000;
cloud1 = cloud1 .* 1000;

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

hgt = (geo + geo_p) ./ 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;

top = 25;

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

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,[top,length(time_model)]);
cloud_fin = reshape(th_avg_cloud,[top,length(time_model)]);
hgt_fin = reshape(th_avg_hgt,[top,length(time_model)]);

th_temp1 = temp1(x_pt-1:x_pt+1,y_pt-1:y_pt+1,1:25,:);
th_cloud1 = cloud1(x_pt-1:x_pt+1,y_pt-1:y_pt+1,1:25,:);

th_avg_temp1 = mean(mean(th_temp1,1),2);
th_avg_cloud1 = mean(mean(th_cloud1,1),2);

temp_fin1 = reshape(th_avg_temp1,[top,length(time_model)]);
cloud_fin1 = reshape(th_avg_cloud1,[top,length(time_model)]);

temp_diff = temp_fin - temp_fin1;

x = zeros(size(hgt_fin));

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

for i = 1:431
    xm(i,:) = master;
end

xm = double(xm);

%%

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

%%

close all

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

subplot(2,1,1)

contourf(double(master),hgts,theta,cntrs)
hold on
[C,h] = contour(double(master),hgts,theta,[250:5:305],'k-','linewidth',2);
clabel(C,h);
contourcmap('jet',[250:305]);
hold on
plot(time_model,zeros(length(time_model))+1800,'k--','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)')

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

subplot(2,1,2)

[Ch,hh] = contour(master,hgts,rh,[70:10:100],'k-','linewidth',1.5);
clabel(Ch,hh);
hold on
plot(time_model,zeros(length(time_model))+1800,'k--','LineWidth',2)
set(gca,'XTickLabel',[]);
set(gca,'xlim',[time_model(1) time_model(end)])
set(gca,'ylim',[min(min(hgt_fin)) min(hgt_fin(25,:))])
ylabel('Height ASL (m)')
datetick('x','mm/dd')
xlabel('Date (Ticks at 00 UTC)')

fl2 = fig_label('b','ul',18,0.003,0.07,0.205,230);
set(fl2,'color','k');

fig_nudge(gca,0,0.1)

%%
figure('Position',[100, 100, 900, 600])

s1 = subplot(3,1,1)

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

contourf(x,hgt_fin,temp_fin,cntrs)
hold on
[C1,h1] = contour(x,hgt_fin,temp_fin,[250:5:305],'k-','linewidth',2);
clabel(C1,h1);
contourcmap('jet',[250:305])
hold on
[C2,h2] = contour(x,hgt_fin,cloud_fin,[0.1 0.1],'b-.','LineWidth',2);
hold on
plot(time_model,zeros(length(time_model))+1800,'k--','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',14,0.0025,0.0525,0.205,230);
set(fl,'color','k');

s2 = subplot(3,1,2)

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

contourf(x,hgt_fin,temp_fin1,cntrs)
hold on
[C1,h1] = contour(x,hgt_fin,temp_fin1,[250:5:305],'k-','linewidth',2);
clabel(C1,h1);
contourcmap('jet',[250:305])
hold on
[C2,h2] = contour(x,hgt_fin,cloud_fin1,[0 0.1],'b-.','LineWidth',2);
hold on
plot(time_model,zeros(length(time_model))+1800,'k--','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(3,1,3)

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

contourf(x,hgt_fin,temp_diff,[-10:1:10],'LineStyle','none');
colormap(s3,redblue)
hold on
[C1, h1] = contour(x,hgt_fin,temp_diff,[-8:1:8],'k-');
clabel(C1,h1,[-8:2:8]);
caxis([-7 7])
set(gca,'xlim',[time_model(1) time_model(end)])
set(gca,'ylim',[min(min(hgt_fin)) 1800])
datetick('x','mm/dd')
xlabel('Date (Ticks at 00 UTC)')

fig_nudge(gca,0,0.1)

fl = fig_label('c','ul',14,0.004,0.1,0.205,100);
set(fl,'color','k');
