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 |
|