function [MLD,VARatMLD]=mld_m(VAR,Depth,threshold,ref);

% Compute MLD using threshold criteria referenced to 10m depth to avoid
% detecting diurnal mixed layer
% June 10, 2016: modified to allow input of threshold
% density threshold (de Boyer Montegut et al 2004) is 0.03 kg/m^3, or temperature threshold is 0.2 C
% bug fix: was using as refdepth always the top line, now MLD>10m
% Interpolating to find depth at which diff is == threshold
% Dec 24, 2016: modified to interpolate only if isempty(max(find(abs(pdenthreshold)<threshold)))~=1

% %test:
% VAR=TEMP;
% Depth=PRES;
% threshold=0.05;

if nargin<4
    ref=10; %reference depth, use 10 by default unless specify as input
end

%%    % Density threshold criteria
    MLD=nan(length(VAR(1,:)),1);
    VARatMLD=nan(length(VAR(1,:)),1);
    % Search for the first level that exceeds the potential density threshold
    for s=1:length(VAR(1,:)); %profile index
        
            % Reference depth (use 10 m instead of the surface value):
            if Depth(1,s)<=ref
                refd=searchclosest_m(Depth(:,s),ref);
            else
                refd=1;
            end
            
         for j = refd:length(VAR(:,1)); %depth level index   
            if abs(VAR(refd,s)-VAR(j,s))>threshold % Density Threshold criteria
                MLD(s,1)=Depth(j,s);
                VARatMLD(s,1)=nanmean([VAR(refd,s) VAR(j,s)]);
                
                if j>=2
                % Interpolate to exactly match the potential density threshold
                % This code was "copied and modified from J Holte
                % This interpolation is particularly good for coarse
                % resolution profiles
                clear pdenseg presseg presinterp pdenthreshold
                presseg = [Depth(j-1,s) Depth(j,s)];
                pdenseg = [VAR(refd,s)-VAR(j-1,s) VAR(refd,s)-VAR(j,s)];
                P = polyfit(presseg,pdenseg,1);
                presinterp = presseg(1):1:presseg(2);
                pdenthreshold = polyval(P,presinterp);
              
                % The potential density threshold MLD value:
                if isempty(max(find(abs(pdenthreshold)<threshold)))~=1
                    MLD(s,1)= presinterp(max(find(abs(pdenthreshold)<threshold)));
                end
                
               end
                
                break
            end
        end
    end
    
% elseif method==2
%     
%     Temp=VAR;
%     %% Temperature threshold criteria
%     MLD_tt=nan(length(Temp(1,:)),1);
%     % Search for the first level that exceeds the temperature threshold
%     for s=1:length(Temp(1,:));
%         for j = 1:length(Temp(:,1));
%             % Reference depth (use 10 m instead of the surface value):
%             if Depth(1,s)<=10
%                 refd=searchclosest(Depth(:,s),10);
%             else
%                 refd=1;
%             end
%             
%             if abs(Temp(refd,s)-Temp(j,s))>threshold % Temperature threshold criteria (de Boyer Montegut, 2004)
%                 MLD_tt(s,1)=Depth(j,s);
%                 break
%             end
%         end
%     end
%     MLD=MLD_tt;
% end
