What means GRho0*z_w(N) in prsgrd32.h?

General scientific issues regarding ROMS

Moderators: arango, robertson

Post Reply

What means GRho0*z_w(N) in prsgrd32.h?

#1 Post by Schimmel » Mon Apr 20, 2009 9:33 am

I try to understand the calculation of baroclinic pressure gradient in prsgrd32.h.
I'm wondering about the following code:

Code: Select all

!  Preliminary step (same for XI- and ETA-components:
      GRho0=1000.0_r8*GRho     <-- What dose the 1000 mean? Why not GRho0=g ?
      DO j=JstrV-1,Jend
        DO i=IstrU-1,Iend
          cff2=0.5_r8*(rho(i,j,N(ng))-rho(i,j,N(ng)-1))*                &
     &         (z_w(i,j,N(ng))-z_r(i,j,N(ng)))*cff1
          P(i,j,N(ng))=GRho0*z_w(i,j,N(ng))+                            &  <-- Why is this added here?
     &                 GRho*(rho(i,j,N(ng))+cff2)*                      &
     &                 (z_w(i,j,N(ng))-z_r(i,j,N(ng)))
        END DO
I found the equation in SHCHEPETKIN AND MCWILLIAMS 2003 (eq. 5.40), but without the addition of GRho0*z_w(i,j,N(ng)).
Why is this term added in the ROMS-Code?
And what dose the term mean?

I think GRho0*z_w(i,j,N(ng)) is the same like g*zeta*1000/rho0, so it represents the pressure of the thin surface elevation water column divided by rho0.
If this is right, why is the density of this water column is estimated with 1000, and not e.g. with rho0?

Why is the term added at all? The second summand includes the water column from the topmost rho-point z_r(N) until zeta=z_w(N), so I see no reason to add the first term.

Thank you.

User avatar
Posts: 185
Joined: Fri Nov 14, 2003 4:57 pm

Re: What means GRho0*z_w(N) in prsgrd32.h?

#2 Post by shchepet » Tue Apr 21, 2009 1:55 am

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


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.

Post Reply