barotropic pres gradient

General scientific issues regarding ROMS

Moderators: arango, robertson

Post Reply
Message
Author
kurapov
Posts: 29
Joined: Fri Sep 05, 2003 4:49 pm
Location: COAS/OSU

barotropic pres gradient

#1 Unread post by kurapov »

I've been trying to do fancy momentum term balance analysis, doing this by post-processing of
ROMS output. I've translated the roms pressure grad term in Matlab and
stumbled upon the following detail. In the code, the barotropic pressure
gradient, e.g., at level N, is apparently twice the value that I would expect.
Since this part of the code has been scrutinized by many people, and since the output
looks fine (surface velocity is close to (g/f) d zeta/ dx), I would presume that there is
a mistake in the way I interpret the code. Can you relieve my pain?

For simplicity, let's assume the grid is Cartesian and square: on_u=const=dy=dx.
Let's assume fluid is homogeneous: rho=const=rho0=1000 (kg/m3).

To my understanding, the pressure gradient term in xi dir., computed by a ROMS subroutine,
should be equal to:

- Hz*dx*dy*g*[d(zeta)/dx] (*)

Comments here:

(1) The baroclinic part of the pressure gradient, the integral, is 0 since rho=const.
(2) In ROMS, terms are scaled by the grid cell area, dx*dy, and the layer hight, Hz,
computed at the appropriate, u-, location.
(3) zeta is used as a generic symbol for SSH; in the code its value is probably Zt_avg1.

Take, e.g., Nonlinear/prsgrd31.h:

Code: Select all

      fac1=0.5_r8*g/rho0
      fac2=1000.0_r8*g/rho0
      fac3=0.25_r8*g/rho0

      DO j=Jstr,Jend
        DO i=IstrU,Iend
          cff1=z_w(i  ,j,N(ng))-z_r(i  ,j,N(ng))+                       &
     &         z_w(i-1,j,N(ng))-z_r(i-1,j,N(ng))
          phix(i)=fac1*(rho(i,j,N(ng))-rho(i-1,j,N(ng)))*cff1
#ifdef RHO_SURF
          phix(i)=phix(i)+                                              &
     &            (fac2+fac1*(rho(i,j,N(ng))+rho(i-1,j,N(ng))))*        &
     &            (z_w(i,j,N(ng))-z_w(i-1,j,N(ng)))
#endif
          ru(i,j,N(ng),nrhs)=-0.5_r8*(Hz(i,j,N(ng))+Hz(i-1,j,N(ng)))*   &
     &                       phix(i)*on_u(i,j)

In our case (RHO_SURF is defined),

phix(i)=(fac2+fac1*2*rho)*(z_w(i,j,N)-z_w(i-1,j,N))

Also,

z_w(i,j,N)=zeta(i,j)
fac2=g
fac1=0.5*g/rho0
rho0=rho

Then,

phix(i)=2*g*(zeta(i,N)-zeta(i-1,N))

and

ru(i,j,N)=-dy*Hzu*2*g*(zeta(i,N)-zeta(i-1,N)) (**)

(where Hzu=0.5*(Hz(i-1,j,N)+Hz(i,j,N)) is the layer depth at u-location

So, apparently, ru is twice as large as the expected value (*).

I look farther to check the next term, Coriolis (in rhs3d.F), and it added to ru
without a factor of 2:

ru=ru-Hz*dx*dy*f*v

I do a similar exercise with prsgrd32.h: same problem, a factor of 2
for the barotropic pressure gradient at level N.

At the same time, I check ROMS output (a coastal upwelling case) and
the Coriolis term, computed using output of surface v, is in balance with the pressure
grad term, computed using output zeta. I have not run my case with diagnostics yet, to
see what term balance output ROMS provides.

Any comments so far about the appearance of 2 in (**) ?

Thanks,

Alex Kurapov

bzhang
Posts: 11
Joined: Wed Mar 26, 2003 9:25 pm
Location: CICS/ESSIC University of Maryland

Re: barotropic pres gradient

#2 Unread post by bzhang »

In ROMS, rho is density anomaly, which is 'true density minus a reference density rho0'. In your case, density and rho0 are both 1000 kg/m3, then rho(i,j,k) is 0.0 anywhere in ROMS. So fac1=0. fac2+fac1=1, NOT 2. Hope I solve the problem.

kurapov
Posts: 29
Joined: Fri Sep 05, 2003 4:49 pm
Location: COAS/OSU

Re: barotropic pres gradient

#3 Unread post by kurapov »

bzhang: your comment makes sense. Thanks!

Post Reply