function [erho,rho,rhoave_3d,pressure] = erho_integral ( ncfile, timestep )
% ERHO_INTEGRAL:  compute potential energy
%
% USAGE:  erho = erho_integral ( ncfile, timestep );
%
% E_rho = integral_-h^zeta (rho-rhoav) g dz
%
% where rhoav is the vertical average rho:
%
% rhoav = 1/h integral_-h^zeta rho dz
% 
% g = 9.81
% 
% rho is in situ density which can be computed from temp and salt 

roms_grid = roms_get_grid ( ncfile, ncfile );

h = nc_varget ( ncfile, 'h' );

start = [timestep  0  0  0];
count = [1        -1 -1 -1];
salt = nc_varget ( ncfile, 'salt', start, count );
temp = nc_varget ( ncfile, 'temp', start, count );

%
% Guess at the pressure.  Is this ok?
pressure = -1 * roms_grid.z_r;
rho = sw_dens ( salt, temp, pressure );

[rho_int, thickness, dz] = roms_zint212 ( rho, roms_grid );
rhoave = rho_int ./ h;

%
% Here, rho is a 3d field, but rhoave is 2D, one entry for
% each separate grid cell.  So replicate up and down the
% water column?
wc_size = [roms_grid.N size(rhoave)];
rhoave_3d = zeros(wc_size);
for j = 1:roms_grid.N
	rhoave_3d(j,:,:) = rhoave;
end

quantity = (rho - rhoave_3d) .* (-1*roms_grid.z_r) * 9.81;

erho = roms_zint212 ( quantity, roms_grid );
return
