I have tried to use roms (Feb. 2007 version)
to reproduce the barotropic linear wave solution in the periodic channel.
I presumed this could be done both in the 2D set-up (#undef SOLVE3D)
and in the 3D set-up (#define SOLVE3D) if unstratified water is used in the
second case. Surprisingly, the 2D and 3D results were qualitatively
different. In the 2D case, the linear wave solution is close to analytical.
In the 3D case, the traveling wave has dissipated quickly even though
I took care to reproduce the case without dissipation
(no surf or bottom stress, AKv=0, AKt=0).
Any comments for why the 3D barotropic wave dissipates are
appreciated. I am trying to find the reason in the 2d-3d coupling
in step2d.F. However, I would also be glad to learn that my 3D set-up
(choice of cpp options and parameters) is not adequate for this test.
Some details are provided below:
==========
Test case:
==========
N-S periodic channel, no advection, no rotation,
no forcing, no dissipation, flat bottom,
and uniform conditions across the channel.
The underlying equations are:
d zeta / d t + H d v / dy = 0
d v / d t + H g d zeta / dy = 0
If zeta0<<H, roms, in an appropriate configuration,
should reproduce this analytical solution:
zeta=zeta0*cos(omega*t-k*y);
u=0;
v=g (k/omega) zeta;
c=(omega/k)=sqrt(g H);
k=2*pi/L, where L is the channel length.
==========
The cpp options used:
==========
#undef SOLVE3D (or #define SOLVE3D)
#define OUT_DOUBLE
#undef UV_ADV
#undef UV_COR
#undef SPLINES /* parabolic splines reconstruction (vert adv.) */
#undef UV_VIS2 /* Laplacian horizontal mixing */
#undef MIX_S_UV
#define UV_LDRAG /* turn ON or OFF linear bottom friction */
#undef UV_QDRAG /* turn ON or OFF quadratic bottom friction */
#undef TS_U3HADVECTION /* define if 3rd-order upstream horiz. advection */
#undef TS_C4VADVECTION /* define if 4th-order centered vertical advection */
#undef TS_DIF2
#undef MIX_S_TS
#ifdef SOLVE3D
# define SALINITY
# define NONLIN_EOS
# define DJ_GRADPS /* Splines density Jacobian (Shchepetkin, 2000) */
#endif
#if defined NL_MODEL && defined SOLVE3D
# undef MY25_MIXING
# define ANA_VMIX
#endif
#define NS_PERIODIC
#define EASTERN_WALL
#define WESTERN_WALL
#if defined SOLVE3D
# define ANA_BSFLUX
# define ANA_BTFLUX
# define ANA_BMFLUX
# define ANA_STFLUX
# define ANA_SSFLUX
#endif
In analytical.F, I set AKv=0, AKt=0 (in SUBROUTINE set_vmix)
In the input file: RDRG=0, VISC2=0
==========
2D case (#undef SOLVE3D):
==========
zeta0=0.1 m/s, H=100 m, L=100 km, dx=dy=2 km, grid 5x50 rho points, DT=8 sec
Works just fine (the ini conditions propagate with correct phase speed and
without dissipation)
==========
3D case (#define SOLVE3D):
==========
Same grid, N=10 layers, DT=240 sec, NDTFAST=30,
ini T=const=10 C, ini S=const=33 psu
The wave propagates with the proper wave speed, but its amplitude decreases
quickly (e-folding scale of 2 days).
barotropic wave in the channel (no forcing, no dissipation)
I see. So the baroclinic mode is effectively doing nothing. I'm not sure, but my guess is that is has to do with the time-filtering of the barotropic mode in the coupled system. I don't know if time-filtering is used in the standalone 2d-code, but I guess not, since it serves to remove instabilities that are caused by the splitting.There is a description of the filter shapes in the paper of Shchepetkin-McWilliams about ROMS, but I don't know what applies to the Rutgers version. Maybe you can try different filters?
You can check if that is the cause of the problem by setting the DU_avg1 and Zt_avg1 weights in a way that they represent the exact values of zeta, ubar and vbar at the last barotropic timestep (every weight zero except the last one, or smthg like that...)
I don't think it's save to use such a filter for anything else then what could be solved with the standalone 2d-version in the first place.
You can check if that is the cause of the problem by setting the DU_avg1 and Zt_avg1 weights in a way that they represent the exact values of zeta, ubar and vbar at the last barotropic timestep (every weight zero except the last one, or smthg like that...)
I don't think it's save to use such a filter for anything else then what could be solved with the standalone 2d-version in the first place.
What is your time step "dt" (for 3D mode) and mod splitting ratio "ndtfast", and maximum
barotropic Courant number "Cu_max" as it is reported by the model standard output?
Also try to estimate wtat is the ratio of
2*pi * dt * wave_phase_speed / wave_length_of_interest ?
where pi=3.141596..., wave_phase_speed is basically sqrt(g*h), and
wave_length_of_interest is period of you cosine wave.
If the above is more than one, then you should expect significant dissipation. This is
normal and the code is designed to do that. If is is smaller than one, but you still have
dissipation, then there is soomething wrong.
Then, assuming that the ratio above is of order of 1, try to play with parameters
in "set_weights". Set p=2, q=4, and play with r=0.1 .... 0.2 (but no more than 0.284)
see what happens.
>
>... I see. So the baroclinic mode is effectively doing nothing.
>
actually barotropic mode is doing pretty much everything in this case.
barotropic Courant number "Cu_max" as it is reported by the model standard output?
Also try to estimate wtat is the ratio of
2*pi * dt * wave_phase_speed / wave_length_of_interest ?
where pi=3.141596..., wave_phase_speed is basically sqrt(g*h), and
wave_length_of_interest is period of you cosine wave.
If the above is more than one, then you should expect significant dissipation. This is
normal and the code is designed to do that. If is is smaller than one, but you still have
dissipation, then there is soomething wrong.
Then, assuming that the ratio above is of order of 1, try to play with parameters
in "set_weights". Set p=2, q=4, and play with r=0.1 .... 0.2 (but no more than 0.284)
see what happens.
>
>... I see. So the baroclinic mode is effectively doing nothing.
>
actually barotropic mode is doing pretty much everything in this case.
Colleagues,
Thank you for your informative and thoughtful replies. I'll check on these and let you know.
Meanwhile, can you answer the simple question? In the 3D case, does rhs3d (specifically, prsgrd) compute the total pressure (- g d zeta/ dx + (g/rho0) int_{z}^{zeta} [d rho / dx] dz') or only its rho-dependent part? In other words, is b/t pressure grad included in the
baroclinic forcing?
Alex Kurapov
Thank you for your informative and thoughtful replies. I'll check on these and let you know.
Meanwhile, can you answer the simple question? In the 3D case, does rhs3d (specifically, prsgrd) compute the total pressure (- g d zeta/ dx + (g/rho0) int_{z}^{zeta} [d rho / dx] dz') or only its rho-dependent part? In other words, is b/t pressure grad included in the
baroclinic forcing?
Alex Kurapov
It should not dissipate, if that number is 0.47. However, you specified the whole length
of your channel as "wavelenght_of_interest". I suspect that the actual wave is shorter
than that, so it may be not 0.47, but is more. As of right now, your signal travels 7.5km
per time step, and it takes about 14 steps to cross the whole channel. That sets limit of
what is resolved and what is filtered.
The rule of thumb is that your Cu_max should be around 0.8. If less, then it means that
you are stepping too slow. There is nothing wrong with that, except that you spend more
CPU time than you actually have to. But the frequency you are interested in MUST BE
RESOLVED by "dt" (meaning "dt" for 3D mode), which means that 2*pi * dt/T < 1
(say 0.5 is OK, "T" is wave period, T=wavelength/c). Because you are actually after
barotropic signal, this boild down to that your mode splitting ratio should be comparable
to the number of grid points needed to accurately describe your sinusoidal wave.
To begin with, chose DT=120 and ndtfast=15 (both half of what you have, hence you
keep the same 8 sec barotropic step) and see what effect it makes.
deally you want to reduce ndtfast even more (try ndtfast=10), but see this post
viewtopic.php?p=645&highlight=#645
and scroll it down to the very bottom.
Then play with p=?,q=?,r=? in set_weights.F and make sure that you averaging with
proper shape. Keep p,q constant, but increase r a little bit and see how it affects
dissipation.
of your channel as "wavelenght_of_interest". I suspect that the actual wave is shorter
than that, so it may be not 0.47, but is more. As of right now, your signal travels 7.5km
per time step, and it takes about 14 steps to cross the whole channel. That sets limit of
what is resolved and what is filtered.
The rule of thumb is that your Cu_max should be around 0.8. If less, then it means that
you are stepping too slow. There is nothing wrong with that, except that you spend more
CPU time than you actually have to. But the frequency you are interested in MUST BE
RESOLVED by "dt" (meaning "dt" for 3D mode), which means that 2*pi * dt/T < 1
(say 0.5 is OK, "T" is wave period, T=wavelength/c). Because you are actually after
barotropic signal, this boild down to that your mode splitting ratio should be comparable
to the number of grid points needed to accurately describe your sinusoidal wave.
To begin with, chose DT=120 and ndtfast=15 (both half of what you have, hence you
keep the same 8 sec barotropic step) and see what effect it makes.
deally you want to reduce ndtfast even more (try ndtfast=10), but see this post
viewtopic.php?p=645&highlight=#645
and scroll it down to the very bottom.
Then play with p=?,q=?,r=? in set_weights.F and make sure that you averaging with
proper shape. Keep p,q constant, but increase r a little bit and see how it affects
dissipation.
Concerning the pressure gradient: I'm looking at prsgrd31.h (probably not good, but maybe for flat ocean bottom it's ok) and there is a term called fac2 that is multiplied by a gradient of z_w(i,j,N). It looks like the shallow-water-part of the prsgrd. z_w(i,j,N) is set in set_depth (called at the last barotropic loop)
and contains Zt_avg1, so this is filtered already. I had difficulties deciphering set_depth, but z_w(i,j,N) represents probably zeta in the baroclinic mode prsgrd calculation. So I think in the comment "Compute surface baroclinic pressure gradient." in prsgrd31.h, the word "baroclinic" wants to say that the calculation is done in the 3d-algorithm, but it also includes the swe-term.
and contains Zt_avg1, so this is filtered already. I had difficulties deciphering set_depth, but z_w(i,j,N) represents probably zeta in the baroclinic mode prsgrd calculation. So I think in the comment "Compute surface baroclinic pressure gradient." in prsgrd31.h, the word "baroclinic" wants to say that the calculation is done in the 3d-algorithm, but it also includes the swe-term.
Sasha and Stef,
Thanks for your valuable comments. I agree, it should really make sense to have b/t time step small enough to resolve the fast b/t wave in time. Here are results of experiments
w/ decreased DT. The number shows the decrease in magnitude over two days (max zeta(t=48h)/ max zeta(t=0)):
DT NDTFAST Decay (over 2 days)
240 30 0.07
120 15 0.15
120 30 0.2
60 30 0.7
It is also nice to learn that with larger time steps, when fast b/t waves are not resolved adequately, those are damped.
Cheers,
Alex Kurapov
Thanks for your valuable comments. I agree, it should really make sense to have b/t time step small enough to resolve the fast b/t wave in time. Here are results of experiments
w/ decreased DT. The number shows the decrease in magnitude over two days (max zeta(t=48h)/ max zeta(t=0)):
DT NDTFAST Decay (over 2 days)
240 30 0.07
120 15 0.15
120 30 0.2
60 30 0.7
It is also nice to learn that with larger time steps, when fast b/t waves are not resolved adequately, those are damped.
Cheers,
Alex Kurapov