function spice = spice_tailleux(SA,CT,p,SAref,CTref) 
% FUNCTION spice = spice_tailleux(SA,CT,p,SAref,CTref,tau_ref)
%
% PURPOSE: compute a spiciness variable that is similar in nature to the
% spiciness variable originally constructed by Jackett and McDougall (1985)
% and subsequently by Flament (2002) and McDougall and Krysik (2005)
%
% Remi Tailleux, University of Reading, 7 January 2020
%--------------------------------------------------------------------------
% Resize the reference CT and SA values to make them arrays same size as
% input arrays
rho00 = 1000.; 
CT_ref = CTref*ones(size(CT));
SA_ref = SAref*ones(size(SA));

CT_zero = 0.*ones(size(CT));
SA_zero = 35.16504*ones(size(SA)); 

% Compute the specific volumes 
rhoSA = gsw_rho_CT_exact(SA,CT_ref,p);
rhoCT = gsw_rho_CT_exact(SA_ref,CT,p); 

rhoSA_zero = gsw_rho_CT_exact(SA_zero,CT_ref,p);
rhoCT_zero = gsw_rho_CT_exact(SA_ref,CT_zero,p); 

rhoSA = rhoSA./rhoSA_zero;
rhoCT = rhoCT./rhoCT_zero; 

% Reference value for density
rho_ref = gsw_rho_CT_exact(SA_ref,CT_ref,p);

% Compute the spiciness to have the value tau_ref at the reference point 
spice = rho00.*log(rhoSA./rhoCT); 
%spice = rho_ref.*log(rhoSA./rhoCT); 

end

