disp('spry_mass');
clear
%************   set mean conditions  ***********
zu=10;
ts=29;%sea temp
t=27;%air temp
dt=ts-t;
P=1000;%air pressure
s=.8;%saturation ratio
sourcestrength=.3;%1=andreas/fairall and 0.1=deLeeuw Hexmax
%***************************************************
%caclualte a few variables, set constants
Q=qsat([t 1000])*s/1e3;%spec hum, kg/kg
tdk=273.6;
Rgas=287.1;
Le=(2.501-.00237*ts)*1e6;%latent heat vaporization
cpa=1004.67;%air specific heat
rhoa=P*100/(Rgas*(t+tdk)*(1+0.61*Q));%air density
visa=1.326e-5*(1+6.542e-3*t+8.301e-6*t*t-4.84e-9*t*t*t);%air kinematic viscosity
cpw=4.2e3;rhow=1e3;%water specific heat.
gamma=241*17.501./(t+241).^2;%Dq_sat/DT=gamma*qsat(T)
beta=1./(1+1e-3*qsat([t 1000])*Le/cpa*gamma);%wet bulb coefficient

%*   setup wind speed bins from 10 to 50 m/s
u=10:50;
nx=length(u);
%nsr=zeros(1,nx);
massflux=zeros(1,nx);r=zeros(1,nx);vf=zeros(1,nx);conc=zeros(1,nx);sa=zeros(1,nx);area=zeros(1,nx);cd10n=zeros(1,nx);nsr=zeros(nx,400);
usr=zeros(1,nx);h=zeros(1,nx);rlev=zeros(1,nx);

%get droplet flux properties as function wind speed
for i=1:nx
        [dr,rr,vf,nsr(i,:),massflux(i),rlev(i),rm(i),vfm(i),conc(i),sa(i),area(i),cd10n(i),usr(i),h(i)]=drop_source_2(sourcestrength, u(i),zu,-1);
        %nsr(i,:)=nsa;
end;

%****************  compute T and radius time constants

ka=2.411e-2*(1+3.309e-3*t-1.441e-6*t^2);%air thermal conductivity
dw=2.11e-5*((t+273)/273)^1.94;%air water vapor diffusivity
Fp=1+.25*(2*vf.*rr/1e6/visa).^.5;%slip factor (P&K)

tauT=rhow*cpw*(rr/1e6).^2/3/ka./Fp;%temperature time constant
tauR=rhow/rhoa*(rr/1e6).^2/(dw)/(1-s)/(qsat([t 1000])*1e-3)/beta./Fp;%radius time constant

%*****************************************************



%*****************************************************
%compute heat/evaporation 
mx=length(rr);
for i=1:nx
    tauf=h(i)./vf;%fall velocity time constant
    exr=max(1/3000,min(3000,tauf./tauR));
    mm=find(exr>0.5);ll=mm(length(mm));rvap(i)=rr(ll);
    mm=find(3*exr>1);ll2=mm(length(mm));
    efr=(1-exp(-exr));
    ext=max(1/3000,min(3000,tauf./tauT));
    eft=(1-exp(-ext));
    rfri=0.27*(1+1.1/(1-s))^.33*efr+exp(-exr);%andreas radius factor
    ffr=1-exp(-3*exr);
    saf(i)=4*pi*dw*sum(Fp(ll2:mx).*nsr(i,ll2:mx).*rr(ll2:mx)/1e6*dr./vf(ll2:mx));%evap coefficient, m/s, as per fairll et al 94 (large droplet part of Andreas)
    sam(i)=4/3*pi*rhow*sum((1-rfri.^3).*nsr(i,:).*(rr/1e6).^3.*dr);%filtered mass flux kg/m^2/s as per Andreas
    fam(i)=4/3*pi*rhow*sum(nsr(i,1:ll).*(rr(1:ll)/1e6).^3*dr);%small droplet part of Andreas
    fat(i)=4/3*pi*rhow*sum(nsr(i,:).*eft.*(rr/1e6).^3*dr);%t-weighted mass flux
end;
wb=3.8e-6*u.^3.4;%white cap fraaction
sfcc=(1-0.27*(1+1/(1-s))^.33);%large droplet (slow) expansion correction factor 
sfcc2=1-0.0214*(1+1/1-s);%small droplet (fast) expansion correction factor
sfcb=h.*rhoa.*qsat([t 1000])*beta*(1-s);%fairall to andreas scaling factor

figure;loglog(u,saf./wb,u,fam*1e3/sfcc2./sfcb./wb,'x',u,sam*1e3/sfcc./sfcb./wb,u,fam*1e3/sfcc2./sfcb./wb+saf./wb,'^',u,sa./wb,'.');xlabel('10-m Wind speed (m/s)');ylabel('Sa, Fairall Velocity Parameter (m/s)');Title('bluline - large drop integral,  x-small drop inte, redline-full integral');
figure;plot(u,rlev,u,rvap);xlabel('10-m Wind speed (m/s)');ylabel('Droplet Radius)');Title('bluline - Toba turbulent levitation, grenline-full evaporation size cutoff');

%figure;loglog(u,cpw*fat*(dt+beta*2.5*qsat([t 1000])*(1-s)),u,rhoa*cpa*1e-3*(dt-.3)*u,'x',u,sam*Le,u,rhoa*Le*1.0e-3*u*1e-3*(.98*qsat([ts-.3 1000])-qsat([t 1000])*s),'o');;xlabel('10-m Wind speed (m/s)');ylabel('Heat Flux (W/m^2)');Title('x-Hs, o-Hl, redline - Ql, bluline - Qs');



%plot figure similar to Andreas Hl versus Ql
uj=[10 12 15 16 18 20 22 25 27 30 32];
ix=length(uj);
for i=1:ix
ij(i)=find(u==uj(i));
end;
figure;loglog(rhoa*Le*1e-3*u(ij)*1e-3*(.98*qsat([ts-0.3 1000])-qsat([t 1000])*s),sam(ij)*Le,'o',[.1 1e3],.1*[.1 1e3],[.1 1e3],[.1 1e3],[.1 1e3],10*[.1 1e3]);axis([100 1000 1e-2 1e3]);
xlabel('Turbulent latent heat flux (W/m^2)');ylabel('Droplet latent heat flux (W/m^2)');Title('Hl vs Ql, windspeed=10 12 15 16 18 20 22 25 27 30 32 m/s');
rfk=[2 4 6 8 10  20  30  50  75   100  150 250   400];
sfk=[3 3 5 9 15 120 350 700 2000 2200 2300 2800 1000]*1e3; %old fairall volume source, mic^3/cm^2/mic/s
%************  Andreas 98 source function
im=find(u==30);ii=im;%select wind speed for comparion
r80=.518.*rr.^.976;
ro=(rr/.518).^(1/.976);
f80=.506./ro.^.024;
zz=fdrop_an(u(ii),r80);%andreas source strength at r assuming droplets have negligible evaporation
%andreas in #/m^2/s/mic
%***************   compare droplet source funcitions at 30 m/s
figure;loglog(rfk,sfk*(1e-6)^3/(1e-2)^2*wb(ii),rr,4/3*pi*(rr/1e6).^3.*nsr(ii,:),rr,zz*4*pi/3.*(rr/1e6).^3);axis([1 1e3 1e-14 1e-8]);;xlabel('Drop Radius  (\mum)');ylabel('Droplet source volume flux (m^3/m^2/s/\mum)');Title('redline - Ql, bluline - Fairall, greenline-wave model, redline-Andreas98');



%***************************   feedback parameterization
wetdep=(1-s)*(1-beta)/gamma;%wet bulb depression
qll=sam*Le;%droplet latent heat flux, no feedback, from full integral
hll=rhoa*Le*1.0e-3*u*1e-3*(.98*qsat([ts-0.3 1000])-qsat([t 1000])*s);%direct turbulent moisture flux, no feedback
hss=rhoa*cpa*1.0e-3*(dt-0.3)*u;%direct turbulent sensible, no feedback (0.3 is coolskin effect)
hsswet=rhoa*cpa*1.0e-3*(dt-0.3+wetdep)*u;%maximum direct turbulent sensible allowed
qqs=cpw*fat*(dt+wetdep);%droplet sensible flux
hs_eps=0.5*rhoa*usr.^3/.4.*(log(h./10)+.4*u./usr);%dissipation heating in droplet layer, factor of 1/2 for energy lost to waves, etc

ch=1e-3*(hss+hll+qqs)./(hss+hll);%'effective' enthalpy bulk transfer coefficient
figure;plot(u,cd10n,u,ch);axis([0 50 0 4e-3]);xlabel('Wind speed (m/s)');ylabel('C_{d10n} & C_{h10n}');


%******************  plot heat flux estimates and droplet with/without feedback
qll0=sourcestrength*0.08./(1+7*exp(-.04*u.^1.33)).*wb.*h.*rhoa*Le*1e-3.*qsat([t 1000])*beta.*(1-s).*(1-.27.*(1+1./(1-s)).^.33);%estimate of qll from parameterization
feed=qll0./(qll0+hsswet+qqs+hs_eps);%feedback factor
dtf=wetdep*feed;%change in air temperature caused by droplet evaporation
%dtf=0;
sp=s+dtf*gamma/(1-beta);%saturation ratio, after feedback efffects
%************   droplet latent heat flux, with feedback
qllp=sourcestrength*0.08./(1+7*exp(-.04*u.^1.33)).*wb.*h.*rhoa*Le*1e-3.*qsat([t 1000])*beta.*(1-sp).*(1-.27.*(1+1./(1-sp)).^.33);
%******************************************************
figure;plot(u,feed);axis([0 50 0 1]);xlabel('Wind speed (m/s)');ylabel('Feedback Coefficient');


figure;loglog(u,qqs,u,hss,'x',u,qll,u,hll,'o',u,qllp,'-^',u,hs_eps,'.');xlabel('10-m Wind speed (m/s)');ylabel('Heat Flux (W/m^2)');Title('x-Hs, o-Hl, redline - Ql, bluline - Qs, triangle-Qll');

for i=1:length(u)
    [hsd, hld, hss_x, hll_x, feed_x, massf,uf]=spray_param(sourcestrength,u(i),zu,ts,t,s,P,-999,-999,-999);
    hsdrp(i)=hsd;
    hldrp(i)=hld;
    hssd(i)=hss_x;
    hlld(i)=hll_x;
    feed_p(i)=feed_x;
    massflux_p(i)=massf;    
    uf_p(i)=uf;
end;
