%{
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
it=input('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 ');
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
break
end
end