## barotropic wave in the channel (no forcing, no dissipation)

### barotropic wave in the channel (no forcing, no dissipation)

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).

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).

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.

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.

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

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.

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.

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