

function plot_sbdart_output()

   datafile = './ex01.out';

   %-- read the SBDART out file.
   fid=fopen(datafile,'r');
   dat = textscan(fid, '%f', 'HeaderLines',3);
   fclose(fid);
   
   %-- reformat to a Matlab matrix
   dat = cell2mat(dat);
   dat = reshape(dat, 8,[]);

   %-- extract paramters.
   wl     = dat(1,:); %wavelength
   topdn  = dat(3,:);
   topup  = dat(4,:);
   topdir = dat(5,:);
   botdn  = dat(6,:);
   botup  = dat(7,:);
   botdir = dat(8,:);
   

   %-- compute integrated broad band fluxes.
   dwl = wl(2:end)-wl(1:end-1);
   Ftopdn  = sum( dwl .* (0.5*(topdn(2:end)+topdn(1:end-1)) ) );
   Ftopup  = sum( dwl .* (0.5*(topup(2:end)+topup(1:end-1)) ) );
   Fbotdn  = sum( dwl .* (0.5*(botdn(2:end)+botdn(1:end-1)) ) );
   Fbotup  = sum( dwl .* (0.5*(botup(2:end)+botup(1:end-1)) ) );
   Fbotdir = sum( dwl .* (0.5*(botdir(2:end)+botdir(1:end-1)) ) );
   disp(sprintf('Ftopdn=%-6.1f, Ftopup=%-6.1f, Fbotdn=%-6.1f, Fbotup=%-6.1f, Fbotdir=%-6.1f',...
                 Ftopdn,Ftopup,Fbotdn,Fbotup,Fbotdir))
   

   %-- plot the SBDART output.
   figure(11)
   set(gca,'FontSize',14)
   plot(wl, [topdn; topup; topdir])
   axis([0,inf,-5,inf]);
   grid on
   legend({'topdn','topup','topdir'})
   title('Spectral radiative fluxes')
   ylabel('W/m^2/{\mu}m');
   xlabel('Wavelength ({\mu}m)');

   %-- plot the SBDART output.
   figure(12)
   set(gca,'FontSize',14)
   plot(wl, [botdn; botup; botdir])
   axis([0,inf,-5,inf]);
   grid on
   legend({'botdn','botup','botdir'})
   title('Spectral radiative fluxes')
   ylabel('W/m^2/{\mu}m');
   xlabel('Wavelength ({\mu}m)');

%{   
   %-- plot the SBDART output.
   figure(11)
   plot(wl, [topdn; topup; topdir; botdn; botup; botdir])
   axis([0,inf,-50,inf]);
   grid on
   legend({'topdn','topup','topdir','botdn','botup','botdir'})
   title('Spectral radiative fluxes')
   ylabel('W/m^2/{\mu}m');
   xlabel('Wavelength ({\mu}m)');


   %-- plot the difference from the contral run
   fid = fopen('control.out','r');
   dat_ctl = textscan(fid,'%f','HeaderLines',3);
   fclose(fid);
   dat_ctl = cell2mat(dat_ctl);
   dat_ctl = reshape(dat_ctl,8,[]);
   diff = dat-dat_ctl;
   figure(10) 
   plot(wl,diff(3:8,:))
   grid on
   legend({'topdn','topup','topdir','botdn','botup','botdir'},'Location','Best')
   title('Differences between the test and the control run')
   ylabel('W/m^2/{\mu}m');
   xlabel('Wavelength ({\mu}m)');
%}
%keyboard
end
