function [badBursts, badCells] = mask_argnCurr_cdf(cdfFileRoot);
% mask_argnCurr_cdf.m  A function to QA/QC Aquadopp velocity data in a netCDF
%                     (.cdf) file.  
%
%    usage:  mask_argnCurr_cdf(cdfFileRoot);
%
%        where:  cdfFileRoot - the name of the .cdf file which contains
%                               Aquadopp velocity data as output by 
%                               AquaPro, surrounded by single quotes with
%                               no file extension 
%                badBursts - a list of bad bursts, which are collected
%                             during instrument deployment and recovery
%                badCells - a list of bad cells, which are located above
%                             MSL and contaminated by sidelobe reflection
%                             errors
%
% Copyright 2004 
% USGS Woods Hole Field Center
% Written by Charlene Sullivan
% csullivan@usgs.gov
% 
% Dependencies:
%   smart_interp.m
%   gmin.m
%   gmax.m
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use of this program is described in:
%
% Sullivan, C.M., Warner, J.C., Martini, M.A., Voulgaris, G., 
% Work, P.A., Haas, K.A., and Hanes, D.H. (2006) 
% South Carolina Coastal Erosion Study Data Report for Observations
% October 2003 - April 2004., USGS Open-File Report 2005-1429.
%
% Program written in Matlab v7.1.0 SP3
% Program ran on PC with Windows XP Professional OS.
%
% "Although this program has been used by the USGS, no warranty, 
% expressed or implied, is made by the USGS or the United States 
% Government as to the accuracy and functioning of the program 
% and related program material nor shall the fact of distribution 
% constitute any such warranty, and no responsibility is assumed 
% by the USGS in connection therewith."
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

% C. Sullivan   06/23/05,   version 1.0
% This mfile applies QA/QC measures including removal of bad data collected
% during deployment and recovery of the instrument, and removal of cells
% (bins) contaminated by sidelobe reflection errors.


more off

version = '1.0';

% Copy cdfFileRoot.cdf to cdfFileRoot-mask.cdf. In the 'mask' version
% all bad data will be masked (set to FILL). Masked data should NOT be
% written to the BBV .nc file. The mask file is an intermediate file,
% between the .cdf file created from AquaPro output, and the .nc file
% created for distribution.
copyfile([cdfFileRoot,'.cdf'],[cdfFileRoot,'-mask.cdf'])

cdf = netcdf([cdfFileRoot,'-mask.cdf'],'w');

% Get the dimensions
theBurstDim = cdf('burst');
nBursts = max(size(theBurstDim));

% Identify bad data collected during instrument deployment and 
% recovery by looking at a plot of pressure. Ask the user to enter the
% first and last good bursts of data. Data up to the first good burst and
% data after the last good burst will be masked.
figure
plot(cdf{'burst'}(:), cdf{'pressure'}(:))
hold on
ylabel('Pressure (m)')
xlabel('Burst number')
disp(' ')
disp('Please view the figure (feel free to zoom!) and look for bad')
disp('bursts collected during instrument deployment and recovery')
answer = input(['Would you like to mask any bad bursts? y/n :  '],'s');
if ~strcmp(answer,'Y') & ~strcmp(answer,'y') & ...
   ~strcmp(answer,'N') & ~strcmp(answer,'n')
        disp(['Valid answers are either yes ("y") or no ("n")'])
        answer = input(['Would you like to mask any bad bursts? y/n :  '],'s');
end
if strcmp(answer,'N') || strcmp(answer,'n')
    first_good = 1;
    last_good = nBursts;
    bad_bursts = [];
    disp(' ')
    disp('Not masking any bursts')
    disp(['All bursts (1:',num2str(nBursts),') remain un-masked in ',...
          cdfFileRoot,'-mask.cdf'])
elseif strcmp(answer,'Y') || strcmp(answer,'y')
    first_good = input('Please view the figure and enter the first GOOD burst:  ');
    if isempty(first_good)
        disp('Valid answers are numeric only')
        first_good = input('Please view the figure and enter the first GOOD burst:  ');
    end
    last_good = input('Please view the figure and enter the last GOOD burst:  ');
    if isempty(last_good)
        disp('Valid answers are numeric only')
        last_good = input('Please view the figure and enter the last GOOD burst:  ');
    end
    if isequal(first_good,1)
        disp(' ')
        disp(['Bursts ',num2str(last_good),':',num2str(nBursts),...
              ' will be masked in ',cdfFileRoot,'-mask.cdf'])
        badBursts = [last_good:nBursts]';
    elseif isequal(last_good,nBursts)
        disp(['Bursts 1:',num2str(first_good),' will be masked in ',...
               cdfFileRoot,'-mask.cdf'])
        badBursts = [1:first_good]';
    else
        disp(['Bursts 1:',num2str(first_good),' and bursts ',num2str(last_good), ...
              ':',num2str(nBursts),' will be masked in ',cdfFileRoot,'-mask.cdf'])
        badBursts = [1:first_good last_good:nBursts]';
    end
end
plot(cdf{'burst'}(1:first_good), cdf{'pressure'}(1:first_good),'r')
plot(cdf{'burst'}(last_good:nBursts), cdf{'pressure'}(last_good:nBursts),'r')
legend('good data','masked data')

% Calculate MSL and the tidal range.
press = cdf{'pressure'}(:);
press(badBursts) = nan;
MSL = gmean(press);
tide = gstd(press);
disp(' ')
disp(['The calculated MSL is ',num2str(MSL), ' meters'])
disp(['The calculated tidal range is ',num2str(tide),' meters'])

% Data within approx. 6% of the surface is contaminated by sidelobe
% reflection errors.  We will mask the data in these cells.
% Reference: http://www.nortekusa.com/principles/Doppler.html
cellHeight = cdf{'height'}(:); %cell height above bottom, m
badCells = find( cellHeight > (0.94*(MSL+(tide))));
disp(' ')
if isempty(badCells)
    disp('No cells found within the specified depth range')
else
    disp(['Cells [',num2str(badCells'),'] are above MSL'])
    disp('Data in these cells will be masked')
end

% Mask data, recalculate min/max values, and
% set NaN's to fill
disp(' ')
disp('Masking data ...')
theVars = var(cdf);
if ~isempty(badBursts) || ~isempty(badCells)
    for i = 1:length(theVars),
    
        if ~strcmp(ncnames(theVars{i}),'time') & ...
           ~strcmp(ncnames(theVars{i}),'time2') & ...
           ~strcmp(ncnames(theVars{i}),'burst') & ...
           ~strcmp(ncnames(theVars{i}),'lat') & ...
           ~strcmp(ncnames(theVars{i}),'lon') & ...
           ~strcmp(ncnames(theVars{i}),'height')
            theFillVal = theVars{i}.FillValue_(:);
            data = theVars{i}(:);
            [nR, nC] = size(data);
            fills = find(data >= theFillVal);
            if ~isempty(fills)
                data(fills) = nan;
            end
            if nC > 1
                if ~isempty(badBursts)
                    data(badBursts,:) = nan;
                end
                if ~isempty(badCells)
                    data(:,badCells) = nan;
                end
            else
                if ~isempty(badBursts)
                    data(badBursts) = nan;
                end
            end
            theVars{i}.minimum = gmin(data(:));
            theVars{i}.maximum = gmax(data(:));
            bads = find(isnan(data));
            data(bads) = theFillVal;
            theVars{i}(:) = data;
        end
    end
else
    disp(' ')
    disp(['Not masking any data in the file ',cdfFileRoot,'-mask.cdf'])
end

% Set attribute WATER_DEPTH to MSL
cdf.WATER_DEPTH = ncfloat(MSL);
cdf.WATER_DEPTH_NOTE = ncchar('Water depth is MSL');

% Update the file history
if ~isempty(badBursts) && ~isempty(badCells)
    hist = cdf.history(:);
    hist_new = ['Bad data collected during instrument deployment and recovery, ',...
                'and data collected above MSL masked by mask_argnCurr_cdf.m V ',...
                version,'; ',hist];
    cdf.history = ncchar(hist_new);
    
elseif ~isempty(badBursts)
    hist = cdf.history(:);
    hist_new = ['Bad data collected during instrument deployment and recovery ',...
                'masked by mask_argnCurr_cdf.m V ',version,'; ',hist];
    cdf.history = ncchar(hist_new);
elseif ~isempty(badCells)
    hist = cdf.history(:);
    hist_new = ['Data collected above MSL ',...
                'masked by mask_argnCurr_cdf.m V ',version,'; ',hist];
     cdf.history = ncchar(hist_new);
end

% Update the file creation date
cdf.CREATION_DATE = ncchar(datestr(now,0));

close(cdf);

disp(' ')
disp('Finished masking data')
disp(' ')
disp('NOTE:  Masking data for sidelobes does not necessarily remove all')
disp('       contaminated cells.  The user should load the file ')
disp(['       ',cdfFileRoot,'-mask.cdf and manually look for and remove if'])
disp('       necessary additional contaminated cells')