function [VAR_interp,t_interp]=tinterp(TIME,VAR,tint);

% Input Variables:
% Depth: matrix with profile depth in each column
% VAR: matrix with VAR data for each profile in columns
% tint: vector with new depth coordinates for which interplated data will
% be computed
% Modified to interpolate over time, and allow for variable Time axis ( a
% matrix where the third column for time changes for each profile
%% Interpolate VAR profile data to 1m depth
s.depth=TIME;
zint=tint;

% Predefine VAR dimensions
VAR_interp=nan(size(s.depth));
[nZ,nP]=size(VAR_interp);

t_interp=nan(length(zint(:,1)),nP);
VAR_interp=nan(length(zint(:,1)),nP);

for p=1:nP
    indd=find(~isnan(s.depth(:,p)));
    x_data=s.depth(indd,p);
    y_data=VAR(indd,p);
    
    if nansum(~isnan(y_data(~isnan(y_data))))>=2 % if there is TS data (e.g. downcast)
    
%     idx = diff([-9.99; x_data]) == 0.0; % Find duplicate x-values
%     y=y_data(~idx);
%     x=x_data(~idx);
%     
%     idy = diff([-9.99; y]) == 0.0; % Find duplicate x-values
%     y=y(~idy);
%     x=x(~idy);
% 
%     [y,indu]=unique(y_data); % 
%     x=x_data(indu);
%     
%     % order data to have it monotonically increasing by pressure
%     [x_mi,indmi]=sort(x);
%     y_mi=y(indmi);
    
%     [~,indu]=unique(x_mi);
%     y_mi=y_mi(indu);
%     x_mi=x_mi(indu);

    y_mi=y_data;
    x_mi=x_data;
    
    VAR_interp(1:length(zint(:,p)),p)=interp1(x_mi(~isnan(y_mi)),y_mi(~isnan(y_mi)),zint(:,p)); %,'nearest','extrap');
    
    t_interp(1:length(zint(:,p)),p)=zint(:,p);
    
%     % find max depth and overide extrap with NaN
%     [maxd, indmaxd]=nanmax(s.depth(:,p));
%     
%     CT_interp(indmaxd:end,p)=NaN;
    else
        continue % skip profile data for which there is no TS (e.g. downcast)
    end
   
end