%  1.  NAME
%
%  adcp_tide - Interpolate tide data to ADCP positions and time
%
%  2.  USAGE
%
%  adcp_tide
%
%  3.  DESCRIPTION
%
%  This Matlab script interpolates first all U and V tide components at
%  ADCP positions for all times. These are then interpolated at ADCP
%  ensemble times. The output file can be used to generate a vector plot
%  of the tidal current at ADCP ensembles, and to detide ADCP currents.
%
%  Tidal data are read from five different input files:
%
%  1. Time in decimal days (suffix .t)
%
%  2. Longitude in degrees (suffix .lon)
%
%  3. Latitude in degrees (suffix .lat)
%
%  4. Eastward tidal component (suffix .u)
%
%  5. Northward tidal component (suffix .v)
%
%  ADCP ensemble times can be read from the output of the getnav utility.
%  See below for the input .
%
%  4.  PARAMETERS
%
%  The following parameters must be edited before running the script:
%
%     year_base
%        Year base corresponding to input decimal days
%
%     tide_file
%        Root name of ASCII files containing tide data. Extensions will
%        be added automatically.
%
%     adcp_file
%        Name of a file with ADCP time,longitude and latitude as output
%        by utility getnav.
%
%     out_file
%        Name of output file. See below for its .
%
%  5.  INPUT FILE FORMAT
%
%  Input data have only one column with one data for each grid point. The
%  grid defines a two dimensional geographical area with associated
%  longitude, latitude, eastward and northward currents values at each
%  point. These data must be listed in the same order in each file. If
%  several times are listed in the time file, current values must repeat
%  themselves in the corresponding files in the same order as well. All
%  of these files can be generated with the xtract(1) utility.
%
%  6.  OUTPUT FILE FORMAT
%
%  The generated file has the following columns:
%
%  1. ADCP time
%
%  2. Longitude
%
%  3. Latitude
%
%  4. Eastward tidal component
%
%  5. Northward tidal component
%
%  7.  SEE ALSO
%
%  xtract(1)
%
%  8.  AUTHORS
%
%  Pierre Jaccard, Geophysical Institute, University of Bergen, 1999
%
       
% Edit the following parameters:       
%------------------------------------------------------------------

year_base = 1996;
tide_file = 'tide2';
adcp_file = 'adcp.nav'; % '../../../../adcp/nav/skag/001/posfix.smn';
out_file  = 'test.tid';

%------------------------------------------------------------------
% No more editing !!!

%--- Read ADCP times and positions:
adcp = tonan(read_asc(adcp_file));
adcp_t = 1; adcp_x = 2; adcp_y = 3;
n_adcp = size(adcp,1);

%--- Read time and positions of tides:
tide_t = tonan(read_asc([ tide_file '.t'   ], 1));
tide_x = tonan(read_asc([ tide_file '.lon' ], 1));
tide_y = tonan(read_asc([ tide_file '.lat' ], 1));
nx = length(tide_x);
nt = length(tide_t);
if nx ~= length(tide_y),
	error('X and Y vectors do not have the same amount of data');
end

%--- Define the grid for adcp positions:
grid{1} = adcp(:,adcp_x);
grid{2} = adcp(:,adcp_y);

%--- U and V must be processed separatly because of memory problems:
uv_files = strvcat( [ tide_file '.u' ], [ tide_file '.v' ]);
uv = NaN * ones(n_adcp, 2, nt);

for f=1:2,
	disp([ 'Reading file ' uv_files(f,:) ]);
	
	%--- Read data
	data = tonan(read_asc(uv_files(f,:), 1));
	data = reshape(data, nx, nt);
	
	%--- Loop through the time profiles:
	for t=1:nt,
		disp([ 'DNMI profile time: ' to_date(year_base, tide_t(t)) ]);
		
		%--- Build a contour array with the following columns:
		%
		%    1. longitude
		%    2. latitude
		%    3. data
		%
		%    Grid points that do not conatin any good value are removed.

		igood = find(~isnan(data(:,t)));
		new = NaN * ones(length(igood), 3);
		new = cat(2, tide_x(igood), tide_y(igood), data(igood,t));
		
		%--- Regrid each profile at ADCP positions:
		tmp = regrid2(new, 1, 2, grid, 'cubic');
		uv(:,f,t) = diag(tmp(:,:,3));
	end
end
clear new;

%--- Regrid at ADCP times:
data = NaN * ones(n_adcp,2);
for n=1:n_adcp,
	data(n,:) = interp1(tide_t, transpose(squeeze(uv(n,:,:))), ...
			adcp(n,adcp_t), 'cubic');
end

%--- Print the results
fid = fopen(out_file, 'w');
for n = 1:n_adcp,
	line = sprintf(' %15.8f %15.8f %15.8f %15.8f %15.8f', ...
			adcp(n,adcp_x), adcp(n,adcp_y), data(n,1), data(n,2), adcp(n,adcp_t));
	line = strrep(line, ' NaN', '1E38');
	fprintf(fid, '%s\n', line);
end
fclose(fid);


	



