function [VAR_interp,Z_interp]=zinterp(Depth,VAR,zint);

% Input Variables:
% Depth: matrix with profile depth in each column
% VAR: matrix with VAR data for each profile in columns
% zint: vector with new depth coordinates for which interplated data will
% be computed

%% Interpolate VAR profile data to 1m depth
s.depth=Depth;
zint=zint(:);

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

Z_interp=nan(nanmax(zint),nP);
VAR_interp=nan(nanmax(zint),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);
    
    VAR_interp(1:length(zint),p)=interp1(x_mi(~isnan(y_mi)),y_mi(~isnan(y_mi)),zint); %,'nearest','extrap');
    
    Z_interp(1:length(zint),p)=zint;
    
%     % 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