function PhiPointHydrostatic=hydrotatic(Phi,z,gamma)
% compute the hydrostatic point water head given the vertical density
% contribution, the head at the top and z-grid
% Phi is at zm(1), z is the cel faces and gamma at the cell centers
% TO 090308


% Copyright 2009 Theo Olsthoorn, TU-Delft and Waternet, without any warranty
% under free software foundation GNU license version 3 or later


rho=gamma; if gamma<1, rho=gamma+1; end  % use density, assume density if gamma<1
zm=0.5*(z(1:end-1)+z(2:end));
p1=(Phi-zm(1))*rho(1);
pm =cumsum([p1;abs(diff(zm)).*(0.5*(rho(1:end-1)+rho(2:end)))]);
PhiPointHydrostatic=pm./rho+zm;
