function nc_4dvar_meld ( ncfile1, ncfile2, t_eps, output_ncfile )
% NC_4DVAR_MELD:  stiches two 4DVAR obs files together.
%
% The output ncfile will basically look like ncfile1 except for the
% record dimension, which should end up being longer in length.  Any
% variables that are in ncfile2 that are not in ncfile1 will not
% be brought over to the output netcdf file.
%
% 
% Requires:  SNCTOOLS
% Requires:  MEXNC



nc1_info = nc_info ( ncfile1 );
nc2_info = nc_info ( ncfile2 );

%
% Ok, figure out how we are going to sort this.  The key will be
% the first time value of each survey, which is stored in survey_time
survey_time1 = nc_varget ( ncfile1, 'survey_time' );
survey_time2 = nc_varget ( ncfile2, 'survey_time' );
survey_time = [survey_time1; survey_time2];
source = [ones(size(survey_time1)); 2*ones(size(survey_time2))];
[sorted_survey_time, I] = sort ( survey_time );
source = source(I);


%
% Sort the Nobs properly
Nobs1 = nc_varget ( ncfile1, 'Nobs' );
Nobs2 = nc_varget ( ncfile2, 'Nobs' );
Nobs = [Nobs1; Nobs2];
Nobs = Nobs(I);

%
% Sort the starting point for each survey
start1 = cumsum ( [0; Nobs1(1:end-1)] );
start2 = cumsum ( [0; Nobs2(1:end-1)] );
start = [start1; start2];
start = start(I);

%
% Any that are too close together?  We need to know how many.
ind = find(diff(sorted_survey_time) < t_eps);
too_close_together = length(ind);


%
% Some may need to be combined.  That's why we keep track of the
% number that are too close together.
num_surveys = length(survey_time) - too_close_together;


%
% Ok, now create the output file.  Mostly we just copy over the input definitions,
% but with some exceptions.
[ncid, status] = mexnc ( 'create', output_ncfile, nc_clobber_mode );
if status, error ( mexnc('strerror', status ) ), end

[status] = mexnc ( 'close', ncid );
if status, error ( mexnc('strerror', status ) ), end

%
% Go thru the dimensions.
for j = 1:length(nc1_info.Dimension)

	if strcmp ( nc1_info.Dimension(j).Name, 'survey' )

		%
		% The survey length will most likely not be the 
		% same as the first file.
		[dimid, status] = nc_add_dimension ( output_ncfile, ...
			nc1_info.Dimension(j).Name, num_surveys );

	elseif nc1_info.Dimension(j).Record_Dimension
		[dimid, status] = nc_add_dimension ( output_ncfile, ...
			nc1_info.Dimension(j).Name, 0 );
	else
		[dimid, status] = nc_add_dimension ( output_ncfile, ...
			nc1_info.Dimension(j).Name, ...
			nc1_info.Dimension(j).Length );
	end
	if status, error ( mexnc('strerror', status ) ), end
end

%
% Go thru the variables.
for j = 1:length(nc1_info.DataSet)
	var = nc1_info.DataSet(j);
	nc_addvar ( output_ncfile, var );

	%
	% Don't copy the survey time, that comes at the end.
	if strcmp ( var.Name, 'survey_time' )
		continue
	end
	%
	% Don't copy Nobs, that comes at the end.
	if strcmp ( var.Name, 'Nobs' )
		continue
	end

	%
	% If it is not a record variable, then transfer the
	% data as well.
	if ~var.IsUnlimitedVariable
		vardata = nc_varget ( ncfile1, var.Name );
		nc_varput ( output_ncfile, var.Name, vardata );
	end
end

%
% Go thru the Global attributes
for j = 1:length(nc1_info.Attribute)
	nc_attput ( output_ncfile, nc_global, ...
		nc1_info.Attribute(j).Name, ...
		nc1_info.Attribute(j).Value );
end





%
% Ok, now lets go thru and get the surveys in order
for j = 1:num_surveys

	current_survey = get_specified_survey ( source(j), start(j), Nobs(j), ncfile1, ncfile2 );

	%
	% Any other surveys that need to be combined with this one?  
	% Jeezuz this is a pain in the ass.
	too_close = find(abs(sorted_survey_time - sorted_survey_time(j)) < t_eps);
	while (length(too_close) > 1)  % we always get one, the current one
		%
		% The 2nd index is the survey we combine with the current survey
		to_be_combined = too_close(2);

		s_tbc = source(to_be_combined);
		st_tbc = start(to_be_combined);
		n_tbc = Nobs(to_be_combined);
		survey_to_be_combined = get_specified_survey ( s_tbc, st_tbc, n_tbc, ncfile1, ncfile2 );

		fmt = 'Combining survey from first file at survey time %.2f length %d with that from 2nd file at survey time %.2f, length %d\n';
		fprintf ( 1, fmt, sorted_survey_time(j), Nobs(j), sorted_survey_time(to_be_combined), n_tbc );

		current_survey = combine_surveys ( current_survey, survey_to_be_combined );

		%
		% remove the survey we just folded into the current survey
		%survey_time(to_be_combined) = [];
		sorted_survey_time(to_be_combined) = [];
		source(to_be_combined) = [];
		start(to_be_combined) = [];

		%
		% adjust the number of observations first
		Nobs(j) = Nobs(j) + Nobs(to_be_combined);
		Nobs(to_be_combined) = [];

		too_close = find(abs(sorted_survey_time - sorted_survey_time(j)) < t_eps);
	end


	fmt = 'Writing survey at time %.2f of length %d\n';
	fprintf ( 1, fmt, sorted_survey_time(j), length(current_survey.obs_time) );

	%
	% Write it to the (j-1) slot of the output file
	nc_add_recs ( output_ncfile, current_survey, 'datum'  );


end

%
% Write the new survey time, Nobs
nc_varput ( output_ncfile, 'survey_time', sorted_survey_time );
nc_varput ( output_ncfile, 'Nobs', Nobs );


%
% Append to the history, saying what we have done.
blurb = sprintf ( 'Merged %s with %s, t_epsilon of %f', ncfile1, ncfile2, t_eps );
nc_addhist ( output_ncfile, blurb );

return;






% COMBINE_SURVEYS
%
% Combine two surveys into one, making sure that it is sorted.
function survey = combine_surveys ( s1, s2 )

	t1 = getfield ( s1, 'obs_time' );
	t2 = getfield ( s2, 'obs_time' );

	[obs_time, I] = sort ( [t1; t2] );

	%
	% go thru each of the fields.
	fnames = fields ( s1 );
	nfields = length(fnames);
	survey = s1;
	for j = 1:nfields
		t1 = getfield ( s1, fnames{j} );
		t2 = getfield ( s2, fnames{j} );

		thefield = cat ( 1, t1, t2 );

		%
		% And finally put the field back into the survey, and
		% make sure it is sorted properly.
		survey = setfield ( survey, fnames{j}, thefield(I,:) );
	end
return



function current_record = get_specified_survey ( source, start, Nobs, ncfile1, ncfile2 )

	if ( source == 1 )
		source_ncfile = ncfile1;
	else
		source_ncfile = ncfile2;
	end

	%
	% Get the records that should constitute the current survey
	current_record = nc_getbuffer ( source_ncfile, start, Nobs );

return
