function [delta,Dter,Dqer] = cool_skin(sal,Tw,rhoa,cpa,Pa, ...
                                        U_star,T_star,Q_star, ...
                                        dlw,dsw,nsw,delta,g,Rgas, ...
                                        CtoK,Qsat_coeff)
% COOL_SKIN: compute the cool-skin parameters.
% COOL_SKIN computes the cool-skin parameters. This code follows 
% the fortran program bulk_v25b.f. For more details, see the cool-skin
% and warm layer paper by Fairall et al (1996), JGR, 101, 1295-1308. 
% All input variables should be vectors (either row or column), except 
% Rgas, CtoK, Qsat_coeff, and g, which can be scalars. Uses some
% functions from CSIRO SEAWATER TOOLBOX.
%
% INPUT:  sal       -  salinity [psu (PSS-78)]
%         Tw        -  water surface temperature [C]
%         rhoa      -  air density [kg/m^3]
%         cpa       -  specific heat capacity of air [J/kg/C] 
%         Pa        -  air pressure [mb]
%         U_star    -  friction velocity including gustiness [m/s]
%         T_star    -  temperature scale [C]
%         Q_star    -  humidity scale [kg/kg]
%         dlw       -  downwelling (INTO water) longwave radiation [W/m^2]
%         dsw       -  measured insolation [W/m^2]
%         nsw       -  net shortwave radiation INTO water [W/m^2]
%         delta     -  cool-skin layer thickness [m]
%         g         -  gravitational constant [m/s^2]
%         Rgas      -  gas constant for dry air [J/kg/K]
%         CtoK      -  conversion factor for deg C to K
%        Qsat_coeff - saturation specific humidity coefficient
%
% OUTPUT: delta     -  cool-skin layer thickness [m]
%         Dter      -  cool-skin temperature difference [C]; positive if 
%                      surface is cooler than bulk (presently no warm skin 
%                      permitted by model)
%         Dqer      -  cool-skin specific humidity difference [kg/kg]
%
% USAGE:  [delta,Dter,Dqer] = cool_skin(sal,Tw,rhoa,cpa,Pa, ...
%                                 U_star,T_star,Q_star, ...
%                                 dlw,dsw,nsw,delta,g,Rgas, ...
%                                 CtoK,Qsat_coeff)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 4/9/99: version 1.2 (contributed by AA)
% 8/5/99: version 2.0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% -> column vectors
sal = sal(:); Tw = Tw(:); rhoa = rhoa(:); cpa = cpa(:); Pa = Pa(:);
U_star = U_star(:); T_star = T_star(:); Q_star = Q_star(:); 
dlw = dlw(:); nsw = nsw(:); delta = delta(:); Rgas = Rgas(:);
CtoK = CtoK(:); Qsat_coeff = Qsat_coeff(:); g = g(:);


size_data = size(Tw);

alpha    = sw_alpha(sal,Tw,0);   % thermal expansion coeff [1/C]
beta_sal = sw_beta(sal,Tw,0);    % saline contraction coeff [1/psu]
cpw      = sw_cp(sal,Tw,0);      % specific heat capacity  [J/kg/C] 
rhow     = sw_dens0(sal,Tw);     % density at atmospheric press [kg/m^3]
viscw    = sw_visc(sal,Tw,0);    % kinematic viscosity of sea-water [m^2/s]
tcondw   = sw_tcond(sal,Tw,0);   % thermal conductivity of sea-water [W/m/K]

% the following are values used for COARE
% alpha    = 2.1e-5*(Tw+3.2).^0.79;% as used for COARE data
% beta_sal = 0.026./sal;           % as used for COARE data
% cpw      = 4000*ones(size(Tw));  % as used for COARE data
% rhow     = 1022*ones(size(Tw));  % as used for COARE data
% viscw    = 1e-6*ones(size(Tw));  % as used for COARE data
% tcondw   = 0.6*ones(size(Tw));   % as used for COARE data

% latent heat of water
Le       = (2.501-0.00237*Tw)*10^6;

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

% a big constant
bigc     = (16.*g.*cpw.*(rhow.*viscw).^3)./(tcondw.^2.*rhoa.^2);

% constant for correction of dq; slope of sat. vap. 
wetc     = 0.622.*Le.*Qs./(Rgas.*(Tw+CtoK).^2);

% compute fluxes out of the ocean (i.e., up = positive)
hsb  = - rhoa.*cpa.*U_star.*T_star;
hlb  = - rhoa.*Le.*U_star.*Q_star;

% net longwave (positive up)
nlw = - lwhf(Tw,dlw,dsw);

% total heat flux out of the water surface
qout = nlw + hsb + hlb;

% compute deltaSc = fc*Sns, see sec. 2.4 (p. 1297-1298) in cool-skin paper
% 3 choices; comment out those that are not used! 
deltaSc  = zeros(size_data);
ipos_nsw = find(nsw > 0);
deltaSc(ipos_nsw) = f_c(delta(ipos_nsw),1).*nsw(ipos_nsw); % Paulson and Simpson (1981)
% deltaSc(ipos_nsw) = f_c(delta(ipos_nsw),2).*nsw(ipos_nsw); % COARE approx. to Paulson
% deltaSc(ipos_nsw) = f_c(delta(ipos_nsw),3).*nsw(ipos_nsw); % Hasse (1971)

qcol = qout - deltaSc;

% initialize
alphaQb   = zeros(size_data);
lamda     = zeros(size_data);
Dter      = zeros(size_data);

ipos_qcol = find(qcol > 0);

% eqn. 17 in cool-skin paper
alphaQb(ipos_qcol) = alpha(ipos_qcol).*qcol(ipos_qcol) + ...
                     sal(ipos_qcol).*beta_sal(ipos_qcol) ...
                     .*hlb(ipos_qcol).*cpw(ipos_qcol)./Le(ipos_qcol);

% eqn. 14 in cool-skin paper
lamda(ipos_qcol)   = 6./(1+(bigc(ipos_qcol).*alphaQb(ipos_qcol) ...
                           ./U_star(ipos_qcol).^4).^0.75).^(1/3);

% eqn. 12 in cool-skin paper
delta(ipos_qcol) = lamda(ipos_qcol).*viscw(ipos_qcol) ...
                   ./(sqrt(rhoa(ipos_qcol)./rhow(ipos_qcol)) ...
                   .*U_star(ipos_qcol));

% eqn. 13 in cool-skin paper
Dter(ipos_qcol)    = qcol(ipos_qcol).*delta(ipos_qcol)./tcondw(ipos_qcol);

Dqer          = wetc.*Dter; 




function fc = f_c(delta,option)
% F_C: computes the absorption coefficient fc. 
% fc=F_C(delta,option) computes the absorption coefficient fc. 
%
% INPUT: delta   -  thickness of cool-skin layer [m]
%        option  -  1  use Paulson and Simpson (1981) data for seawater; 
%                      See also p. 1298 of Fairall et al (1996) JGR, 101, 
%                      cool-skin and warm-layer paper. 
%    
%                   2  use approximation to Paulson as given in 
%                      p. 1298 of Fairall et al (1996) JGR, 101, cool-skin
%                      and warm-layer paper. 
%
%                   3  use fc = const. = 0.19, as suggested by Hasse (1971).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 8/5/99: version 1.2 (contributed by AA)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if option == 1
  % Use Paulson and Simpson data

  % Wavelength bands for the coefficients [um]
  % 1) 0.2-0.6
  % 2) 0.6-0.9
  % 3) 0.9-1.2
  % 4) 1.2-1.5
  % 5) 1.5-1.8
  % 6) 1.8-2.1
  % 7) 2.1-2.4
  % 8) 2.4-2.7
  % 9) 2.7-3.0

  % F_i is the amplitude
  F_i    = [0.237 0.360 0.179 0.087 0.080 0.0246 0.025 0.007 0.0004];
  F_i1   = repmat(F_i,length(delta),1);

  % Gam_i is the absorption length [m]
  Gam_i  = [34.8 2.27 3.15e-2 5.48e-3 8.32e-4 1.26e-4 3.13e-4 7.82e-5 1.44e-5]; 
  Gam_i1 = repmat(Gam_i,length(delta),1);

  delta1 = repmat(delta,1,length(Gam_i));

  % fc is the absorption in the cool-skin layer of thickness delta
  fc = sum(F_i1.*(1-(Gam_i1./delta1).*(1-exp(-delta1./Gam_i1))), 2);

elseif option == 2
  % use COARE approximation to Paulson and Simpson data

  fc = 0.137+11.*delta-(6.6e-5./delta).*(1-exp(-delta/8e-4));

elseif option == 3
  % use Hasse simple approximation

  fc = 0.19;

end