% script to compare heat flux as a function of time between simulations

close all
clear all

scrz=get(0,'screensize');

%% add a bunch of paths to make sure netcdf handling is easier and to point at relevant scripts/data
addpath('/uufs/chpc.utah.edu/common/home/u0559900/MATLAB_SCRIPTS/'); %we call some scripts from here
addpath('/uufs/chpc.utah.edu/common/home/u0559900/MATLAB_SCRIPTS/mexcdf/mexnc/') % netcdf
addpath('/uufs/chpc.utah.edu/common/home/u0559900/MATLAB_SCRIPTS/mexcdf/snctools'); %netcdf

%% catalog of runs/names/locations
%inflow controls runs
runs(1).paths='/uufs/chpc.utah.edu/common/home/horel-group2/nlareau/les_output/production_runs/control_runs/mountains/output/';
runs(1).names='inflow_mountains';
runs(2).paths='/uufs/chpc.utah.edu/common/home/horel-group2/nlareau/les_output/production_runs/control_runs/plane/output/';
runs(2).names='inflow_plane';

%linear profile
runs(3).paths='/uufs/chpc.utah.edu/common/home/horel-group2/nlareau/les_output/production_runs/linear_mountains/output/';
runs(3).names='linear_mountains';
runs(4).paths='/uufs/chpc.utah.edu/common/home/horel-group2/nlareau/les_output/production_runs/linear_plane/output/';
runs(4).names='linear_plane';

%three layer
runs(5).paths='/uufs/chpc.utah.edu/common/home/horel-group2/nlareau/les_output/production_runs/three_layer_mountains/output/';
runs(5).names='three_layer_mountains';
runs(6).paths='/uufs/chpc.utah.edu/common/home/horel-group2/nlareau/les_output/production_runs/three_layer_plane/output/';
runs(6).names='three_layer_plane';

%capping layer
runs(7).paths='/uufs/chpc.utah.edu/common/home/horel-group2/nlareau/les_output/production_runs/capping_mountains/output/';
runs(7).names='capping_mountains';
runs(8).paths='/uufs/chpc.utah.edu/common/home/horel-group2/nlareau/les_output/production_runs/capping_plane/output/';
runs(8).names='capping_plane';

%shallow
runs(9).paths='/uufs/chpc.utah.edu/common/home/horel-group2/nlareau/les_output/production_runs/shallow_mountains/output/';
runs(9).names='shallow_mountains';
runs(10).paths='/uufs/chpc.utah.edu/common/home/horel-group2/nlareau/les_output/production_runs/shallow_plane/output/';
runs(10).names='shallow_plane';

%data dir
ddir='/uufs/chpc.utah.edu/common/home/horel-group2/nlareau/matlab_files/';

labeler={'','','a) Linear','e) Linear','b) Multi-Layer','f) Multi-Layer','c) Capping','g) Capping','d)Shallow','h) Shallow'};


%% initialize deficity matrix and x grid

deficit=zeros(10,20*60);
surfdef=zeros(10,20*60);
hovmol=zeros(10,20*60,2039);
dx=50;%50 m grid spacing
xvec=0:dx:(2038*50);
xmat=repmat(xvec,[50,1,1200]);xmat=permute(xmat,[3 1 2]);


%% set up z grid to interp to
zvec=10:10:600;
zmat=repmat(zvec',[1 2039]);
zmat3=repmat(zvec', [1 2039 1200]); zmat3=permute(zmat3,[3 1 2]);

xnew=repmat(xvec,[length(zvec) 1]);
xmat3=repmat(xvec,[length(zvec) 1 1200]); xmat3=permute(xmat3,[3 1 2]);
%% load terrain

hf='/uufs/chpc.utah.edu/common/home/horel-group2/nlareau/les_output/production_runs/linear_mountains/output/wrfout_d01_0001-01-01_00:00:00';
mountain=double(nc_varget(hf,'HGT'));
mountain=squeeze(mountain(1,1,:));
hf='/uufs/chpc.utah.edu/common/home/horel-group2/nlareau/les_output/production_runs/linear_plane/output/wrfout_d01_0001-01-01_00:00:00';
plane=double(nc_varget(hf,'HGT'));
plane=squeeze(plane(1,1,:));

%% generate ground mask to blank extrapolated data
mmask=zeros(size(zmat));
for xx=1:length(xvec);
    hidx=find(zmat(:,xx)<mountain(xx));
    mmask(hidx,xx)=1;
end

pmask=zeros(size(zmat));
for xx=1:length(xvec);
    hidx=find(zmat(:,xx)<plane(xx));
    pmask(hidx,xx)=1;
end

%% Loop over each model run

HFX=zeros(10,75);
TKE=zeros(10,75);
%uideal=linspace(.5,25,length(tvec));

runner=[3 4 5 6 7 8 9 10];
mt1=900;%600
mt2=1100;%1440
cdef=zeros(10,1200,2039);
sdef=zeros(10,1200,2039);%surface def
Fr=zeros(10,1200,2039);
Frmean=zeros(10,1200);
umean=zeros(10,1200);
Bdef=zeros(10,1200);
for mm=runner;
    %load the .mat files
    namer=strcat('Zgrid_',runs(mm).names,'.mat');
    load(namer);
    %%
    hidx=find(zvec==500);
    tref=t(:,hidx,1);
    trefmat=repmat(tref,[1 size(t,2) size(t,3)]);
    def=trefmat-t;% temperature deficit at each level
    def(def<0)=0;
    def(:,:,1:600)=0;
    def(:,:,1440:end)=0;
    def(isnan(def))=0;
    sdef(mm,:,:)=squeeze(nanmean(def(:,1:10,:),2));
    cdef(mm,:,:)=squeeze(trapz(10*def(:,1:hidx,:),2));% integrate the deficit over the column
    gwave=sqrt((9.81./tref(1)).*squeeze(cdef(mm,:,:)));
    umax=squeeze(max(u(:,1:hidx+2,:),[],2));
    ubar=squeeze(nanmean(u(:,1:hidx,:),2));
    [shear,~,~]=gradient(u,10);
    
    %%
    [maxshear,Shidx]=(max(shear(:,1:hidx,:),[],2));
    maxshear=squeeze(maxshear);
    Shidx=squeeze(Shidx);
    MMshear=nanmean(maxshear(:,mt1:mt2),2);
    MMShidx=nanmean(Shidx(:,mt1:mt2),2);
    
    Nsq=(9.81./tref(1)).*squeeze(cdef(mm,:,:))./500;
    Ri=Nsq./(maxshear.^2);
    Fr=umax./gwave;
    Fr(isinf(Fr))=nan;
    %Frmean(mm,:)=nanmin(Fr(:,mt1:mt2),[],2);
    Frmean(mm,:)=nanmean(Fr(:,mt1:mt2),2);
    umean(mm,:)=nanmean(umax(:,mt1:mt2),2);
    Bdef(mm,:)=squeeze(trapz(50.*gwave(:,mt1:mt2),2));
    
    %%
    if mm==3;
        scrz=get(0,'screensize');
        f=figure(100+mm);clf;
        set(f,'renderer','painters','color','w','pos',[scrz(1) scrz(2) scrz(3)*.8 scrz(4).*5]);
        xvec2=xvec./10^3;
        subplot(3,3,[3 6 9]);
        colormap(flipud(gray(20)));
        %pcolor(xvec2,tvec,squeeze(cdef(3,:,:)));shading flat;
        pcolor(xvec2,tvec,gwave(:,:).^2);shading flat;
        caxis([0 100])
        %pcolor(xvec2,tvec,squeeze(sdef(3,:,:)));shading flat;
        xlim([25 77]);
        %subplot(3,3,[3 6 9]);
        hold on;
        uideal=linspace(.5,25,1200);
        [~,tidx]=min(abs(uideal-5));%find(tvec==5);
        plot(28,tvec(tidx),'>','markersize',8,'markerfacecolor','g');
        [~,tidx]=min(abs(uideal-10));%find(tvec==5);
        plot(28,tvec(tidx),'>','markersize',8,'markerfacecolor','g');
        [~,tidx]=min(abs(uideal-15));%find(tvec==5);
        plot(28,tvec(tidx),'>','markersize',8,'markerfacecolor','g');
        set(gca,'ydir','rev','box','on','layer','top','yaxislocation','right');
        freezeColors;
        cbh=colorbar('location','south');
        %cpos=get(cbh,'pos');
        %ylabel(cbh,'K m');
        xlabel('Dist. [km]','fontsize',10,'fontweight','bold');
        %set(cbh,'pos',[cpos(1)+.0751 cpos(2)+.25 cpos(3) .4.*cpos(4)]);
        cbfreeze(cbh);
        set(gca,'fontsize',10,'fontweight','bold')
        ylabel('Time [hrs]','fontsize',10,'fontweight','bold')
        text(27,1,'d)','fontsize',10,'fontweight','bold');

        subplot(3,3,1:2);
        load('goodcmap');
        colormap(cmap);
        [~,tidx]=min(abs(uideal-5));%find(tvec==5);
        pcolor(xvec2,zvec,squeeze(mean(u((tidx-10:tidx+10),:,:),1)));
        shading flat;
        caxis([-5 25]);
        hold on;
        contour(xvec2,zvec,squeeze(mean(t((tidx-10:tidx+10),:,:),1)),0:1:500,'-k');
        ha=area(xvec2,mountain);
        set(ha,'facecolor',[.7 .7 .7],'linewidth',2);
        ylim([0 600]);
        xlim([20 82]);
        box on;
        set(gca,'layer','top','xticklabe',[]);
        ylabel('Height [m]','fontsize',10,'fontweight','bold');
        set(gca,'fontsize',10,'fontweight','bold')
        text(21,520,'a)','fontsize',10,'fontweight','bold');
        
        
        
        subplot(3,3,4:5);
        [~,tidx]=min(abs(uideal-10));%find(tvec==5);
        pcolor(xvec2,zvec,squeeze(mean(u((tidx-10:tidx+10),:,:),1)));
        shading flat;
        caxis([-5 25]);
        hold on;
        contour(xvec2,zvec,squeeze(mean(t((tidx-10:tidx+10),:,:),1)),0:1:500,'-k');
        ha=area(xvec2,mountain);
        set(ha,'facecolor',[.7 .7 .7],'linewidth',2);
        ylim([0 600]);
        xlim([20 82]);
        box on;
        set(gca,'layer','top','xticklabe',[]);
        ylabel('Height [m]','fontsize',10,'fontweight','bold');
        set(gca,'fontsize',10,'fontweight','bold')
        cbh=colorbar;
        cpos=get(cbh,'pos');
        set(cbh,'pos',[cpos(1)+.06 cpos(2)+.0001 cpos(3)*.5 cpos(4)]);
        ylabel(cbh,'m s^{-1}');
        text(21,520,'b)','fontsize',10,'fontweight','bold');
        
        subplot(3,3,7:8);
        [~,tidx]=min(abs(uideal-15));
        pcolor(xvec2,zvec,squeeze(mean(u((tidx-10:tidx+10),:,:),1)));
        shading flat;
        caxis([-5 25]);
        
        hold on;
        [xmat,zmat]=meshgrid(xvec2,zvec);
        unow=squeeze(mean(u((tidx-10:tidx+10),:,:),1));
        wnow=squeeze(mean(w((tidx-10:tidx+10),:,:),1));
        spr=2;
        spr2=30;
        %qm=quiver(xmat(1:spr:end,1:spr2:end),zmat(1:spr:end,1:spr2:end),.025*unow(1:spr:end,1:spr2:end),wnow(1:spr:end,1:spr2:end),.3);
        %set(qm,'linewidth',2)
        contour(xvec2,zvec,squeeze(mean(t((tidx-10:tidx+10),:,:),1)),0:1:500,'-k');
        ha=area(xvec2,mountain);
        set(ha,'facecolor',[.7 .7 .7],'linewidth',2);
        ylim([0 600]);
        xlim([20 82]);
        box on;
        set(gca,'layer','top');
        set(gca,'fontsize',10,'fontweight','bold')
        
        xlabel('Dist. [km]','fontsize',10,'fontweight','bold');
        ylabel('Height [m]','fontsize',10,'fontweight','bold');
        text(21,520,'c)','fontsize',10,'fontweight','bold');
    end
    %%
    %#### CAP BH
    scrz=get(0,'screensize');
    if mm==3
        panels=[ 1 5 2 6 3 7 4 8];
        f=figure(13);clf;
        set(f,'renderer','painters','color','w','pos',[scrz(1) scrz(2) scrz(3)*.6 scrz(4)]);
        ha13=subplot1(2, 4, 'Gap', [.01 .03], 'XTickL', 'None', 'YTickL', 'None');
        
    end
    figure(13);
    hold on;
    %subplot1(panels(mm-2))
    axes(ha13((mm-2)))
    BH=squeeze(gwave.^2);
    BH(isnan(BH))=5;
    colormap([(jet(20));(pink(5))]);
    
    pcolor(xvec2,tvec,BH);shading flat;caxis([0 140]);
    hold on;
    set(gca,'ydir','rev');
    xlim([22 78]);
    box on;
    if (mm-2)==8;
        cbh=colorbar;
        cpos=get(cbh,'pos');
        cpos(1)=cpos(1)+.065;
        cpos(2)=cpos(2)+.1;
        cpos(4)=.5.*cpos(4);
        set(cbh,'pos',cpos);
    end
    set(gca,'layer','top','linewidth',2,'fontsize',10,'fontweight','bold');
    if ismember(mm,[4 6 8 10])
        
        xlabel('Dist [km]','fontsize',10,'fontweight','bold');
    else
        set(gca,'xticklabel',[]);
    end
    if ismember(mm,[3 4])
        ylabel('Time [hrs]','fontsize',10,'fontweight','bold');
    else
        set(gca,'yticklabel',[]);
    end
    th=text(24,1,labeler(mm),'fontsize',10,'fontweight','bold','color','w');
    
    

    %%
    %#### CAP FROUDE NUMBER
    scrz=get(0,'screensize');
    if mm==3
        panels=[ 1 5 2 6 3 7 4 8];
        f=figure(14);clf;set(f,'renderer','painters','color','w','pos',[scrz(1) scrz(2) scrz(3)*.6 scrz(4)]);
        ha=subplot1(2, 4, 'Gap', [.01 .03], 'XTickL', 'None', 'YTickL', 'None');
        
    end
    figure(14);
    %subplot1(panels(mm-2))
    axes(ha((mm-2)))
    FRnow=squeeze(Fr);
    FRnow(isnan(FRnow))=5;
    colormap((flipud(jet(20))));
    
    pcolor(xvec2,tvec,FRnow);shading flat;caxis([0 5]);
    hold on;
    set(gca,'ydir','rev');
    xlim([22 78]);
    box on;
    if (mm-2)==8;
        cbh=colorbar;
        cpos=get(cbh,'pos');
        cpos(1)=cpos(1)+.065;
        cpos(2)=cpos(2)+.1;
        cpos(4)=.5.*cpos(4);
        set(cbh,'pos',cpos);
    end
    set(gca,'layer','top','linewidth',2,'fontsize',10,'fontweight','bold');
    if ismember(mm,[4 6 8 10])
        
        xlabel('Dist [km]','fontsize',10,'fontweight','bold');
    else
        set(gca,'xticklabel',[]);
    end
    if ismember(mm,[3 4])
        ylabel('Time [hrs]','fontsize',10,'fontweight','bold');
    else
        set(gca,'yticklabel',[]);
    end
    th=text(24,1,labeler(mm),'fontsize',10,'fontweight','bold','color','w');
    
end


return
%% close all
clrs={'r','--r','g','--g','b','--b','k','--k','m','--m','c','--c'};
close all
f1=figure(1);clf; set(f1,'color','w','pos',.8.*scrz);
a=0;
hold on;
Bdef(Bdef<.4*10^5)=nan;
for mm=[3:8]%runner
    a=a+1;
    plot(tvec,Frmean(mm,:),clrs{a});
    ylim([0 5]);
end
lh=legend(runs(runner).names);
set(lh,'interpreter','none');
box on;
grid on;
set(gca,'layer','top')
return

%% rate of buoyancy reduction


[dbdt,~]=gradient(Bdef,60);
f1=figure(2);clf; set(f1,'color','w','pos',.8.*scrz);
a=0;
hold on;
Bdef(Bdef<.5*10^5)=nan;
Bdef(isnan(Bdef))=.5*10^5;

[bb,cc]=butter(5,1/25);
for mm=runner
    a=a+1;
    bsmooth=filtfilt(bb,cc,Bdef(mm,:));
    subplot(2,1,1);hold on;
    plot(tvec,bsmooth,clrs{a});
    subplot(2,1,2);hold on;
    dbdt=gradient(bsmooth(:));
    plot(tvec,dbdt,clrs{a});
    ylim([-1200 100])
end
lh=legend(runs(runner).names);
set(lh,'interpreter','none');
box on;
grid on;
set(gca,'layer','top')


%%
clrs={'r','--r','g','--g','b','--b','k','--k','m','--m','c','--c'};
f1=figure(2);clf; set(f1,'color','w','pos',.8.*scrz);
a=0;
hold on;
for mm=runner
    a=a+1;
    plot(umean(mm,1:5:end),Bdef(mm,1:5:end),clrs{a},'linewidth',2);
    %ylim([0 3]);
end
lh=legend(runs(runner).names);
set(lh,'interpreter','none');
box on;
grid on;
set(gca,'layer','top')

%%
clrs={'r','-.r','b','-.b','k','-.k'};
f1=figure(3);clf; set(f1,'color','w','pos',.8.*scrz);
a=0;
hold on;
umean(umean==0)=nan;
for mm=runner
    a=a+1;
    plot(tvec,(Bdef(mm,:)),clrs{a},'linewidth',2);
end
lh=legend(runs(runner).names);
set(lh,'interpreter','none');
box on;
grid on;
set(gca,'layer','top')