function [dr,rr,vfrr,nsr,massflux,rlev,rm,vfm,massconc,sa,area,cd10n,ustr,h]=drop_source(sourcestrength,Uz,zu,ust)
disp('drop_source');

%Droplet Model code	cwf 6/01/02

%**********************Inputs
%Uz=wind speed (m/s) @ height of zu (m)
%zu=height of wind speed (m)
%ust=u* (m/s) if specified as <0, then ustr is computed from charnock;
	
%********************** computed values
%ustr=u* (m/s)
%z0=roughness (m)
%h=sig waveheight/2 (m)
%droplet fall velocity

%********************** returned variables
%rr             droplet radius
%nsr            dF/dr  #droplets producesd at effective height h/2 in #/m^2/micron/sec
%massflux       total liquid water mass flux g/m^2/sec
%rm             mass-weighted mean radius of droplet production spectrum

	
%********************* physical constants
visa=1.5e-5;%viscosity of air (m^2/s)
visw=1e-6;%viscosity of water (m^2/s)
kon=0.4;%van Karman constant
grav=9.83;%acceleration of gravity
sig_water=7.4e-5;%surface tension of seawater(m^3/s^2)

%*******************  adjustable parameters
%dll=10            'other' length scale associated with spectral transfer of breaking energy
%waveslope=0.1      slope of typical breaking wave
%Tp_frac=0.01       fraction of breaking dissipation available to make droplets
%gust=0.6           fraction of 10-m wind speed realized as gusts at breaking wave top


%***********   set radius variable
dr=2;%set radius increment
rr=2:dr:800;%droplet radii in micrometers
for i=1:length(rr)
    vfrr(i)=wf2(rr(i)/1e6);%fall velocity
end;
mx=length(vfrr);
%*********************************
 

%******************  determine u*

if ust>0
    ustr=ust;
    z0=zu/exp(kon*Uz/ustr);
    
else
    charn=0.018;
    ustr=.04*Uz;
    for i=1:10
        z0=charn*ustr^2/grav;
        cd10n=min(3.5e-3,(kon/log(10/z0)).^2);%cd is not allowed to exceed 3.5e-3
        z0=10./exp(kon./sqrt(cd10n));
        cdzn=(kon/log(zu/z0)).^2;
        ustr=Uz.*sqrt(cdzn);
    end;
    %z0=zu./exp(kon./sqrt(cd10n));
end;
cd10n=(kon./log(10/z0)).^2;
U10=ustr*log(10/z0)/kon;

%************** equilibrium wave properties
a=.018;
b=.729;
twave=b*U10;%                     wave period (s)
hwave=a*U10.^2;%                   wave significant height (m)
lwave=grav/2/pi*twave^2;%       wave length (m)
cwave=grav/2/pi*twave;%         wave phase speed (m/s)

h=hwave/2;%wave top height above mean sea level
ud=ustr/kon.*log(.01*h./z0);%mean wind speed at near the wave top
ustw=sqrt(1.15/1e3)*ustr;%ustar in the water

p=.1*ustw^2*cwave*10;% energy input by breaking of waves (m/s)^3, right behind the breaker, cwave is wave phase speed
%factor of 0.1 is fraction of waves breaking; factor of 10 concentrates breaking energy to 1/10 of the wave
dl=visw/p^.3333;%kolmogorov length scale (m)
dll=10;%other 'appropriate' length scale (m)

gust=0.6;%fraction of 10-m wind speed realized as gusts at breaking wave top
sigud=gust*U10;%assume very large velocity gusts at wave top
waveslope=0.1;%typical slope of a breaking wave
velf=(ud-vfrr/waveslope)/sigud;%normalize gaussian parameter
nlf=(1+erf(velf))/2;%probability of exceeding critical velocity [erf(-inf)=-1; erf(0)=0;erf(inf)=1]

%****************   set fraction of wave breaking dissipation
%Tp_frac=0.6;%fit old Fairall and Andreas models
%Tp_frac=0.06;%fit de Leeuw HEXMAX 22 m/s spectrum
Tp_frac=0.9*sourcestrength;%midway between those two; manifested at a height h/2; referred to as 'f' in the writeup

nd=Tp_frac*3/4/pi/sig_water;

%*****************  source function
nsr=nd./(rr/1e6).^2.*nlf.*p/dll.*exp(-9/4*1.55*(pi*dl./rr/1e-6).^1.33)*1e-6;%#droplets/m^2/sec per micron radius increment
%**********************************

vsr=nsr*4/3*pi.*rr.^3*1e-4;%volume spectrum, micron^2/cm^2/sec.  The factor 1e4 converts from /m^2 to /cm^2 
vol=sum(vsr*dr)*(1e-4)^3;%integrate spectrum to get total vol production cm^3/cm^2/s, 1e-4 converts from micron to cm

massflux=vol*1.00*(1e2)^2;%mass flux g/m^2/s, 1.00 is density water g/cm^3, factor of 1e2 converts from cm to m
rm=(sum(nsr.*rr.^3*dr)/sum(dr*nsr))^0.333;%mass weighted mean radius, micrometers
vfm=(sum(vfrr.*nsr.*rr.^3*dr)/sum(dr*nsr.*rr.^3));%mass weighted mean fall velocity, m/s
massconc=4/3*pi*sum(nsr.*(rr/1e6).^3*dr./vfrr)*1.00e6;%mass concentration below h, gm/m^3
Fp=1+.25*(2*vfrr.*rr/1e6/visa).^.5;%slip factor
mm=find(vfrr./(.4*ustr)<0.5);
rlev=rr(mm(length(mm)));
ll=10;sa=4*pi*2e-5*sum(Fp(ll:mx).*nsr(ll:mx).*rr(ll:mx)/1e6*dr./vfrr(ll:mx));%evap coefficient, m/s
ll=1;area=4*pi*sum(nsr(ll:mx).*(rr(ll:mx)/1e6).^2*dr./vfrr(ll:mx)).*h;%total suspended droplet area/ocean area in the droplet layer
%figure;loglog(rr,nsr*4/3*pi.*rr.^3*1e-4);xlabel('Droplet radius (\mum)');ylabel('Volume Source Strength (# mic^2/cm^2/s)');
%title('Source strength Fairall 2002');axis([1 1e3 1e-1 1e6]);%
