diff --git a/SETS/RiverTomCrouteinf.c b/SETS/RiverTomCrouteinf.c
new file mode 100755
index 00000000..4566b9a3
--- /dev/null
+++ b/SETS/RiverTomCrouteinf.c
@@ -0,0 +1,217 @@
+#include "math.h"
+#include "mdoodz.h"
+#include "stdbool.h"
+#include "stdlib.h"
+
+/*--------------------------------------------------------------------------------------------------------------------*/
+/*------------------------------------------------------ M-Doodz -----------------------------------------------------*/
+/*--------------------------------------------------------------------------------------------------------------------*/
+
+double SetSurfaceZCoord(MdoodzInput *instance, double x_coord) {
+  double TopoLevel = 0.0;
+  // if (fabs(x_coord) < instance->model.surf_Winc/2.0) TopoLevel = -500.0/instance->scaling.L;
+  return TopoLevel;
+}
+
+/*--------------------------------------------------------------------------------------------------------------------*/
+/*------------------------------------------------------ M-Doodz -----------------------------------------------------*/
+/*--------------------------------------------------------------------------------------------------------------------*/
+
+int SetPhase(MdoodzInput *instance, Coordinates coordinates) {
+  const double H_crust  = instance->model.user1 / instance->scaling.L;
+  const double H_uppercrust  = instance->model.user2 / instance->scaling.L;
+    int phase = 2;
+    
+    if (coordinates.z> -H_crust) phase = 1;
+    if (coordinates.z> -H_uppercrust) phase = 0;
+  
+  return phase;
+}
+
+/*--------------------------------------------------------------------------------------------------------------------*/
+/*------------------------------------------------------ M-Doodz -----------------------------------------------------*/
+/*--------------------------------------------------------------------------------------------------------------------*/
+
+int SetDualPhase(MdoodzInput *input, Coordinates coordinate, int phase) {
+    
+    int    il, dual_phase = phase;
+    double Lx = input->model.xmax - input->model.xmin;
+    double Lz = input->model.zmax - input->model.zmin;
+    double Ax, Az;
+
+    // Set checkerboard for phase 0
+    Ax = cos( 50.0*M_PI*coordinate.x / Lx  );
+    Az = sin( 80.0*M_PI*coordinate.z / Lz  );
+    if ( ( (Az<0.0 && Ax<0.0) || (Az>0.0 && Ax>0.0) ) && dual_phase==0 ) {
+        dual_phase += input->model.Nb_phases;
+    }
+
+    if ( ( (Az<0.0 && Ax<0.0) || (Az>0.0 && Ax>0.0) ) && dual_phase==1 ) {
+        dual_phase += input->model.Nb_phases;
+    }
+    
+  return dual_phase;
+}
+
+/*--------------------------------------------------------------------------------------------------------------------*/
+/*------------------------------------------------------ M-Doodz -----------------------------------------------------*/
+/*--------------------------------------------------------------------------------------------------------------------*/
+
+double SetTemperature(MdoodzInput *instance, Coordinates coordinates) {
+  const double lithosphereThickness = instance->model.user1 / instance->scaling.L;
+  const double surfaceTemperature   = 273.15 / instance->scaling.T;
+  const double mantleTemperature    = (instance->model.user0 + 273.15) / instance->scaling.T;
+  const double particleTemperature  = ((mantleTemperature - surfaceTemperature) / lithosphereThickness) * (-coordinates.z) + surfaceTemperature;
+  if (particleTemperature > mantleTemperature) {
+    return mantleTemperature;
+  } else {
+    return particleTemperature;
+  }
+}
+
+/*--------------------------------------------------------------------------------------------------------------------*/
+/*------------------------------------------------------ M-Doodz -----------------------------------------------------*/
+/*--------------------------------------------------------------------------------------------------------------------*/
+
+double SetHorizontalVelocity(MdoodzInput *instance, Coordinates coordinates) {
+  return -coordinates.x * instance->model.bkg_strain_rate;
+}
+
+/*--------------------------------------------------------------------------------------------------------------------*/
+/*------------------------------------------------------ M-Doodz -----------------------------------------------------*/
+/*--------------------------------------------------------------------------------------------------------------------*/
+
+double SetVerticalVelocity(MdoodzInput *instance, Coordinates coordinates) {
+  return coordinates.z * instance->model.bkg_strain_rate;
+}
+
+/*--------------------------------------------------------------------------------------------------------------------*/
+/*------------------------------------------------------ M-Doodz -----------------------------------------------------*/
+/*--------------------------------------------------------------------------------------------------------------------*/
+
+SetBC SetBCVx(MdoodzInput *instance, POSITION position, Coordinates coordinates) {
+  SetBC bc;
+  const double Lx    = instance->model.xmax - instance->model.xmin;
+  const double V_tot =  Lx * instance->model.bkg_strain_rate; 
+  const double x = coordinates.x, z = coordinates.z;
+
+  // Evaluate velocity of W and E boundaries
+  const double VxW =  0.5*V_tot;
+  const double VxE = -0.5*V_tot;
+  
+  // Assign BC values
+  if (position == N || position == S || position == NW || position == SW || position == NE || position == SE) {
+    bc.type  = constant_shear_stress;
+    bc.value = 0;
+  } else if (position == W) {
+    bc.type  = constant_velocity;
+    bc.value = VxW;
+  } else if (position == E) {
+    bc.type  = constant_velocity;
+    bc.value = VxE;    
+  } else {
+    bc.type  = inside;
+    bc.value = 0.0;  
+  }
+  return bc;
+}
+
+/*--------------------------------------------------------------------------------------------------------------------*/
+/*------------------------------------------------------ M-Doodz -----------------------------------------------------*/
+/*--------------------------------------------------------------------------------------------------------------------*/
+
+SetBC SetBCVz(MdoodzInput *instance, POSITION position, Coordinates coordinates) {
+  SetBC bc;
+  const double Lx = instance->model.xmax - instance->model.xmin;
+  double V_tot    =  Lx * instance->model.bkg_strain_rate; // |VxW| + |VxE|
+ 
+  // Compute bottom boundary velocity and account for erosion
+  const double surf_Winc = instance->model.surf_Winc;
+  const double surf_Vinc = instance->model.surf_Vinc;
+  const double Vz_corr   = surf_Winc*surf_Vinc / Lx;
+  const double VzS = (-0.5 * V_tot * (instance->topo_height->west - instance->model.zmin) / Lx) -(0.5 * V_tot * (instance->topo_height->east - instance->model.zmin) / Lx) + Vz_corr;
+
+  // Set boundary nodes types and values
+  if (position == W || position == SW || position == NW ) {
+    bc.type  = constant_shear_stress;
+    bc.value = 0.0;
+  } else if ( position == E || position == SE || position == NE) {
+    bc.type  = constant_shear_stress;
+    bc.value = 0.0;
+  } else if (position == S || position == N) {
+    bc.type  = constant_velocity;
+    bc.value = VzS;
+  } else {
+    bc.type  = inside;
+    bc.value = 0.0;
+  }
+  return bc;
+}
+
+/*--------------------------------------------------------------------------------------------------------------------*/
+/*------------------------------------------------------ M-Doodz -----------------------------------------------------*/
+/*--------------------------------------------------------------------------------------------------------------------*/
+
+SetBC SetBCT(MdoodzInput *instance, POSITION position, double particleTemperature) {
+  SetBC     bc;
+  double surface_temperature =          zeroC  / instance->scaling.T;
+  double mantle_temperature  = (1330. + zeroC) / instance->scaling.T;
+  if (position == S) {
+    bc.type  = constant_temperature;
+    bc.value = mantle_temperature;
+  }
+  if (position == free_surface || position == N) {
+    bc.type  = constant_temperature;
+    bc.value = surface_temperature;
+  } 
+  if (position == W || position == E) {
+    bc.type  = constant_heatflux;
+    bc.value = 0.;
+  }
+  return bc;
+}
+
+/*--------------------------------------------------------------------------------------------------------------------*/
+/*------------------------------------------------------ M-Doodz -----------------------------------------------------*/
+/*--------------------------------------------------------------------------------------------------------------------*/
+
+void AddCrazyConductivity(MdoodzInput *input) {
+  int               *asthenospherePhases = (int *) malloc(sizeof(int));
+  CrazyConductivity *crazyConductivity   = (CrazyConductivity *) malloc(sizeof(CrazyConductivity));
+  asthenospherePhases[0]                 = 3;
+  crazyConductivity->phases              = asthenospherePhases;
+  crazyConductivity->nPhases             = 1;
+  crazyConductivity->multiplier          = 1000;
+  input->crazyConductivity               = crazyConductivity;
+}
+
+/*--------------------------------------------------------------------------------------------------------------------*/
+/*------------------------------------------------------ M-Doodz -----------------------------------------------------*/
+/*--------------------------------------------------------------------------------------------------------------------*/
+
+int main() {
+  MdoodzSetup setup = {
+          .BuildInitialTopography = &(BuildInitialTopography_ff){
+                  .SetSurfaceZCoord = SetSurfaceZCoord,
+          },
+          .SetParticles = &(SetParticles_ff){
+                  .SetPhase              = SetPhase,
+                  .SetTemperature        = SetTemperature,
+                  .SetHorizontalVelocity = SetHorizontalVelocity,
+                  .SetVerticalVelocity   = SetVerticalVelocity,
+                  .SetDualPhase          = SetDualPhase,
+          },
+          .SetBCs = &(SetBCs_ff){
+                  .SetBCVx    = SetBCVx, // manuel
+                  .SetBCVz    = SetBCVz,
+                  .SetBCT     = SetBCT,
+          },
+          .MutateInput = AddCrazyConductivity,
+
+  };
+  RunMDOODZ("RiverTomCrouteinf.txt", &setup);
+}
+
+/*--------------------------------------------------------------------------------------------------------------------*/
+/*------------------------------------------------------ M-Doodz -----------------------------------------------------*/
+/*--------------------------------------------------------------------------------------------------------------------*/
diff --git a/SETS/RiverTomCrouteinf.txt b/SETS/RiverTomCrouteinf.txt
new file mode 100644
index 00000000..554cbd70
--- /dev/null
+++ b/SETS/RiverTomCrouteinf.txt
@@ -0,0 +1,129 @@
+/**** RESTART ****/
+istep = 00020
+irestart = 0
+
+/**** OUTPUT FILES ****/
+writer         = 1
+writer_step    = 50
+writer_markers = 0
+
+/**** SCALES ****/ Scaling parameters for
+eta = 1.e+22      /viscosity
+L   = 1.e4        /length
+V   = 1.0e-9      /velocity
+T   = 1.e+02      /temperature
+ 
+/**** SPACE****/
+Nx      = 251 
+Nz      = 81
+xmin    = -250.000000e3
+xmax    =  250.000000e3
+zmin    = -150.000000e3
+zmax    =   10.000000e3 
+
+/**** TIME ****/
+Nt          = 100000
+dt          = 2.63e10
+constant_dt = 0
+Courant     = 0.25
+RK          = 4
+
+/**** SWITCHES ****/
+mechanical        = 1
+pure_shear_ALE    = -1  /1: box stretched at constant strain rate; -1 box of constant size
+elastic           = 1
+thermal           = 1
+free_surface      = 1
+free_surface_stab = 1.0
+initial_cooling   = 1
+subgrid_diffusion = 2
+shear_heating     = 1
+adiab_heating     = 0
+track_T_P_x_z     = 1 / initial P and T field record
+conserv_interp    = 0
+
+/**** SURFACE PROCESSES ****/
+surface_processes  = 5 /0 = none; 1 = diffusion only; 2 = fills instantaneously the basin up to the base level;  3 = diffusion + source term
+surf_diff          = 6e-6 /topography diffusion coefficient
+surf_Winc          = 1e3
+surf_Vinc          = 5e-9
+
+/**** SETUP DEPENDANT ****/
+bkg_strain_rate = 7.5e-16 / Background strain rate [s-1] 
+user0           = 1330.0    / Mantle temperature    
+user1           = 60.0e3    / Crust thickness    
+user2           = 20.0e3    /Upper crust Thickness
+
+
+/**** GRAVITY ****/
+gx = 0.0000
+gz = -9.81
+
+/**** PHASE PROPERTIES ****/
+Nb_phases = 3 
+
+ID            = 0       / upper Crust
+rho           = 2800.00 / ref. density
+G             = 3e10    / shear modulus
+alp           = 30.0e-6 / thermal expansivity
+bet           = 1e-11   / compressibility
+Cv            = 1050    / heat capacity
+k             = 2.0     / thermal conductivity
+Qr            = 1e-7  / radiogenic heat production
+C             = 5e7     / cohesion
+phi           = 30.00   / friction angle
+eta_vp        = 1.333e21    / viscoplastic viscosity
+pwlv          = 10      / power law viscosity -> database flow_laws.c
+density_model = 3
+
+ID            = 1       / Lower Crust
+rho           = 2800.00 / ref. density
+G             = 3e10    / shear modulus
+alp           = 30.0e-6 / thermal expansivity
+bet           = 1e-11   / compressibility
+Cv            = 1050    / heat capacity
+k             = 2.0     / thermal conductivity
+Qr            = 1e-7  / radiogenic heat production
+C             = 5e7     / cohesion
+phi           = 30.00   / friction angle
+eta_vp        = 1.333e21    / viscoplastic viscosity
+pwlv          = 10      / power law viscosity -> database flow_laws.c
+density_model = 3
+
+ID            = 2       / Mantle
+rho           = 3260.00 / ref. density
+G             = 3e10    / shear modulus
+alp           = 30.0e-6 / thermal expansivity
+bet           = 1e-11   / compressibility
+Cv            = 1050    / heat capacity
+k             = 3.0     / thermal conductivity
+Qr            = 1.0e-10 / radiogenic heat production
+C             = 5e7     / cohesion
+phi           = 30.00   / friction angle
+eta_vp        = 1.333e21    / viscoplastic viscosity
+pwlv          = 40      / power law viscosity -> database flow_laws.c
+linv          = 40      / linear creep viscosity -> database flow_laws.c
+expv          = 40      / exponential creep viscosity -> database flow_laws.c
+density_model = 3
+
+/**** PARTICLES ****/
+Nx_part  = 4
+Nz_part  = 4
+min_part_cell = 16
+
+/**** LINEAR SOLVER ****/
+lin_abs_mom = 5.0e-8
+lin_abs_div = 5.0e-8
+lin_rel_mom = 5.0e-8
+lin_rel_div = 5.0e-8
+penalty     = 1e3
+
+/**** NON-LINEAR ITERATIONS ****/
+Newton         = 1
+Newton2Picard  = 1
+nit_max        = 20 // for high resolution use from 10 to 20
+line_search    = 0
+nonlin_abs_mom = 5.0e-8
+nonlin_abs_div = 5.0e-8
+min_eta        = 5.0e18
+max_eta        = 1.0e25