function SP=data_merge(DD,kdim,iiview,ii_concat,ii_Jday_ax);
% Call SP=data_merge(DD,kdim,iiview,ii_concat,ii_Jday_ax););
%
% DD            = structure containing data structures to join or string containing data file root;
% kdim		    = dimension over which to join datasets;
% iiview	    = 1 to call view5 on the output SP (default is 0);
% ii_concat     = 0 to interpolate across gaps and non-uniform spacings, 1 to just concatenate;
% ii_Jday_ax    = 1 to convert kdim axis to time in Julian days;
%
% See also data_join, data_interp, data_join_mult
%
%
% COPYRIGHT (c) 2008
% 
% Applied Research Laboratories, The University of Texas at Austin
% (ARL:UT)
% All rights reserved.
% 
% This material may be reproduced by or for the U.S. Government
% only. This notice must appear in all copies of this file and its
% derivatives. Further distribution outside the U.S. Government is
% prohibited without the explicit written authorization from the
% Director of the Environmental Sciences Laboratory, ARL:UT.
% 
% THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND,
% EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO THE
% IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
% PURPOSE, OR NONINFRINGEMENT.
% 
% IN NO EVENT WILL ARL:UT BE LIABLE FOR ANY DAMAGES WHATSOEVER ARISING
% OUT OF THE USE OR INABILITY TO USE THIS SOFTWARE, EVEN IF ARL:UT HAS
% BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.
if nargin<4, ii_concat=0; end
if nargin<5, ii_Jday_ax=0; end
%
if nargin<1, DD='*'; end
if ischar(DD),
    ii_slash=findstr(DD,'\');
    if isempty(ii_slash), ii_slash=0; end
    data_dir=DD(1:ii_slash(end));
    mat_files=dir([DD,'*.mat']); nfile=length(mat_files);
    if nfile==0, disp('No files found that match input string'); SP=[]; return; end
    for jfile=1:nfile,
        disp([int2str(jfile),': ',mat_files(jfile).name])
    end
    iif=str2num(input('Enter file numbers to join ([]=all): ','s'));
    if isempty(iif), iif=[1:nfile]; end
    DD=cell(length(iif),1);
    for jfile=1:length(iif);
        str=['xx=load(''',data_dir,mat_files(iif(jfile)).name,''');']; eval(str);
        fields=fieldnames(xx); str=['DD{jfile}=xx.',fields{1},';']; eval(str);
    end
end
ndata=length(DD); ndim=DD{1}.ndim; 
if nargin<2 | isempty(kdim),
    for jdata=1:ndata,
        str=sprintf('    %s   ',DD{jdata}.labels{1:DD{jdata}.ndim}); 
        disp(['DATA ',int2str(jdata),' axis labels: ',str]);
        for jdim=1:DD{jdata}.ndim,
            lims=DD{jdata}.scales{jdim}([1,end]);
            str=['   Axis limits = ',sprintf('[%1.2f %1.2f] nval = %i',[lims(:)',DD{jdata}.dims(jdim)])]; 
            disp(str);
        end
    end
    kdim=       str2num(input(['Enter dimension number over which to join: '],'s'));
end
if nargin<3,
    iiview=     str2num(input(['          Automatically view output? (1=yes,0=no,[]=yes)   : '],'s'));
end
if isempty(iiview), iiview=1; end
% Interpolate if not same size:
ax_data=[]; delta=zeros(1,ndata);
disp(['Merging over dimension labeled ',DD{1}.labels{kdim}]);
dims_max=DD{1}.dims;
label=lower(DD{1}.labels{kdim});
if ~isempty(findstr(label,'time')) & ~isempty(findstr(label,'day')), ii_time=1; else, ii_time=0; end
for jdata=1:ndata,
    if DD{jdata}.ndim~=ndim, error('All data structures must have same number of dimensions'); end
	label=lower(DD{jdata}.labels{kdim}); scales=DD{jdata}.scales{kdim}(:);
    if ii_time,
        ii=findstr(label,'day'); Jday(jdata)=str2num(label(ii+4:end));
        scales=scales + Jday(jdata)*24;
    end
    ax_data=[ax_data; scales]; 
    indices{jdata}=[length(ax_data)-length(scales)+1:length(ax_data)];
    if DD{jdata}.dims(kdim)>1,
        delta(jdata)=mean(diff(DD{jdata}.scales{kdim}));
    else
        delta(jdata)=1;
    end
    dims_max=max(dims_max,DD{jdata}.dims);
end
if ii_time,
    [Jday0,jdata0]=min(Jday); ax_data=ax_data - Jday0*24;
else
    jdata0=1;
end
ax_sort=sort(ax_data);
delta_avg=median(delta);
if ii_concat==0,
    scale_new=[ax_sort(1):delta_avg:ax_sort(end)];
else
    scale_new=ax_sort; 
end
ldim_new=length(scale_new);
% Compute permutation so that kdim is in first dimension:
dim_order=[1:ndim]; temp=dim_order(kdim); dim_order(kdim)=dim_order(1); dim_order(1)=temp;
SP=DD{jdata0}; SP.dims=dims_max; SP.dims(kdim)=ldim_new; 
SP.scales{kdim}=scale_new(:)';
if isfield(SP,'scales2'), SP.scales2=SP.scales; end
data2=repmat(NaN,SP.dims(dim_order));
if ii_concat==0,
	for jdata=1:ndata,
		scales=ax_data(indices{jdata}); npt=length(scales);
        if npt>0,
            ii_scale_new=find(scale_new>=scales(1) & scale_new<=scales(end));
            data=permute(DD{jdata}.data,dim_order);
            dims=DD{jdata}.dims(dim_order);
            dims_out=[length(ii_scale_new), dims(2:end)];
            indx_full{1}=ii_scale_new;
            for jdim=2:ndim, indx_full{jdim}=[1:dims(jdim)]; end
            if npt>1,
                data2(indx_full{:})=reshape(interp1(scales,data,scale_new(ii_scale_new)),dims_out);
            else
                data2(indx_full{:})=reshape(data,dims_out);
            end
        end
        SP.data_files{jdata}=DD{jdata}.title;
	end
else
	for jdata=1:ndata,
		scales=ax_data(indices{jdata}); npt=length(scales);
        data=permute(DD{jdata}.data,dim_order);
        dims=DD{jdata}.dims(dim_order);
        dims_out=[dims(1), dims(2:end)];
        indx_full{1}=indices{jdata};
        for jdim=2:ndim, indx_full{jdim}=[1:dims(jdim)]; end
        data2(indx_full{:})=reshape(data,dims_out);
        SP.data_files{jdata}=DD{jdata}.title;
	end
end
SP.data=permute(data2,dim_order);
if ldim_new>0 & ii_time & ii_Jday_ax,
    SP.scales{kdim}=Jday0 + SP.scales{kdim}/24;
    SP.labels{kdim}='Time - Julian Day';
    if isfield(SP,'scales2'), SP.scales2{kdim}=SP.scales{kdim}; SP.labels2{kdim}=SP.labels{kdim}; end
end
if iiview, view5(SP); end
