function r=eosben07_rho(p,th,s)
%EOSBEN07_RHO Density.
%   EOSBEN07_RHO(P, TH, S) computes the seawater density [kg m-3] at
%   pressure, P [dbar], potential temperature (referenced to 0 dbar), TH
%   [degC], and salinity, S [g kg-1].

eosben07_const

r=(a11+th.*(a12+a14*th+a15*s)+s.*(a13+a16*s)+p.*(b11+b12*th+b13*s)) ...
./(a21+th.*(a22+a24*th+a25*s)+s.*(a23+a26*s)+p.*(b21+b22*th+b23*s));
