clear
close all
% Program to estimate PM10 from AOD observation
% Defines thermodynamics profiles
z=[0:10:5000];
mod='night',  % day or night
if strcmp(mod,'day')
 RH0=30;
 T0=20+273.16;
 GRAD=-1;
 q=10;
 T=T0+GRAD*z/100;
elseif strcmp(mod,'night')
  RH0=95;
  T0=20+273.16;
  GRAD=3;
  GRAD0=-0.65;
  Zinv=200;
  T=T0+GRAD*z/100;
  ind=find(z>Zinv);
  T(ind)=T(ind(1))+GRAD0*z(ind)/100;
end
p=barometric(1013,z,T);
for k=1:length(T)
 E(k)=wilg(T(k)-273.16,100,p(k));
end 
q=RH0/100*E(1)/p(1)*622;
e=q*p/622;
RH=e./E*100;
RH(RH>100)=100;
% Find the PBL Heigh -condensation level
if strcmp(mod,'day')
    ind=find(RH==100);
    Zpbl=z(ind(1));
    % aerosol profiles above the PBL top
    H=100;
    EXTsurf=0.1;
elseif strcmp(mod,'night')
    [v,ind]=min(RH);
    Zpbl=0;
    H=100;
    EXTsurf=0.5;
end    
plot(RH,z)
% Aerosol higroscopic growh
GF=0.58;
EXT0=EXTsurf/1000;
EXT=EXT0*(1+RH/30).^GF;
ind=find(z>Zpbl);
EXT(ind)=EXT(ind(1))*exp(-(z(ind)-Zpbl)/H);
figure
% Adding Gaussian aerosol layer in free troposhere
EXTdust=EXT0/2*exp(-1e-5*(z-4000).^2);
%EXT=EXT+EXTdust;
plot(EXT,z)
xlabel('Extinction coefficient [1/km]')
ylabel('Altitude [m]')
% computing AOD
AOD=trapz(z,EXT);
lambda=0.5;
% computing aserosol single scattering properties from MIE theory
ro=2000;
reff=0.2;
[S1,S2,Q]=mie(2*pi*reff/lambda,ref_sulfate(lambda),10);

if strcmp(mod,'night')
 ind=find(EXT<EXT(1)/10);
 Zpbl=z(ind(1));
end 
EXTmean=AOD/Zpbl;
EXTsurf=EXTmean
% Estimation PM10 based on mean extinction in the PBL
PM10=EXTsurf*4/3*reff*ro/Q*1000
% Estimation PM10 based surface extinction
PM10=EXT(1)*4/3*reff*ro/Q*1000

reff=reff*0.8;
[S1,S2,Q]=mie(2*pi*reff/lambda,ref_sulfate(lambda),10);
PM10=EXTsurf*4/3*reff*ro/Q*1000