17.1. Some important points about Large-Eddy Simulations (LES)

  • Momentum and tracer advection schemes :

Many advection schemes are now implented in CROCO (i.e. section 4.3 of the model documentation) leading to one recurrent question: which scheme is the most apropriate for my configuration? Unfortunately, no advection scheme is perfect for all applications.

In eddy-resolving ocean simulations significant gradients due to detached eddies or upwelling can be found. More particularly, in coastal eddy-resolving configurations, salinity may vary from river concentration to ocean concentration within a few kilometres in horizontal leading to the formation of even more stronger gradients. In such cases (if your solution has strong gradients, shocks or propagating fronts), it is recommend to use a total variation bounded ( or a monotonicity-preserving scheme (section

An exemple of numerical artefacs related to the advection schemes: Gibbs phenomenon

Numerical experiments with non-monotonic scheme show artificial oscillations of the solution near regions of sharp gradients (figure below).


Non-monotonic vertical advection schemes (Akima for TS, Spline for UV)

On the contrary, TVD and WENO5 schemes enable sharper shock predictions and as they preserve monotonicity they do not generate spurious oscillations in the solution (figure below).


Monotonic or quasi-monotonic vertical advection schemes (WENO5 for TS, TVD for UVW)

Recommended advection schemes for LES :

CPP options of Momentum Advection


Activate 5th-order WENOZ quasi-monotone lateral advection scheme for UV


Activate 5th-order WENOZ quasi-monotone vertical advection scheme for UV


Activate 5th-order WENOZ quasi-monotone lateral advection scheme for W (in NBQ simulation)


Activate 5th-order WENOZ quasi-monotone vertical advection scheme for W (in NBQ simulation)



Activate Total Variation Diminushing lateral advection scheme for UV


Activate Total Variation Diminushing vertical advection scheme for UV


Activate Total Variation Diminushing lateral advection scheme for W (in NBQ simulation)


Activate Total Variation Diminushing vertical advection scheme for W (in NBQ simulation)

CPP options of Tracer advection


Activate 5th-order WENOZ quasi-monotone lateral tracer advection scheme


Activate 5th-order WENOZ quasi-monotone vertical tracer advection scheme

  • Turbulence schemes : MILES & LES approaches

In LES, direct transfer ends at the lowest scale resolved, and subgrid dissipation of energy is accomplished by implicit mixing of advection schemes, as well as by explicit parametrization provided by turbulent closure schemes. The choices of advection schemes and/or turbulent closure schemes are thus critical to represent correctly the turbulent energy cascade.

Small scales tend to be more isotropic and homogeneous than the large ones, thus LES requires 3D turbulent closure schemes. Two options of 3D turbulent closure schemes are available in CROCO : A generic two-equation turbulence closure model called Generalized Length Scale (GLS) scheme & a Smagorinsky model (i.e. model documentation). An alternative approach is monotonically integrated LES (MILES). In MILES, the dissipative nature of monotonic advection schemes is exploited to provide an implicit model of turbulence.

Related CPP options (for users):


Activate 3D Generic Length Scale scheme


Activate 3D Smagorinsky SGS model

  • Options of Bottom boundary layer

In coastal seas, the bottom mixed layers may occupy a considerable fraction of the water depth. In contrast, bottom mixed layers in ocean bassins cover only a small portion of the total depth of several thousands of meters. Moreover the strong dissipation of kinetic energy generated by the bed friction can be enhanced in shallow water. Hence the parametrization of the bottom boundary layer dynamic is particularly imortant in coastal large eddy simulations. Some new parametrization options are under development in CROCO to potentially improve the representation of the bottom boundary layers.

BSTRESS_FAST allows solving the bottom friction term of the momentum equations at the fast time step (using part of the code structure inherited from the Non-Boussinesq solver). It avoids reducing the slow-mode (baroclinic) time step for cases with high bottom friction or/and high near bottom vertical resolution. This is not yet a default option as it needs further evaluation in various configurations.

NBQ_FREESLIP imposes a free-slip boundary condition on the bottom (the normal component of the fluid velocity field is set to zero at the bottom level but the tangential component is unrestricted). This is not a default option, by default a no-slip condition is imposed on the bottom. Further evaluation in various configurations is needed.

Related CPP options (for users):


Activate free-slip boundary condition on the bottom


solve the bottom friction term of the momentum equations at the fast time step

17.2. KH_INST Test Case

  1. Create a configuration directory:
    mkdir ~/CONFIGS/KH_INST
  2. Copy the input files for compilation from croco sources:
    cp ~/croco/croco/OCEAN/cppdefs.h .
    cp ~/croco/croco/OCEAN/param.h .
    cp ~/croco/croco/OCEAN/jobcomp .
  3. Edit cppdefs.h for using KH_INST case
    # define KH_INST
    # undef REGIONAL

    Explore the CPP options selected for KH_INST case and undef MPI

    after #elif defined KH_INST
    # undef MPI

    You can check the KH_INST settings in param.h.

  4. Edit the compilation script jobcomp:
    # set source, compilation and run directories
    # determine operating system
    # compiler options
    # set MPI directories if needed
    MPIDIR=$(dirname $(dirname $(which $MPIF90) ))
    MPILIB="-L$MPIDIR/lib -lmpi -limf -lm"
    # set NETCDF directories
    # Use :
    #-lnetcdf           : version netcdf-3.6.3                --
    #-lnetcdff -lnetcdf : version netcdf-4.1.2                --
    #-lnetcdff          : version netcdf-fortran-4.2-gfortran --
    #NETCDFLIB="-L/usr/local/lib -lnetcdf"
    NETCDFLIB=$(nf-config --flibs)
    NETCDFINC=-I$(nf-config --includedir)
  5. Compile the model:
    ./jobcomp > jobcomp.log

    If compilation is successful, you should have a croco executable in your directory.

    You will also find a Compile directory containing the model source files:

    • .F files: original model source files that have been copied from ~/croco/croco/OCEAN

    • _.f files: pre-compiled files in which only parts defined by cpp-keys are kept

    • .o object files

  6. Copy the namelist input file for KH_INST case:
    cp ~/croco/croco/TEST_CASES/croco.in.KH_INST croco.in

    Eventually edit it.

  7. Run the model:
    ./croco croco.in > croco.out

    If your run is successful you should obtain the following files:

    khinst_rst.nc # restart file
    khinst_his.nc # instantaneous output file
  8. Have a look at the results:
    ncview khinst_his.nc
  9. Test: some questions:

    • What is the impact of the relaxation of the non-hydrostatic hypothesis?

    • What are the impacts of the advection schemes?

    • What is the impact of adding a turbulent scheme?

17.3. Set up your own NBQ configuration

  • In cppdefs.h you should activate

  • NBQ : activate the non-Boussinesq and non-hydrostatic kernel

/* Non-Boussinesq */
# define NBQ
  • To set up adapted time steps to your NBQ configuration (dt & NDTFAST in croco.in file), you can activate in cppdefs_dev.h

  • DIAG_CFL : activate diagnostics of the CFL criteria

    # define DIAG_CFL

If DIAG_CFL is defined, at each NINFO during the run, CFL criteria are indicated in your output file :

  • INT_3DADV : Slow (baroclinic) mode CFL criterion. This parameter depends on your mesh grid size and your ocean current intensity (time-varying diagnostic). It should be inferior to approximately 1 (depending on the advection scheme, i.e. section 4.2.5).

  • EXT_GWAVES : CFL criterion based on the barotropic wave speed. It should be inferior to 0.89 (i.e. section 4.2.5).

  • NBQ_HADV : CFL criterion based on the pseudo-acoustic wave speed. This parameter should be inferior to 1.7.

  • Compile your model

  • Edit croco.in file, add the following line

    time_stepping_nbq: NDTNBQ    CSOUND_NBQ    VISC2_NBQ
                      1          "5xsqrt(gHmax)"         1.e-2
    • NDTNBQ : irrelevant parameter

    • CSOUND_NBQ : Pseudo-acoustic waves speed. This parameter should be at least superior to five times the barotropic wave speed (sqrt(gHmax)) in your domain and inferior or equal to the acoustic waves speed (1500 m/s).

    • VISC2_NBQ : Bulk Viscosity (i.e. section 1.4)

  • Run your simulation.

If it’s blow-up, change your time steps to respect the CFL criteria (increase NDTFAST such as NBQ_HADV<1.7). Relaxing the hydrostatic hypothesis, change the dynamic of the small-scale processes and thus potentially the intensity of the small-scale currents leading to more drastic baroclinic CFL conditions. So if it’s still blow up, reduce the baroclinic time step (dt).


  • In which cases, do I need to activate the NBQ Precise option?

Two versions of CROCO-NBQ are currently available: NBQ_PERF & NBQ_PRECISE. NBQ_PERF solve the compressible and non-hydrostatic Navier-Stokes equations and conserve precisely the volume. This version is the most efficient in terms of computational time. NBQ_PRECISE solve also the compressible and non-hydrostatic Navier-Stokes equations and conserve precisely the mass. However, this version is more time consuming (to conserve precisely the mass, an update of the sigma vertical grid at each fast time step is needed). By default, the NBQ_PERF option is defined in the cppdefs_dev.h file. In regional or coastal configurations (resolution ranges of 50-300m), no significant differences in terms of oceanic dynamics have been observed so far. However, specific care is needed when surface waves are explicitly represented. In such configurations, the fluctuations of the vertical grid are more pronounced and NBQ_PRECISE considerably improve the representation of the surface waves. At resolution ranges of 1m, further investigations are needed.

Related CPP options (for users):


The most efficient version in terms of computational time


The most precise version in terms of mass conservation

  • Options of open boundaries conditions

An Orlanski radiation condition (OBC_NBQORLANSKI) has been applied to the internal mode velocities, temperature, and salinity at the open boundaries. Whereas barotropic and acoustic waves are radiated through the boundary using the methods of characteristics. It is recommended to use a sponge layer to deal with strong nonlinearities (as for example to avoid reflexion of solitary waves at the lateral boundaries).

Related CPP options (for users):




Radiative conditions


interior/bdy forcing/nudging


interior/bdy forcing/nudging


bdy forcing/nudging

17.5. Appendix : some words on CROCO-NBQ kernel

In CROCO-NBQ, the ”fast mode” includes in addition to the external (barotropic) mode, the pseudo-acoustic mode that allows computation of the nonhydrostatic pressure within a non-Boussinesq approach (Auclair et al., 2018). A two-level time-splitting kernel is thus conserved, but the fast time step integrates a 3D-compressible flow. Hence, acoustic waves or “pseudo-acoustic” waves have indeed been re-introduced to avoid Boussinesq-degeneracy which inevitably leads to a 3D Poisson-system in non-hydrostatic Boussinesq methods and to reduce computational costs. As long as “pseudo-acoustic” waves remain faster than the fastest physical processes in the domain, their phase-velocity can artificially be slowed down rendering unphysical high-frequency processes associated with bulk compressibility but preserving a coherent slow non-hydrostatic dynamics with a softening of the CFL criterion. More details are given on http://poc.omp.obs-mip.fr/auclair/WOcean.fr/SNH/Pub/Tutorials/CROCO/Html_maps/Croco2018_map.html.

Auclair, F., Bordois, L., Dossmann, Y., Duhaut, T., Paci, A., Ulses, C., Nguyen, C., 2018. A non-hydrostatic non-Boussinesq algorithm for free-surface ocean modelling. Ocean Modelling 132, 12–29. https://doi.org/10.1016/j.ocemod.2018.07.011

Related CPP options (for developers):


The equation of motion for vertical velocity is solved implicitly in the vertical direction.


The semi-implicit theta method is used to reduce the numerical dissipation iduced by the implicitation of the vertical velocity equation in the vertical direction (i.e. Fringer et al. 2006)


Prognostic the grid evolution


Classical fourth-order Adams-Moulton (AM4) time-stepping method


Forward-Backward time-stepping method


Perfect conservation of mass (undef NBQ_MASS : perfect conservation of volume)


The sigma vertical grid is updated at each fast time step to reflect the newly solved elevations (as the free surface is now explicitly resolved at each fast time step).


The sigma vertical grid is updated only at each slow time step (reduce the computational time).


Trick to change the name of a variable in the equation of mass conservation