MITgcm_contrib/arnaud_matlab/calc_LCL.m

%
%     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
    
1	%
2	% function [LCL] = calc_LCL(L0,T0,Q0)
3	%
4	% Computes the pressure level LCL where
5	% a parcel lifted adiabtically from L0
6	% reaches saturation.
7	%
8	% L0 = initial level of teh parcel (1<L0<40)
9	% T0 = parcel temperature (in K) at L = L0
10	% Q0 = parcel specific humidity (in kg/kg) at L = L0
11	%
12	% NB: LCL is an interval, meaning the actual
13	% LCL is in between the calculated pressure levels.
14	%
15	% Algorithm: raise the parcel by DP along a dry adiabat;
16	% computes its new vapor pressure ef and temperature Tf.
17	% loop until ef>= esat
18
19	function [LCL] = calc_LCL(L0,T0,Q0)
20
21
22	path('/u/u0/czaja/MATLAB',path);
23
24	% Constants
25	mb = 100;
26	EPS = 0.622; %Rv/Rd
27
28	% Pressure grid
29	DP = 25; %model vertical resolution (in mb)
30	P = [987.5:-25:0]; %model pressure levels (in mb)
31	Pref = 1000; %ref. pressure (in mb) but not used...
32
33	% Dry adiabat
34	TDA = dry_adiabat(T0); %dry adiabatic temperature profile
35
36	% Loop on vertical levels
37	w = Q0 / (1-Q0); %mixing ratio
38	ef = P(L0) * w / (w + EPS); %vapor pressure in mb
39	Tf = T0;
40	[esat,xx] = saturation_humidity(T0,Pref);
41	if ef >= esat,
42	LCL = L0;
43	return;
44	else
45	L = L0;
46	while ef<esat
47	L = L + 1;
48	ef = ef - DP * w / (w+EPS);
49	Tf = TDA(L);
50	[esat,xx] = saturation_humidity(Tf,Pref);
51	end
52	LCL = [L-1 L];
53	end
54