23. CROCO/MUSTANG tutorial & tips

23.1. Get to know the CROCO/MUSTANG coupling

Read the documentation on CROCO/MUSTANG (MUSTANG Sediment model)


The MUSTANG learning curve is a steep one. Understanding the documentation strongly benefits from reading the code itself.


In this tutorial $croco refers to the main directory of CROCO source code

23.2. Run a test case

  • Choose a test case (Sediment test cases)

  • Copy the various configuration files that you need

    cp -r $croco/TEST_CASES/ ./TEST_CASES
    cp $croco/OCEAN/cppdefs.h .
    cp $croco/OCEAN/param.h .
    cp $croco/OCEAN/Makefile .
    cp $croco/OCEAN/jobcomp .
  • Modify your jobcomp to point to the location of your CROCO source code

  • Edit the cppdefs.h file, e.g.:

    # define DUNE
    # define MUSTANG

Make sure MUSTANG is activated. For some test cases SEDIMENT (USGS sediment model) is activated by default in cppdefs.h.

23.3. Create your own configuration

  1. First choice: V1 or V2 ?

    If you need bedload - you dont’t have the choice:

    # define key_MUSTANG_V2

    What MUSTANG V2 has to offer:

    • Bedload

    • A new conceptual model for sediment mixture erosion

    • A new model to compute porosity


  2. Modify the param.h file

    Define the number of substances ntrc_subs Define the number of layers ksdmin,*ksdmax*

  3. Create your SUBSTANCE & MUSTANG input files

    Write down the name and location of your files in the croco.in file

    Create the substance file by copying parasubstance_MUSTANG_full_example.txt. Keep only the sections that matter (e.g. don’t keep a listing of sand and gravel parameters if your run only includes muds). Remember that the number of variables should correspond to ntrc_subs in param.h

    Create a user-defined MUSTANG namelist by copying the default one (paraMUSTANG_default.txt). Keep only the parameters that matter for your configuration. If MUSTANG does not find a parameter in the user-defined namelist file, it will use the value defined in the default namelist file.

  4. Modify the cppdefs.h file

    Choose what you want to model with the main CPP keys:

    • Without special CPP key, the model is morpho-static. The seabed evolution does not impact the bathymetry seen by the ocean model. If you want do to morphodynamics run:

      # define MORPHODYN

      This should also activate automatically MORPHODYN_MUSTANG_byHYDRO in the cppdefs_dev.

      Plus, you will have to put l_morphocoupl=.true. in paraMUSTANG*.txt

    • Sand transported in suspension in 3D (no CPP key needed) : that might be very cost effective for regional scale modelling (i.e. if your CROCO time step is large compared to the time step neeeded to guarantee the stability of the explicit settling scheme). Sand transported in suspension in pseudo 2D

      # define key_sand2D
      # define MUSTANG_CORFLUX


    If you want to add a source of sand (e.g. rivers) with the pseudo-2D scheme, it has not been tested yet. Most probably your discharge will only be a fraction of what you wanted. You will need to either adjust the concentration or to modify step3D_t.F in the following section to sum up the water column fluxes in the bottom layer:

    ! Apply point sources for river runoff simulations

    • Read wave files:

      # define WAVE_OFFLINE

    Activates the reading of wave data (this is an existing CROCO CPP option). If combined with #define MUSTANG, it reads significant wave height, wave period, wave direction and bottom orbital velocity. Then the wave-induced bottom shear stress is computed in sed_MUSTANG_CROCO.F90. Note that the significant wave height (or wave amplitude) has to be given as for now but is not used to compute the bed shear stress.

    Header of an example wave file:

    wwv_time = UNLIMITED ; // (2586 currently)
    eta_rho = 623 ;
    xi_rho = 821 ;
    double wwv_time(wwv_time) ;
    double hs(wwv_time, eta_rho, xi_rho) ;
            hs:_FillValue = -32767. ;
    double t01(wwv_time, eta_rho, xi_rho) ;
            t01:_FillValue = -32767. ;
    double dir(wwv_time, eta_rho, xi_rho) ;
            dir:_FillValue = -32767. ;
    double ubr(wwv_time, eta_rho, xi_rho) ;
            ubr:_FillValue = -32767. ;
    • Read netcdf files for solid discharge in river:

      #  define PSOURCE_NCFILE
      #  define PSOURCE_NCFILE_TS

    It reads the concentration values in get_psource_ts.F

    Header of an example source file:

    qbar_time = 7676 ;
    n_qbar = 6 ;
    runoffname_StrLen = 30 ;
    temp_src_time = 8037 ;
    salt_src_time = 8037 ;
    MUD1_src_time = 7676 ;
    double qbar_time(qbar_time) ;
            qbar_time:long_name = "runoff time" ;
            qbar_time:units = "days" ;
            qbar_time:cycle_length = 0 ;
            qbar_time:long_units = "days since 1900-01-01" ;
    double Qbar(n_qbar, qbar_time) ;
            Qbar:long_name = "runoff discharge" ;
            Qbar:units = "m3.s-1" ;
    char runoff_name(n_qbar, runoffname_StrLen) ;
    double temp_src_time(temp_src_time) ;
            temp_src_time:cycle_length = 0 ;
            temp_src_time:long_units = "days since 1900-01-01" ;
    double salt_src_time(salt_src_time) ;
            salt_src_time:cycle_length = 0 ;
            salt_src_time:long_units = "days since 1900-01-01" ;
    double temp_src(n_qbar, temp_src_time) ;
            temp_src:long_name = "runoff temperature" ;
            temp_src:units = "Degrees Celcius" ;
    double salt_src(n_qbar, salt_src_time) ;
            salt_src:long_name = "runoff salinity" ;
            salt_src:units = "psu" ;
    double MUD1_src_time(MUD1_src_time) ;
            MUD1_src_time:long_name = "runoff time" ;
            MUD1_src_time:units = "days" ;
            MUD1_src_time:long_units = "days since 1900-01-01" ;
    double MUD1_src(n_qbar, MUD1_src_time) ;
  6. Initial conditions for the sediment cover

    There are mainly 2 options:

    • Uniform sediment cover

      In paraMUSTANG*.txt:

      l_unised = .true.       ! boolean set to true for a uniform bottom initialization
      fileinised =  './Init.nc' ! File for initialisation (if l_unised is False)
      hseduni = 0.03           ! initial uniform sediment thickness(m)
      cseduni= 1500.0          ! initial sediment concentration
      csed_mud_ini = 550.0      ! mud concentration into initial sediment (if =0. ==> csed_mud_ini=cfreshmud)
      ksmiuni = 1             ! lower grid cell indices in the sediment
      ksmauni = 10             ! upper grid cell indices in the sediment

      And then, the fraction of each sediment variable in the seafloor is defined with cini_sed_n() in parasubsance_MUSTANG.txt

    • Read the sediment cover from a netcdf file or restart from a RESTART file

      In paraMUSTANG*.txt:

      l_repsed=.true.        ! boolean set to .true. if sedimentary variables are initialized from a previous run
      filrepsed='./repsed.nc' ! file from which the model is initialized for the continuation of a previous run

      The netcdf file needs to have the concentration values under the names NAME_sed, with NAME corresponding to the names defined in the SUBSTANCE input files. The number of vertical levels (ksmi, ksma) and the layer thickness (DZS) also need to be defined. The file structure is similar to the RESTART netcdf file, and filerepsed can be used to restart from a CROCO RESTART file.

      Header of an example sediment cover file:

      ni = 821 ;
      nj = 623 ;
      time = UNLIMITED ; // (1 currently)
      level = 10 ;
      double latitude(nj, ni) ;
      double longitude(nj, ni) ;
      double time(time) ;
      double level(level) ;
      double ksmi(time, nj, ni) ;
      double ksma(time, nj, ni) ;
      double DZS(time, level, nj, ni) ;
      double temp_sed(time, level, nj, ni) ;
      double salt_sed(time, level, nj, ni) ;
      double GRAV_sed(time, level, nj, ni) ;
      double SAND_sed(time, level, nj, ni) ;
      double MUD1_sed(time, level, nj, ni) ;

    Alternatively there is a 3rd option possible. If l_repsed=.false. and l_unised=.false., you can specify the filename of your sediment cover dataset (fileinised), but then it is up to you to write yourself the piece of code to read it in initMUSTANG.F90 in the subroutine MUSTANG_sedinit.

    How to prescribe the concentration for the initialisation :

    • Uniform sediment cover. If you use a uniform sediment cover, the initial fraction of each sediment class is read in parasubastance_MUSTANG.txt. Then the concentration of each sediment class is a fraction of cseduni defined in paraMUSTANG.txt (i.e. cv_sed(iv)=cini_sed_n(iv) x cseduni). However, since you prescribe cseduni, it is not necessarily similar to what the model total concentration should be for the same sediment mixture, unless you used the same porosity model as in MUSTANG to compute cseduni.

      With MUSTANG V2, after initialisation, the sediment concentration is adjusted in every layers to match the model porosity law. Hence the initial mass is not preserved, but the bed height and the sediment class fractions are.

      With MUSTANG V1, by default the sediment concentration is not adjusted. In this case, what will happen is that the first time erosion happens, the very first deposit could have a very different porosity to the initial state, and induce an abrupt bed height change. You can select l_init_hsed=.true. to bypass this issue. While adjusting the sediment concentration, it will also adjust the sediment height to conserve the initial mass.


      With MUSTANG V2 we recommend using l_init_hsed=.false. since the subroutine associated with this boolean uses the porosity model of V1.

    • RESTART. If you use l_repsed=.true., l_init_hsed is not even read. In V1, the sediment concentrations that you specify will not be overwritten. It means that you have to start with concentrations that follow the porosity law of the model. In V2, concentrations are overwritten in all layers after computing the porosity for the sediment mixture. In this cas you can specify concentration that are just a fraction of an arbitrary constant total sediment concentration.


    In version 1, you can impose no sediment in a grid cell by imposing ksmi=1 and ksma=0. This could be useful to define reefs for instance. In Version 2 you need at least one sediment layer everywhere. The first layer is never eroded, but is needed to manage the small sediment mass that can be left in the layer just above. To avoid potential issue when computing concentrations for very thin layers, thin layers are merged with underlying layers. Therefore, when initializing sediment concentrations make sure to have at least one layer everywhere.