function theta = potemp80(s,t,p,pref)

%function theta = potemp80(s,t,p, [pref])
%
%-----------------------------------------------------------------------------
%       Local potential temperature at pref
%
%     Input:
%            s  -  Salinity in PSU-78
%            t  -  Temperature in IPTS-68
%            p  -  Pressure in decibars
%            pref - Reference pressure (Optional default = 0)
%     Output:
%            theta -Potential Temperature
%
%
%     References:   using atg80 procedure for adiabadic lapse rate
%                   Fofonoff,N.,Deep-Sea Res.,24,489-491
%
%
%     Checkvalues:
%
%---------------------------------------------------------------------
%
%
%  S. Chiswell 1991
% rev. 28-07-95 W. Erasmi

if nargin<4, pref=0; end;
if nargin < 3,
   error('Not enough input arguments.');
elseif any(size(s)~=size(t)) | any(size(t)~=size(p)),
   error('S, T and P must have the same size.');
end

[ns,ms]=size(s);

s = s(:);
t = t(:);
p = p(:);
h = pref - p;
xk = h .* atg80(s,t,p);

t = t + 0.5 .* xk;
q = xk;
p = p + 0.5 .* h;

xk = h .* atg80(s,t,p);
t = t + 0.29289322 .* (xk-q);
q = 0.58578644 .* xk + 0.121320344 .* q;

xk = h .* atg80(s,t,p);
t = t + 1.707106781 .* (xk-q);
q = 3.414213562 .* xk - 4.121320344 .* q;
p = p + 0.5 .* h;

xk = h .* atg80(s,t,p);
theta = t + (xk - 2.0 .* q)/6.0;

theta = reshape(theta,ns,ms);
