Yes, the problem here is due to the fact that the pressure gradient scheme described
in the 2003 paper WAS NEVER FULLY implemented in Rutgers ROMS.
That paper introduces the notion of adiabatic differences and harmonic averaging of
adiabatic differences as a method to enforce non-oscillatory stratification of cubic
polynomial interpolants. Note that just monotonic decrease of "in situ" density with
vertical coordinate z (pointing upward) does not necessarily guarantees that the water
column is positively stratified: in fact, most of density variation in the ocean occurs
due to what I would call "bulk compressibility" -- imagine than T,S are spatially
uniform and equal to T=3.8 degrees, S=34.5 PSU. This would still yield about 90% of
the observed variation of "in situ" density. However, in can be shown (Dukowicz, 2001)
that "bulk compressibility" does not create any contribution to pressure gradient and,
in fact, is dynamically passive.
The 2003 paper argues that pressure gradient should not be computed from "in situ"
density, but rather splits EOS into two parts to facilitate computing of adiabatic
differencing, see Sec. 7 there. Once this is properly done, pressure gradient force
should vanish identically when T,S are spatially uniform. Yes, this requires a special
form of EOS.
The PGF errors associated with neglecting this issue (i.e. computing PGF from straight
from "in situ" density vs. splitting EOS) are shown on Fig. 16 in 2003 paper.
The dynamical errors due to entirely neglecting seawater compressibility (missing
thermobaric effect) are shown on Figs. 17 and 18.
Rutgers ROMS instead took the route of bringing "latest and greatest", and supposedly
more accurate versions of EOS, not realizing that it leads to the admission of much
cruder errors due to sigma-problem than the improvement of the physical accuracy of
EOS relative to the Jackett and McDougal, 1995 version.
Another related dilemma is a long-term, 40+ years, misunderstanding of the Boussinesq
approximation within the oceanic modeling community leading to the belief that the more
accurate EOS is used the better it is for your model. This is slowly being sorted out
starting with the year 2010, but overall it already took 100 years (counting from the
1903 lecture noted by Boussinesq) to realize that if the Boussinesq approximation is
desired it should be applied to the entire system of equations including EOS.
In fact, most of the effort to improve EOS starting with the original UNESCO standard
of early 198x was undertaken pursuing two directions: (i) accuracy of representation
of "in situ" density under the conditions of extreme pressure, and (ii) to ensure
internal thermodynamic consistency among the derived properties resulting from EOS
(after work of Feistel, 2003). While both are important goals, neither of them leads
to any improvement in a Boussinesq model: accuracy of "in situ" density is irrelevant,
since in some places you have is as "in situ", in other places it is replaced by rho0.
Thermodynamics is grossly simplified by the Boussinesq approach to start with, reducing
it to just conservation. So the only relevant properties needed from EOS are translations
of spatial derivatives of T,S into "in situ" adiabatic derivatives of density (derivatives
under constant pressure equal to the actual local pressure). This opens the possibility
to significantly simplify EOS and yet, have a very physically accurate Boussinesq model.