clear all
close all

run = 'nlcd2011-dec31';

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

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');
u     = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run,'/wrfout_d03_2011-01-01_00:00:00'),'U');
v     = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run,'/wrfout_d03_2011-01-01_00:00:00'),'V');
w     = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run,'/wrfout_d03_2011-01-01_00:00:00'),'W');

lat   = double(ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run,'/wrfout_d03_2011-01-01_00:00:00'),'XLAT'));
lon   = double(ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run,'/wrfout_d03_2011-01-01_00:00:00'),'XLONG'));
hgt   = double(ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run,'/wrfout_d03_2011-01-01_00:00:00'),'HGT'));

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_u = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run,'/wrfout_d03_',int2str(year(d)),'-',mon,'-',day,'_',hour,':00:00'),'U');
      
      u = cat(4,u,temp_u);
      
      temp_v = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run,'/wrfout_d03_',int2str(year(d)),'-',mon,'-',day,'_',hour,':00:00'),'V');
      
      v = cat(4,v,temp_v);
      
      temp_w = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run,'/wrfout_d03_',int2str(year(d)),'-',mon,'-',day,'_',hour,':00:00'),'W');
      
      w = cat(4,w,temp_w);
      
      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).

pot_temp = temp + 300;

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

u(:,:,:,1:6) = [];
v(:,:,:,1:6) = [];
w(:,:,:,1:6) = [];

geo_hgt = (geo + geo_p) ./ 9.8;

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

%% Get east-west cross section arrays ready.

max_lev = 41;
ew_min  = 215;
ew_max  = 235;
y_pt    = 300;

ns_min  = 290;
ns_max  = 340;
x_pt    = 225;

x_min = 210;   % Valley domain boundaries (find min height).
x_max = 240;
y_min = 290;
y_max = 340;

figure(1)
ax = usamap([lat(x_min,y_min),lat(x_max,y_max)],[lon(x_min,y_max),lon(x_max,y_min)])
contourm(lat(:,:,1),lon(:,:,1),hgt(:,:,1),[1000:100:2500],'k-')
hold on
time_height_x = plotm(lat(ew_min:ew_max,y_pt),lon(ew_min:ew_max,y_pt),'b-')
hold on
time_height_y = plotm(lat(x_pt,ns_min:ns_max),lon(x_pt,ns_min:ns_max),'r-')

hgt_avg = mean(geo_hgt,4);
temp_avg = mean(pot_temp,4);
u_avg = mean(u,4);
v_avg = mean(v,4);

slc_ew_hgt = reshape(geo_hgt(ew_min:ew_max,y_pt,1:max_lev,1),[(ew_max-ew_min+1),max_lev])'; 
slc_ns_hgt = reshape(geo_hgt(x_pt,ns_min:ns_max,1:max_lev,1),[(ns_max-ns_min+1),max_lev])';

x1 = zeros(size(slc_ew_hgt(:,:)));
x2 = zeros(size(slc_ns_hgt(:,:)));

for i = 1:max_lev
   x_ew(i,:) = x1(i,:) + [1:1:(ew_max-ew_min+1)];
   x_ns(i,:) = x2(i,:) + [1:1:(ns_max-ns_min+1)];
end

%% Average temp cross sections.

avg_ew_hgt = reshape(hgt_avg(ew_min:ew_max,y_pt,1:max_lev),[(ew_max-ew_min+1),max_lev])';
avg_ew_temp = reshape(temp_avg(ew_min:ew_max,y_pt,1:max_lev),[(ew_max-ew_min+1),max_lev])';

ew_temp_06 = squeeze(pot_temp(ew_min:ew_max,y_pt,1:max_lev,78))';
ew_temp_12 = squeeze(pot_temp(ew_min:ew_max,y_pt,1:max_lev,133))';

ew_u_06 = squeeze(u(ew_min:ew_max,y_pt,1:max_lev,78))';
ew_u_12 = squeeze(u(ew_min:ew_max,y_pt,1:max_lev,133))';

ns_temp_06 = squeeze(pot_temp(x_pt,ns_min:ns_max,1:max_lev,55))';
ns_temp_12 = squeeze(pot_temp(x_pt,ns_min:ns_max,1:max_lev,133))';

ns_v_06 = squeeze(v(x_pt,ns_min:ns_max,1:max_lev,55))';
ns_v_12 = squeeze(v(x_pt,ns_min:ns_max,1:max_lev,133))';

avg_ns_hgt = reshape(geo_hgt(x_pt,ns_min:ns_max,1:max_lev,1),[(ns_max-ns_min+1),max_lev])';
avg_ns_temp = reshape(temp_avg(x_pt,ns_min:ns_max,1:max_lev,1),[(ns_max-ns_min+1),max_lev])'-273.14;

%%
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 800])

s1 = subplot(2,1,1)

contourf(x_ew,avg_ew_hgt,ew_u_06,[-3:0.25:3],'LineStyle','none');
colormap(redblue);
caxis([-3 3]);
hold on;
contour(x_ew,avg_ew_hgt,ew_temp_06,cntrs,'k-','LineWidth',1);
hold on;
[C, h] = contour(x_ew,avg_ew_hgt,ew_temp_06,[250:5:305],'k-','LineWidth',2);
clabel(C,h);
hold on;
plot(x_ew(1,:),avg_ew_hgt(1,:),'k-');
axis([1 (ew_max-ew_min) min(min(slc_ew_hgt(:))) 1600]);
set(gca,'XTick',[0:5:50])
set(gca,'XTickLabel','')
set(gca,'YTick',[1300:50:1600]);


fl = fig_label('a','ul',16,0.015,0.04,0.6,25);
set(fl,'color','k');


%%

close all

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

s1 = subplot(2,1,2)

contourf(x_ew,avg_ew_hgt,ew_u_06,[-3:0.25:3],'LineStyle','none');
colormap(redblue);
caxis([-3 3]);
hold on;
contour(x_ew,avg_ew_hgt,ew_temp_06,cntrs,'k-','LineWidth',1);
hold on;
[C, h] = contour(x_ew,avg_ew_hgt,ew_temp_06,[250:5:305],'k-','LineWidth',2);
clabel(C,h);
hold on;
plot(x_ew(1,:),avg_ew_hgt(1,:),'k-');
axis([1 (ew_max-ew_min) min(min(slc_ew_hgt(:))) 1600]);
set(gca,'XTick',[0:5:50])
set(gca,'XTickLabel',[0:5:50].*1300)
set(gca,'YTick',[1300:50:1600]);

xlabel('Distance (m)')
ylabel('Height ASL (m)')

fl = fig_label('b','ul',16,0.015,0.04,0.6,25);
set(fl,'color','k');

c = colorbar;
hL = title(c,'m s^-^1','FontSize',10);
set(c,'position',[0.92 0.13 0.015 0.75]) 

%%
close all

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

s1 = subplot(2,2,2)

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

contourf(x_ns,avg_ns_hgt,ns_v_06,[-10:1:10],'LineStyle','none')
colormap(redblue)
caxis([-10 10])
hold on;
contour(x_ns,avg_ns_hgt,ns_temp_06,cntrs,'k-','LineWidth',1);
hold on
[C, h] = contour(x_ns,avg_ns_hgt,ns_temp_06,[250:5:305],'k-','LineWidth',2)
clabel(C,h)
hold on
plot(x_ns(1,:),avg_ns_hgt(1,:),'k-')
axis([1 (ns_max-ns_min+1) min(min(slc_ns_hgt(:))) 2400]);
set(gca,'XTickLabel','')
set(gca,'YTick',[1400:200:2200])
set(gca,'YTickLabel','')

fl = fig_label('b','ul',16,0.015,0.04,3,90);
set(fl,'color','k');

fig_nudge(gca,0,-0.04)

s2 = subplot(2,2,4)

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

contourf(x_ns,avg_ns_hgt,ns_v_12,[-10:1:10],'LineStyle','none')
colormap(redblue)
caxis([-10 10])
hold on;
contour(x_ns,avg_ns_hgt,ns_temp_12,cntrs,'k-','LineWidth',1);
hold on
[C, h] = contour(x_ns,avg_ns_hgt,ns_temp_12,[250:5:305],'k-','LineWidth',2)
clabel(C,h)
hold on
plot(x_ns(1,:),avg_ns_hgt(1,:),'k-')
axis([1 (ns_max-ns_min+1) min(min(slc_ns_hgt(:))) 2400]);
set(gca,'YTickLabel','')
set(gca,'XTick',[0:10:50])
set(gca,'XTickLabel',[0:10:50].*1300)
xlabel('Distance (m)')

fl = fig_label('d','ul',16,0.015,0.04,3,90);
set(fl,'color','k');

c = colorbar;
hL = title(c,'m s^-^1','FontSize',10);
set(c,'position',[0.92 0.13 0.015 0.8425]) 

%%

close all

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

s1 = subplot(2,2,1)

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

contourf(x_ns,avg_ns_hgt,ns_v_06,[-10:1:10],'LineStyle','none')
colormap(redblue)
caxis([-10 10])
hold on;
contour(x_ns,avg_ns_hgt,ns_temp_06,cntrs,'k-','LineWidth',1);
hold on
[C, h] = contour(x_ns,avg_ns_hgt,ns_temp_06,[250:5:305],'k-','LineWidth',2)
clabel(C,h)
hold on
plot(x_ns(1,:),avg_ns_hgt(1,:),'k-')
axis([1 (ns_max-ns_min+1) min(min(slc_ns_hgt(:))) 2400]);
set(gca,'XTickLabel','')
set(gca,'YTick',[1400:200:2200])
ylabel('Height ASL (m)')

fl = fig_label('a','ul',16,0.015,0.04,3,90);
set(fl,'color','k');

fig_nudge(gca,0,-0.04)

s2 = subplot(2,2,3)

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

contourf(x_ns,avg_ns_hgt,ns_v_12,[-10:1:10],'LineStyle','none')
colormap(redblue)
caxis([-10 10])
hold on;
contour(x_ns,avg_ns_hgt,ns_temp_12,cntrs,'k-','LineWidth',1);
hold on
[C, h] = contour(x_ns,avg_ns_hgt,ns_temp_12,[250:5:305],'k-','LineWidth',2)
clabel(C,h)
hold on
plot(x_ns(1,:),avg_ns_hgt(1,:),'k-')
axis([1 (ns_max-ns_min+1) min(min(slc_ns_hgt(:))) 2400]);
set(gca,'YTick',[1400:200:2200])
set(gca,'XTick',[0:10:50])
set(gca,'XTickLabel',[0:10:50].*1300)
xlabel('Distance (m)')
ylabel('Height ASL (m)')

fl = fig_label('c','ul',16,0.015,0.04,3,90);
set(fl,'color','k');


%% Make a GIF.

figure('Position', [100, 100, 2000, 1000]);

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

%%

set(0,'defaultfigureposition',[100 100 1500 1000])

for i = 1:9
   
   figure(i)
   set(i,'defaultfigureposition',[100 100 1500 1000])
   
   slc_ns_v  = reshape(v(x_pt,ns_min:ns_max,1:max_lev,i),[51,max_lev])';
   slc_ns_temp  = reshape(pot_temp(x_pt,ns_min:ns_max,1:max_lev,i),[51,max_lev])';
   
   contourf(x_ns,slc_ns_hgt(:,:),slc_ns_v(:,:),[-8:0.5:8],'edgecolor','none');
   contourcmap('jet',[-10:0.5:10])
   colorbar;
   hold on
   contour(x_ns,slc_ns_hgt(:,:),slc_ns_temp(:,:),cntrs,'k-')
   hold on
   [C,h] = contour(x_ns,slc_ns_hgt(:,:),slc_ns_temp(:,:),[250:5:305],'k-','LineWidth',2)
   clabel(C,h)
   axis([1 (ns_max-ns_min+1) min(min(slc_ns_hgt(:))) 2900]);
   title1 = {'Model Potential Temperature (Contours) and v Wind Component (Shading),',datestr(time_model(i))}
   title(title1)
   
   saveas(i,strcat('/uufs/chpc.utah.edu/common/home/u0875930/public_html/paper_figs/gif_nlcd/image00',int2str(i)),'png')
   
   close(i)
end

for i = 10:99
   
   figure(i)
   set(i,'defaultfigureposition',[100 100 1500 1000])
   
   slc_ns_v  = reshape(v(x_pt,ns_min:ns_max,1:max_lev,i),[51,max_lev])';
   slc_ns_temp  = reshape(pot_temp(x_pt,ns_min:ns_max,1:max_lev,i),[51,max_lev])';
   
   contourf(x_ns,slc_ns_hgt(:,:),slc_ns_v(:,:),[-10:0.5:10],'edgecolor','none');
   contourcmap('jet',[-10:0.5:10])
   colorbar;
   hold on
   contour(x_ns,slc_ns_hgt(:,:),slc_ns_temp(:,:),cntrs,'k-')
   hold on
   [C,h] = contour(x_ns,slc_ns_hgt(:,:),slc_ns_temp(:,:),[250:5:305],'k-','LineWidth',2)
   clabel(C,h)
   axis([1 (ns_max-ns_min+1) min(min(slc_ns_hgt(:))) 2900]);
   title1 = {'Model Potential Temperature (Contours) and v Wind Component (Shading),',datestr(time_model(i))}
   title(title1)
   
   saveas(i,strcat('/uufs/chpc.utah.edu/common/home/u0875930/public_html/paper_figs/gif_nlcd/image0',int2str(i)),'png')
   
   close(i)
end

for i = 100:length(time_model)
   
   figure(i)
   set(i,'defaultfigureposition',[100 100 1500 1000])
   
   slc_ns_v  = reshape(v(x_pt,ns_min:ns_max,1:max_lev,i),[51,max_lev])';
   slc_ns_temp  = reshape(pot_temp(x_pt,ns_min:ns_max,1:max_lev,i),[51,max_lev])';
   
   contourf(x_ns,slc_ns_hgt(:,:),slc_ns_v(:,:),[-10:0.5:10],'edgecolor','none');
   contourcmap('jet',[-10:0.5:10])
   colorbar;
   hold on
   contour(x_ns,slc_ns_hgt(:,:),slc_ns_temp(:,:),cntrs,'k-')
   hold on
   [C,h] = contour(x_ns,slc_ns_hgt(:,:),slc_ns_temp(:,:),[250:5:305],'k-','LineWidth',2)
   clabel(C,h)
   axis([1 (ns_max-ns_min+1) min(min(slc_ns_hgt(:))) 2900]);
   title1 = {'Model Potential Temperature (Contours) and v Wind Component (Shading),',datestr(time_model(i))}
   title(title1)
   
   saveas(i,strcat('/uufs/chpc.utah.edu/common/home/u0875930/public_html/paper_figs/gif_nlcd/image',int2str(i)),'png')
   
   close(i)
end

%%

set(0,'defaultfigureposition',[100 100 1500 1000])

for i = 1:9
   
   figure(i)
   set(i,'defaultfigureposition',[100 100 1500 1000])
   
   slc_ew_u  = reshape(u(ew_min:ew_max,y_pt,1:max_lev,i),[21,max_lev])';
   slc_ew_temp  = reshape(pot_temp(ew_min:ew_max,y_pt,1:max_lev,i),[21,max_lev])';
   
   contourf(x_ew,slc_ew_hgt(:,:),slc_ew_u(:,:),[-8:0.5:8],'edgecolor','none');
   colormap(redblue)
   caxis([-10 10])
   colorbar;
   hold on
   contour(x_ew,slc_ew_hgt(:,:),slc_ew_temp(:,:),cntrs,'k-')
   hold on
   [C,h] = contour(x_ew,slc_ew_hgt(:,:),slc_ew_temp(:,:),[250:5:305],'k-','LineWidth',2)
   clabel(C,h)
   axis([1 (ew_max-ew_min+1) min(min(slc_ew_hgt(:))) 2000]);
   title1 = {'Model Potential Temperature (Contours) and u Wind Component (Shading),',datestr(time_model(i))}
   title(title1)
   
   saveas(i,strcat('/uufs/chpc.utah.edu/common/home/u0875930/public_html/paper_figs/gif_nlcd_ew/image00',int2str(i)),'png')
   
   close(i)
end

for i = 10:99
   
   figure(i)
   set(i,'defaultfigureposition',[100 100 1500 1000])
   
   slc_ew_u  = reshape(u(ew_min:ew_max,y_pt,1:max_lev,i),[21,max_lev])';
   slc_ew_temp  = reshape(pot_temp(ew_min:ew_max,y_pt,1:max_lev,i),[21,max_lev])';
   
   contourf(x_ew,slc_ew_hgt(:,:),slc_ew_u(:,:),[-8:0.5:8],'edgecolor','none');
   colormap(redblue)
   caxis([-10 10])
   colorbar;
   hold on
   contour(x_ew,slc_ew_hgt(:,:),slc_ew_temp(:,:),cntrs,'k-')
   hold on
   [C,h] = contour(x_ew,slc_ew_hgt(:,:),slc_ew_temp(:,:),[250:5:305],'k-','LineWidth',2)
   clabel(C,h)
   axis([1 (ew_max-ew_min+1) min(min(slc_ew_hgt(:))) 2000]);
   title1 = {'Model Potential Temperature (Contours) and u Wind Component (Shading),',datestr(time_model(i))}
   title(title1)
   
   saveas(i,strcat('/uufs/chpc.utah.edu/common/home/u0875930/public_html/paper_figs/gif_nlcd_ew/image0',int2str(i)),'png')
   
   close(i)
end

for i = 100:length(time_model)
   
   figure(i)
   set(i,'defaultfigureposition',[100 100 1500 1000])
   
   slc_ew_u  = reshape(u(ew_min:ew_max,y_pt,1:max_lev,i),[21,max_lev])';
   slc_ew_temp  = reshape(pot_temp(ew_min:ew_max,y_pt,1:max_lev,i),[21,max_lev])';
   
   contourf(x_ew,slc_ew_hgt(:,:),slc_ew_u(:,:),[-8:0.5:8],'edgecolor','none');
   colormap(redblue)
   caxis([-10 10])
   colorbar;
   hold on
   contour(x_ew,slc_ew_hgt(:,:),slc_ew_temp(:,:),cntrs,'k-')
   hold on
   [C,h] = contour(x_ew,slc_ew_hgt(:,:),slc_ew_temp(:,:),[250:5:305],'k-','LineWidth',2)
   clabel(C,h)
   axis([1 (ew_max-ew_min+1) min(min(slc_ew_hgt(:))) 2000]);
   title1 = {'Model Potential Temperature (Contours) and u Wind Component (Shading),',datestr(time_model(i))}
   title(title1)
   
   saveas(i,strcat('/uufs/chpc.utah.edu/common/home/u0875930/public_html/paper_figs/gif_nlcd_ew/image',int2str(i)),'png')
   
   close(i)
end