KELVIN case typo
KELVIN case typo
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
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
Kelvin wave bug
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?
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?
- arango
- Site Admin
- Posts: 1361
- Joined: Wed Feb 26, 2003 4:41 pm
- Location: DMCS, Rutgers University
- Contact:
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:
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:
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).
Code: Select all
U = c*exp(-f*y/c)*cos(omega*time)
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))
clarification on KELVIN case fix
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
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
yp problem--clarification
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
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
kelvin case
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!)
Southern Hemisphere Kelvin
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.
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.
Who is Kevin?!?!
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!
Best,
Nicky
just kidding...Rachael.
Thanks, Hernan, for accepting the apology that I sent by email!
Best,
Nicky
just kidding...Rachael.