function [tau_rain,heat_rain] = rain_flux(Ta,Pa,rh,rain,Ts,sal,u,zu)

% RAIN_FLUX: computes heat and momentum flux due to rain.

% RAIN_FLUX computes heat flux and momentum flux due to rain. This 

% code follows the Fortran program bulk_v25b.f. For more details, 

% see Fairall et al. (1996), JGR, 101, 3751-3752.

%

% INPUT:  Ta        -  air temperature          [C]

%         Pa        -  air pressure             [mb]

%         rh        -  relative humidity        [%]

%         rain      -  rain rate                [mm/hr]

%         Ts        -  sea surface temperature  [C]

%         sal       -  salinity                 [psu (PSS-78)]

%         u         -  wind speed               [m/s]

%         zu        -  wind measurement height  [m]

%

% OUTPUT: tau_rain  -  momentum flux of rainfall            [N/m^2]

%         heat_rain -  heat flux of rainfall (OUT of ocean) [W/m^2]

%

% USAGE: [tau_rain,heat_rain] = RAIN_FLUX(Ta,Pa,rh,rain,Ts,sal,u,zu) 

  

% NOTE: All input variables should be vectors (either row or column), zu

%       may also be a fixed scalar. Output variables are column vectors.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 4/3/99: version 1.2 (contributed by AA)

% 8/5/99: version 2.0

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% -> column vectors

Ta = Ta(:); Pa = Pa(:); rh = rh(:); rain = rain(:); Ts = Ts(:);

sal = sal(:); u = u(:); zu = zu(:);



% get constants

as_consts;

cpa = cp;





o61  = 1/eps_air - 1;              % ~0.61 (moisture correction for temp.)

Qa   = 0.01*rh.*qsat(Ta,Pa);       % specific humidity of air [kg/kg]

T    = Ta + CtoK;                  % C -> K

Tv   = T.*(1 + o61*Qa);            % air virtual temperature

rhoa = (100*Pa)./(gas_const_R*Tv); % air density

Le   = (2.501-0.00237*Ts)*1e6;     % latent heat of vaporization at Ts

Qs   = Qsat_coeff*qsat(Ts,Pa);     % saturation specific humidity





% compute heat flux of rainfall OUT of ocean

dwat   = 2.11e-5*(T./CtoK).^1.94;    % water vapour diffusivity

dtmp   = (1+3.309e-3*Ta-1.44e-6*Ta.^2)*0.02411./(rhoa.*cpa); % heat diffusivity

dqs_dt = Qa.*Le./(gas_const_R*T.^2); % Clausius-Clapeyron

alfac  = 1./(1+0.622*(dqs_dt.*Le.*dwat)./(cpa.*dtmp)); % wet bulb factor

cpw    = sw_cp(sal,Ts,0);            % heat capacity of sea water 



heat_rain = rain.*alfac.*cpw.*((Ts-Ta)+(Qs-Qa).*Le./cpa)/3600;



% compute momentum flux of rainfall

[cd10,u10] = cdntc(u,zu,Ta);% use Smith's formula to compute wind speed at 10m

tau_rain   = rain.*u10/3600;



