several doubts for my first application

[mariafattorini]

Hello all,
I am setting my first simple application and I have need to know some things:

1. is Akv the classic eddy viscosity or is it VISC?

2. is there a condition to choice the size of grid cells and the size of time step to maintain numerical stability, eg. a condition on the Courant number?

3. What is the difference between Istr and IstrR ( and similarly between Iend and IendR) ?

4. If I change some value in the ana_*.h file (eg. if I change the value of Uwind in the ana_winds.h),
have I to compiler again the model?

Thank you very much for any help,
Maria

[shchepet]

1. is Akv the classic eddy viscosity or is it VISC?

No. VISC is horizontal viscosity, sometimes known as "eddy viscosity", basically parameterization of horizontal momentum exchange by interaction of eddies smaller than horizontal grid size. In practice it is chosen most of the time to ensure smoothness and integrity of the solution on the grid scale, but at the same time be as small as possible to preserve physically important processes -- most of the time the goal is to simulate high-Reinolds flows. For this reason it is also common practice to reduce VISC every time when grid is refined, or to set it proportional to local dx, on highly curvilinear grids, where dx changes significantly depending on location.

Akv is vertical viscosity. It is always physically defined (computed by a vertical mixing parameterization scheme -- GLS, MY, or KPP). Its primary role is to distribute wind stress throughout vertical column (without Akv only wind stress would apply only on the uppermost grid box, which is physically incorrect; similarly near the bottom -- the bottom drag causes appearance of bottom boundary layer which tends to slow down flow within some finite thickness near the bottom).

In contrast with VISC, there is no need to have Akv for numerical reasons.

2. is there a condition to choice the size of grid cells and the size of time step to maintain numerical stability, eg. a condition on the Courant number?

Yes, there are stability limits imposed by phase speed of surface and internal waves, Coriolis and advection terms.

For surface gravity waves dt * sqrt[gh(1/dx^2+1/dy^2)] should be less that 1, or strictly speaking, less than certain number of order of 1, which depends on specific time stepping algorithm (say 0.87 for
Forward-Backward spep, so 1.1 for LF-AM3); here g=9.81 is acceleration of gravity; h is local depth; and
dt is barotropic time step, that is baroclinic dt / ndtfast (baroclinic dt and ndtfast are set in roms.in file); Advection imposes limit dt*u/dx is less than 0.5 or so; and internal gravity wave speed also imposes limit like dt*c/dx less than 1 or so -- again specific numbers depend on specific time stepping algorithms of varian of ROMS kernel used, but in any case, all of them are designed with anticipation that typical internal wave speed is greater than advection speed, physically more restrictive process receives different numerical treatment. It is usually hard to have an a-priori estimate phase speed of internal gravity wave, but common values are 2...3 m/sec.

Vertical viscosity/diffusivity is treated implicitly, so they do not impose any extra restriction.

Horizontal viscosity/diffusivity are always explicit, so in principle they may restrict time step, but their coefficient are always chosen small enough for advective CFL be always more restrictive than viscous (recall that if these two are equally restrictive, the advective scheme essentially became equivalent to first order upstream scheme, which is highly viscous).

3. What is the difference between Istr and IstrR ( and similarly between Iend and IendR) ?


Istr,Iend,Jstr,Jend are known as primary tile bounds. They are designed to split internal portion of computational domain into a set of non-overlapping rectangles. Here internal means all the RHO-points where time stepping of all RHO-type variables is performed. This excludes physical boundary points and, ghost points associated with periodic bounary conditions and MPIexchange zones (if MPI is used).

IstrU and JstrU are lower bounds for U-and V-points. They are characterized by having one point less near WESTERN and SOUTHERN physical boundaries. Thus, IstrU=Istr+1 if the tile is adjacent to the WESTERN boundary; and IstrU=Istr elsewhere. This is because on Arakawa C-grid U-points are placed half-way between RHO-points, so if you imagine a closed rectangular basin, the lateral boundaries pass exactly through points of velocity component normal to the boundary, while nearest RHO-points end up being half-interval inside. As the result, the number of internal U-points in XI-direction is less than RHO-points by 1.

IstrR,IendR,JstrR,JendR are "auxiliary bounds" which are the same as Istr,Iend,Jstr,Jend, but they also include physical boundaries; periodic and MPI ghost pointrs are still excluded. Thus, initially R-ed bounds are always set to non-R-ed, and then it is checked whether the subdomain has a segment physical boundary; if so, the ranges are extended, i.e., at the WESTERN boundary IstrR=Istr-1; ad ESTERN IendR=Iend-1.

In the case of periodic conditions in either direction the difference between R-ed and non-R-ed bounds disapears, and so does the difference between IstrU and Istr; and JstrV and Jstr.

Very rarely and typically limited to few analytical routines IstrR,IendR,JstrR,JendR are designed to include periodic and MPI ghost points in addition to physical boundaries. They are called "extended bounds" in this case.

Note: similarly to their non-R-ed prototypes, index ranges IstrR:IendR,JstrR:JendR are also non-overlapping.


Also note: the definition of WESTERN_ ESTERN_ something in the code is purely logical as it would be that model XI-direction (the first FORTRAN dimension, index "i") goes from west to east; ETA-direction (second dimension, index "j") goes from south to north. If the grid is curvilinear and rotated, the WESTERN_EDGE may end up pointing to the geographical south, but it is still called WESTERN_EDGE in the code.

A useful illustration related to this matter may be found in

http://people.arsc.edu/~kate/ROMS/whole_grid.png

for the whole grid, and

http://people.arsc.edu/~kate/ROMS/Istr.png

for tiling.
4. If I change some value in the ana_*.h file (eg. if I change the value of Uwind in the ana_winds.h),
have I to compiler again the model?

Yes, you have to recompile every time once you change a .F or .h file.

[mariafattorini]

Hello Shchepet and thank you very much for your clear explanations.

I have another question:

I tried to simulate a tank with flat bathymetry and a hump in the middle of the base of the tank without forcing (not waves, nor tide, no wind, no inflow, ..) and with zero initial conditions, but the result is not a zero velocity field.
The velocity are very some as 10^-13 or so.

Why does the model create some small velocity? is it a numerical question?

And is the model able to shock capture?

Thank you very much again,
Maria

[mariafattorini]

Hello all,

Studying the roms, I have found the other following doubts:

1. What is the difference between R0 (density for state equation) and RHO0 (mean density)?
2. Usually what are common values of mixing coefficients (both vertical, Akt and Akv, and horizontal, VISC2 and TNU2)?
3. There are many scheme about vertical mixing: my study concerns a submerged downward cold sea water jet that comes from a tube of diameter=1.9m and located at 12m under the free surface.I have no idea on what is the better vertical mixing in my case. Can someone give me any suggestion?


Many thanks for any help,
Maria

[kate]

What scale are you modeling this at? It could be that ROMS is completely inappropriate for what you are trying to do.

[mariafattorini]

Hello Kate,

I did roms model with a dimension grid about 10m and discharge features given by a "near field" model (done with a other software), as suggested by jivica in the other post.

I have an other question:
in ana_psource, how is it distinguished a sink from a source?

Thank you very much,
Maria