KELVIN case typo

Bug reports, work arounds and fixes

Moderators: arango, robertson

Post Reply
Message
Author
rmueller
Posts: 26
Joined: Thu May 04, 2006 4:15 am
Location: Earth & Space Research

KELVIN case typo

#1 Unread post by rmueller »

This may be fixed in ROMS 3.X, but there is a bug in the KELVIN case for ROMS 2.2.

The Kelvin wave solution should be of the form
U = sqrt(g/H)*exp(-f*y/sqrt(g*H))*cos(omega*time).

The KELVIN case in Analytical.F file has the form of:
U = sqrt(g/H)*exp(-f*y/sqrt(g/H))*cos(omega*time).

I'm still a new-be with FORTRAN, so perhaps I'm off the mark on this? Just a heads up to others who may be using this. And again...note that I'm referring to ROMS 2.2!

Salud,
Rachael

jcwarner
Posts: 1173
Joined: Wed Dec 31, 2003 6:16 pm
Location: USGS, USA

#2 Unread post by jcwarner »

it appears to be incorrect as well in the current version.
Go ahead and submit this as a bug report.
In the new version it is in Functionals/ana_m2obc.h

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

#3 Unread post by kurapov »

Re: U = sqrt(g/H)*exp(-f*y/sqrt(g*H))*cos(omega*time).

Anybody looking into that, check dimensions. Is U a velocity (m/s)?
sqrt(g/H) has units of (1/s). Does the wave aplitude have to be sqrt(g*H)?

rmueller
Posts: 26
Joined: Thu May 04, 2006 4:15 am
Location: Earth & Space Research

Kelvin wave bug

#4 Unread post by rmueller »

Thanks for mentioning this Alexander.

The dimensions actually do work out because the equation I posted is incomplete. There is an amplitude scaling in the code that provides the 'm' dimension.

I really should have written:
U = eta_o * sqrt(g/H) * exp(-f*y/sqrt(g*H))*cos(omega*time), for f>0

In the current iteration, the scaling is 1m.

Do I need to do anything else to "submit a bug report"? Or suffice that it's posted here?

User avatar
arango
Site Admin
Posts: 1347
Joined: Wed Feb 26, 2003 4:41 pm
Location: DMCS, Rutgers University
Contact:

#5 Unread post by arango »

Yes, good catch. Kelvin waves are nondispersive gravity waves moving at the speed of c=SQRT(g*h). If you derive the dispersion relationship from the linear shallow water equations, you will get:

Code: Select all

U = c*exp(-f*y/c)*cos(omega*time)
Notice that c needs to be used inside and outside of the exponential. This renders the exponential nondimensional. So your equation above is partially correct. As Alex suggested, it is always a good idea to check the dimensions. This test was contributed by an user and I didn't checked it carefully.

The correction is very easy in subroutine ana_m2obc. You just need to change the value of cff in two places:

Code: Select all

cff=SQRT(g*GRID(ng)%h(Istr-1,j))

...

cff=1.0_r8/SQRT(g*GRID(ng)%h(Iend,j))
for the western and eastern boundary, respectively. By the way, notice that the commented code for elevation in subroutine ana_initial is correct. It is using SQRT(g*h).

rmueller
Posts: 26
Joined: Thu May 04, 2006 4:15 am
Location: Earth & Space Research

clarification on KELVIN case fix

#6 Unread post by rmueller »

Thanks for the commentary Kevin.

Just to prevent further confusion in interpreting your representation and my earlier comment, the equation can be written in two ways:

U = Uo*exp(-f*y/c)*cos(omega*time), for f>0

where
Uo = c*eta_0/H, c = sqrt(g*H)

OR

Uo = eta_0*sqrt(g/H)

I think there may be another bug here too, but I haven't tried to run the KELVIN case as is (and don't have time to right now). In my case, the "GRID(ng)%yp(Istr-1,j)" values are null. This means that instead of an off-shore (y) decay there is a constant wave forcing in the y-direction. I am not using an analytical grid, which may be the cause, but I don't see any assignment of yp for the KELVIN case that would make these value non-zero. Again...I'm new to FORTRAN here, so this could be rookie mistake; but if anyone runs this case it may be worth double checking that the yp values are actually populated with distance values.

ciao,
Rachael

rmueller
Posts: 26
Joined: Thu May 04, 2006 4:15 am
Location: Earth & Space Research

yp problem--clarification

#7 Unread post by rmueller »

In the aforementioned reference to "GRID(ng)%yp(Istr-1,j)" values being null, I forgot to clarify that my y_psi(x_psi) values in my non-analytical grid are non-zero. If I'm understanding all this correctly, my grid's y_psi(x_psi) values should populate the yp(xp) variable in mod_grid.F. It doesn't seem to be happening because my print out of yp(xp) is null.

As mentioned above, however, this may just be a problem with using an non-analytical grid, analytical boundary forcing, ROMS 2.2, and/or (key point) being a novice user. I just figure it's worth mentioning 'cause it's the kind of bug that could easily slide under the radar.

-R

johnluick

kelvin case

#8 Unread post by johnluick »

Since the topic is being discussed, perhaps it is worth mentioning that the kelvin case boundary conditions presently only work in the northern hemisphere. Down here where f changes sign the exponential blows up. (I added an "ABS" to the exponential.) Also as I recall the combination of COS and SIN functions in the time dependence in fsobc/m2obc is only suitable for the NH (ensure that geostrophy is obeyed - the expression in m2obc may require a change of sign). I have since erased my mental and computer memory so I don't remember what I did but it will be pretty obvious. Come to think of it, there is also a bunch of junk in the kelvin case relating to the non-driven open boundary which can be heavily pruned. (No need to be a control freak, just let it radiate out!)

rmueller
Posts: 26
Joined: Thu May 04, 2006 4:15 am
Location: Earth & Space Research

Southern Hemisphere Kelvin

#9 Unread post by rmueller »

Good idea to describe the more general case that I was conveniently avoiding. :-)

I'm working in the southern hemisphere too and had to convince myself of the signs. My coast is also aligned N-S rather than E-W (hence, u=0). The full expression for my case is a simplified version of:

v = eta_0 * sqrt(g/H) * exp(i*(l*y - omega*t))*exp(-abs(f)*x/sqrt(g*H))
v = sqrt(g/H) * eta(x,y,t)

The sign on the velocity actually works out in my derivation through the definition that -1*abs(f) = f, in the Southern hemisphere. Also, since c = (+/-) sqrt(g*H), we choose the sign of c to ensure a decaying wave that supports the use of abs(f) instead of f, where f<0. The point being that the sign of v (and u too, I believe) is the same as in the northern hemisphere case. Unless, of course, I made a mistake or am talking about a different sign than you are discussing in your post.

Thanks for mentioning the excessive detail on boundary handling. I am probably being a control freak too!!! I'll take a closer look.

rmueller
Posts: 26
Joined: Thu May 04, 2006 4:15 am
Location: Earth & Space Research

Who is Kevin?!?!

#10 Unread post by rmueller »

To prevent any further assault on his name, I just wanted to publicly acknowledge that the "Kevin" mentioned in an earlier post is NOT a nickname for Hernan! I hope my error doesn't inspire new users to do the same. It's just a mistake that I made to publicly embarrass myself for calling Hernan by the wrong name. Not something that I'd recommend. It defies rule #1: Don't call people by the wrong name, especially if they are helping you out.

Thanks, Hernan, for accepting the apology that I sent by email!
:oops:

Best,
Nicky

just kidding...Rachael.

Post Reply