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