The difference between the 2003 paper and what you see in the code is explained by the fact

that equations in the paper are written assuming that density "rho" is the whole density (meaning

that it is about 1030. kg/m^3), while the code is written assuming that density "rho" is density

perturbation (with at least 1000 subtracted). As the result, in the code there is a need to add

the 1000 (or so) back wherever is needed -- basically just the barotropic contribution of the

total pressure.

Splitting density into bulk part in the code is motivated by the need to avoid roundoff errors:

keeping density as "whole" would degrade the accuracy of computation of pressure gradient by at

least two (actually close to 3) decimal places, provided that everything else is equivalent.

[Note that in the actual pressure gradient routine rho(i,j,k) (which is anomaly) is multiplied by

z_r(i,j,k), which is fairly large in comparison with free-surface perturbation, but 1000*g/rho0

is multipled only by free surface "zeta". As the result, cancellation of large terms is done by

hand, and does not compromise roundoff errors.]

There are also different ways of dealing with it, regarding what value is subtracted from the real

ocean density. Rutgers ROMS traditionally subtracts 1000. As the result, fortran variable "rho"

can sometimes (at least at surface) be interpreted as "sigma_t" in the oceanographic tradition,

and you have appearance of the expression like GRho0= 1000.*g/rho0. Recal that in a Boussinesq

model density appears only as "boyancy", that is in combitation -g*rho/rho0, so the combination

GRho0 + Grho*rho(i,j,k)

is merely translation from anomaly to full density with multiplication by g/rho0 at the same time.

UCLA and Agrif ROMS subtract rho0 instead of 1000. As a consequence of this choice, the above

expression becomes just

g + Grho*rho(i,j,k)

and "rho" stored in fotran array is no longer interpretable as "sigma_t", but in fact it is

much smaller, depending on the choice of rho0, which can be selected by the user.

Furthermore, in-sity density in the ocean is not just

rho0 + rho_perturbation

but one can actually define it as

rho0 + rho_(z) + rho_perturbation

where rho_(z) is a universal profile, basically close to a linear function, which can

be universally defined. This is because speed of sound in the ocean (a natural measure of

compressibility) does not change very much -- 1480...1540 m/c -- relatively to its mean

value, and because variations of T and S tend to be smaller at depth where compressibility

effects matters most. As the result, rho_perturbation is now defined not just relatively

to a constant reference value rho0, but relatively to the profile rho0+rho_(z), which means

that rho_perturbation is now much much smaller. At the same time one can prove that rho_(z)

is dynamically passive, even thought rho_(z) is a nonlinear function of z (pressure). This

is done in Appendix B in the 2003 paper in a somewhat crude way, and in a more refined way

in

http://www.atmos.ucla.edu/~alex/ROMS/EOSArticleRev1.pdf
the latter is motivated by not just reducing the roundoff and pressure gradient errors,

but also by the desire to reduce some errors associated with the Boussinesq approximation

as well.