function [p_dens]=potential_density(p,t,s,p_ref)
%
% function [dens]=potential_density(p,t,s,p_ref)
% Note that p, t and s are in situ values.
% p	m*1 column vector of standard pressures in decibars
% t	m*n matrix of temperatures (NaN used as flags) in celsius
% s	m*n matrix of salinities (NaN used as flags) in parts per
%                  thousand
% p_ref   scalar reference pressure in decibars
% p_dens  m*n matrix of sigma, the potential density anomoly, i.e.,
% potential density in kg/m^3 - 1000.

%     Copyright J. V. Mansbridge, CSIRO, march 16 1992


[m,n]=size(t);
p_matrix = p*ones(1, n);
thet = theta(p_matrix,t,s,p_ref);
p_ref_vec = p_ref*ones(m, 1);
p_dens = density( p_ref_vec, thet, s);
