function [theta, thetae, q, rh, rhi, qsat, qsati] = thermo2_td(t,p,td) % thermo_td generates thermodynamic variables from t(oC) p(mb) td(oC) vectors % output [theta, thetae, q, rh, rhi, qsat, qsati] = thermo_td(t,p,td) % in K, K, g/kg, %, %, g/kg, g/kg % Note: Assumes constant specific heat capacity % % changed rsat denominator from p to p-es (3 may 07) % Problem with NaN^-K, so added section to cope with NaNs - May 2007 % % Realised thetae calculation is incorrect (17 Nov 2009) % Changed to used formulae of Bolton (1980), for pseudo-equivalent p.t. % correction of rhi (24 July 2017) %convert t,p,td to SI units tk = t + 273.15; %K p = p*100; %Pa tdk = t +273.15; %K %calculate the potential temperature from temperature array p0=1000*100; %reference pressure Pa R=287; %gas constant cp=1004; %specific heat wrt pressure K= R/cp; theta = NaN*ones(size(p)); a = NaN*ones(size(p)); e1 = find(~isnan(p)); a(e1) = ((p(e1)./p0).^(-K)); theta(e1) = tk(e1).*a(e1); %calculate equivalent potential temperature Lv = 2.5e6; %latent heat of vapourisation eps = 0.622; %Rd/Rv es0 = 0.61e3; %reference saturation vapour pressure Pa tk0 = 273.15; %reference temperature [es esi] = thermo_es(t); %...wrt to water and ice (mb) es = es*100; esi = esi*100; %Pa rsat = eps*es./(p-es); %saturation mixing ratio rsati = eps*esi./(p-esi); %sat mixing ration wrt ice qsat = rsat./(rsat+1); %sat specific humidity qsati = rsati./(rsati+1); %sat specific humidity wrt ice %thetae_wrong = theta.*exp((Lv.*rsat)./(cp.*tk)); %equivalent potential temp % This is incorrect, missing the humidity contribution. %calculate r= mass water vapour/mass dry air & % q= mass water vapour/mass moist air [e ei] = thermo_es(td); %e in mb e = e*100; ei= ei*100; %in Pa r = (eps*e)./(p-e); q = r./(1+r); q = q*1e3; r = r*1e3; qsat = qsat*1e3; qsati = qsati*1e3; % calculate relative humidity w.r.t. water and ice rh = e./es*100; rhi = e./esi*100; % calculate pseudo-equivalent potential temperature, from Bolton, Mon Wea Rev, 1980 % r = is g/kg % Firstly calculate Temp at LCL, note e must be in mb. Tlcl = 2840./(3.5*log(tk)-log(e/100)-4.805) + 55; % eqn 21, Bolton thetae = theta.*exp(((3.376./Tlcl)-0.00254).*r.*(1+0.81*r*1e-3)); % eqn 38, Bolton