SELFE - FAQ

The FAQ below applies to both serial and parallel versions.

1. What are the requirements in generating a horizontal grid?
2. How do I set up a vertical grid?
3. My results show time step dependency.
4. My results show platform/compiler/CPU dependency.
5. Run crashed with a fort.11 error message "QUICKSEARCH: no intersecting edge....".
6. How to do a tidal simulation with SELFE?
7. How do I compile GOTM libraries?
8. How do I choose the flags in lqk.gr3?
9. Why are my river and plume too fresh?
10. I'd like to add some passive tracer transport to SELFE. How can I interface my own code to it?
11. How to impose river discharge if the depth is negative there?
12. How do I set up sflux/ using NARR (North America only) files in the Utility Library for atmospheric model (nws=2)?
13. How do I prepare my own netcdf files for atmospheric model (nws=2)?
14. Why does my non-hydrostatic run blow up?
15. I have large velocity at open boundary?
16. My run crashed. How can I find out why?

1. How do I set up a horizontal grid?

The two important dimensionless numbers to keep in mind when generating a horizontal grid for SELFE are:

1. Dimensionless wavelength: make sure that each wavelength is resolved by at least 20 grid points; this is usually an easy task since shallow-water waves have long wavelength;
2. CFL number=dt*sqrt(g*h)/dx>0.2; otherwise the default linear ELM inside becomes diffusive. If possible, CFL>0.5 is even better. This is satisfied with either a large dt or a small dx. For a given dt, it gives a maximum dx. If for some reason dt has to be small (e.g. tsunami simulations), you may be able to bypass this condition by turning off advection in deeper depths if the advection is negligible there.

go back

2. How do I set up a vertical grid?

You should set it up according to your application. The first question to ask is whether a "pure S" model is sufficient for your problem. We found it's adequate for most problems (99.9%) except when strong stratification encounters steep bottom slope (e.g., strongly stratified Columbia River plume flows over a shelf break with slope in excess of 1/20). In the S part of vgrid.in, adjust the constants h_c, theta_b, theta_f to resolve appropriate boundary layers. Note that h_c also determines the boundary between S-zone and sigma-zone (where traditional sigma-coordinates are used because S-transformation is invalid when h < h_c). If (in the unlikely case that) a SZ model is needed, h_s may be chosen as the max. depth of the shelf, and make sure there is a smooth transition between S and Z coordinates. You can get a preview of some sample z-coordinates at various depths by running the code with the pre-processor flag (ipre in param.in) turned on.

go back

3. My results show time ste dependency.

ELCIRC and SELFE use Eulerian-Lagrangina method for advection, and therefore have an operational time step range; that is, for a given grid resolution, reducing time step dt beyond certain value dt1 (e.g. when CFL<0.2) will not reduce the numerical error (see Fig.5 in the SELFE paper). On the other hand, when dt is increased beyond certain value (dt2), the truncation errors will dominate. Within the operational range (dt1<=dt<=dt2) the results should be relatively insensitive to the time step.

For convergence study, you need to keep CFL number fixed, i.e. reducing dx and dt simultaneously.

go back

4. My results show platform/compiler/CPU dependency.

Due to some intricate differences between different compilers/platforms some minor differences in results are expected. Bit-by-bit match of results using different numbers of CPUs is impossible also. But the differences should be small and should stablize over time iteration.

Related to this, some tools (e.g. MD5 checksum on binary outputs) won't give identical results. You shoud use the post-processing tools included in the package (e.g. read_output* etc) to compare results.

go back

5. Run crashed with error message "QUICKSEARCH: no intersecting edge...".

This is related to the Eulerian-Lagrangian method (ELM) used in SELFE. ELM requires backtracking from sides/nodes and in the case of non-hydrostatic model, elements. In a nutshell, fluid particle paths are being followed in 3D and backwards in time. During the process we need to determine the containing element for a given point; "point-in-a-polygon" is performed for this task, which suffers from the well known underflow problem. To avoid this, we nudge the starting point of the backtracking (sides/nodes) from its original location towards the element center; the nudging factor is specified in param.in as "btrack_noise" (or in newer versions, "btrack_nudge"). Proper choice of this parameter usually avoids the underflow problem.

The fatal error is therefore usually related to problems in the velocity field, e.g. very large or noisy velocity locally. So the first thing you'd check is to look at the animation of the velocity field right before the crash to spot this type of problems, and if so, fix them by adjusting some parameters. For example, use of nadv=0 in some region may lead to instability etc. Reducing 'rmaxvel' (maximum u or v allowed) may also help.

If the velocity field looks smooth, the problem is then related to the underflow problem. In the newer versions (v3.1dc and up) you should not need to adjust "btrack_nudge" etc anymore. Contact the development team if you still cannot figure out the problem.

go back

6. How to do a tidal simulation with SELFE?

The easiest way is to use 2D SELFE.

SELFE 3D formulation assumes that the bottom boundary layer (BBL) is well resolved (i.e. the velocity at the top of the bottom cell must be inside the BBL), and therefore it'll be a problem if you use only a few vertical levels in a barotropic run. In other words, the process of calibrating tidal amplitude/phases before a baroclinic run is different in SELFE from most other models. You have two options:

1. Use as many levels as you'd use in a baroclinic run and make sure the BBL is resolved by adjusting h_c, theta_b and theta_f in vgrid.in; then use a non-zero drag coefficient or roughness;
2. Use a few levels with Cd=0 in deep depths (e.g., h > 50m) in recognition that the BBL is only a small portion of the total depth. Use a non-zero drag in shallow depths.

go back

7. How do I compile GOTM libraries?

Due to conflict of licenses, we cannot distribute GOTM library with SELFE (nor should users re-distribute it). Therefore you first need to download it from gotm.net.

Here is an excerpt for gotm.net (you need netcdf library installed on your system):

• setenv NETCDFHOME /usr/local/netcdf/
• setenv NETCDFINC $NETCDFHOME/include
• setenv NETCDFLIBNAME $NETCDFHOME/lib/libnetcdf.a
• setenv GOTMDIR /home/users/yinglong/GOTM/gotm-3.2.5
• setenv FORTRAN_COMPILER IFORT
• cd src and make. The compiled libraries (including libturbulence_prod.a) are in /lib/IFORT/. Don't forget to update the path in Makefile of SELFE to point to the location of these libraries.

Note: the library currently has problems in shallow waters and may crash with NaN sometimes. If this happens try itur=3 option.

go back

8. How do I choose the flags in lqk.gr3?

The flag of "2" indicates quadratic interpolation for S,T, and this usually does not work well in the open sea region. Therefore you'd use "2" for estuaries only, and use "1" for open sea. Also consider using upwind or TVD option for transport as they are mass conservative.

go back

9. Why are my river and plume are too fresh.

This is most likely to be the result of applying freshwater discharge too close to the mouth of the estuary, i.e., the length of the river is not long enough.

go back

10. I'd like to add some passive tracer transport to SELFE. How can I interface my own code to it?

There are 3 interface points in the serial version; you can find them by searching for "user-defined tracer part". (1) In module global, define the total # of tracers you have in ntracers; (2) Define your own initial condition for each tracer at prism centers; (3) Define 3 arrays: source/sink (bdy_frc), surface (flx_sf) and bottom (flx_bt) b.c. Also don't forget to add the output flags, transport method (upwind or TVD) and horizontal b.c. etc in param.in. The actual call looks like: do_transport_tvd(1,up_tvd,tvd_mid2,flimiter2,ntracers). Consult the source code inline comments for more details.

go back

11. How to impose river discharge if the depth is negative there?

In many estuaries the depths (with reference to MSL) are often negative. In SELFE, the total depth (i.e., depth + elevation) must be positive at all times on open boundary nodes. A simple workaround is to artificially set a positive depth (say, 5m) along a small segment near the river boundaries (or add a small segment of "channel" with 5m beyond the boundary; usually ~100m will do). The depth can change sharply from 5m to negative values. Because of the imposed discharge, the inundation algorithm inside SELFE will quickly push water downstream and set up a surface slope automatically.

An alternative is to use an initial condition for elevation (elev.ic).

go back

12. How do I set up sflux/ using NARR (North America only) files in the Utility Library for atmospheric model (nws=2)?

1. Download NARR files for your application period. Each NARR file covers ~ 1 day (e.g., narr_air.1981_01_28.nc is for Jan. 28, 1981).
2. In your run directory, mkdir sflux and inside it, create symbolic links to the NARR files. e.g., if you run starts from June 10, 2004 and ends June 20, 2004, then
   
   sflux_air_1.001.nc --> narr_air.2004_06_10.nc
   
   sflux_air_1.002.nc --> narr_air.2004_06_11.nc
   
   ...
   
   sflux_air_1.011.nc --> narr_air.2004_06_20.nc
   
   sflux_air_1.012.nc --> narr_air.2004_06_21.nc (extra day to account for time zone difference).
   
   Similarly for sflux_rad_*.nc and sflux_prc_*.nc. The number "1" after "air_" denotes first data set used; you can use up to 2 sets in SELFE (which combines them with some given weights) but we recommend using only 1 set.
3. In sflux, copy the file sflux_inputs.txt and edit it to reflect the start time. The field "utc_start" is hours behind UTC.

That's it!

go back

13. How do I prepare my own netcdf files for atmospheric model (nws=2)?

You should start by studying the existing netcdf inputs available in the SELFE Utility Library . The format can be seen in the FORTRAN code called gen_atmos.f90 there, which can also be used as a template to generate your own inputs. Some very useful matlab scripts (readnc*.m) can be found in the source code bundle (under utility/), and you can clearly see all the steps it takes to generate netcdf files yourself.

A few notes on the sflux code:

• The maximum numbers of input files and time steps (used in the atmospheric forcing and not to be confued with the time step in main SELFE) are set in the module netcdf_io (you can search for it in the sflux*.f90 code) as max_files and max_times. If your application requires larger values simply increase them in the module.
• The grids used in atmospheric forcings (usually structured) can be different from hgrid.gr3, but must encompass the latter for successful interpolation, which is done in SELFE at run time (both in space and in time).

Also you may find some useful information from other users at the mailing list archive.

go back

14. Why does my non-hydrostatic run blow up?

Here are some tips for a successful non-hydrostatic run, particularly for highly non-hydrostatic problems.

• Make sure the vertical resolution is not much larger than horizontal resolution (i.e. dz <= dx), otherwise the iterative solver (pre-conditioned CG) may not converge. This shouldn't be a problem for most practical applications.
• to get the non-hydrostatic effects right you may need a lot more resolution in regions of interest than the hydrostatic simulations.

go back

15. I have large velocity at open boundary

In hydrostatic models the incoming velocity should be specified at open boundary. Over-specification i.e. elevation+velocity b.c. there usually avoids the problem.

In reality it's often difficult to find velocity b.c., and so a useful way is to use 1-way nesting approach. Run SELFE (in 3D or 2D, barotropic or baroclinic modes; obviously 2D model is the cheapest) for a larger or same grid, with indvel=1, and then use the scripts in utility/OneWayNestScripts/ (in the source code bundle) to generate uv3D.th (follow the instructions in README) together with ifltype=4 in bctides.in. This approach is especially useful when the boundary is away from region of interest.

go back

16. My run crashed. How can I find out why?

The most useful way is to look at animation of the surface velocity to find out potential "trouble" regions. It's best to let the model output until the last step just before the crash (as found in mirror.out) - adjust nspool and ihfskip in param.in to achieve this. You can usually find out the reason easily this way (e.g. boundary condition; wrong atmospheric forcing etc).

go back