function rho = roms_eos(T, S, Z)

%
% ROMS_EOS:  Computes ROMS "in situ" density anomaly
%
% rho = roms_eos(T, S, Z)
%
% This function computes "in situ" density using ROMS nonlinear equation
% of state for seawater. It assumes no pressure variation along
% geopotential surfaces, that is, depth (meters; negative) and pressure
% (dbar; assumed negative here) are interchangeable. 
%
% On Input:
%
%    T           Potential temperature (3D array, Celsius)
%    S           Salinity (3D array)
%    Z           Depth (3D array, meter, negative)
%
% On Output:
%
%    rho         "In situ" density anomaly (3D array, km/m3)
%
% Check Values: (T=3 C, S=35.5, Z=-5000 m).
%
%     rho   = 50.3639165364       (kg/m3)
%     den1  = 1028.2845117925     (kg/m3)
%     bulk  = 23786.056026320     (Pa)
%
% Note: Salinity does not have physical units. Check the following forum
%       post for details:
%
%       https://www.myroms.org/forum/viewtopic.php?f=14&t=294
%
%  Reference:
%
%  Jackett, D. R. and T. J. McDougall, 1995, Minimal Adjustment of
%    Hydrostatic Profiles to Achieve Static Stability, J. of Atmos.
%    and Oceanic Techn., vol. 12, pp. 381-389.
%

% svn $Id: roms_eos.m 895 2018-02-11 23:15:37Z arango $
%=========================================================================%
%  Copyright (c) 2002-2018 The ROMS/TOMS Group                            %
%    Licensed under a MIT/X style license                                 %
%    See License_ROMS.txt                           Hernan G. Arango      %
%=========================================================================%

%--------------------------------------------------------------------------
% Set equation of state expansion coefficients.
%--------------------------------------------------------------------------

A00=+19092.56D0;  A01=+209.8925D0;   A02=-3.041638D0;   A03=-1.852732D-3;
A04=-1.361629D-5; B00=+104.4077D0;   B01=-6.500517D0;   B02=+0.1553190D0;
B03=+2.326469D-4; D00=-5.587545D0;   D01=+0.7390729D0;  D02=-1.909078D-2;
E00=+4.721788D-1; E01=+1.028859D-2;  E02=-2.512549D-4;  E03=-5.939910D-7;
F00=-1.571896D-2; F01=-2.598241D-4;  F02=+7.267926D-6;  G00=+2.042967D-3;
G01=+1.045941D-5; G02=-5.782165D-10; G03=+1.296821D-7;  H00=-2.595994D-7;
H01=-1.248266D-9; H02=-3.508914D-9;  Q00=+999.842594D0; Q01=+6.793952D-2;
Q02=-9.095290D-3; Q03=+1.001685D-4;  Q04=-1.120083D-6;  Q05=+6.536332D-9;
U00=+0.824493D0;  U01=-4.08990D-3;   U02=+7.64380D-5;   U03=-8.24670D-7;
U04=+5.38750D-9;  V00=-5.72466D-3;   V01=+1.02270D-4;   V02=-1.65460D-6;
W00=+4.8314D-4; 

grav = 9.81;

%--------------------------------------------------------------------------
% Compute density (kg/m3) at standard one atmosphere pressure.
%--------------------------------------------------------------------------

ind = find(T < -2.0);
if (~isempty(ind)),
  T(ind) = -2.0;                 % lower temperature valid minumum value
end

ind = find(S < 0.0);
if (~isempty(ind)),
  S(ind) = 0.0;                  % lower salinity valid minimum value
end

clear ind

sqrtS = sqrt(S);

den1 = Q00 + Q01.*T + Q02.*T.^2 + Q03.*T.^3 + Q04.*T.^4 + Q05.*T.^5 +   ...
       U00.*S + U01.*S.*T + U02.*S.*T.^2 + U03.*S.*T.^3 + U04.*S.*T.^4 +...
       V00.*S.*sqrtS + V01.*S.*sqrtS.*T + V02.*S.*sqrtS.*T.^2 +         ...
       W00.*S.^2;

%--------------------------------------------------------------------------
% Compute secant bulk modulus (bulk = K0 - K1*z + K2*z*z).
%--------------------------------------------------------------------------

K0 = A00 + A01.*T + A02.*T.^2 + A03.*T.^3 + A04.*T.^4 +                 ...
     B00.*S + B01.*S.*T + B02.*S.*T.^2 + B03.*S.*T.^3 +                 ...
     D00.*S.*sqrtS + D01.*S.*sqrtS.*T + D02.*S.*sqrtS.*T.^2;

K1 = E00 + E01.*T + E02.*T.^2 + E03.*T.^3 +                             ...
     F00.*S + F01.*S.*T + F02.*S.*T.^2 +                                ...
     G00.*S.*sqrtS;

K2 = G01 + G02.*T + G03.*T.^2 +                                         ...
     H00.*S + H01.*S.*T + H02.*S.*T.^2;

bulk = K0 - K1.*Z + K2.*Z.^2;

%--------------------------------------------------------------------------
% Compute "in situ" density anomaly (kg/m3).
%--------------------------------------------------------------------------

rho = (den1.*bulk) ./ (bulk + 0.1.*Z);

rho = rho - 1000;                                       % density anomaly

return

