%
%     function [LCL] = calc_LCL(L0,T0,Q0)
%
% Computes the pressure level LCL where
% a parcel lifted adiabtically from L0
% reaches saturation. 
%
% L0 = initial level of teh parcel (1<L0<40) 
% T0 = parcel temperature (in K) at L = L0
% Q0 = parcel specific humidity (in kg/kg) at L = L0
%
% NB: LCL is an interval, meaning the actual
%     LCL is in between the calculated pressure levels.
%
% Algorithm: raise the parcel by DP along a dry adiabat;
%            computes its new vapor pressure ef and temperature Tf.
%            loop until ef>= esat

     function [LCL] = calc_LCL(L0,T0,Q0)


path('/u/u0/czaja/MATLAB',path);

% Constants
  mb = 100;
  EPS = 0.622; %Rv/Rd

% Pressure grid
  DP = 25; %model vertical resolution (in mb)
  P = [987.5:-25:0]; %model pressure levels (in mb)
  Pref = 1000; %ref. pressure (in mb) but not used...

% Dry adiabat
  TDA = dry_adiabat(T0); %dry adiabatic temperature profile
  
% Loop on vertical levels
  w = Q0 / (1-Q0); %mixing ratio
  ef = P(L0) * w / (w + EPS); %vapor pressure in mb
  Tf = T0;
  [esat,xx] = saturation_humidity(T0,Pref);
  if ef >= esat, 
    LCL = L0; 
    return; 
  else
    L = L0;
    while ef<esat
      L = L + 1;
      ef = ef - DP * w / (w+EPS);
      Tf = TDA(L);  
      [esat,xx] = saturation_humidity(Tf,Pref);
    end
    LCL = [L-1 L];
  end