function [varargout] = getnc2(varargin);
%
%  [var, lat, lon, level, time] = getnc2(filin, varn, lims, lev, tim);
%      for 4D variables (multi-level)
%  [var, lat, lon, time] = getnc2(filin, varn, lims, lev, tim);
%      for 3d variables (single level)
%  You can output any or all of the above variables.
%
%  Inputs include:
%  filin           = name of input netcdf file
%  varn (optional) = name of variable to retrieve.  If it is the 
%                    only variable in the file, then it need not be input.
%  lims (optional) = [minlon maxlon minlat maxlat] of output data.  
%                    Default is [0 360 -90 90];
%  lev  (optional) = level to output (may be multiple, staggered). 
%                    Default is lev = 1.  lev = 'all' returns all levels
%  tim  (optional) = time index to output (may be multiple, staggered)
%                    Default is tim = 1.  tim = 'all' returns all levels
%
%  If you do not input one or more of the input arguements (and hence
%  wish to use the default values), then the names of the inputs must
%  match the above input arguements (e.g., don't name your level vector 
%  'level', name it 'lev').  This may be fixed some time, but I don't 
%  currently care to do it.
%
%  Dan Vimont, May 2000
%
% 26 Jan 2005, Camille Li 
% Fix for case where desired variable is specified, but there
% is more than one variable in the netcdf file and the desired
% variable is not the first one
% 
% 8 Mar 2006, Camille Li
% Modify to make work for 2D [lat,lon] fields as well as 
% 3D [tim,lat,lon] and 4D [tim,lev,lat,lon].

if nargchk(1, 5, nargin);
  error('Check input arguments - must have between one and five')
end

%  Get filename
filin = varargin{1};

%  Set defaults
lev = 'all';
tim = 'all';
lims = [0 360 -90 90];
varn = 'all';
inputname={'varn';'lims';'lev ';'tim '};

%  Define inputs
if nargin ~= 5;
  for i = 2:(nargin);
    eval([inputname{i-1} ' = varargin{i};']);
  end
else
  varn = varargin{2};
  lims = varargin{3}; 
  lev = varargin{4};
  tim = varargin{5};
end

if isstr(tim); tim = lower(tim); end;

%  Open NetCDF files

nc = netcdf(filin, 'nowrite');
if isempty(nc); error('Incorrect filename'); end;

  varnames = get_varn(nc);

  if strcmp(varn, 'all');
    if length(varnames) > 1;
      error('Too many NetCDF variables - must input varn');
    end
    varn = name(varnames{1});    
  end

% fix for case where desired variable is specified, but there
% is more than one variable in the netcdf file and the desired
% variable is not the first one
	for i1=1:length(varnames),
		temp=name(varnames{i1});	
		if isequal(temp,varn);
  			ncvarnames = varnames{i1};
		end
	end
  vardims = dim(ncvarnames);

  if (length(vardims)+1 == nargout);  %  This saves some time
    time = nc{name(vardims{1}), 1}(:);
  end

  if length(vardims) == 3;  %  3-D var: tim, lat, lon

    lattot = nc{name(vardims{2}), 1}(:);
    lontot = nc{name(vardims{3}), 1}(:);
    [xk, yk] = keep_var(lims, lontot, lattot); %%%%%%%%%%%%%%%%
if isempty(lontot)==0,
    if strcmp(tim, 'all');
      varout = nc{varn, 1}(:,yk,xk);
    elseif isempty(find(diff(tim) ~= 1));
      varout = nc{varn, 1}(tim,yk,xk);
    else
      varout = NaN*ones(length(tim), length(yk), length(xk));
      for i = 1:length(tim);
        varout(i,:,:) = nc{varn, 1}(tim(i), yk, xk);
      end
    end
else
    if strcmp(tim, 'all');
      varout = nc{varn, 1}(:,:,:);
    elseif isempty(find(diff(tim) ~= 1));
      varout = nc{varn, 1}(tim,:,:);
%    else
%      varout = NaN*ones(length(tim), length(yk), length(xk));
%      for i = 1:length(tim);
%        varout(i,:,:) = nc{varn, 1}(tim(i), yk, xk);
%      end
    end
end


  elseif length(vardims) == 2;  %  2-D var: lat, lon

    lattot = nc{name(vardims{1}), 1}(:);
    lontot = nc{name(vardims{2}), 1}(:);
    [xk, yk] = keep_var(lims, lontot, lattot);
	if isempty(lontot)==0,
      varout = nc{varn, 1}(yk,xk);
	else
      varout = nc{varn, 1}(:,:,:);
	end
    
  elseif length(vardims) == 4;  %  4-D var: tim, lev, lat, lon
    levtot = nc{name(vardims{2}), 1}(:);
	if lev==1; levtot=1; end
    levind=[1:1:length(levtot)];
    if ~strcmp(lev, 'all'); levind = lev; end
    lattot = nc{name(vardims{3}), 1}(:);
    lontot = nc{name(vardims{4}), 1}(:);
    [xk, yk] = keep_var(lims, lontot, lattot);

if isempty(lontot)==0,
    if strcmp(tim, 'all');
		if strcmp(lev, 'all');
			varout = nc{varn, 1}(:,:,yk,xk);
      	elseif isempty(find(diff(levind) ~= 1));
			varout = nc{varn, 1}(:,levind,yk,xk);
      	else
			varout = [];
			for i = 1:length(tim);
	  		vartem = NaN*ones(1, length(levind), length(yk), length(xk));
	  		for j = 1:length(levind);
	    	  vartem(1,j,:,:) = nc{varn, 1}(tim(i), levind(j), yk, xk);
	  		end
	  		varout = [varout; vartem];
			end
      	end
    elseif isempty(find(diff(tim) ~= 1));
      	if strcmp(lev, 'all');
			varout = nc{varn, 1}(tim,:,yk,xk);
      	elseif isempty(find(diff(levind) ~= 1));
			varout = nc{varn, 1}(tim,levind,yk,xk);
      	else
			varout = NaN*ones(length(tim), length(levind), length(yk), ...
			 	length(xk));
			for i = 1:length(tim);
	  			for j = 1:length(levind);
	    	  	  varout(i,j,:,:) = nc{varn, 1}(tim(i), levind(j), yk, xk);
	  			end
			end
      	end
    else
      	if strcmp(lev, 'all');
			varout = NaN*ones(length(tim), length(levtot), length(yk), ...
			  length(xk));
			for i = 1:length(tim);
	  			varout(i,:,:,:) = nc{varn, 1}(tim(i),:,yk,xk);
	end  
      elseif isempty(find(diff(levind) ~= 1));
	varout = NaN*ones(length(tim), length(levind), length(yk), ...
			  length(xk));
	for i = 1:length(tim);
	  varout(i,:,:,:) = nc{varn, 1}(tim(i),levind,yk,xk);
	end
      else
	varout = NaN*ones(length(tim), length(levind), length(yk), ...
			  length(xk));
	for i = 1:length(tim);
	  for j = 1:length(levind);
	    varout(i,j,:,:) = nc{varn, 1}(tim(i), levind(j), yk, xk);
	  end
	end
      end
    end  %  If strcmp(tim, 'all')
else
    if strcmp(tim, 'all');
      	if strcmp(lev, 'all');
			varout = nc{varn, 1}(:,:,:,:);
			size(varout)
      	elseif isempty(find(diff(levind) ~= 1));
			varout = nc{varn, 1}(:,levind,:,:);
%      else
%	varout = [];
%	for i = 1:length(tim);
%	  vartem = NaN*ones(1, length(levind), length(yk), length(xk));
%	  for j = 1:length(levind);
%	    vartem(1,j,:,:) = nc{varn, 1}(tim(i), levind(j), yk, xk);
%	  end
%	  varout = [varout; vartem];
%	end
    	end
    elseif isempty(find(diff(tim) ~= 1));
      	if strcmp(lev, 'all');
			varout = nc{varn, 1}(tim,:,:,:);
      	elseif isempty(find(diff(levind) ~= 1));
			varout = nc{varn, 1}(tim,levind,:,:);
%      else
%	varout = NaN*ones(length(tim), length(levind), length(yk), ...
%			  length(xk));
%	for i = 1:length(tim);
%	  for j = 1:length(levind);
%	    varout(i,j,:,:) = nc{varn, 1}(tim(i), levind(j), yk, xk);
%	  end
%	end
      	end
    %else
    %	if strcmp(lev, 'all');
	%		varout = NaN*ones(length(tim), length(levtot), length(yk), ...
	%		  length(xk));
	%		for i = 1:length(tim);
	%  		  varout(i,:,:,:) = nc{varn, 1}(tim(i),:,yk,xk);
	%		end  
    %  	elseif isempty(find(diff(levind) ~= 1));
	%		varout = NaN*ones(length(tim), length(levind), length(yk), ...
	%		  length(xk));
	%		for i = 1:length(tim);
	%  		  varout(i,:,:,:) = nc{varn, 1}(tim(i),levind,yk,xk);
	%		end
    %  	else
	%		varout = NaN*ones(length(tim), length(levind), length(yk), ...
	%		  length(xk));
	%		for i = 1:length(tim);
	%  			for j = 1:length(levind);
	%    		  varout(i,j,:,:) = nc{varn, 1}(tim(i), levind(j), yk, xk);
	%  			end
	%		end
    %  	end
    end  %  If strcmp(tim, 'all')
end
  end  %  If ndims = 4
  
  mv = nc{varn}.missing_value(:); if isempty(mv); mv = NaN; end;
  ao = nc{varn}.add_offset(:); if isempty(ao); ao = 0; end;
  sf = nc{varn}.scale_factor(:); if isempty(sf); sf = 1; end;
  mv = sf*mv+ao;
  
  if ~isnan(mv);
    varout(varout == mv) = NaN;
  end
  
nc = close(nc);


%  Deal out output variables

if nargout >= 1; varargout{1} = varout; end
if nargout >= 2; varargout{2} = lattot(yk); end;
if nargout >= 3; varargout{3} = lontot(xk); end;

if length(vardims) == 3;  %  3-D var: tim, lat, lon
  if nargout == 4;
    if ~isstr(tim); 
      varargout{4} = time(tim); 
    else
      varargout{4} = time;
    end
  end
elseif length(vardims) == 4;  %  4-D var: time, lev, lat, lon
  if nargout >= 4; varargout{4} = levtot(levind); end;
  if nargout == 5;
    if ~isstr(tim); 
      varargout{5} = time(tim); 
    else
      varargout{5} = time;
    end
  end
end

