function oms_makeSeabassProducts(inFilename)
% oms_makeSeabassProducts - Shell for creating product files for HF OMS data
%
% Reads data from a HF data file and the daily avg data file. Derives Lw for
% the HF data based on the daily avg K (attenuation coef) and writes the data
% to a *.prod.nc file
%
% Use as: oms_makeSeabassProducts(infile)
%
% Inputs: inFilename = String name of the netcdf file containing the high frequency 
%                  oasis mooring spectoradiometer data. This function will
%                  requires access to daily averaged data also. If 
%                  outfile =  'ep1.nc', then it will look for 'ep1.avg.nc'.
%
% Requires: NEAR

% Brian Schlining
% 13 Dec 1999
% 15 Mar 2000; Updated to handle new MBARI netCDF conventions

%==================
% Open the HF file
%==================
[p filename ext] = fileparts(inFilename);
infile = netcdf(inFilename, 'nowrite');
if isempty(infile)			% Exit if all the files can't be opened
   fprintf(1, '\nError: Unable to open %s\n', inFilename)
   return;
end
t0     = nctime(infile);

%=========================
% Open the daily avg file
%=========================
avgFilename = [p filesep filename '.avg.nc'];
avgfile = netcdf(avgFilename, 'nowrite');
if isempty(avgfile) 
   fprintf(1, '\nError: Unable to open %s\n', avgFilename)
   return
end

%==============================================
% Create a quality control flag file if needed
%==============================================
qcFilename = [p filesep filename '_qc' ext];
qcfile = netcdf(qcFilename);
OK = 1;
if ~isempty(qcfile)
   tqc = nctime(qcfile); 
   if length(tqc) ~= length(t0); % Check to see if the QC flag is he correct size
      ncquit(qcfile);
      OK = 0;
   end
   
else 
   OK = 0;
end

if ~OK % If qc file is not current or doesn't exist create a new one
   fprintf(1, 'Performing quality control on source data...');
   ep_qcflag(infile);
   fprintf(1,'completed\n');
   qcfile = netcdf(qcFilename);
end

%========================
% Create new netcdf file
%========================
% Create the products file
outFilename = [p filesep filename '.prod.nc'];
fprintf(1, 'Creating %s\n', outFilename);
outfile = netcdf(outFilename, 'clobber');

% Add global attributes
outfile.conventions = 'MBARI/timeSeries/mooring/spectroradiometer';
outfile.creationDate = datestr(now, 1);
outfile.lastModified = datestr(now, 1);
outfile.mooring = infile.mooring(:);
outfile.description = ['Radiometer data products derived from ' filename ext];
[f p e] = fileparts(char(ncnames(infile)));
outfile.keywords = 'spectroradiometer, attenuation, water-leaving radiance';
outfile.sourceFile = ncchar([f filesep p e]);
outfile.dataComment0 = ['This archive contains products derived from data in ' p e];

% Add dimensions
outfile('latitude')            = 1;                     % Dimension
outfile{'latitude'}            = ncdouble('latitude');  % Specify Co-ordinate variable dimensions
outfile{'latitude'}.long_name  = 'Latitude';
outfile{'latitude'}.units      = 'degrees_north';
outfile{'latitude'}(:)         = infile{'latitude'}(:);

outfile('longitude')           = 1;                     % Dimension
outfile{'longitude'}           = ncdouble('longitude'); % Specify Co-ordinate variable dimensions
outfile{'longitude'}.long_name = 'Longitude';
outfile{'longitude'}.units     = 'degrees_east';
outfile{'longitude'}(:)        = infile{'longitude'}(:);

outfile('depth')               = 1;                     % Dimension
outfile{'depth'}               = ncdouble('depth'); % Specify Co-ordinate variable dimensions
outfile{'depth'}.long_name     = 'Depth';
outfile{'depth'}.units         = 'meters';
outfile{'depth'}.positive      = 'down';
outfile{'depth'}(:)            = 0;

copy(infile{'time'}, outfile, 1, 1); 

%===========================================
% Find out what Es wavelengths are measured 
%===========================================
S         = oms_getinfo(avgfile);
good      = strmatch('es',S.symbol);
names_Es  = S.name(good);
lambda_Es = S.wavelength(good);

% Find out what Lu wavelengths are measured by the ~20m instruments
% If no 20m instument is found use whatever one is closest
i         = strmatch('lu',S.symbol);
depths    = unique(S.instrumentDepth(i));
depth     = depths(near(depths, 20));
j         = find(S.instrumentDepth == depth);
good      = intersect(j, i);
names_Lu  = S.name(good);
lambda_Lu = S.wavelength(good);

% Get the wavelengths that are common to Es and Lu[~20m]
[lambda iEs iLu] = intersect(lambda_Es, lambda_Lu);

% loop through each wavelength and derive Lw for that wavelength
lon    = infile{'longitude'}(:);
offset = gmtoffset(lon);
tin    = t0 + offset;  % Get the HF time (local)
tavg   = nctime(avgfile) + offset; % Get the daily average time (local)
n      = noon(tavg,0)./24;		  % Get the time of noon (local)
tNoon  = floor(tavg) + n;
def    = ones(size(tin))*-999;
count = 0;
for i = lambda
   count = count + 1;
   
   % Calculate a daily average K
   k = oms_k1(avgfile, i);
   v = var(infile, char(names_Lu{iLu(count)})); % Get the HF variable object for Lu(lambda)
   
   % Create a new netcdf variable in the output file to hold Lw
   lamN = num2str(i);
   nLw = ['Lw' lamN];
   outfile{nLw} = ncfloat('time', 'depth', 'latitude', 'longitude');
   outfile{nLw}.long_name = ['Water-Leaving radiance at ' lamN ' nm'];
   outfile{nLw}.units = 'uwatts/cm2/nm/str';
   outfile{nLw}.instrumentDepth = 0;
   outfile{nLw}.symbol = 'Lw';
   outfile{nLw}.wavelength = i;
   outfile{nLw}.source = [char(ncnames(v)) ' at ' num2str(v.instrumentDepth(:)) ' m'];
   outfile{nLw}.FillValue_ = -999;
   outfile{nLw}.missing_value = -999;
   outfile{nLw}(:) = def;
   
   nLwn = ['Lwn' lamN];
   outfile{nLwn} = ncfloat('time', 'depth', 'latitude', 'longitude');
   outfile{nLwn}.long_name = ['Solar-normalized water-Leaving radiance at ' lamN ' nm'];
   outfile{nLwn}.units = 'uwatts/cm2/nm/str';
   outfile{nLwn}.instrumentDepth = 0;
   outfile{nLwn}.symbol = 'Lwn';
   outfile{nLwn}.wavelength = i;
   outfile{nLwn}.source = [char(ncnames(v)) ' at ' num2str(v.instrumentDepth(:)) ' m'];
   outfile{nLwn}.FillValue_ = -999;
   outfile{nLwn}.missing_value = -999;
   outfile{nLwn}(:) = def;
   
   nRrs = ['Rrs' lamN];
   outfile{nRrs} = ncfloat('time', 'depth', 'latitude', 'longitude');
   outfile{nRrs}.long_name = ['Remote sensing reflectance at ' lamN ' nm'];
   outfile{nRrs}.units = '1/str';
   outfile{nRrs}.instrumentDepth = 0;
   outfile{nRrs}.symbol = 'Rrs';
   outfile{nRrs}.wavelength = i;
   outfile{nRrs}.source = [char(ncnames(v)) ' at ' num2str(v.instrumentDepth(:)) ' m'];
   outfile{nRrs}.FillValue_ = -999;
   outfile{nRrs}.missing_value = -999;
   outfile{nRrs}(:) = def;
   
   %====================================================
   % Apply the daily average K to the HF data to get Lw
   %====================================================
   for d = 1:length(tavg) 
      good = find((tin > tNoon(d) - 0.5) & (tin <= tNoon(d) + 0.5));
      if ~isempty(good)
         dat = v(good);
         bad = find(dat == v.FillValue_(:));
         dat(bad) = NaN;
         Lw = calc_lw(dat, k(d), v.instrumentDepth(:));
         if ~isempty(Lw) 
            bad = find(isnan(Lw));
            Lw(bad) = -999;
            outfile{nLw}(good) = Lw;
            outfile{nLwn}(good) = lw2lwn(Lw, t0(good), infile{'latitude'}(:), infile{'longitude'}(:), i);
         end
      end
      
   end
   
end

ncquit(outfile);  % Have to close it to update the file

%=========================================
% Calculate Rrs and apply quality control
%=========================================
outfile = netcdf(outFilename, 'write');   % Open it back up to get Lw
EsIndex      = strmatch('es',S.symbol);
for i = lambda
   
   lamN = num2str(i);
   nRrs = ['Rrs' lamN];
   nLwn = ['Lwn' lamN];
   nLw  = ['Lw' lamN];
   nQc  = ['qc' lamN];
   flag = qcfile{nQc}(:);
   bad  = find(flag ~= 0);
   good = find(S.wavelength(EsIndex) == i);  
   
   
   if ~isempty(good)
      % Calc Rrs
      Rrs = outfile{nLw}(:)./infile{S.name{EsIndex(good)}}(:);
      Rrs(find(Rrs < 0)) = -999;
      
   end
   
   
   Lwn = outfile{nLwn}(:);
   Lwn(find(Lwn < 0)) = -999;
   Lw  = outfile{nLw}(:);
   Lw(find(Lw < 0)) = -999;
   
   Rrs(bad) = -999;
   Lwn(bad) = -999;
   Lw(bad)  = -999;
   
   
   outfile{nRrs}(:) = Rrs;
   outfile{nLwn}(:) = Lwn;
   outfile{nLw}(:)  = Lw;
   
end

%===============
% Add OC2V4 Chl 
%===============
%ncquit(outfile);  % Have to close it to update the file
%outfile = netcdf(outFilename, 'write');   % Open it back up to get Lw

%nChl = ['CHL'];
%outfile{nChl} = ncfloat('time', 'depth', 'latitude', 'longitude');
%outfile{nChl}.long_name = ['OC2V4 Chlorophyll from Rrs(490)/Rrs(555)'];
%outfile{nChl}.units = 'mg/m^3';
%outfile{nChl}.instrumentDepth = 0;
%outfile{nChl}.symbol = 'CHL';
%outfile{nChl}.wavelength = i;
%outfile{nChl}.source = [char(ncnames(v)) ' at ' num2str(v.instrumentDepth(:)) ' m'];
%outfile{nChl}.FillValue_ = -999;
%outfile{nChl}.missing_value = -999;
%outfile{nChl}(:) = def;

%Chl = oc2v4(outfile{'Rrs490'}(:), outfile{'Rrs555'}(:));
%Chl = real(Chl);
%bad = find(Chl < 0 | isnan(Chl) | isinf(Chl);
%Chl(bad) = -999;

ncquit(outfile);
ncquit(qcfile);
ncquit(avgfile);
ncquit(infile);
