function [mlh] = mlh_dropsonde(drop, land, method)
% Calculates Mixing Layer Height (MLH) from dropsonde measurements.
%
% INPUTS:
%   drop   - Dropsonde measurements (e.g., temperature, pressure).
%   land   - Flag for land-based or water-based measurement (e.g., 0 or 1).
%   method - Method used for MLH calculation (e.g., parcel,thetav,qv,calvin) 
%
% OUTPUTS:
%   mlh    - Mixed Layer Height


% -- User defined variables and options --
deltas_Land = 1.0; % [K]
deltau_Land = 0.5; % [K]
dthetar_Land = 4; % overshooting threshold; [K/km]

deltas_Ocean = 0.2; % [K]
deltau_Ocean = 0.15; % [K]
dthetar_Ocean = 0.5; % overshooting threshold; [K/km]

if land==0
    deltas=deltas_Land;
    deltau=deltau_Land;
    dthetar=dthetar_Land;
end
if land ==1
    deltas=deltas_Ocean;
    deltau=deltau_Ocean;
    dthetar=dthetar_Ocean;
end

% -- Data Check --
fldn = fieldnames(drop);
badsonde = false;
%disp(size(drop.Temp))
%disp(size(drop.Time))
%disp(size(drop.Altitude))
badsonde = badsonde || (sum(~isnan(drop.Temp)) < 1);
badsonde = badsonde || max(drop.Altitude) < 1000;
badsonde = badsonde || min(drop.Altitude) > 200;
badsonde = badsonde || (max(drop.Press <= 200));
badsonde = badsonde || range(drop.Temp((drop.Time - drop.Time(1)) < 10)) > 30+273.15; 
badsonde = badsonde || max(drop.Temp) > 50 + 273.15;
badsonde = badsonde || min(drop.Temp) < -90 + 273.15;

if ~badsonde
    % Remove above 4000 m
    ia = drop.Altitude < 4000;
    for i = 1:length(fldn)
        if isnumeric(drop.(fldn{i}))
            drop.(fldn{i}) = drop.(fldn{i})(ia);
        end
    end

    % Subsample to 5mb grid
    press_ss = (floor(min(drop.Press)/5)*5:5:ceil(max(drop.Press)/5)*5)';
    level_ss = discretize(drop.Press,press_ss);
    [~,iss,~]=unique(level_ss,'last'); 

    for i = 1:length(fldn)
        if isnumeric(drop.(fldn{i}))
            drop_ss.(fldn{i}) = drop.(fldn{i})(iss);
        else
            drop_ss.(fldn{i}) = drop.(fldn{i}); % carry the header info
        end
    end

    % Determine Boundary Layer Stability Using Stability Threshold
    dTheta52 = drop.Theta(end-4) - drop.Theta(end-1);

    if dTheta52<-1*deltas 
        bltype = -1; % Convective Boundary Layer'\
        fprintf('Convective Boundary Layer \n');
    elseif dTheta52>=deltas
        bltype = 0; % Stable Boundary Layer
        fprintf('Stable Boundary Layer \n');
    else 
        bltype = 1;% Neutral Boundary Layer
        fprintf('Neutral Boundary Layer \n');

    end

    % -- Rearrange the varaibles from bottom to top of profile
    Altitude_ss = flip(drop_ss.Altitude); 
    Theta_ss = flip(drop_ss.Theta);
    Temp_ss = flip(drop_ss.Temp);
    Press_ss = flip(drop_ss.Press) .* 100; %convert to Pa 
    RH_ss = flip(drop_ss.RH);
    Thetav_ss = flip(drop_ss.Thetav);
    Uwnd_ss = flip(drop_ss.Uwnd);
    Vwnd_ss = flip(drop_ss.Vwnd);
    qv_ss = flip(drop_ss.qv); %unitless
    maxVL=size(Temp_ss,2);    

    % -- Compute the Mixed Layer Height -- 
    switch method
        case 'parcel'
            % Parcel method of Liu and Lang (2010)
            if any(bltype)
                % Convective or Neutral Boundary Layers
                deltaTheta1 = drop_ss.Theta - drop_ss.Theta(end);
                k = find(deltaTheta1 >= deltau,1,'last');
                dTheta = diff(drop_ss.Theta)./diff(drop_ss.Altitude./1000);
                kr = find(dTheta(1:k-1)>=dthetar,1,'last');
                mlh = drop_ss.Altitude(kr);

            else
                % Stable Boundary Layers
                dTheta = diff(drop_ss.Theta)./diff(drop_ss.Altitude./1000);
                mlh=NaN;
            end 
        case 'theta'
            mlh = mixed_layer_height_model(Thetav_ss,Altitude_ss,'thetav');
        case 'ri'
            mlh = mixed_layer_height_model_Ri(Thetav_ss,Uwnd_ss,Vwnd_ss,Altitude_ss);
        case 'qv'    
            mlh = mixed_layer_height_model(qv_ss,Altitude_ss,'qv');
    end
else
    mlh = NaN;
    fprintf(' Badsonde  \n') 


end






















