clear all
close all

%%

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

master = master';

run  = 'jan2011_alldoms';

days = [1    2    3    4    5    6    7    8    9   ];%[10   11   12   13   14   15   16   17];
year = [2011 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    1   ];%[2    2    2    2    2    2    2    2];

temp  = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group3/horel_data/cfoster/wrf_daq_runs/',run,'/wrfout_d03_2011-01-01_00:00:00'),'T');
geo   = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group3/horel_data/cfoster/wrf_daq_runs/',run,'/wrfout_d03_2011-01-01_00:00:00'),'PH');
geo_p = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group3/horel_data/cfoster/wrf_daq_runs/',run,'/wrfout_d03_2011-01-01_00:00:00'),'PHB');
cloud = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group3/horel_data/cfoster/wrf_daq_runs/',run,'/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-group3/horel_data/cfoster/wrf_daq_runs/',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-group3/horel_data/cfoster/wrf_daq_runs/',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-group3/horel_data/cfoster/wrf_daq_runs/',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-group3/horel_data/cfoster/wrf_daq_runs/',run,'/wrfout_d03_',int2str(year(d)),'-',mon,'-',day,'_',hour,':00:00'),'QCLOUD');
      
      cloud = cat(4,cloud,temp_cloud);   
       
      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;

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 = 216;
y_pt = 317;

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

x = zeros(size(hgt_fin));

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

%%

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, 1200, 700])

s1 = subplot(2,1,1)

p = get(s1,'position');
p(4) = p(4)*1.2;
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
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)')

fig_nudge(gca,0,-0.05)

s2 = subplot(2,1,2)

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

contourf(x,hgt_fin,cloud_fin,cntrs)
hold on
[C2,h2] = contour(x,hgt_fin,cloud_fin,[0:0.1:0.3],'k-');
hold on
plot(time_model,zeros(length(time_model))+1800,'k--','LineWidth',2)
clabel(C2,h2);
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)')
datetick('x','mm/dd')
xlabel('Date (Ticks at 00 UTC)')

%%

close all

figure('Position',[100, 100, 1200, 300])

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])
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)')
ylabel('Height ASL (m)')

%%

