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.
