% Function to determine the boundary layer height 
% based on Richardson number 

function mlh = mixed_layer_height_model_Ri(thv,u,v,z,variable)

  ind = find(~isnan(u) & ~isnan(v)); % Exclude NaN values
  thve = thv(ind);
  ze = z(ind);
  ue = u(ind);
  ve = v(ind);
  Nz = numel(ze);
  if Nz < 2
      %error('Insufficient data points for computation.');
      mlh = NaN;
      return;
  end
  z_interp=0.5*(ze(1:Nz-1)+ ze(2:Nz));
  thv_interp=0.5*(thve(1:Nz-1)+ thve(2:Nz));
  g = 9.80665;

  %disp(ze)
  %disp(ue)
  %disp(ve)
  % Compute Rib 
  Rib = zeros([numel(ze)-1,1]);
  for k = 1:numel(Rib)
    du = ue(k+1) - ue(1);
    dv = ve(k+1) - ve(1);
    Rib(k) = g*(thve(k+1)-thve(1))*(ze(k+1)-ze(1))/thv_interp(k)/(du^2+dv^2);
  end

  % find PBLH as the height where Rib = 0.2
  validRib = Rib(isfinite(Rib)); % Exclude NaN and Inf (AEOLUS)  values from Rib

  if isempty(validRib)
      mlh = NaN; % All values in Rib are NaN
  else
      %disp(validRib);
      %disp(z_interp(~isnan(Rib)));
      hi = interp1(validRib, z_interp(isfinite(Rib)), 0.2);
      mlh = max(hi);
  end




end






















