function main()

   %--read aerosol optical properties
   aerFile = ['~/sbdart/sbdart2.4/meto624_2012/ARM2012/',...
              'sgpaerosolbe1turnC1.c1.20070101.000000.cdf']; 

   [aer] = read_aerosol_data( aerFile );

   %--sampling the aerosol profile data, 
   %  since the max dim. for "zbaer" in SBDART is 65.
   %  length of aer.height = 177.
   htIdx = 1:5:length(aer.height);

   %--averget the single scattering albedo and asymmetry factor vertically.
   %  for simplicity we don't input profiles for these two parameters here.
   ssa = nanmean(aer.ssaprof);
   asy = nanmean(aer.asyprof);

   %-- get day of year and GMT time.
   idayV = ceil( aer.time./3600./24 ); %aer.time is seconds since 2007/1/1, 0:0:0
   gmtV = (aer.time - ( (idayV-1)*24*3600 ) )./3600;  %!time in hours GMT

   %-- calculate solar zenith angle and determine day time cases.
   %   For this exercise, we consider solar radiatoin only, so
   %   skip the night time cases.
   %   The SBDART is a plane-parallel model. don't work very well in very
   %   low sun. so we set the solar zenith angle cosine limit to 0.15.
   zenV = solar_zenith( idayV, gmtV, aer.lat, aer.lon );
   tmIdx = find( cos(zenV)>0.15 );
   

keyboard

   %-- setup the SBDART input parameters
   %
   isat=0;         %!user defined wavelength range with no filter function
   wlinf=0.25;     %!starting wavelength of 0.25 um
   wlsup=4.0;      %!ending wavelength of 4.0 um
   wlinc=-0.01;    %!wavelength increment is dL/L=.01
   nf=3;           %!use the more accurate MODTRAN3 solar spectrum

   alat=aer.lat;   %!latitude on Earth
   alon=aer.lon;   %!longitude on Earth

   idatm=6;        %!atmosphere: =3 midLatWinter, =6 is US standard

   isalb=6;        %!=6 vegetation albedo model.

   iaer=5;         %! =5 user defined spectral dependence of BLA 

   zout=[0,100];,  %!get output at the surface and top of atmosphere
   iout=10;        %!output mode for TOA and surface spectral fluxes
   nstr=4;         %!number of discrete ordinates in both hemispheres in DISORT


tic;
   %-- sbdart output to a file
   outFileName = [ './ex04.out' ];
   fout = fopen( outFileName, 'w');

   for it = 1:10:length(tmIdx) %To save time, run every 10th cases only.
      
      itm = tmIdx(it);

      iday = idayV(itm);
      time = gmtV(itm);


      %-- set up the aerosol parameters.
      %
      zbaer  = aer.height(htIdx);       %altitude grid for custom aerosol in (km)
      dbaer  = aer.extprof(htIdx, itm); %aerosol density profile.
      wlbaer = [.355, .500];                       %wavelengthes at which aerosol properties are specified.
      qbaer  = [aer.aod500(itm), aer.aod355(itm)]; %aerosol optical depth
      wbaer  = ones(size(wlbaer)) * ssa(itm);      %aerosol single scattering albedo
      gbaer  = ones(size(wlbaer)) * asy(itm);      %aerosol asymmetry factor
      abaer  = aer.angst(itm);                     %aerosol Angstrom exponent.

      %-- if any of following parameters is missing,
      %   no calculation is done.
      if any( isnan( [dbaer', qbaer, wbaer, gbaer, abaer] ) )
         continue;
      end


      disp(sprintf('itm=%-d, iday=%-d, gmt=%-6.3f',itm, iday, time))

      %--write the INPUT file
      %         
      fid = fopen( 'INPUT', 'w' );
         
      fprintf( fid, '&INPUT\n' );
      fprintf( fid, 'ISAT   = %-d,\n',    isat  ); %!user defined wavelength range with no filter function
      fprintf( fid, 'WLINF  = %-6.3f,\n', wlinf );
      fprintf( fid, 'WLSUP  = %-6.3f,\n', wlsup );
      fprintf( fid, 'WLINC  = %-6.3f,\n', wlinc );
      fprintf( fid, 'NF     = %-d,\n',    nf    );

      fprintf( fid, 'IDAY   = %-d,\n',    iday  );
      fprintf( fid, 'TIME   = %-9.5f,\n', time  );
      fprintf( fid, 'ALAT   = %-9.5f,\n', alat  );
      fprintf( fid, 'ALON   = %-9.5f,\n', alon  );

      fprintf( fid, 'IDATM  = %-d,\n',    idatm );
     
      fprintf( fid, 'ISALB  = %-d,\n',    isalb );

      fprintf( fid, 'IAER   = %-d,\n',    iaer  );
      fprintf( fid, 'ZBAER  = '); fprintf( fid, '%-6.3f, ', zbaer ); fprintf( fid, '\n');
      fprintf( fid, 'DBAER  = '); fprintf( fid, '%-7.4f, ', dbaer ); fprintf( fid, '\n');
      fprintf( fid, 'WLBAER = '); fprintf( fid, '%-6.3f, ', wlbaer); fprintf( fid, '\n');
      fprintf( fid, 'QBAER  = '); fprintf( fid, '%-9.5f, ', qbaer ); fprintf( fid, '\n');
      fprintf( fid, 'WBAER  = '); fprintf( fid, '%-10.4E,', wbaer ); fprintf( fid, '\n');
      fprintf( fid, 'GBAER  = '); fprintf( fid, '%-10.4E,', gbaer ); fprintf( fid, '\n');
      fprintf( fid, 'ABAER  = %-10.4E,\n', abaer );

      fprintf( fid, 'ZOUT   = '); fprintf( fid, '%-6.2f, ', zout ); fprintf( fid, '\n');
      fprintf( fid, 'IOUT   = %-d\n', iout );
      fprintf( fid, 'NSTR   = %-d\n', nstr );
      fprintf( fid, '/' );
      
      fclose(fid); %fid = fopen( 'INPUT', 'w' );



      %-- Run the sbdart model with the above given INPUT file.
      %
      [stus, outFluxStr] = system('./sbdart');
      if ( stus ~= 0 ) 
         disp('--- calling sbdart failed !')
         keyboard
         return
      end
      outF = sscanf(outFluxStr, '%f %f %f %e %e %e %e %e %e');

     
      %-- output foramt:
      %   one output record per run, integrated over wavelength.
      %   output quantities are, (integrations by trapezoid rule)
      %
      %   time zen iaer WLINF,WLSUP,FFEW,TOPDN,TOPUP,TOPDIR,BOTDN,BOTUP,BOTDIR
      %
      %     time = seconds since 2007/1/1 0:0:0
      %     zen  = solar zenith angle
      %     iaer = iaer value 
      %     WLINF = lower wavelength limit             (microns)
      %     WLSUP = upper wavelength limit             (microns)
      %     FFEW  = filter function equivalent width   (microns)
      %     TOPDN = total downward flux at ZOUT(2) km  (w/m2)
      %     TOPUP = total upward flux at ZOUT(2) km    (w/m2)
      %     TOPDIR= direct downward flux at ZOUT(2) km (w/m2)
      %     BOTDN = total downward flux at ZOUT(1) km  (w/m2)
      %     BOTUP = total upward flux at  ZOUT(1) km   (w/m2)
      %     BOTDIR= direct downward flux at ZOUT(1) km (w/m2)
      %
      fprintf( fout, ...
               ['%-d %-6.2f %-d   ', '%7.3f %7.3f %7.3f %12.4E %12.4E %12.4E %12.4E %12.4E %12.4E\n'],...
                aer.time(itm), zenV(itm), iaer, outF );


      
   end %for itm = 1:length(aer.time)

   fclose(fout);      
toc;


keyboard
end



%----------------------------------------------------------------------
function [aer] = read_aerosol_data( fileName )
%----------------------------------------------------------------------

   %fileName = ['~/sbdart/sbdart2.4/meto624_2012/ARM2012/',...
   %            'sgpaerosolbe1turnC1.c1.20070101.000000.nc']; 
   ncid = netcdf.open(fileName,'NC_NOWRITE');


   %-- get time, location and dimension values
   %
   varID = netcdf.inqVarID(ncid,'lat'); %in deg.
   lat   = netcdf.getVar(ncid,varID);

   varID = netcdf.inqVarID(ncid,'lon'); %in deg.
   lon   = netcdf.getVar(ncid,varID);

   varID = netcdf.inqVarID(ncid,'alt'); %meters above sea level
   alt   = netcdf.getVar(ncid,varID);

   varID = netcdf.inqVarID(ncid,'time'); % seconds since 2007/1/1 0:0:0
   time  = netcdf.getVar(ncid,varID);

   varID  = netcdf.inqVarID(ncid,'height'); %in (km)
   height = netcdf.getVar(ncid,varID);



   %-- get variables
   %
   varID  = netcdf.inqVarID(ncid,'be_aod_500'); %unitless, [nTimes]
   aod500 = netcdf.getVar(ncid,varID);

   varID  = netcdf.inqVarID(ncid,'be_aod_355'); %unitless, [nTimes]
   aod355 = netcdf.getVar(ncid,varID);

   varID = netcdf.inqVarID(ncid,'be_angst_exp'); %unitless, [nTimes]
   angst = netcdf.getVar(ncid,varID);

   varID   = netcdf.inqVarID(ncid,'extinction_profile'); %unitless, [nTimes]
   extprof = netcdf.getVar(ncid,varID);

   varID   = netcdf.inqVarID(ncid,'single_scattering_albedo'); %unitless, [nHeight,nTimes]
   ssaprof = netcdf.getVar(ncid,varID);

   varID   = netcdf.inqVarID(ncid,'asymmetry_parameter'); %unitless, [nHeight,nTimes]
   asyprof = netcdf.getVar(ncid,varID);

   %varID  = netcdf.inqVarID(ncid,'rh'); %in (%), [nHeight,nTimes]
   %rhprof = netcdf.getVar(ncid,varID);

   
   %-- deal with missing values
   %
   MISG = netcdf.getAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'missing-value');
   MISG = str2num(MISG);

   aod500(  aod500  ==MISG ) = NaN;
   aod355(  aod355  ==MISG ) = NaN;
   angst(   angst   ==MISG ) = NaN;
   extprof( extprof ==MISG ) = NaN;
   ssaprof( ssaprof ==MISG ) = NaN;
   asyprof( asyprof ==MISG ) = NaN;
   %rhprof(  rhprof  ==MISG ) = NaN;
   
   %-- deal with unphsysical values.
   extprof(extprof<0)=0;
   ssaprof(ssaprof<0)=nan;
   asyprof(asyprof<0)=nan;
   

   %-- group into a structure.
   aer.lat     = lat;
   aer.lon     = lon;
   aer.alt     = alt;
   aer.time    = time;
   aer.height  = height;
   aer.aod500  = aod500;
   aer.aod355  = aod355;
   aer.angst   = angst;
   aer.extprof = extprof;
   aer.ssaprof = ssaprof;
   aer.asyprof = asyprof;


   %-- close the netcdf file.
   netcdf.close(ncid);


%keyboard
end




function [zen] = solar_zenith( iday, gmt, lat, lon )

   dtor = pi/180;

   iday = mod(iday-1,365)+1;

   Dangle = 2 * pi * ( iday-1 ) ./ 365;
   
   rsun = 1.00011 + 0.034221*cos(  Dangle) + 0.001280*sin(  Dangle) +...
                    0.000719*cos(2*Dangle) + 0.000077*sin(2*Dangle);
            
   decl = 0.006918 - 0.399912*cos(  Dangle) + 0.070257*sin(  Dangle) -...
                     0.006758*cos(2*Dangle) + 0.000907*sin(2*Dangle) -...
                     0.002697*cos(3*Dangle) + 0.001480*sin(3*Dangle);

   eqtime = ( 0.000075 + 0.001868*cos(  Dangle) - 0.032077*sin(  Dangle) -...
                         0.014615*cos(2*Dangle) - 0.040890*sin(2*Dangle) ) *12/pi;


   sunlon = -15*(gmt-12+eqtime./60);
      
   t0 = (90-lat )*dtor;
   t1 = (90-decl)*dtor;
   p0 = lon*dtor;
   p1 = sunlon*dtor;
   zz = cos(t0).*cos(t1) + sin(t0).*sin(t1).*cos(p1-p0);
   zen = acos(zz)./dtor; %in deg.

end








