function [MLD_n2]=mld_bmax(DENSITY,DEPTH)

% close all
addpath(genpath('/glade/u/home/magdac/MATLABhelp_new'));

%%  Purpose:
% Compute MLD as in Large et al. 1997, ~ maxima in buoyancy frequency
% (page 2427)

% DEPTH units='meters';
% DENSITY units-'kg/m^3'; % Potential density, p_ref=0 dbar 

DENSITY_10m=DENSITY(:);
DEPTH_10m=DEPTH(:);
% 
% DENSITY_10m=D.RHO(:,p); % To test code
% DEPTH_10m=D.DEPTH(:,p);

[nd,np]=size(DENSITY_10m);

% Compute local buoyancy profile
%% Compute N2 for all stations

np=1;
N2=nan(nd,np);
dRHOdz=nan(nd,np);
drhodz=nan(nd,np);
MLD_n2=nan(np,1);


% for n=1:np
% display(['Profile # ',num2str(n),'/',num2str(np)]);
for n=1;

%     indnan=find(isnan(DEPTH_10m(:,n)),1,'First');
%     nd=indnan-1;
    
    dpth=DEPTH_10m(1:nd,n);
    sigma=DENSITY_10m(1:nd,n);
    
    % Find surface value (some profiles may not have acquired
    % measurements at the surface):
    ind_surf=searchclosest_m(dpth,0);
    
    sd=50;
    if dpth(ind_surf)>sd % if shallowest depth is > 20m,
        clear ind_surf     % consider there is no surface value, skip profile
        display(['shallowest depth is greater than ',num2str(sd),' m']);
        MLD_n2(n)=NaN;
        break
    end
    
    % Compute drho./dz for each depth against the surface
    drdz=nan(nd,1);
    for d=ind_surf+1:nd
        drdz(d)=(sigma(d)-sigma(ind_surf))./(dpth(d)-dpth(ind_surf));
    end
    dRHOdz=drdz;
    
    % Compute drho./dz locally
    drhodz=nan(nd,1);
    for d=ind_surf+1:nd
        drhodz(d)=(sigma(d)-sigma(d-1))./(dpth(d)-dpth(d-1));
    end
    
    % Compute MLD, i.e.,
    % 1) Search for the MAXIMUM in dRHOdz,
    % 2) find the depth at which the maximum value is equal to drhodz
    [PKS,LOCS]=nanmax(dRHOdz); % findpeaks(dRHOdz(:,n),depth); % findpeaks.m requires MATLAB signal processing toolbox
    
    % find intersection between max value and drhodz curve
    t=dpth; y=drhodz; y2=ones(length(dpth),1)*PKS(1);
    [zout,yout] = intersections(t,y,t,y2,1);
    
    if isempty(zout)~=1
        [MLD_n2(n),indp]=nanmin(zout);
    elseif length(zout)>=2 && zout(indp)<20 % if the first value occurs at the surface,
        [MLD_n2(n),indp]=nanmin(zout(2:end)); % use the subsequent value
    end
    
    % Compute MLD using j Holte's code:
    %    [mldinfo]=findmld_m(pres,temp,sal,mldindex,yesplot,floatnumber);
    
end