function [argout1, argout2, argout3, argout4, argout5, argout6] = getcusta2(...
    mode, get_ranges, get_dset, get_vars, get_dset_stride, ...
    get_num_urls, get_georange, get_variables, get_archive, ...
    whichurl)
%
% these are used only locally
%
global urlinfolist urlinfotext pcheckbool old_var old_dset old_ranges depths dectime
%
% explanation:  This function is intended to be used to access any 
% number of 3-D datasets with irregular spacing in the vertical,
% e.g., the Levitus datasets.
%
% get_ranges = [WestLon EastLon   -- units: decimal degrees
%           SouthLat NorthLat -- units: decimal degrees
%           MinDepth MaxDepth -- units: meters
%           MinTime  MaxTime]  -- units: decimal years
%
% Mode may be 'cat' (Catalogue data holdings), 'datasize' (estimate
% size of data request, or 'get' (Get Data).

% The preceding empty line is important.
% $Log$

% $Id$
% klee 08/25/99
% It's modified from getcubic from released 3.0.2. Now it calls multiple files one
% at a time and rearrage and combine same-named variable into one big one, as -k
% in loaddods, but then the url is too long to handle w/ -k.  
% It's also capable of multiple time and depth names. 
% "pixchek" now is an option, since one location for one mooring station, no subset 
% is necessary.
% Added urlinfotext in addtion to urlinfolist, for text variables return.

if exist(get_archive) == 2
  eval(get_archive)
  % fill empty time string 
  if isempty(TimeName)
    TimeName = 'browse_time';
  end
else
  dodsmsg(['Problem reading dataset metadata ' get_archive '.m'])
  return
end

if nargin >= 1
  % initialize return arguments
  argout1 = [];
  argout2 = [];
  argout3 = [];
  argout4 = [];
  argout5 = [];
  argout6 = [];
  
  if strcmp(mode,'cat')
    % save these values
    old_ranges = get_ranges;
    old_dset = get_dset;
    old_var = get_vars;
    %
    %  Need to get the vector of times with values for middle of each day that is
    %  within the range of start and end days. This portion of the program
    %  assumes that there is one field per day exactly in the middle of the day.
    %
    
    % Before call, make sure that the range requested overlaps with the data range.
    
    if (get_ranges(4,1) > TimeRange(2)) | (get_ranges(4,2) < TimeRange(1))
      dodsmsg('Requested range does not overlap data. Respecify')
      return
    end
    eval(['[dectime, urlinfolist, urlinfotext, url, n, pcheckbool] = ', Cat_m_File, ...
	  '(CatalogServer, TimeRange, get_archive, get_ranges, get_variables);'])
    % catch rude catalogue servers that do not return args if empty!
    % D. B. 98/04/13
    if ~(exist('urlinfolist') == 1)
      urlinfolist = [];
      dectime = [];
      argout5 = 0;
      return
    else
      argout5 = n;
    end
    %EACH ROW oF URLINFOLIST IS ORGANIZED LIKE SO
    % [lons lats depth_index depths...]   % modefied from getcubic, since dectimes 
                                          % is global already        klee 08/17/99
    % with the idea that at some point start and end indices will not
    % equal each other ..
    % one argument for each range (x,y,z,t) accepted by browser
    if ~isempty(urlinfolist)
      if ~all(isnan(urlinfolist(:,1))), argout1 = urlinfolist(:,1); end
      if ~all(isnan(urlinfolist(:,2))), argout2 = urlinfolist(:,2); end
      if ~all(isnan(urlinfolist(:,4))), argout3 = urlinfolist(:,4); end
    end
    if ~isempty(dectime), argout4 = dectime; end
    if ~isempty(url), argout6 = url; end
    % this is a kludge to determine the number of URLs required ...
    if isempty(argout5)
      if any(length(argout1) > 0 | length(argout2) > 0 | ...
	    length(argout3) > 0 | length(argout4) > 0)
	%argout5 = max(length(argout1),1) * max(length(argout2),1) * ...
	%    max(length(argout3),1) * max(length(argout4),1);
        % We need only one loaddods call for all the files.
        argout5 = 1;
      else
	argout5 = 0;
      end
    end
    return
    
  elseif strcmp(mode,'datasize')
    %  see if dataset request has changed at all get inventory again.

    if isempty(old_var) | isempty(old_dset) | isempty(old_ranges)
      [x,y,z,t,n,URL] = getcusta2('cat', get_ranges, get_dset, get_vars, ...
	  get_dset_stride, get_num_urls, get_georange, get_variables, ...
	  get_archive);
      get_num_urls = n;
    end
    
    if size(get_vars) ~= size(old_var) | (get_dset ~= old_dset) | ...
	  any(any(get_ranges ~= old_ranges))
      [x,y,z,t,n,URL] = getcusta2('cat',get_ranges, get_dset, get_vars, ...
	  get_dset_stride, get_num_urls, get_georange, get_variables, ...
	  get_archive);
      get_num_urls = n;
    end
    if get_num_urls > 0
      %  The browser always sets the longitude range to go from -180 to 180.  Need
      %  to find out how the data set specifies the longitude range: -180 to 180 or
      %  0 to 360 and set up start and end longitudes accordingly. If the range
      %  selected in the browser is between 0 and 180 it does not matter.
      LL = LonRange(1);
      RL = LonRange(2);
      LonDirection = 1;
      if LL > RL % Data goes from east to west in the data set. Swap LL and RL.
	temp = LonRange(1);
	LL = LonRange(2);
	RL = temp;
	LonDirection = -1;
      end
      startLon = get_ranges(1,1);
      endLon = get_ranges(1,2);
      if endLon < startLon 
	startLon = startLon - 360.0;
      end
      % Make sure that the west end of the range is to the 
      % east of the west end of the data set.
      if startLon < LL
	temp = (1 + floor((LL - startLon) / 360.0)) * 360;
	startLon = startLon + temp;
	endLon = endLon + temp;
      end
      if abs(get_ranges(1,1) - get_ranges(1,2)) >= 360
	if LL < RL
	  startLon = LL;
	  endLon = RL;
	else
	  startLon = RL;
	  endLon = LL;
	end
      else % this is still incorrect if the constraints match up! dab 98/04/14
	% Now examine the east end of the data set.
	if startLon < RL
	  if endLon <= RL
	  else
	    if (endLon-360) > LL
	      %   Set up to make A SECOND CALL for the same URL
	      endLon(2) = endLon - 360;
	      endLon(1) = RL;
	      startLon(2) = LL;
	    else
	      endLon = RL;
	    end
	  end
        else
	  startLon = startLon - 360;
	  endLon = endLon - 360;
	  if endLon <= LL
	    dodsmsg('Longitude range does not overlap the data set. ')
	    return
	  else
	    startLon = LL;
	    % East end of range east of west end of data set, 
	    %now check it against east end.
	    if endLon > RL
	      endLon = RL;
	    end	      
	  end
	end
      end
    else
      % num urls is zero; return
      return
    end
    if ~all(size(startLon) == size(endLon))
      dodsmsg('problem with longitude transformation')
    else
      for i = 1:length(startLon)
	sLon = startLon(i);
	eLon = endLon(i);
	% run a script to transform longitude to pixel number
        if pcheckbool == 1
	  pixcheck
        else
          StartColumn = 0; EndColumn = 0; StartRow = 0; EndRow = 0;
        end
	% accumulate volumes if two calculations were necessary
	if isempty(argout1)
	  argout1 = length(StartRow:get_dset_stride:EndRow) * ...
	      length(StartColumn:get_dset_stride:EndColumn) * ...
	      length(get_vars) * 8 / 1e6 * get_num_urls;
	else
	  argout1 = argout1 + ...
	      length(StartRow:get_dset_stride:EndRow) * ...
	      length(StartColumn:get_dset_stride:EndColumn) * ...
	      length(get_vars) * 8 / 1e6 * get_num_urls;
	end
	VERSION = version; VERSION = str2num(VERSION(1));
	if VERSION < 5
	  % number of depths will be reflected in number of URLs
	else
	  % THIS IS DISABLED UNTIL JAMES FIXES LOADDODS TO HANDLE
	  % CUBIC MATRICES  -- dbyrne 98/11/24
	  %	    if ~isempty(depths)
	  %	      d = find(depths >= get_ranges(3,1) & depths <= get_ranges(3,2));
	  %	      argout1 = argout1 * length(d);
	  %	    end
	end
	argout2 = get_num_urls;
      end
    end
    return
    
  elseif strcmp(mode,'get')
    
    if get_num_urls > 0 % we have a catalogue -- we can just issue the request
      %  The browser always sets the longitude range to go from -180 to 180.  Need
      %  to find out how the data set specifies the longitude range: -180 to 180 or
      %  0 to 360 and set up start and end longitudes accordingly. If the range
      %  selected in the browser is between 0 and 180 it does not matter.
      LL = LonRange(1);
      RL = LonRange(2);
      LonDirection = 1;
      if LL > RL % Data goes from east to west in the data set. Swap LL and RL.
	temp = LonRange(1);
	LL = LonRange(2);
	RL = temp;
	LonDirection = -1;
      end
      startLon = get_ranges(1,1);
      endLon = get_ranges(1,2);
      if endLon < startLon 
	startLon = startLon - 360.0;
      end
      % Make sure that the west end of the range is to the 
      % east of the west end of the data set.
      if startLon < LL
	temp = (1 + floor((LL - startLon) / 360.0)) * 360;
	startLon = startLon + temp;
	endLon = endLon + temp;
      end
      if abs(get_ranges(1,1) - get_ranges(1,2)) >= 360
	if LL < RL
	  startLon = LL;
	  endLon = RL;
	else
	  startLon = RL;
	  endLon = LL;
	end
      else % this is still incorrect if the constraints match up! dab 98/04/14
	% Now examine the east end of the data set.
	if startLon < RL
	  if endLon <= RL
          else
	    if (endLon-360) > LL
	      %   Set up to make A SECOND CALL for the same URL
	      endLon(2) = endLon - 360;
	      endLon(1) = RL;
	      startLon(2) = LL;
	    else
	      endLon = RL;
	    end
	  end
	else
	  startLon = startLon - 360;
	  endLon = endLon - 360;
	  if endLon <= LL
	    dodsmsg(['Longitude range does not overlap the data set. ', ...
	      'No data request posted.'])
	    return
	  else
	    startLon = LL;
	    % East end of range east of west end of data set, 
	    %now check it against east end.
	    if endLon > RL
	      endLon = RL;
	    end	      
	  end
	end
      end
      %  end % of /* if abs(ranges(1,1) - ranges(1,2)) >= 360 */
    end

    if ~all(size(startLon) == size(endLon))
      dodsmsg('problem with longitude transformation')
    else

      % construct lists from which to make URLs
      whichvariables = [];
      for i = 1:size(get_variables,1)
	for j = 1:size(SelectableVariables,1)
	  if strcmp(deblank(get_variables(i,:)), deblank(SelectableVariables(j,:)))
	    whichvariables = [whichvariables j];
	  end
	end
      end
      if ~isempty(whichvariables)
	variablelist = DodsName(whichvariables,:); 
	%% take into account that there may be mutliple servers
	%if size(Server,1) > 1
	%  serverlist = Server(whichvariables,:); 
	%else
	  serverlist = Server; 
	%end
      else
	dodsmsg('Getcusta2: No variables were selected!')
      end
      VERSION = version; VERSION = str2num(VERSION(1));
      % make up a time variable to pass back for those datasets
      % that do not have a time variable.
      if exist('urlinfolist') == 1
	browse_time = dectime(whichurl);
      end
      % loop through as many longitude spans as necessary
      for k = 1:length(startLon)
	sLon = startLon(k);
	eLon = endLon(k);
        if pcheckbool == 1
	  pixcheck
        else
          StartColumn = 0; EndColumn = 0; StartRow = 0; EndRow = 0;
        end
	if size(DodsName,1) ~= size(SelectableVariables,1)
	  dodsmsg([ 'There is a problem in the archive.m file: ', ...
		'the number of DodsNames and SelectableVariables does not', ...
		'match']) 
	  return
	end
	
	% construct and evaluate URL(s)
        if size(serverlist,1) == 1
	  % there will be only one URL.  Initialize it.
	  iURL = '';
	  % return the URL to be displayed by the browser. D.B. 98/4/5

          % Pass whole array of urlinfolist in with serversize == 1, 
          % since we can now access multiple servers and concatenate
          % variables with the same name in loaddods.  klee 08/17/99
	  eval([ 'iURL =  ', URL_m_File, '(get_archive, urlinfolist, urlinfotext, ', ...
		'StartColumn, EndColumn, StartRow, EndRow, ', ...
		'get_dset_stride, get_ranges, variablelist,', ...
		'serverlist, whichurl);'])

	  if ~isempty(iURL)
	    argout5 = iURL;
	    % because of v4 name length bug and the fact we have no info
	    % about the ORDER of returning variables, we might have to 
	    % evaluate some more URLs here:
            snum = findstr(iURL, 'http');
            if size(snum,2) < 2
              if VERSION < 5
                for i = 1:size(variablelist,1)
		  s(i) = size(dods_ddt(dods_dbk(variablelist(i,:))),2);
	        end
	        if any(s > 19)  
                  loaddods('-T',iURL);
		else
		  loaddods(iURL);
		end
              else
                loaddods(iURL);
              end
	    else                
              tmpvarlist = str2mat(TimeName,LongitudeName,LatitudeName, ...
	                           DepthName, variablelist);
              timesize = size(TimeName,1);
              depthat = size(tmpvarlist,1) - size(variablelist,1);
              for k = 1:size(tmpvarlist,1)
                eval(['tmptmp', deblank(tmpvarlist(k,:)), ' = [];'])
              end
disp(size(snum,2))
              for i = 1:size(snum,2)
disp(i)
                %tmpurl = ['''',iURL(snum(i):snum(i+1)-2),''''];
                if i < size(snum,2)
                  tmpurl = [iURL(snum(i):snum(i+1)-2)];
                else
                  tmpurl = [iURL(snum(i):length(iURL))];
                end
                for k = 1:size(tmpvarlist,1)
                  eval([deblank(tmpvarlist(k,:)) ' = '''';'])
                  %eval(['tmptmp', deblank(tmpvarlist(k,:)), ' = '''';'])
                end
                if VERSION < 5
                  for i = 1:size(variablelist,1)
		    s(i) = size(dods_ddt(dods_dbk(variablelist(i,:))),2);
	          end
                  if any(s > 19)
                    loaddods('-T',tmpurl)
                  else
                    loaddods(tmpurl);
                  end
                else
		  loaddods(tmpurl);
                end
                thesize = []; 
                % tmpvar is the current var, tmptmpvar is all of the files
                for k = 1:size(tmpvarlist,1)
                  tmpvar = deblank(tmpvarlist(k,:));
                  eval(['tmp', tmpvar, '= '''';'])
                  if exist(tmpvar)
                    eval(['tmp', tmpvar, '=', tmpvar, ';'])
                  end
                  thesize = [thesize, size(eval(['tmp', tmpvar]),1)];
                end
                bsize = max(thesize);
                for k = 1:size(tmpvarlist,1)
                  tmpvar = deblank(tmpvarlist(k,:));
                  if isempty(eval(['tmp', tmpvar])) 
                    eval(['tmp', tmpvar, '= nan*ones(bsize,1);'])
                  else
                    if k > timesize & k <= depthat
                      eval(['tmp', '= tmp', tmpvar, ';'])
                      eval(['tmp', tmpvar, '= tmp*ones(bsize,1);'])
                    end
                  end
                  eval(['tmptmp', tmpvar, '=[tmptmp', tmpvar, '; tmp', tmpvar, '];'])
                end
	      end % end of loop through multiple file calls
              % put tmptmpvar back to var
              tmptmpdy = invjul(tmptmptime,1);
              for k = 1:size(tmpvarlist,1)
                tmpvar = deblank(tmpvarlist(k,:));
                if strcmp(tmpvar, 'decimal_year')
                  eval([tmpvar ' = tmptmpdy; clear tmptmpdy; clear tmptmpdecimal_year;'])
                else
                  eval([tmpvar '=tmptmp', tmpvar, '; clear tmptmp', tmpvar, '; clear tmp', tmpvar])
                end
              end
            end % end of if size(snum,2) <= 2
	  end % end of if ~isempty(iURL) 
        end % end of if serverlist size == 1
	
disp('end of loaddods calls');

	% now make up a list of all returned variables
	% (above list plus time, longitude, latitude, depth),
	% evaluate their sizes, scale them and convert nulls to NaNs.
	% this has to be a temporary variable list -- 'variablelist'
	% is still in use
	if isempty(LongitudeName)
	  LongitudeName = 'LonVector';
	end
	if isempty(LatitudeName)
	  LatitudeName = 'LatVector';
	end
	tmpvariablelist = str2mat(TimeName,LongitudeName,LatitudeName, ...
	    DepthName, variablelist);
	%tmpwhichvariables = [1 2 3 4 whichvariables+4];
        defaultholder = str2mat('TimeName','LongitudeName','LatitudeName','DepthName');
        for i = 1:4
          tmpvname = eval(deblank(defaultholder(i,:)));
          if ~isempty(tmpvname)
            tmpwhichvariables = [tmpwhichvariables, ...
                                 reshape((i*ones(size(tmpvname,1),1)),1,size(tmpvname,1))];
          else
            tmpwhichvariables = [tmpwhichvariables, i];
          end
        end
        tmpsize = size(tmpwhichvariables,2);
        tmpwhichvariables = [tmpwhichvariables, whichvariables+4];
	for i = 1:size(tmpvariablelist,1)
	  tmpname = dods_ddt(dods_dbk(tmpvariablelist(i,:)));
	  if VERSION < 5 % fix v4 bug
	    if size(tmpname,2) > 19
	      tmpname = tmpname(1:19);
	    end
	  end
          if ~isempty(deblank(tmpname))
	    if exist(tmpname) == 1
	      % make a temporary variable to work with
	      tmpout = eval(tmpname);
              clear(deblank(tmpvariablelist(i,:)));
	      sizes(i,:) = size(tmpout);
	      if exist('DataNull') == 1
		if ~isnan(DataNull(tmpwhichvariables(i)))
		  j = find(tmpout == DataNull(tmpwhichvariables(i)));
		  % this caused a bug in v4 so using conditional:
		  if ~isempty(tmpout) & ~isempty(j)
		    tmpout(j) = NaN*j;
		  end
		else
		  % do nothing because null values are already NaNs
		end
	      else
		% if we have no defined null, do nothing
		% we used to do this!: (bad, bad)
		% j = find(tmpout < -1.0e+33);
		% tmpout(j) = NaN*j;
	      end
	      % stack data in a column and scale it if possible
	      if exist('DataScale') == 1
		if size(DataScale,1) ~= size(SelectableVariables,1)+4
		  dodsmsg(sprintf('%s\n%s', 'Data cannot be auto-scaled.', ...
		  ['Error is in ' get_archive '.m']))
		  argout1 = [argout1; tmpout(:)];
		else
		  if ~isnan(DataScale(tmpwhichvariables(i),1))
		    argout1 = [argout1; DataScale(tmpwhichvariables(i),1) + ...
			  tmpout(:)*DataScale(tmpwhichvariables(i),2)];
		  else % scale exists but is a nan -- don't scale
		    argout1 = [argout1; tmpout(:)];
		  end
		end
	      else % there is no DataScale
		argout1 = [argout1; tmpout(:)];
	      end
	    else % the supposedly returned variable DOES NOT EXIST 
	      % in this workspace.  Put in an empty placeholder.
	      if ~isempty(tmpname)
		eval([tmpname '= [];'])
	      end
	      sizes(i,:) = [0 0];
	    end
	  end
disp(i);
        end
        clear tmpout;

	% now convert variablelist to user-friendly names
	%names = str2mat('Time','Longitude','Latitude', 'Depth', ...
	%    SelectableVariables(whichvariables,:));
        defaultnames = str2mat('Time','Longitude','Latitude','Depth');
        tmptimename = eval(defaultholder(1,:));
        if ~isempty(tmptimename)
          if size(tmptimename,1) == 1, names = [deblank(defaultnames(1,:))];
          else, names = [upper(tmptimename(1,1)), tmptimename(1,2:length(deblank(tmptimename(1,:))))];
                for i = 2:size(tmptimename)
                  names = str2mat(names, [upper(tmptimename(i,1)), ...
                          tmptimename(i,2:length(deblank(tmptimename(i,:))))]);
                end
          end
        end
        for i = 2:4
          tmph = eval(defaultholder(i,:));
          if ~isempty(tmph)
            if size(tmph,1) == 1, names = str2mat(names, deblank(defaultnames(i,:))); 
            else, for i = 1:size(tmph,1)
                    names = str2mat(names, [upper(tmph(i,1)), tmph(i,2:length(deblank(tmph(i,:))))]);
                  end
            end
          else
            names = str2mat(names, '');
          end
        end
        names = str2mat(names, SelectableVariables(whichvariables,:));
	% ------------------------ start producing output arguments
	if length(startLon) > 1
	  % we have to unpack the variables
	  total_size = 0;
	  for j = 1:size(sizes,1);
	    string = sprintf('tmp%i%4.4i = %s;',k, j, ...
		[ 'reshape(argout1((1:sizes(j,1)*sizes(j,2))+',...
		  'total_size),sizes(j,1),sizes(j,2))']);
	    eval(string)
	    total_size = sizes(j,1)*sizes(j,2)+total_size;
	  end
	  argout1 = [];
	  if k == 1
	    URL = argout5;
	  else
	    URL = [URL ' ' argout5];
	    for j = 1:size(sizes,1)
	      % if it's not Time or Latitude, 'combine' it
	      if ~strcmp(deblank(names(j,:)),'Time') & ...
		    ~strcmp(deblank(names(j,:)),'Latitude')
		eval(sprintf('tmp%i%4.4i=combine(tmp%i%4.4i,tmp%i%4.4i);', ...
		    k, j, k-1, j, k, j))
	      end
	      eval(sprintf('sizes(j,:) = size(tmp%i%4.4i);',k, j))
	      eval(sprintf('argout1 = [argout1; tmp%i%4.4i(:)];',k,j))
	    end % end of loop through j = num_args 
	  end
	else
	  URL = argout5;
	end % end of if length(startLon) > 1
      end % end of 'for' loop through length of startLon
      argout2 = sizes;
      argout3 = names;
      argout4 = tmpwhichvariables;
      argout5 = URL;
    end
    return
  end % strcmp(mode) 
end % if nargin >= 1
