function [N2_mld,MLD_n2]=n2atmld(DENSITY,DEPTH)

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

%%  Purpose:
% Adapted from mld_bmax.m to compute N


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

DENSITY_10m=DENSITY(:);
DEPTH_10m=DEPTH(:);
% 
% DENSITY_10m=rho; % To test code
% DEPTH_10m=dpth;


[nd,np]=size(DENSITY_10m);

% Compute local buoyancy profile
%% Compute N^2 for all stations
% density jump at the base of the mixed layer
np=1;
dRHOdz=nan(nd,np);
drhodz=nan(nd,np); 
MLD_n2=nan(np,1);
N2_mld=nan(np,1);

g=9.81; %m/s^2;
rho_ref=1000; % but maybe should use the midpoint for eachlayer?

% 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);
    % rho_ref=nan(nd,1);
    for d=ind_surf+1:nd
        drdz(d)=(sigma(d)-sigma(ind_surf))./(dpth(d)-dpth(ind_surf));
        % rho_ref(d)=(sigma(d)+sigma(ind_surf))./2; %reference density: average density between end points of layer
    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);
        N2_mld(n)=(g./rho_ref).*yout(indp);
    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
        N2_mld(n)=(g./rho_ref).*yout(indp); 
    end
    
    % Compute MLD using j Holte's code:
    %    [mldinfo]=findmld_m(pres,temp,sal,mldindex,yesplot,floatnumber);
    
%     % Compute N2
%     N2=(g./rho_ref).*drhodz;
end