function [DENSITY,DEPTH,SA,CT]=density_m(TEMP,SAL,PRES,LON,LAT,refpres);

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PURPOSE: Computes DENSITY from in situ temperature and salinity,
% using Thermodynamic functions from TEOS-10 (http://www.TEOS-10.org)
% 
% INPUT: matrices with profile data in columns
%     TEMP = in-situ temperature (ITS-90)                         [ deg C ]
%     SAL = Practical Salinity (PSS-78)
%     PRES = sea pressure [dbar], (i.e. absolute pressure - 10.1325 dbar)
%     LON = longitude in decimal degrees, [ 0 ... +360 ] or [ -180 ... +180 ]
%     LAT = latitude in decimal degrees north, [ -90 ... +90 ] 
%     refpres = reference pressure to compute potential density
%
% OUTPUT: 
%     DENSITY = in-situ density in [kg/m^3]
%     DEPTH = depth [m]
%     SA = ABSOLUTE SALINITY, [g/kg ]
%     CT = CONSERVATIVE TEMPERATURE (ITS-90)                      [ deg C ]
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PREDEFINE MATRICES for OUTPUT VARS
[nz,np]=size(PRES);
Z=nan(nz,np);
SA=nan(nz,np);
CT=nan(nz,np);
in_ocean=nan(nz,np);
DENSITY=nan(nz,np);

% Compute vars fro each profile/column
for i=1:length(LAT)
    % Compute Depth from pressure
    Z(:,i)=gsw_z_from_p(PRES(:,i),LAT(i)); % Thermodynamic functions from TEOS-10
    
    % Compute Absolute Salinity (from practical salinity)
    [SA(:,i), in_ocean(:,i)] = gsw_SA_from_SP(SAL(:,i),PRES(:,i),LON(i),LAT(i));
    
    % Compute Conservative Temperature from in-situ temperature and Absolute Salinity
    CT(:,i)=gsw_CT_from_t(SA(:,i),TEMP(:,i),PRES(:,i));
    
    % Compute Density
    % refpres=0; % in dbar, reference pressure to compute potential density
    DENSITY(:,i)=gsw_rho(SA(:,i),CT(:,i),refpres);
end

DEPTH=-Z; clear Z
