% Function to determine the boundary layer height 
% using model PM2.5 concentrations data 
% method: find the level with minimum of local vertical derivative along profile 
% variable: a string of the variable used to find PBLH. Use any one of below: 
%           'thetav': virtual potential temperature 
%           'pm25':   PM2.5 concentration 


function mlh = mixed_layer_height_model(f,z,variable)

  ind = find(~isnan(f));
  fe = f(ind);
  ze = z(ind); 

   % Compute derivative
   ddz = diff(fe,1)./diff(ze,1);
   Nz = numel(fe);
   fe_interp=0.5*(fe(1:Nz-1)+ fe(2:Nz));
   z_interp=0.5*(ze(1:Nz-1)+ ze(2:Nz));

if strcmp(variable,'pm25')
  % find highest local minimum gradient 
  if numel(ddz)>=3 
    [lm,inds] = findpeaks(-1*ddz);
    k = find(lm > 2.5e-3);
    if ~isempty(k)
      ind = inds(k(end));
      if ddz(ind) < -3.5e-3
      % find the highest level where gradient is smaller than a threshold 
        ind2 = find(ddz<-3.0e-3);
        mlh = z_interp(ind2(end));
      else 
        mlh = z_interp(ind);
      end
    else 
      [mm,ind]=min(ddz);
      if mm<-3.5e-3
        ind2 = find(ddz<-3.0e-3);
        mlh = z_interp(ind2(end));
      else
        mlh = z_interp(ind);
      end 
    end 
  else
       mlh = NaN; 
                
  end 
 

elseif strcmp(variable,'qv')
   % find minimum gradient
   [mm,ind]=min(ddz);
   %disp(ind)
   mlh = z_interp(ind);

elseif strcmp(variable,'thetav') 
  % find lowest hevel with vertical gradient > a threshold 
  inds = find(ddz>1.3e-3); 
  mlh = z_interp(inds(1)); 
end 


