matlab计算EIS，LTS，LCL，potential te

function [EIS,LTS] = EIS_Height(T0,p0,T700)

%Inputs: T0 is surface temp in Kelvin,p0 is surface pressure,i.e.,SLP

%    T700 is air temperature at 700 hPa in Kelvin

%    Output: EIS in Kelvin

%    EIS is a normalized measure of lower tropospheric stability acccounting for

%    temperature-dependence of the moist adiabat.

%    created by jianbida

load constants.mat

% Interpolate to 850 hPa

T850 = (T0+T700)./2;

t850=T850-273.15;

% Lhvap  = 2470000 ;

Lhvap  = 2500000-2323.*(T850-273.15) ;

% Assume 80% relative humidity to compute LCL, appropriate for marine boundary layer

LCL = lifting_condensation_level(T0,80,p0);

% Lower Tropospheric Stability (theta700 - theta0)

LTS = potential_temperature(T700,700)-potential_temperature(T0,p0);

es = 6.11*exp(17.3*T850/(T850+237.3)) ;% surface saturation vapor pressure;

qs  = 0.622*es./(850-es);% mixing ratio in surface

%  Gammam  = -dtheta/dz is the rate of potential temperature decrease along the moist adiabat

%  in K / m

Gammam=(g./cp.*(1.0-(1.0 + Lhvap.*qs./Rd./T850)./(1.0 + Lhvap.^2.*qs./cp./Rv./T850.^2)));

%  Assume exponential decrease of pressure with scale height given by surface temperature

z700 = (Rd.*T0./g).*log(p0/700);

EIS=LTS -Gammam.*(z700 - LCL);

end

function [theta] = potential_temperature(T,p)

% Compute potential temperature for an air parcel.

%  Input:  T is temperature in Kelvin

%          p is pressure in mb or hPa

%  Output: potential temperature in Kelvin.

load constants.mat

theta=T.*(ps./p).^kappa;

end

function [LCL] =lifting_condensation_level(T0,RH,p0)

% Compute the Lifiting Condensation Level (LCL) for a given temperature,pressure and relative humidity

%  Inputs:  T0 is surface temperature in Kelvin,p0 is the surface pressure

%          RH is relative humidity (dimensionless)

%  Output: LCL in meters

load constants.mat

T=T0-273.15;

% g=9.8;

% Rd=287;

Td=243.04.*(log(RH/100)+((17.625.*T)./(243.04+T)))./(17.625-log(RH/100)-((17.625.*T)./(243.04+T))); % dewpoint temperature

P=p0./(((T-Td)./223.15+1).^3.5);  % LCL pressure

TL=T0./((T-Td)./223.15+1);      % LCL temperature

LCL=(Rd.*((T0))./g).*log(p0./P); %LCL height

end