function [data, long_name, units, data_std] = mcgetavg(file, var_name, corner, edges, navg);
%       function [data, long_name, units, data_std] = 
%			mcgetavg(file, var_name, corner, edges, navg);
%
%       Loads data from a netCDF datafile.  If corner and edges
%       are omitted, all data is retrieved.
%
%       data = variable requested by 'name'
%		long_name = the netCDF long name
%		units 	= the units
%       file 	= netCDF data file name and path
%       var_name = variable name to load
%       corner	= starting index, record dimension must be first
%       edges	= length of hyperslab, record dimension must be first
%		navg	= number of points to average together, if desired
% 				averaging will have problems for variables with
%				more than two dimensions
%		data_std = standard deviation if navg is used
%
%       WARNING: this function assumes that corners and edges are correct.
%       the unlimited index always comes first
%
%       Written by Marinna Martini, U.S. Geological Survey
%
mexcdf('SETOPTS',0);

cdfid = mexcdf('OPEN', file, 'NOWRITE');
if cdfid < 0, error(['Cannot open file ',file]); end

%       get global info
[ndims, nvars, natts, recdim, rcode] = mexcdf('INQUIRE', cdfid);
varid = mexcdf('VARID', cdfid, var_name);
if varid < 0, rcode = varid; 
	mexcdf('CLOSE', cdfid);
	error([var_name, ' does not exist in ', file]);	
end

if exist('corner')~=1 | exist('edges')~=1,
	% find out how many dims this variable has
	[junk, type, nvdims, vdims, nvatts, rcode] = mexcdf('VARINQ', cdfid, varid);
    corner = zeros(nvdims,1);
    edges = ones(nvdims,1);
	% set the dimensions
	for i=1:nvdims,
        [junk, edges(i), rcode] = mexcdf('DIMINQ', cdfid, i-1);
    end	
end

fill = mexcdf('ATTGET',cdfid,var_name,'_FillValue');

% grab the whole data set!
if exist('navg')==1,
	nrec = edges(1);
	navgrec = floor((nrec-corner(1))/navg);
	edges(1) = navgrec;
	if length(edges) == 1,
		data = zeros([edges 1]);
		data_std = data;
	else
		data = zeros(edges);
		data_std = data;
	end
	edges(1) = navg;
	for i=1:navgrec,
		corner(1)=(i-1)*navg+1;
		dummy = mexcdf('VARGET', cdfid, var_name, corner, edges);
		flags = find(dummy == fill);
		dummy(flags) = ones(size(flags)).*NaN;
		keyboard
                data(i,:) = mean(dummy');
		data_std(i,:) = std(dummy');
	end
	data=data';
	data_std=data_std';
else
	[data, rcode] = mexcdf('VARGET', cdfid, var_name, corner, edges);
	data_std = [];
end

long_name = mexcdf('ATTGET',cdfid,var_name,'long_name');
units = mexcdf('ATTGET',cdfid,var_name,'units');

mexcdf('CLOSE', cdfid);

