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	edhill	1.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