function [S]=salinity(Sguess, T, P, rho, deltaS, eps, iter);

%
% SALINITY:  Compute salinity from T-S relationship.
%
% [S]=salinity(Sguess, T, Z, rho, iter)
%
% Given a surface of constant density (rho), holding pressure (P)
% constant, this function uses a Newton-Raphson iterative scheme
% to compute salinity (S) as a function of potential temperature
% from a first salinity estimate (Sguess):
%
%      d(S)/d(T) = - [d(rho)/d(T)] / [d(rho)/d(S)]
%                = - d(rho)/d(T) d(S)/d(rho)
%
% This function requires a seawater equation of state (EOS):
%
%            den = eos(S, T, P)
%
% ROMS equation of state for seawater is used here.  Recall that
% ROMS EOS assumes no pressure variation along geopotentials, so
% pressure (dbar) and depth (meters; negative) are interchangeable.
%
% On Input:
%
%    Sguess      First guess salinity (from reference T-S profile)
%    T           Potential temperature (Celsius; from model or data)
%    P           Pressure (dbar) or depth (meters, negative).
%    rho         Density anomaly (kg/m3; from model or data)
%    deltaS      Salinity increment (0.005; OPTIONAL)
%    iter        Number or Netwton-Raphoson iterations (40; OPTIONAL)
%    eps         Convergence criteria (0.001; OPTIONAL)
%
% On Output:
%
%    S           Salinity
%

% svn $Id$
%===========================================================================%
%  Copyright (c) 2002-2009 The ROMS/TOMS Group                              %
%    Licensed under a MIT/X style license                                   %
%    See License_ROMS.txt                           Hernan G. Arango        %
%===========================================================================%

% Initialize.

if (nargin < 7),
  iter=40;
end,

if (nargin < 6),
  eps=0.001;
end,

if (nargin < 5),
  deltaS = 0.005;
end,

S=Sguess;

epsilon=eps.*S;

% Compute salinity, iterate until convergence.

for i=1:iter,
  f=rho-roms_eos(S, T, P);
  df=-0.5.*(roms_eos(S+deltaS, T, P) - roms_eos(S-deltaS, T, P))./deltaS;
  Sold=S;
  S=S-f./df;
  if (abs(S-Sold) < epsilon),
    break
  end,
end,

return
