function ep_qcflag(infile);
% EP_QCFLAG - Create a quality control file for eqpac data products
%
% Use as: ep_qcflag(infile)
%
% Inputs: infile = an eqpac OMS netcdf object
%

% Brian Schlining
% 17 Aug 2000

filename = name(infile);
[p f e] = fileparts(filename);
outFilename = [p filesep f '_qc' e];

% Create a new netcdf file
outfile = netcdf(outFilename, 'clobber');

if ~isempty(outfile)
   OK = 1;
   fprintf(1, 'Reading %s\n', char(ncnames(infile)));
   
   %=======================================
   % Grab the data and products
   %=======================================
   X = ep_cmplNc(infile);
   flag = zeros(length(infile{'time'}(:)), length(X.lambda));
  
   %========================
   % Create new netcdf file
   %========================
   fprintf(1, 'Writing %s\n', char(ncnames(outfile)));
   outfile.conventions = 'MBARI/timeSeries/mooring/spectroradiometer';
   outfile.creationDate = datestr(now, 1);
   outfile.lastModified = datestr(now, 1);
   outfile.mooring = infile.mooring(:);
   outfile.description = ['Quality control flags for ' f e];
   
   project = infile.project(:);
   if isempty(project)
      project = 'null';
   end
   outfile.project = project;
   
   waterDepth = infile.waterDepth(:);
   if isempty(waterDepth)
      waterDepth = 'null';
   end
   outfile.waterDepth = waterDepth;
   
   outfile.keywords = 'quality control, spectroradiometers';
   outfile.dataComment0 = 'QC flags are bit based. 0 = good 1 = bad'; 
   outfile.dataComment1 = 'The first (least-signifigant) bit is for the Lu[20m].';
   outfile.dataComment2 = 'The second bit is for the Ed[20m].';
   outfile.dataComment3 = 'The Third bit is for the Es[0m].';
   outfile.dataComment4 = 'The Fourth bit indicates an unspecified problem';
   
   % Copy Dimensions
   d    = dim(infile);
   names = ncnames(d);
   
   for i = 1:length(d)
      if isrecdim(d{i})
         fprintf(1,'Copying unlimited-dimension "%s"\n', names{i});
         outfile(names{i}) = 0;
      else
         if ~strcmp(names{i},'bankDepth')
            fprintf(1,'Copying dimension "%s"\n', names{i});
            copy(d{i}, outfile);
         end
      end
   end
   
   % Coy the important co-ordinate variables
   copy(infile{'time'}, outfile, 1, 1);
   copy(infile{'latitude'}, outfile, 1, 1);
   copy(infile{'longitude'}, outfile, 1, 1);
   def = zeros(size(X.t));
   
   for Lambda = X.lambda
      
      qcName = ['qc' num2str(Lambda)];
      outfile{qcName} = ncint('time', 'longitude', 'latitude');
      outfile{qcName}.long_name = ['Quality control flag for ' num2str(Lambda) ' nm channels'];
      outfile{qcName}.units = 'bit = 1 is bad. bits: 1= Lu[20m], 2 = Ed[20m], 3 = Es';
      outfile{qcName}.symbol = 'flag';
      outfile{qcName}.wavelength = Lambda;
      outfile{qcName}.FillValue_ = -999;
      outfile{qcName}.missing_value = -999;
      outfile{qcName}(:) = def;
      
   end

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %========================
   % Apply quality controls
   %========================
   
   %====================================================
   % Calculate Frouins clear-sky surface radiance (Et).
   % Flag Es > Et*1.15 as bad
   %====================================================
   Ef = frouin(nctime(infile), infile{'latitude'}(:), infile{'longitude'}(:));
   Ef = Ef(:, 1:5);
   
   % Verify Es < Ef*1.15
   Ef2     = Ef*1.15;
   EsDiff  = Ef2 - X.Es;
   bad     = find(EsDiff < 0);
   X.Es(bad) = NaN;
   
   %==========================================================
   % Locate Night time data and get rid of it.
   % This removes data where the solar zenith angle is below the 
   % horizon.
   %==========================================================
   sunRiseHour = sunrise(X.t, infile{'latitude'}(:), infile{'longitude'}(:)) + gmtoffset(infile{'longitude'}(:));
   sunSetHour  = sunset(X.t, infile{'latitude'}(:), infile{'longitude'}(:)) + gmtoffset(infile{'longitude'}(:));
   [localYear localMonth localDay localHour localMinute localSecond] = datevec(X.t + gmtoffset(infile{'longitude'}(:)));
   localHour   = localHour + (localMinute + localSecond/60)/60;
   riseCheck   = localHour - sunRiseHour;   % Negative if before sunrise
   setCheck    = sunSetHour - localHour;     % Negative if after sunset
   bad         = find(riseCheck < 0 & setCheck < 0 );	% Flag data before sunrise/ after sunset
   X.Es(bad, :)  = NaN;
   X.Ed(bad, :)  = NaN;
   X.Lu(bad, :)  = NaN;
   
   %==============================================================
   % Calculate fouling ratio. Flag ratios greater than 0.1 as bad
   %==============================================================
   % Not implemented Ed683 and Lu683 values are not trustworthy.
   % Probably due to instrument self-shading and mooring shadows
   
   %============================================================
   % Locate irradiance attenuation less than that of pure water
   % Flag it as bad
   %============================================================
   % Kw - Pure water based on Morel
   Kw = repmat(morel_k(0, X.lambda), length(X.t), 1);
   Kw(:,5)=0.06; %Cheating for lambda 555
   KDiff = X.K - Kw;  % Negative mean Ke is less than pure water
   bad = find(KDiff < 0);
   X.Es(bad) = NaN;
   X.Ed(bad) = NaN;
   
   %==================================================
   % Flag bad areas that are 'known' to be unreliable.
   % This is determined by eyeballing the data
   %==================================================
   badStart = NaN;
   badEnd = NaN;
   switch lower(infile.mooring(:))
   case 'ep1'
      badStart = [0 729543.920724017 729963.9221371198 730274.922211625 datenum('07 Jan 2000') datenum('06 Dec 2000')];
      badEnd   = [datenum('01 Nov 1997') 729708.916940933 730062.9169316513 730313.922136382 datenum('22 Oct 2000') datenum('07 Jun 2001')];
      
   case 'ep2'
      badStart = [0;datenum('11 Nov 1998')];
      badEnd   = [datenum('01 Nov 1998');datenum('11 Nov 1999')];
   end
   
   for i = 1:length(badStart)
      bad = find(X.t >= badStart(i) & X.t <= badEnd(i));
      switch lower(infile.mooring(:))
   	case 'ep1'
         Essave=find(X.t(bad)>=datenum('03 Jun 1997')&X.t(bad)<=datenum('12 Nov 1997')|X.t(bad)>=datenum('24 Jul 1998')&X.t(bad)<=datenum('20 Oct 1998')|...
            X.t(bad)>=datenum('04 Jun 1999')&X.t(bad)<=datenum('09 Jul 1999')|X.t(bad)>=datenum('06 Dec 2000')&X.t(bad)<=datenum('07 Jun 2001'));%Shutter closed 06 Dec 2001/this saves 0m Ed
   	X.Lu(bad, :)= NaN;%Shutter closed 06 Dec 2001 ep1/this saves 0m Ed
   	X.Ed(bad, :)=NaN;%Shutter closed 06 Dec 2001 ep1/this saves 0m Ed
   	bad(Essave)=[];%Shutter closed 06 Dec 2001 ep1/this saves 0m Ed
		flag(bad, :) = NaN;
		case 'ep2'
   	flag(bad, :) = NaN;
  		end;
   end
   
   %=========================================================
   % Use OC2 CHL algorythmsn to flag bad areas of data.
   % When the coeffiecint of variation betwen calculated chl 
   % is too large flag it bad
   %=========================================================
   
   bad = find(X.oc_cov > 0.4);
   switch lower(infile.mooring(:))
   case 'ep1'
      Essave=find(X.t(bad)>=datenum('03 Jun 1997')&X.t(bad)<=datenum('12 Nov 1997')|X.t(bad)>=datenum('24 Jul 1998')&X.t(bad)<=datenum('20 Oct 1998')|...
         X.t(bad)>=datenum('04 Jun 1999')&X.t(bad)<=datenum('09 Jul 1999')|X.t(bad)>=datenum('06 Dec 2000')&X.t(bad)<=datenum('07 Jun 2001'));%Shutter closed 06 Dec 2001/this saves 0m Ed
   X.Lu(bad, :)= NaN;%Shutter closed 06 Dec 2001 ep1/this saves 0m Ed
   X.Ed(bad, :)=NaN;%Shutter closed 06 Dec 2001 ep1/this saves 0m Ed
   bad(Essave)=[];%Shutter closed 06 Dec 2001 ep1/this saves 0m Ed
	flag(bad, :) = NaN;
	case 'ep2'
   flag(bad, :) = NaN;
   end;
   

   for i = 1:length(X.lambda)
      
      buf = zeros(size(X.t));
      bad = find(isnan(X.Lu(:,i)));       % bitset 1 if Lu is bad
      if ~isempty(bad)
         buf(bad) = bitset(buf(bad),1);
      end
      bad = find(isnan(X.Ed(:,i)));       % bitset 2 if Ed is bad
      if ~isempty(bad)
         buf(bad) = bitset(buf(bad),2);
      end
      bad = find(isnan(X.Es(:,i)));       % bitset 3 if Es is bad
      if ~isempty(bad)
         buf(bad) = bitset(buf(bad),3);
      end
      bad = find(isnan(flag(:,i)));       % bitset 4 for unspecified problems
      if ~isempty(bad)
         buf(bad) = bitset(buf(bad),4);
      end
      
      qcName = ['qc' num2str(X.lambda(i))];  % Store it in the netcdf file
      outfile{qcName}(:) = buf;
      
   end
   
   
   ncquit(outfile);
else
   OK = 0;
end