function [ga, pe] = sw_gpan(S,T,P,LAT)

% SW_GPAN    Geopotential anomaly
%=========================================================================
% SW_GPAN  $Revision: 1.3 $  $Date: 1994/10/10 05:01:00 $
%          Copyright (C) CSIRO, Phil Morgan 1992.
%
% USAGE:  [gpan, pe]= sw_gpan(S,T,P,LAT)
%
% DESCRIPTION:
%   Geopotential Anomaly calculated as the integral of svan from the
%   the sea surface to the bottom.  Thus RELATIVE TO SEA SURFACE.
%
% INPUT:  (all must have same dimensions)
%   S = salinity    [psu      (PSS-78)]
%   T = temperature [degree C (ITP-68)]
%   P = Pressure    [db]
%   LAT = latitude
%       (P may have dims 1x1, mx1, 1xn or mxn for S(mxn) )
%
% OUTPUT:
%  gpan = Geopotential Anomaly  [m^3 kg^-1 Pa == m^2 s^-2 == J kg^-1]
%  pe   = anomalie d'energie potentielle en kg/s^2 * 10
%
% AUTHOR:  Phil Morgan 92-11-05  (morgan@ml.csiro.au)
%
% DISCLAIMER:
%   This software is provided "as is" without warranty of any kind.  
%   See the file sw_copy.m for conditions of use and licence.
%
% REFERENCE: S. Pond & G.Pickard  2nd Edition 1986
%            Introductory Dynamical Oceanogrpahy
%            Pergamon Press Sydney.  ISBN 0-08-028728-X
%
% Note that older literature may use units of "dynamic decimeter' for above.
%
% Adapted method from Pond and Pickard (p76) to calc gpan rel to sea 
% surface whereas P&P calculated relative to the deepest common depth.
%=========================================================================

%
% CALLER: general purpose
% CALLEE: sw_svan.m 

%----------------------
% CHECK INPUT ARGUMENTS
%----------------------
if nargin ~= 4
   error('sw_gpan.m: Must pass 4 parameters')
end %if

% CHECK S,T,P dimensions and verify consistent
[ms,ns] = size(S);
[mt,nt] = size(T);
[mp,np] = size(P);

  
% CHECK THAT S & T HAVE SAME SHAPE
if (ms~=mt) | (ns~=nt)
   error('check_stp: S & T must have same dimensions')
end %if

% CHECK OPTIONAL SHAPES FOR P
if     mp==1  & np==1      % P is a scalar.  Fill to size of S
   P = P(1)*ones(ms,ns);
elseif np==ns & mp==1      % P is row vector with same cols as S
   P = P( ones(1,ms), : ); %   Copy down each column.
elseif mp==ms & np==1      % P is column vector
   P = P( :, ones(1,ns) ); %   Copy across each row
elseif mp==ms & np==ns     % PR is a matrix size(S)
   % shape ok 
else
   error('check_stp: P has wrong dimensions')
end %if
[mp,np] = size(P);
 

  
% IF ALL ROW VECTORS ARE PASSED THEN LET US PRESERVE SHAPE ON RETURN.
Transpose = 0;
if mp == 1  % row vector
   P       =  P(:);
   T       =  T(:);
   S       =  S(:);   

   Transpose = 1;
end %if
%***check_stp

%------
% BEGIN
%------
db2Pascal  = 1e4;
[m,n]      = size(P);
svan       = sw_svan(S,T,P);
mean_svan  = 0.5*(svan(2:m,:) + svan(1:m-1,:) );

if n==1
   top = svan(1,1).*P(1,1)*db2Pascal;
else
   top = svan(1,:).*P(1,:)*db2Pascal;
end %if

%press_diff = diff(P);

delta_ga   = (mean_svan.*diff(P))*db2Pascal;
ga         = cumsum([ top; delta_ga]);

% Rajoute par Y. Gouriou en se basant sur la routine de Millard
%  pe[ref] = hdy[ref] * (ctd[ref].pres + ctd[ind].pres)*0.5F / 
%		          (float)grav( (double)ctd[ref].pres,lat ); */
%  pe[ref] =  (float) ( fabs( (double)( pres[ref] - pres[ind]) )
%                     * (anvs[ref]*pres[ref] + anvs[ind]*pres[ind])
%                     * 0.5e-5 / (float)grav( (double)pres[ref],lat ));
mean_P     = 0.5*( P(2:m,:) + P(1:m-1,:) );
delta_pe   = delta_ga .* mean_P ./  sw_g( LAT, -sw_dpth(P(2:end),LAT));
pe         = cumsum([ top/sw_g( LAT, 0); delta_pe]); 

if Transpose
   ga = ga';
   pe = pe';
end %if

return
%--------------------------------------------------------------------