
MUSTANG tips
------------

Get to know the CROCO/MUSTANG coupling
++++++++++++++++++++++++++++++++++++++

Read the documentation on CROCO/MUSTANG (:ref:`CROCO_MUSTANG`)

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


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

Run a test case
++++++++++++++++

* Choose a test case (:ref:`sediment`)

* Copy the various configuration files that you need 
  ::

         cp -r $croco/MUSTANG/MUSTANG_NAMELIST/ ./MUSTANG_NAMELIST
         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.


Create your own configuration
++++++++++++++++++++++++++++++

 
#. 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 

     http://www.ifremer.fr/docmars/html/doc_MUSTANG/doc.MUSTANG.process.html

#. Modify the param.h file

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

#. 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. 

#. 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

     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 
     needed to guarantee the stability of the explicit settling scheme). 
     Sand transported in suspension in pseudo 2D 
     ::

           # define key_sand2D
           # define MUSTANG_CORFLUX
           
     .. warning ::

        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
           !----------------------------------------------------------



#. CROCO/MUSTANG CPP keys
   
   * 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: 
     ::

                dimensions:
            wwv_time = UNLIMITED ; // (2586 currently)
            eta_rho = 623 ;
            xi_rho = 821 ;
                variables:
            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: 
     ::

                dimensions:
            qbar_time = 7676 ;
            n_qbar = 6 ;
            runoffname_StrLen = 30 ;
            temp_src_time = 8037 ;
            salt_src_time = 8037 ;
            MUD1_src_time = 7676 ;
                variables:
            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) ;


#. 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: 
     ::

                    dimensions:
                ni = 821 ;
                nj = 623 ;
                time = UNLIMITED ; // (1 currently)
                    level = 10 ;
                    variables:
                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. 

  .. note :: 
     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.

  .. warning :: 
     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.  
