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
several doubts for my first application

 Posts: 52
 Joined: Tue Mar 03, 2009 2:39 pm
 Location: C.N.R.  LaMMA
Re: several doubts for my first application
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 highReinolds 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
ForwardBackward spep, so 1.1 for LFAM3); 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 apriori 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 nonoverlapping rectangles. Here internal means all the RHOpoints where time stepping of all RHOtype 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 Uand Vpoints. 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 Cgrid Upoints are placed halfway between RHOpoints, so if you imagine a closed rectangular basin, the lateral boundaries pass exactly through points of velocity component normal to the boundary, while nearest RHOpoints end up being halfinterval inside. As the result, the number of internal Upoints in XIdirection is less than RHOpoints 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 Red bounds are always set to nonRed, 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=Istr1; ad ESTERN IendR=Iend1.
In the case of periodic conditions in either direction the difference between Red and nonRed 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 nonRed prototypes, index ranges IstrR:IendR,JstrR:JendR are also nonoverlapping.
Also note: the definition of WESTERN_ ESTERN_ something in the code is purely logical as it would be that model XIdirection (the first FORTRAN dimension, index "i") goes from west to east; ETAdirection (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.
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 highReinolds 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
ForwardBackward spep, so 1.1 for LFAM3); 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 apriori 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 nonoverlapping rectangles. Here internal means all the RHOpoints where time stepping of all RHOtype 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 Uand Vpoints. 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 Cgrid Upoints are placed halfway between RHOpoints, so if you imagine a closed rectangular basin, the lateral boundaries pass exactly through points of velocity component normal to the boundary, while nearest RHOpoints end up being halfinterval inside. As the result, the number of internal Upoints in XIdirection is less than RHOpoints 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 Red bounds are always set to nonRed, 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=Istr1; ad ESTERN IendR=Iend1.
In the case of periodic conditions in either direction the difference between Red and nonRed 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 nonRed prototypes, index ranges IstrR:IendR,JstrR:JendR are also nonoverlapping.
Also note: the definition of WESTERN_ ESTERN_ something in the code is purely logical as it would be that model XIdirection (the first FORTRAN dimension, index "i") goes from west to east; ETAdirection (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.

 Posts: 52
 Joined: Tue Mar 03, 2009 2:39 pm
 Location: C.N.R.  LaMMA
Re: several doubts for my first application
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
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

 Posts: 52
 Joined: Tue Mar 03, 2009 2:39 pm
 Location: C.N.R.  LaMMA
Re: several doubts for my first application
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
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
Re: several doubts for my first application
What scale are you modeling this at? It could be that ROMS is completely inappropriate for what you are trying to do.

 Posts: 52
 Joined: Tue Mar 03, 2009 2:39 pm
 Location: C.N.R.  LaMMA
Re: several doubts for my first application
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
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