function [idx_zref, idx_MLD, ML_mask] = GetMLD(depth, sigma0, MLD_criteria, cast_direction)
% GetMLD - Calculate mixed layer depth from glider profiles
%
% This function calculates the mixed layer depth (MLD) using density 
% threshold methods outlined in Montegut et al. 2004:
% https://doi.org/10.1029/2004JC002378
%
% SYNTAX:
%   [idx_zref, idx_MLD] = GetMLD(depth, sigma0, MLD_criteria)
%
% INPUTS:
%   depth        - 2D array (n×m) of depth values [m]
%   sigma0       - 2D array (n×m) of potential density anomaly [kg/m³]
%   MLD_criteria - String specifying method: 'suga' or 'thompson'
%   cast_direction - calculate MLD for the upcast or downcast direction:
%   'upcast'/'downcast'
%
% OUTPUTS:
%   idx_zref     - Linear indices of reference depth (zref) for each profile
%   idx_MLD      - Linear indices of mixed layer depth for each profile
%   ML_mask      - Logical mask: true for mixed layer, false below MLD
%
%
% METHODS:
%   'suga'     - Suga et al., 2004: zref = 10m, Δσ₀ = 0.125 kg/m³
%   'thompson' - Thomson and Fine 2003: zref = 2.5m, Δσ₀ = 0.03 kg/m³
%   'upper10'  - Upper 10 m
%   'jones'    - Jones et al., 2014: zref = 0, Δσ₀ = 0.03 kg/m³
% EXAMPLE:
%   [idx_zref, idx_MLD] = GetMLD(depth, sigma0, 'suga');
%   MLD_depths = depth(idx_MLD);
%
% REFERENCES:
%   - Montegut et al. 2004, JGR, doi:10.1029/2004JC002378
%   - Suga et al. 2004
%   - Thomson and Fine 2003

% Author: Ben Werb, bwerb@mbari.org
% Date: 6/30/25

% Input validation
if ~ischar(MLD_criteria) && ~isstring(MLD_criteria)
    error('MLD_criteria must be a string: ''suga'' or ''thompson''');
end

% Parse cast direction
switch lower(cast_direction)
    case 'upcast'
        % Upcast: depth decreasing
        ddz = floor([nan(1,width(depth)); diff(depth)]);
        icast = ddz < 0;
        depth(~icast) = NaN; % NaN out the downcast
        sigma0(~icast) = NaN;

    case 'downcast'
        % Downcast: depth increasing
        ddz = floor([nan(1,width(depth)); diff(depth)]);
        icast = ddz > 0;
        depth(~icast) = NaN; % NaN out the downcast
        sigma0(~icast) = NaN;
        
    case 'all'
        % Do nothing

    otherwise
        error("cast_direction must be 'upcast' or 'downcast'");
end

% Set parameters based on chosen method
switch lower(MLD_criteria)
    case 'suga'
        zref = 10;          % Reference depth [m]
        MLD_threshold = 0.125;  % Potential density threshold [kg/m³]
        
    case 'thompson'
        zref = 2.5;         % Reference depth [m]  
        MLD_threshold = 0.03;   % Density threshold [kg/m³]

    case 'upper10'
        zref = 10;         % Reference depth [m]  
        MLD_threshold = 0;      % Potential density threshold [kg/m³]

    case 'jones'
        zref = 0.1;         % Reference depth [m]  
        MLD_threshold = 0.03;      % Potential density threshold [kg/m³]

    otherwise
        error('Unknown MLD_criteria: ''%s''. Use ''suga'' or ''thompson''.', MLD_criteria);
end

% Find indices closest to reference depth for each profile
[~, idx_zref_subscript] = min(abs(depth - zref), [], 1);

% Convert subscript indices to linear indices
idx_zref = sub2ind(size(depth), idx_zref_subscript, 1:length(idx_zref_subscript));

% Get density values at reference depth
zref_sigma0 = sigma0(idx_zref);

% Calculate target density (reference density + threshold)
target_sigma0 = zref_sigma0 + MLD_threshold;

% Find mixed layer depth as depth where density equals target density
% This finds the first depth where density exceeds the threshold
[~, idx_MLD_subscript] = min(abs(sigma0 - target_sigma0), [], 1);

% Convert subscript indices to linear indices  
idx_MLD = sub2ind(size(sigma0), idx_MLD_subscript, 1:length(idx_MLD_subscript));

% Create mixed layer mask (true for mixed layer, false below)
ML_mask = sigma0 > target_sigma0;
end