%{
file=input('Enter subdirectory name (default = output)','s');
if isempty(file)
   file='output';
end
%}
file = 'output';

load(strcat(file,filesep,'time.out'))
load(strcat(file,filesep,'profile.out'))
load(strcat(file,filesep,'cldhov.out'))
load(strcat(file,filesep,'thov.out'))
load(strcat(file,filesep,'qhov.out'))
load(strcat(file,filesep,'rhhov.out'))
load(strcat(file,filesep,'mhov.out'))
load(strcat(file,filesep,'mdhov.out'))
load(strcat(file,filesep,'mphov.out'))
load(strcat(file,filesep,'buoy.out'))
p=-profile(:,1);
time1=time(:,1);
timem=time(:,12);
tout=mean(time(:,13));
tk=profile(:,2)+273.15;
qm=0.001.*profile(:,3);
tv=tk.*(1+0.608.*qm);
aa=size(p);
a=aa(1,1);
z(1)=0.0;
for i=2:a,
    z(i)=z(i-1)+(287./9.81).*0.5.*(tv(i)+tv(i-1)).*log(p(i-1)./p(i));
end
zkm=0.001.*z;
i=1;
while i == 1
    clc
    dispstr = 'Menu choices: \n \n Evolution: \n   1) Precipitation and Evaporation \n   2) Surface air and sea surface temperatures \n   3) 500 mb T and q \n   4) TOA shortwave and longwave \n   5) Potential intensity   \n  \n  Vertical Profiles:  \n   6) Temperature \n   7) Specific humidity \n   8) Buoyancy \n   9) Relative humidity \n  10) Heating rates  \n  11) Mass fluxes \n  12) Entrainment and detrainment  \n  13) Moist static energy  \n  14) Cloud fraction  \n  15) Cloud condensed water  \n \n Time-height sections: \n  16) Buoyancy \n  17) Cloud fraction \n  18) T diff \n  19) Q diff \n  20) RH diff \n  21) Upward mass flux \n  22) Penetrative mass flux \n  23) Unsaturated downdraft mass flux \n  24) Net mass flux     \n \n  0) exit \n ';
    fprintf(1,dispstr);
    it=input('Your choice? ');
    if isempty(it)
        
        pause(0.05); % This is a total hack to suppress some weird behavior where it interprets blank input as 1.
    elseif ~isempty(it)
        
        figure(1)
        if it == 1
            precip=time(:,2);
            evap=time(:,3);
            h=plot(timem,precip,'-',timem,evap,'--');
            set(h,'linewidth',2)
            set(gca,'fontweight','bold')
            xlabel('Time (days)','fontweight','bold')
            ylabel('Precipitation and Evaporation (mm/day)','fontweight','bold')
            legend('Precipitation','Evaporation',0)
        elseif it == 2
            sst=time(:,5);
            ta=time(:,4);
            h=plot(time1,sst,'-',time1,ta,'--');
            set(gca,'fontweight','bold')
            set(h,'linewidth',2)
            xlabel('Time (days)','fontweight','bold')
            ylabel('SST and Ta (C)','fontweight','bold')
            legend('SST','Air temperature',0)
        elseif it == 3
            t500=time(:,6);
            q500=time(:,7);
            h=plot(time1,t500,'-',time1,q500,'--');
            set(h,'linewidth',2)
            set(gca,'fontweight','bold')
            xlabel('Time (days)','fontweight','bold')
            ylabel('500 mb T (C) and q (g/Kg)','fontweight','bold')
            legend('T','q',0)
        elseif it == 4
            sw=time(:,8);
            lw=time(:,9);
            h=plot(time1,sw,'-',time1,lw,'--');
            set(h,'linewidth',2)
            set(gca,'fontweight','bold')
            xlabel('Time (days)','fontweight','bold')
            ylabel('TOA Shortwave and Longwave (W/m^2)','fontweight','bold')
            legend('SW','LW',0)
        elseif it == 5
            vmax=time(:,10);
            h=plot(time1,vmax);
            set(h,'linewidth',2)
            set(gca,'fontweight','bold')
            xlabel('Time (days)','fontweight','bold')
            ylabel('V_{max} (m/s)','fontweight','bold')
        elseif it == 51
            pmin=time(:,11);
            h=plot(time1,pmin);
            set(h,'linewidth',2)
            set(gca,'fontweight','bold')
            xlabel('Time (days)','fontweight','bold')
            ylabel('p_{min} (mb)','fontweight','bold')
        elseif it == 6
            tpro=profile(:,2);
            h=plot(tpro,p);
            set(h,'linewidth',2)
            set(gca,'fontweight','bold')
            u=get(gca,'YTick');
            u1=transpose(-u);
            v=cellstr(num2str(u1));
            set(gca, 'YTickLabel',v)
            xlabel('Temperature (C)','fontweight','bold')
            ylabel('Pressure (mb)','fontweight','bold')
        elseif it == 7
            qpro=profile(:,3);
            h=plot(qpro,p);
            set(h,'linewidth',2)
            set(gca,'fontweight','bold')
            u=get(gca,'YTick');
            u1=transpose(-u);
            v=cellstr(num2str(u1));
            set(gca, 'YTickLabel',v)
            xlabel('Specific humidity (gm/Kg)','fontweight','bold')
            ylabel('Pressure (mb)','fontweight','bold')
        elseif it == 8
            buoy1=profile(:,4);
            h=plot(buoy1,p);
            set(h,'linewidth',2)
            set(gca,'fontweight','bold')
            u=get(gca,'YTick');
            u1=transpose(-u);
            v=cellstr(num2str(u1));
            set(gca, 'YTickLabel',v)
            xlabel('Buoyancy (C)','fontweight','bold')
            ylabel('Pressure (mb)','fontweight','bold')
        elseif it ==9
            rh=100.*profile(:,5);
            h=plot(rh,p);
            set(h,'linewidth',2)
            set(gca,'fontweight','bold')
            u=get(gca,'YTick');
            u1=transpose(-u);
            v=cellstr(num2str(u1));
            set(gca, 'YTickLabel',v)
            xlabel('Relative humidity (%)','fontweight','bold')
            ylabel('Pressure (mb)','fontweight','bold')
        elseif it ==10
            cq=profile(:,6);
            rq=profile(:,7);
            tq=profile(:,8);
            aq=profile(:,9);
            h=plot(cq,p,rq,p,tq,p,aq,p);
            set(h,'linewidth',2)
            set(gca,'fontweight','bold')
            u=get(gca,'YTick');
            u1=transpose(-u);
            v=cellstr(num2str(u1));
            set(gca, 'YTickLabel',v)
            xlabel('Heating (C/day)','fontweight','bold')
            ylabel('Pressure (mb)','fontweight','bold')
            legend('Convective','Radiative','Turbulent', 'Large-scale adiabatic',0)
        elseif it == 11
            mu=profile(:,10);
            mpen=profile(:,11);
            md=profile(:,12);
            h=plot(mu,p,mpen,p,md,p);
            set(h,'linewidth',2)
            set(gca,'fontweight','bold')
            u=get(gca,'YTick');
            u1=transpose(-u);
            v=cellstr(num2str(u1));
            set(gca, 'YTickLabel',v)
            xlabel('Mass flux (10^{-3} Kg m^{-2} s^{-1})','fontweight','bold')
            ylabel('Pressure (mb)','fontweight','bold')
            legend('Net upward','Penetrative','Unsaturated',0)
        elseif it == 12
            ent=profile(:,13);
            det=profile(:,14);
            h=plot(ent,p,det,p);
            set(h,'linewidth',2)
            set(gca,'fontweight','bold')
            u=get(gca,'YTick');
            u1=transpose(-u);
            v=cellstr(num2str(u1));
            set(gca, 'YTickLabel',v)
            xlabel('Entrainment and detrainment (10^{-3} Kg m^{-2} s^{-1})','fontweight','bold')
            ylabel('Pressure (mb)','fontweight','bold')
            legend('Entrainment','Detrainment',0)
        elseif it == 13
            mse=profile(:,15);
            h=plot(mse,p);
            set(h,'linewidth',2)
            set(gca,'fontweight','bold')
            u=get(gca,'YTick');
            u1=transpose(-u);
            v=cellstr(num2str(u1));
            set(gca, 'YTickLabel',v)
            xlabel('Moist static energy (K)','fontweight','bold')
            ylabel('Pressure (mb)','fontweight','bold')
        elseif it == 14
            clf=profile(:,16);
            h=plot(clf,p);
            set(h,'linewidth',2)
            set(gca,'fontweight','bold')
            u=get(gca,'YTick');
            u1=transpose(-u);
            v=cellstr(num2str(u1));
            set(gca, 'YTickLabel',v)
            xlabel('Cloud fraction','fontweight','bold')
            ylabel('Pressure (mb)','fontweight','bold')
        elseif it == 15
            clq=profile(:,17);
            h=plot(clq,p);
            set(h,'linewidth',2)
            set(gca,'fontweight','bold')
            u=get(gca,'YTick');
            u1=transpose(-u);
            v=cellstr(num2str(u1));
            set(gca, 'YTickLabel',v)
            xlabel('Cloud condensed water (g/Kg)','fontweight','bold')
            ylabel('Pressure (mb)','fontweight','bold')
        elseif it == 16
            buoy=max(buoy,-0.1);
            [ntemp,mtemp]=size(buoy);
            if ntemp > mtemp
                buoy=transpose(buoy);
            end
            %cfrac=transpose(buoy);
            cfrac=buoy;
            contourf(timem,zkm,cfrac)
            %shading interp
            xlabel('Time (days)','fontweight','bold')
            %ylabel('Pressure (mb)','fontweight','bold')
            ylabel('Altitude (km)','fontweight','bold')
            title('Lifted parcel buoyancy (K)','fontweight','bold')
            set(gca,'fontweight','bold')
            u=get(gca,'YTick');
            u1=transpose(-u);
            v=cellstr(num2str(u1));
            %set(gca, 'YTickLabel',v)
            set(gca,'YLim',[0 16])
            colorbar
            pause
            tht1=cfrac(3:35,:);
            thtm=mean(max(tht1,0));
            h=plot(timem,thtm);
            set(h,'linewidth',3)
            set(gca,'fontweight','bold')
            xlabel('Time (days)','fontweight','bold')
            ylabel('Vertically averaged buoyancy','fontweight','bold')
        elseif it == 17
            cfrac=transpose(cldhov);
            contourf(timem,zkm,cfrac)
            xlabel('Time (days)','fontweight','bold')
            %ylabel('Pressure (mb)','fontweight','bold')
            ylabel('Altitude (km)','fontweight','bold')
            title('Cloud Fraction','fontweight','bold')
            set(gca,'fontweight','bold')
            u=get(gca,'YTick');
            u1=transpose(-u);
            v=cellstr(num2str(u1));
            %set(gca, 'YTickLabel',v)
            set(gca,'YLim',[0 16])
            colorbar
            pause
            tht1=cfrac(3:35,:);
            thtm=mean(tht1);
            h=plot(timem,thtm);
            set(h,'linewidth',3)
            set(gca,'fontweight','bold')
            xlabel('Time (days)','fontweight','bold')
            ylabel('Vertically averaged cloud fraction','fontweight','bold')
        elseif it == 18
            tht=transpose(thov);
            [m,n]=size(tht);
            ttemp=zeros(m,n);
            for j=1:n,
                ttemp(:,j)=tht(:,j)-tht(:,1);
            end
            tht=ttemp;
            contourf(timem,p,tht)
            xlabel('Time (days)','fontweight','bold')
            ylabel('Pressure (mb)','fontweight','bold')
            title('Temperature perturbation from initial state (K)','fontweight','bold')
            set(gca,'fontweight','bold')
            u=get(gca,'YTick');
            u1=transpose(-u);
            v=cellstr(num2str(u1));
            set(gca, 'YTickLabel',v)
            colorbar
            pause
            tht1=tht(3:35,:);
            thtm=mean(tht1);
            h=plot(timem,thtm);
            set(h,'linewidth',3)
            set(gca,'fontweight','bold')
            xlabel('Time (days)','fontweight','bold')
            ylabel('Vertically averaged temperature perturbation from initial state (K)','fontweight','bold')
        elseif it == 19
            tht=transpose(qhov);
            [m,n]=size(tht);
            ttemp=zeros(m,n);
            for j=1:n,
                ttemp(:,j)=tht(:,j)-tht(:,n);
            end
            tht=ttemp;
            contourf(timem,p,tht)
            xlabel('Time (days)','fontweight','bold')
            ylabel('Pressure (mb)','fontweight','bold')
            title(' Specific humidity perturbation from initial state (g/Kg)','fontweight','bold')
            set(gca,'fontweight','bold')
            u=get(gca,'YTick');
            u1=transpose(-u);
            v=cellstr(num2str(u1));
            set(gca, 'YTickLabel',v)
            colorbar
            pause
            tht1=tht(3:35,:);
            thtm=mean(tht1);
            h=plot(timem,thtm);
            set(h,'linewidth',3)
            set(gca,'fontweight','bold')
            xlabel('Time (days)','fontweight','bold')
            ylabel('Vertically averaged specific humidity perturbation from initial state (g/Kg)','fontweight','bold')
        elseif it == 20
            tht=transpose(rhhov);
            [m,n]=size(tht);
            ttemp=zeros(m,n);
            for j=1:n,
                ttemp(:,j)=tht(:,j)-tht(:,n);
            end
            tht=ttemp;
            contourf(timem,p,tht)
            xlabel('Time (days)','fontweight','bold')
            ylabel('Pressure (mb)','fontweight','bold')
            title('Relative humidity perturbation from initial state (%)','fontweight','bold')
            set(gca,'fontweight','bold')
            u=get(gca,'YTick');
            u1=transpose(-u);
            v=cellstr(num2str(u1));
            set(gca, 'YTickLabel',v)
            colorbar
            pause
            tht1=tht(3:35,:);
            thtm=mean(tht1);
            h=plot(timem,thtm);
            set(h,'linewidth',3)
            set(gca,'fontweight','bold')
            xlabel('Time (days)','fontweight','bold')
            ylabel('Vertically averaged relative humidity perturbation from initial state (%)','fontweight','bold')
        elseif it == 21
            tht=transpose(mhov);
            contourf(timem,p,tht)
            xlabel('Time (days)','fontweight','bold')
            ylabel('Pressure (mb)','fontweight','bold')
            title('Upward convective upward mass flux (10^{-3} Kg m^{-2} s^{-1})','fontweight','bold')
            set(gca,'fontweight','bold')
            u=get(gca,'YTick');
            u1=transpose(-u);
            v=cellstr(num2str(u1));
            set(gca, 'YTickLabel',v)
            colorbar
            pause
            tht1=tht(3:35,:);
            thtm=mean(tht1);
            h=plot(timem,thtm);
            set(h,'linewidth',3)
            set(gca,'fontweight','bold')
            xlabel('Time (days)','fontweight','bold')
            ylabel('Vertically averaged upward convective mass flux','fontweight','bold')
        elseif it == 22
            tht=transpose(mdhov);
            contourf(timem,p,tht)
            xlabel('Time (days)','fontweight','bold')
            ylabel('Pressure (mb)','fontweight','bold')
            title('Penetrative downdraft mass flux (10^{-3} Kg m^{-2} s^{-1})','fontweight','bold')
            set(gca,'fontweight','bold')
            u=get(gca,'YTick');
            u1=transpose(-u);
            v=cellstr(num2str(u1));
            set(gca, 'YTickLabel',v)
            colorbar
            pause
            tht1=tht(3:35,:);
            thtm=mean(tht1);
            h=plot(timem,thtm);
            set(h,'linewidth',3)
            set(gca,'fontweight','bold')
            xlabel('Time (days)','fontweight','bold')
            ylabel('Vertically averaged penetrative downdraft mass flux','fontweight','bold')
            %  %
        elseif it == 23
            tht=transpose(mphov);
            contourf(timem,p,tht)
            xlabel('Time (days)','fontweight','bold')
            ylabel('Pressure (mb)','fontweight','bold')
            title('Unsaturated downdraft mass flux (10^{-3} Kg m^{-2} s^{-1})','fontweight','bold')
            set(gca,'fontweight','bold')
            u=get(gca,'YTick');
            u1=transpose(-u);
            v=cellstr(num2str(u1));
            set(gca, 'YTickLabel',v)
            colorbar
            pause
            tht1=tht(3:35,:);
            thtm=mean(tht1);
            h=plot(timem,thtm);
            set(h,'linewidth',3)
            set(gca,'fontweight','bold')
            xlabel('Time (days)','fontweight','bold')
            ylabel('Vertically averaged unsaturated downdraft mass flux','fontweight','bold')
        elseif it == 24
            tht=transpose((mhov+mdhov+mphov));
            contourf(timem,p,tht)
            xlabel('Time (days)','fontweight','bold')
            ylabel('Pressure (mb)','fontweight','bold')
            title('Net convective upward mass flux (10^{-3} Kg m^{-2} s^{-1})','fontweight','bold')
            set(gca,'fontweight','bold')
            u=get(gca,'YTick');
            u1=transpose(-u);
            v=cellstr(num2str(u1));
            set(gca, 'YTickLabel',v)
            colorbar
            pause
            tht1=tht(3:35,:);
            thtm=mean(tht1);
            h=plot(timem,thtm);
            set(h,'linewidth',3)
            set(gca,'fontweight','bold')
            xlabel('Time (days)','fontweight','bold')
            ylabel('Vertically averaged net convective upward mass flux','fontweight','bold')
        elseif it == 0
            i = 0;
        end
        
        figure(gcf)
        drawnow;
        pause(0.05);
    end
end

