/* ReSCAL - Predefined models
 *
 * Copyright (C) 2011
 *
 * Author: Olivier Rozier <rozier@ipgp.fr>
 *
 * This file is part of ReSCAL.
 *
 * ReSCAL is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
 */


// USER'S GUIDE
// ------------
//
// Basically a model is a set of doublet or cell transitions.
// Here we describe all the functions that help define those transitions.
//
// trans(direction, cel1_start, cel2_start, cel1_end, cel2_end, rate):
//   defines a transition of oriented doublet [cel1_start, cel2_start] -> [cel1_end, cel2_end] occuring at the specified rate.
//   Possible values for directions are EST_OUEST, NORD_SUD, VERTICAL, HORIZONTAL, ISOTROPE.
//
// trans_ref( ref,  directions,  cel1_start,  cel2_start,  cel1_end,  cel2_end,  rate):
//   defines a transition of oriented doublet [cel1_start, cel2_start] -> [cel1_end, cel2_end] occuring at the specified rate,
//   with a unique reference number 'ref' (arbitrary positive integer) which identifies the transition.
//
// trans_cel( cel_start,  cel_end,  rate):
//   defines a cell transition [cel_start] -> [cel_end] occuring at the specified rate.
//
// trans_cel( ref,  cel_start,  cel_end,  rate):
//   defines a cell transition [cel_start] -> [cel_end] occuring at the specified rate, with a unique reference number 'ref'
//
// trans_link( ref1,  ref2,  ncel, prob):
//   creates a probability chain rule from doublet transition 'ref1' to doublet transition 'ref2'.
//   The occurrence of transition 'ref1' generates the occurence of transition 'ref2' with probability 'prob' in its surrounding area.
//   The value 'ncel' (1 or 2) specifies the position of the cell which is linking both transitions. Its state must be the same in the ending and starting doublets of the first and second transitions, respectively.
//   Chain of aligned doublets is not supported yet.
//
// trans_regul( ref,  callback_reg):
//   regulates the global rate of the doublet transition 'ref' by a callback function 'callback_reg'.
//
// trans_cel_regul( ref,  callback_reg):
//   regulates the global rate of the cell transition 'ref' by a callback function 'callback_reg'.
//
// trans_check( ref,  callback_chk,  ncel, time_mode):
//   locally controls the doublet transition 'ref' by a callback function 'callback_chk'. The transition is blocked if a condition is not met.
//   The value 'ncel' (1 or 2) specifies the position of the cell to be tested in the starting doublet of the transition.
//   Possible values for 'time_mode' are TIME_EVOL, TIME_CORR and TIME_NO_EVOL. It indicates whether the transition is taken into account for time evolution (TIME_EVOL) or not (TIME_NO_EVOL).
//   It is highly recommended to use the value TIME_NO_EVOL, except with check_grad_vel() function where the simulated time is automatically corrected if time_mode=TIME_CORR.
//
// trans_check_color( ref,  callback_chk,  ncel, time_mode, lambda_col):
//   declares that the transition 'ref' may have a different rate depending on the color of a cell.
//   The only possible value for 'callback_chk' is check_color.
//   The value 'ncel' (1 or 2) specifies the position of the cell to be tested in the starting doublet of the transition.
//   Possible values for 'time_mode' are TIME_EVOL and TIME_NO_EVOL. It is highly recommended to use the value TIME_NO_EVOL here.
//   The argument 'lambda_col' is the rate of the transition for the colored cells.
//
// trans_type( ref, type):
//   forces the type of the doublet transition 'ref'.
//   The only possible value for 'type' is TR_TRANSPORT. Generally, the transport transitions are automatically detected, unless the state of the moving cell changes.
//


#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "defs.h"
#include "macros.h"
#include "param.h"
#include "transitions.h"
#include "space.h"
#ifdef LGCA
#include "lgca.h"
#endif

extern int Ncel[];    // nombre de cellules par type
extern int LNS;       //largeur nord-sud
extern int ava_trans;        //transitions d'avalanches
extern float ava_delay;
extern int boundary;  //boundary conditions
extern int pbc_mode; //periodic boundary conditions

char *model_name = NULL; //name of model

#if defined(MODELE_DUN)
float lambda_E, lambda_C, lambda_D, lambda_G, lambda_J; //rates
float lambda_T = -1;
float lambda_I = 1.0;
#ifdef CELL_COLOR
float lambda_E_col = -1;
float lambda_C_col = -1;
float lambda_T_col = -1;
float lambda_D_col = -1;
#endif
#elif defined(MODELE_AVA)
float lambda_I = 0.1;
#else
float lambda_I, lambda_D, lambda_L, lambda_C, lambda_R, lambda_E_v, lambda_E_h;  //rates
#endif
float lambda_A = 1.0;
float lambda_A_unstable = 0;
#ifdef MODELE_LIFE
float lambda_B = 1.0;
float lambda_D = 1.0;
extern int B_min, B_max, S_min, S_max;
#endif
float gravite;  // intensite de la gravite
float tr_v, tr_h, diff_h, debit;     //termes de transports verticaux/horizontaux, diffusion horizontale et debit

#if defined(MODELE_CMB) || defined(MODELE_RIV)
float flux_eo, flux_ns;  // intensite flux est-ouest, flux nord-sud
#else
float flux;                        // intensite flux horizontal
#endif

#if defined(LGCA) && defined(MODELE_CMB)
float prob_asv_D, prob_asv_L, prob_asv_C;  //asservissement des transitions vis-a-vis du flux
#endif

#if defined(MODELE_DUN)
float coef_injection = 1.0; //regulation de l'injection de grain
#endif

//callbacks for the regulation of transitions
#ifdef MODELE_DUN
Callback_regul callback_bord_dun; //automatic reinjection of sand grains
Callback_regul callback_injection_coef; //regulation of the injection of sand grains
#endif
#ifdef ALTI
Callback_regul callback_ava_trans; //regulation of the transitions of avalanches
#endif

//callbacks for the control of transitions
#ifdef ALTI
Callback_check check_alti;
Callback_check check_grain_seul;
Callback_check check_slope;  //checks the slope in the east-west direction (for the erosion)
Callback_check check_ava; //checks the value of the steepest slope in all the directions (for the avalanches)
#endif
#ifdef CELL_COLOR
Callback_check check_color; //checks the color of a cell
#endif
#ifdef LGCA
Callback_check check_mvt_bas;
Callback_check check_mvt_haut;
Callback_check check_mvt_est;
Callback_check check_mvt_ouest;
Callback_check check_mvt_est_et_bas;
Callback_check check_no_mvt_est_et_bas;
Callback_check check_mvt_mean_bas;
Callback_check check_mvt_mean_est;
Callback_check check_no_mvt_mean_bas;
Callback_check check_mvt_EB;
Callback_check check_grad_vel;
#ifdef CELL_COLOR
Callback_check check_grad_vel_color;
#endif
#endif
#ifdef MODELE_LIFE
Callback_check check_birth;
Callback_check check_death;
#endif



void params_modele()
{
#ifdef MODELE_DUN
  //lambdas
  parametre("Lambda_E", "erosion rate", &lambda_E, PARAM_FLOAT);
  parametre("Lambda_C", "deposition rate", &lambda_C, PARAM_FLOAT);
  parametre("Lambda_T", "transport rate", &lambda_T, PARAM_FLOAT);
  parametre("Lambda_D", "diffusion rate", &lambda_D, PARAM_FLOAT);
  parametre("Lambda_G", "gravitation rate", &lambda_G, PARAM_FLOAT);
  parametre("Lambda_I", "injection rate - optional", &lambda_I, PARAM_FLOAT);
#ifdef AVALANCHES
  parametre("Lambda_A", "avalanches rate", &lambda_A, PARAM_FLOAT);
  parametre("Lambda_A_unstable", "unstable avalanches rate", &lambda_A_unstable, PARAM_FLOAT);
#endif
#ifdef CELL_COLOR
  parametre("Lambda_E_col", "erosion rate for white grains (Lambda_E by default)", &lambda_E_col, PARAM_FLOAT);
  parametre("Lambda_C_col", "deposition rate for white grains (Lambda_C by default)", &lambda_C_col, PARAM_FLOAT);
  parametre("Lambda_T_col", "transport rate for white grains (Lambda_T by default)", &lambda_T_col, PARAM_FLOAT);
  parametre("Lambda_D_col", "diffusion rate for white grains (Lambda_D by default)", &lambda_D_col, PARAM_FLOAT);
#endif
  //parametre("Lambda_D", "non-utilisé", &lambda_D, PARAM_FLOAT);
  //parametre("Lambda_J", "non-utilisé", &lambda_J, PARAM_FLOAT);
#endif

#ifdef MODELE_AVA
  //lambdas
  parametre("Lambda_I", "injection rate - optional", &lambda_I, PARAM_FLOAT);
  parametre("Lambda_A", "avalanches rate", &lambda_A, PARAM_FLOAT);
  parametre("Lambda_A_unstable", "unstable avalanches rate", &lambda_A_unstable, PARAM_FLOAT);
#endif //MODELE_AVA

#ifdef MODELE_RIV
  //lambdas
  parametre("Lambda_I", "coefficient d'infiltration", &lambda_I, PARAM_FLOAT);
  parametre("Lambda_D", "coefficient de diffusion boue/eau", &lambda_D, PARAM_FLOAT);
  parametre("Lambda_L", "coefficient de dissolution", &lambda_L, PARAM_FLOAT);
  parametre("Lambda_E_v", "coefficient d'erosion verticale par transport", &lambda_E_v, PARAM_FLOAT);
  parametre("Lambda_E_h", "coefficient d'erosion horizontale par transport", &lambda_E_h, PARAM_FLOAT);
  parametre("Lambda_C", "coefficient de cristallisation", &lambda_C, PARAM_FLOAT);
  parametre("Gravite", "coefficient de gravite boue/eau", &gravite, PARAM_FLOAT);
  parametre("Flux_eo", "intensite du flux est-ouest", &flux_eo, PARAM_FLOAT);
  parametre("Flux_ns", "intensite du flux nord-sud", &flux_ns, PARAM_FLOAT);
  parametre("Debit", "debit a la source", &debit, PARAM_FLOAT);
#endif

#ifdef MODELE_CHE
  //lambdas
  parametre("Lambda_I", "coefficient d'infiltration", &lambda_I, PARAM_FLOAT);
  parametre("Lambda_D", "coefficient de diffusion boue/eau", &lambda_D, PARAM_FLOAT);
  parametre("Lambda_L", "coefficient de dissolution", &lambda_L, PARAM_FLOAT);
  parametre("Lambda_C", "coefficient de cristallisation", &lambda_C, PARAM_FLOAT);
  parametre("Gravite", "coefficient de gravite boue/eau", &gravite, PARAM_FLOAT);
  parametre("Flux", "transport horizontal boue/eau", &flux, PARAM_FLOAT);
  parametre("Tr_v", "transport vertical eau/air", &tr_v, PARAM_FLOAT);
  parametre("Tr_h", "transport horizontal eau/air", &tr_h, PARAM_FLOAT);
  parametre("Diff_h", "diffusion horizontale eau/air", &diff_h, PARAM_FLOAT);
  parametre("Debit", "debit de l'eau en entree", &debit, PARAM_FLOAT);
#endif

#ifdef MODELE_CMB
  //lambdas
  parametre("Lambda_I", "coefficient d'infiltration", &lambda_I, PARAM_FLOAT);
  parametre("Lambda_D", "coefficient de diffusion", &lambda_D, PARAM_FLOAT);
  parametre("Lambda_L", "coefficient de dissolution", &lambda_L, PARAM_FLOAT);
  parametre("Lambda_C", "coefficient de cristallisation", &lambda_C, PARAM_FLOAT);
  parametre("Gravite", "coefficient de gravite", &gravite, PARAM_FLOAT);
  parametre("Flux_eo", "intensite du flux est-ouest", &flux_eo, PARAM_FLOAT);
  parametre("Flux_ns", "intensite du flux nord-sud", &flux_ns, PARAM_FLOAT);
#endif

#ifdef MODELE_CRI
  //lambdas
  parametre("Lambda_I", "coefficient d'infiltration", &lambda_I, PARAM_FLOAT);
  parametre("Lambda_D", "coefficient de diffusion", &lambda_D, PARAM_FLOAT);
  parametre("Lambda_L", "coefficient de dissolution", &lambda_L, PARAM_FLOAT);
  parametre("Lambda_C", "coefficient de cristallisation", &lambda_C, PARAM_FLOAT);
  parametre("Gravite", "coefficient de gravite", &gravite, PARAM_FLOAT);
  parametre("Flux", "intensite du flux est-ouest", &flux, PARAM_FLOAT);
#endif

#ifdef MODELE_ICB
  //lambdas
  //parametre("lambda_P", "coefficient de compaction", &lambda_P, PARAM_FLOAT);
  parametre("Lambda_R", "coefficient de refroidissement", &lambda_R, PARAM_FLOAT);
  parametre("Lambda_D", "coefficient de diffusion", &lambda_D, PARAM_FLOAT);
  parametre("Lambda_L", "coefficient de dissolution", &lambda_L, PARAM_FLOAT);
  parametre("Lambda_C", "coefficient de cristallisation", &lambda_C, PARAM_FLOAT);
  parametre("Gravite", "ntensite de la gravite", &gravite, PARAM_FLOAT);
  parametre("Flux", "intensite du flux horizontal", &flux, PARAM_FLOAT);
#endif

#ifdef MODELE_DIF
  //lambdas
  parametre("Lambda_D", "coefficient de diffusion", &lambda_D, PARAM_FLOAT);
  parametre("Debit", "debit a la source", &debit, PARAM_FLOAT);
#endif

#ifdef MODELE_LIFE
  //lambdas
  parametre("Lambda_B", "birth rate", &lambda_B, PARAM_FLOAT);
  parametre("Lambda_D", "death rate", &lambda_D, PARAM_FLOAT);
  parametre("B_min", "min. number of neighbours for birth transitions", &B_min, PARAM_INT);
  parametre("B_max", "max. number of neighbours for birth transitions", &B_max, PARAM_INT);
  parametre("S_min", "min. number of neighbours that prevents death transitions", &S_min, PARAM_INT);
  parametre("S_max", "max. number of neighbours that prevents death transitions", &S_max, PARAM_INT);
#endif
}

void init_modele()
{
  //params_modele();

  if (model_name && strcmp(model_name, MOD_NAME)){
    ErrPrintf("ERROR: Wrong model name %s in parameter file. Check model in defs.h and rebuild all.\n", model_name);
    exit(-2);
  }

#if defined(MODELE_CMB) || defined(MODELE_RIV)

  double flux_e, flux_o, flux_n, flux_s;

  if (flux_eo >= 0){
    flux_o = 0;
    flux_e = flux_eo;
  }
  else{
    flux_o = -flux_eo;
    flux_e = 0;
  }

  if (flux_ns >= 0){
    flux_n = 0;
    //flux_n = flux_ns;//0; //ortho
    flux_s = flux_ns;
  }
  else{
    flux_n = -flux_ns;
    flux_s = 0;
  }
#else
  double flux1, flux2;

  if (flux >= 0){
    flux1 = 0;
    flux2 = flux;
  }
  else{
    flux1 = -flux;
    flux2 = 0;
  }
#endif

  gravite += 1.0; //= 1.05;

  //printf("ReSCAL - %s model\n",MOD_NAME);

#if defined(TRANSITION_BORD) && defined(CYCLAGE_HOR)
  ErrPrintf("ERROR: do not activate CYCLAGE_HOR and TRANSITION_BORD simultaneously !\n");

  exit(-10);
#endif


#ifdef MODELE_DUN

/* Parameters:
  lambda_E   	erosion
  lambda_C    deposition
  lambda_T    transport
  lambda_D	  diffusion (NS)
  lambda_G	  gravity
  lambda_I    injection
  lambda_A    avalanches
  lambda_J    not used
  _______________________________________________

  7 lambdas + 1 not used

*/

  char tmode = TIME_CORR; //time correction (default value for check_grad_vel)

  /***** erosion *****/
  trans_ref(8,  EST_OUEST, EAUC, GR,  EAUC, GRJ, lambda_E );
  trans_ref(7,  VERTICAL,  EAUC, GR,  GRJ,  EAUC, lambda_E/100000.0 );
  trans_type(7, TR_TRANSPORT);
#ifdef CELL_COLOR
  //erosion of colored grains
  if (lambda_E_col >= 0){
    tmode = TIME_NO_EVOL; // time correction not possible here !
    trans_check_color(8, check_color, 2, tmode, lambda_E_col);
    trans_check_color(7, check_color, 2, tmode, lambda_E_col/100000.0);
  }
#endif

#ifdef CGV
  //control of the erosion by shear stress
  Callback_check *pchk = check_grad_vel;
#ifdef CELL_COLOR
  extern float grdvc_max_col;
  if (grdvc_max_col >= 0){
    LogPrintf("check_grad_vel_color() attached to the erosion transitions\n");
    pchk = check_grad_vel_color;
  }
#endif //CELL_COLOR
  trans_check(8, pchk, 2, tmode);
  trans_check(7, pchk, 2, tmode);
#endif //CGV

  /***** gravity *****/
  trans(  VERTICAL,  GRJ, EAUC, EAUC, GRJ, lambda_G );

  /***** transport *****/
  if (lambda_T<0) lambda_T = 3/*6*/*lambda_C; //default value of lambda_T (with correlations of transitions)
  trans_ref(10,  VERTICAL,  EAUC, GRJ,  GRJ,  EAUC, 0 );
  //trans_check(10, check_slope, 2, TIME_NO_EVOL);
  //trans_ref(11,  EST_OUEST,  GRJ,  EAUC,  EAUC, GRJ, 6/*6*/*lambda_C );
  trans_ref(11,  EST_OUEST,  GRJ,  EAUC,  EAUC, GRJ, lambda_T );
#ifdef CELL_COLOR
  //transport of colored grains
  if (lambda_T_col >= 0){
    trans_check_color(11, check_color, 1, TIME_NO_EVOL, lambda_T_col);
  }
#endif

  /***** coupled transitions for saltation *****/
  trans_link(8, 10, 2, 0.5/*0.5*/);
  trans_link(10, 11, 1, 0.5/*0.5*/);
  trans_link(7, 11, 1, 1.0);

  /***** deposition *****/
  //trans_ref(9,  EST_OUEST, GRJ, EAUC, EAUC, GR, 2*lambda_C );
  //TODO: remove coefficient 2 -> implies to update all parameter files (we must warn users !)
  trans_ref(9,  EST_OUEST, GRJ, EAUC, GR, EAUC, lambda_C/*2*lambda_C*/ );
  trans_ref(16,  EST_OUEST, GRJ, GR, GR, GR, 10*lambda_C/*20*lambda_C*/ );
  trans_ref(17,  EST_OUEST, GRJ, DUM, GR, DUM, 10*lambda_C/*20*lambda_C*/ );
  trans(  EST_OUEST, GR, GRJ, GR, GR, 1000 );
#ifdef CELL_COLOR
  //deposition of colored grains
  if (lambda_C_col >= 0){
    trans_check_color(9, check_color, 1, TIME_NO_EVOL, lambda_C_col);
    trans_check_color(16, check_color, 1, TIME_NO_EVOL, 10*lambda_C_col);
    trans_check_color(17, check_color, 1, TIME_NO_EVOL, 10*lambda_C_col);
  }
#endif

  /***** gravity *****/
  trans(VERTICAL, GR, EAUC, EAUC,  GR, lambda_G);

  /***** diffusion *****/
  if (LNS > 1){
    trans_ref(12, NORD_SUD,   GR, EAUC, EAUC,  GR, lambda_D );
    trans_ref(13, NORD_SUD,   EAUC,  GR, GR, EAUC, lambda_D );
#ifdef CELL_COLOR
  //diffusion of colored grains
  if (lambda_D_col >= 0){
    trans_check_color(12, check_color, 1, tmode, lambda_D_col);
    trans_check_color(13, check_color, 2, tmode, lambda_D_col);
  }
#endif
#ifdef CGV
    //control of the diffusion by shear stress
    trans_check(12, check_grad_vel, 1, tmode);
    trans_check(13, check_grad_vel, 2, tmode);
#endif //CGV
  }

  /***** injection *****/
  //injection of sand grains (for the computation of Qsat and Lsat)
  //trans(VERTICAL, IN,  EAUC, IN, GRJ, lambda_I/*0.05*/);

  /***** output *****/
  //output of sand
  if (Ncel[OUT]>0)
    trans(VERTICAL, GR, OUT, EAUC, OUT, 1000);
  //trans(VERTICAL, GRJ, OUT, EAUC, OUT, 1000);

#ifdef AVALANCHES
  /***** avalanches *****/
  if (ava_trans){
    /*if (LNS > 1){
      trans_ref(14, HORIZONTAL, GR, EAUC, EAUC, GR, lambda_A);
      trans_ref(15, HORIZONTAL, EAUC, GR, GR, EAUC, lambda_A);
    }
    else{
      trans_ref(14, EST_OUEST, GR, EAUC, EAUC, GR, lambda_A);
      trans_ref(15, EST_OUEST, EAUC, GR, GR, EAUC, lambda_A);
    }*/
    trans_ref(14, HORIZONTAL, GR, EAUC, EAUC, GR, lambda_A);
    trans_ref(15, HORIZONTAL, EAUC, GR, GR, EAUC, lambda_A);
    trans_check(14, check_ava, 1, TIME_NO_EVOL);
    trans_check(15, check_ava, 2, TIME_NO_EVOL);
    if (ava_delay){
      trans_regul(14, callback_ava_trans);
      trans_regul(15, callback_ava_trans);
    }
  }
#endif
//#ifdef TRANSITION_BORD
  if (!pbc_mode){
    /***** boundary transitions *****/
    LogPrintf("transitions sur le bord\n");

    //trans(EST_OUEST,  BORD, GR,  BORD, GRJ, 0.01);
    //trans(EST_OUEST,  BORD, EAUC,  BORD, GRJ, 0.0002);
    trans(EST_OUEST,  GRJ,  BORD, EAUC, BORD, 20);
    trans(EST_OUEST,  GR,  BORD, EAUC, BORD, 20);

    /***** reinjection *****/
    if ((boundary==BC_REINJECTION) && (Ncel[IN]>0)){
      //trans_ref(2, EST_OUEST, IN,  EAUC, IN, GR, 0); // value 2 sets the (arbitrary) reference number of the transition
      trans_ref(2, VERTICAL, IN,  EAUC, IN, GR, lambda_I/*0.5*/); // value 2 sets the (arbitrary) reference number of the transition

      //reinjection of sand grains
      LogPrintf("reinjection de grains\n");
      trans_regul(2, callback_bord_dun);

      //temporal regulation of the reinjection rate
      //trans_regul(2, callback_injection_coef);
    }
  }
//#endif //TRANSITION_BORD

#endif //MODELE_DUN


#ifdef MODELE_AVA
  trans(VERTICAL, GR, AIR, AIR,  GR, 1000);
// source de grains
  trans(VERTICAL, IN, AIR, IN, GR, lambda_I);
// sortie de grains
  //trans(VERTICAL, GR, OUT, AIR, OUT, 1000);
  trans(VERTICAL, GR, OUT, AIR, OUT, 0.1);

  if (ava_trans){
// avalanches
#if 1 //1 pour conserver le couplage des transitions
    trans_ref(0, HORIZONTAL, GR, AIR, AIR, GR, lambda_A);
    trans_ref(1, HORIZONTAL, AIR, GR, GR, AIR,  lambda_A);
    trans_check(0, check_ava, 1, TIME_NO_EVOL);
    trans_check(1, check_ava, 2, TIME_NO_EVOL);
    if (ava_delay){
      trans_regul(0, callback_ava_trans);
      trans_regul(1, callback_ava_trans);
    }
#else //decouplage
    trans_ref(0, EST_OUEST, GR, AIR, AIR, GR, lambda_A);
    trans_ref(1, EST_OUEST, AIR, GR, GR, AIR,  lambda_A);
    trans_ref(2, NORD_SUD, GR, AIR, AIR, GR, lambda_A);
    trans_ref(3, NORD_SUD, AIR, GR, GR, AIR,  lambda_A);
    trans_check(0, check_ava, 1, TIME_NO_EVOL);
    trans_check(1, check_ava, 2, TIME_NO_EVOL);
    trans_check(2, check_ava, 1, TIME_NO_EVOL);
    trans_check(3, check_ava, 2, TIME_NO_EVOL);
#endif
  }
#endif //MODELE_AVA


#ifdef MODELE_CMB
#ifdef LGCA
// influence sur les transitions du flux genere par le module collisions
// si coef > 1, action positive (transistion favorisee)
// si coef < 1, action negative (transistion inhibee)

float coef_D = 30.0;    //3.0; //10     //diffusion
float coef_L = 10.0; //10.0; //2        //dissolution
float coef_C = 0.1; //0.1; //0.5        //cristallisation

// calcul de l'asservissement
#define CalcProbAsv(C) (C > 1) ? 1 - (1/C) : 1 - C
  prob_asv_D = CalcProbAsv(coef_D);
  prob_asv_L = CalcProbAsv(coef_L);
  prob_asv_C = CalcProbAsv(coef_C);

  if (coef_D > 1.0) lambda_D *= coef_D;
  if (coef_L > 1.0) lambda_L *= coef_L;
  if (coef_C > 1.0) lambda_C *= coef_C;

  LogPrintf("coef_D = %f\tlambda_D = %f\tprob_asv_D = %f\n", coef_D, lambda_D, prob_asv_D);
  LogPrintf("coef_L = %f\tlambda_L = %f\tprob_asv_L = %f\n", coef_L, lambda_L, prob_asv_L);
  LogPrintf("coef_C = %f\tlambda_C = %f\tprob_asv_C = %f\n", coef_C, lambda_C, prob_asv_C);

#ifdef SC_L
  if (coef_L != 1.0) LogPrintf("condition simple sur la dissolution\n");
#endif

#ifdef SC_C
  if (coef_C != 1.0) LogPrintf("condition simple sur la cristallisation\n");
#endif


#endif //LGCA

  //dissolution
  trans(  ISOTROPE,  PLUS, MOINS,  ZERO,  ZERO, lambda_L );
  trans(  ISOTROPE,  PLUS, MOINS,  PLUS,  ZERO, lambda_L );

  //cristallisation
  trans(  ISOTROPE,  PLUS,  ZERO,  PLUS,  PLUS, lambda_C );
  trans(  ISOTROPE,  PLUS,  ZERO,  PLUS, MOINS, lambda_C*2 );

  //diffusion
  trans(  VERTICAL, MOINS,  ZERO,  ZERO, MOINS, lambda_D*gravite);
  trans(  VERTICAL,  ZERO, MOINS, MOINS,  ZERO, lambda_D/gravite);
  //trans(  VERTICAL, MOINS,  ZERO,  ZERO, MOINS, (lambda_D + flux_n)*gravite); //ortho
  //trans(  VERTICAL,  ZERO, MOINS, MOINS,  ZERO, (lambda_D + flux_s)/gravite); //ortho
  if (flux_eo || flux_ns){
    trans(EST_OUEST, MOINS,  ZERO,  ZERO, MOINS, lambda_D+flux_o );
    trans(NORD_SUD, MOINS,  ZERO,  ZERO, MOINS, lambda_D+flux_n );
  }
  else{
    trans(HORIZONTAL, MOINS,  ZERO,  ZERO, MOINS, lambda_D );
  }
  if (flux_eo || flux_ns){
    trans(EST_OUEST,  ZERO, MOINS, MOINS,  ZERO, lambda_D+flux_e );
    trans(NORD_SUD,  ZERO, MOINS, MOINS,  ZERO, lambda_D+flux_s );
  }
  else{
    trans(HORIZONTAL,  ZERO, MOINS, MOINS,  ZERO, lambda_D );
  }
/*
  //etat PIERRE
  trans(  ISOTROPE,PIERRE,  ZERO,PIERRE,  PLUS, lambda_C );
  trans(  ISOTROPE,PIERRE,  ZERO,PIERRE, MOINS, lambda_C*2 );

  //PIERRE = source de ZERO
  trans(  ISOTROPE,PIERRE, MOINS,PIERRE,  ZERO, 0.001 );
*/
  //porosite PLUS/ZERO
  //trans(  ISOTROPE,  ZERO,  PLUS,  PLUS,  ZERO, 0.001 );

#ifdef TRANSITION_BORD
  LogPrintf("transitions sur le bord\n");
  //trans(  VERTICAL,  BORD,  ZERO,  BORD, MOINS, lambda_D*gravite );

  //modelisation fracture
  //trans(  VERTICAL,  BORD,  MOINS,  BORD, ZERO, 0.01 );
  //trans(  VERTICAL,  ZERO,  BORD,  MOINS, BORD, lambda_D/gravite );
  //trans(  VERTICAL,  BORD,  MOINS,  BORD, ZERO, 0.1/*lambda_D/gravite*/ );
  //trans(  VERTICAL,  ZERO,  BORD,  MOINS, BORD, 0.01 );

  trans(  ISOTROPE,  BORD,  MOINS, BORD,  ZERO, 0.1/*lambda_D/gravite*/ );

  //trans(  VERTICAL,  BORD,  MOINS, BORD,  ZERO, 0.1/*lambda_D/gravite*/ );
  //trans(  VERTICAL, MOINS,  BORD,  ZERO,  BORD, 0.1 );
  //trans(HORIZONTAL,  BORD,  MOINS, BORD,  ZERO, 0.1 );
  //trans(HORIZONTAL, MOINS,  BORD,  ZERO,  BORD, 0.1 );

  //porosite PLUS/ZERO aux bords
  //trans(  ISOTROPE,  BORD,  PLUS,  BORD,  ZERO, 0.0001 );
#endif

#endif //MODELE_CMB


#ifdef MODELE_ICB

  trans(  VERTICAL,  PLUS, MOINS,  ZERO,  ZERO, lambda_L );
  trans(  VERTICAL,  PLUS, MOINS,  PLUS,  ZERO, lambda_L*2 );
  trans(  VERTICAL,  PLUS,  ZERO,  PLUS,  PLUS, lambda_C );
  trans(  VERTICAL,  PLUS,  ZERO,  PLUS, MOINS, lambda_C*3 );
  trans(  VERTICAL, MOINS,  PLUS,  ZERO,  ZERO, lambda_L );
  trans(  VERTICAL, MOINS,  PLUS,  ZERO,  PLUS, lambda_L*2 );
  trans(  VERTICAL, MOINS,  ZERO,  ZERO, MOINS, lambda_D/gravite );
  trans(  VERTICAL,  ZERO,  PLUS,  PLUS,  PLUS, lambda_C );
  trans(  VERTICAL,  ZERO,  PLUS, MOINS,  PLUS, lambda_C*3 );
  trans(  VERTICAL,  ZERO, MOINS, MOINS,  ZERO, lambda_D*gravite );
  trans(HORIZONTAL,  PLUS, MOINS,  ZERO,  ZERO, lambda_L );
  trans(HORIZONTAL,  PLUS, MOINS,  PLUS,  ZERO, lambda_L*2 );
  trans(HORIZONTAL,  PLUS,  ZERO,  PLUS,  PLUS, lambda_C );
  trans(HORIZONTAL,  PLUS,  ZERO,  PLUS, MOINS, lambda_C*3 );
  trans(HORIZONTAL, MOINS,  PLUS,  ZERO,  ZERO, lambda_L );
  trans(HORIZONTAL, MOINS,  PLUS,  ZERO,  PLUS, lambda_L*2 );
  trans(HORIZONTAL, MOINS,  ZERO,  ZERO, MOINS, lambda_D+flux1 );
  trans(HORIZONTAL,  ZERO,  PLUS,  PLUS,  PLUS, lambda_C );
  trans(HORIZONTAL,  ZERO,  PLUS, MOINS,  PLUS, lambda_C*3 );
  trans(HORIZONTAL,  ZERO, MOINS, MOINS,  ZERO, lambda_D+flux2 );

  //LogPrintf("version avec refroidissement\n");
  LogPrintf("lambda_R = %f\n",lambda_R);
  trans_cel(MOINS, ZERO, lambda_R);

#endif

#ifdef MODELE_RIV

  double coef_v_L = 1;//0.01;
  double coef_h_L = 1;//100; //dans le sens du flux seulement
  double coef_h_C = 1;//0.01;
  double coef_h_C2 = 1;//0.1;
  double coef_C = 1;//0.01
  double coef_C_BT = 1;//10;//1
  double coef_D_BT = 10.0;

  double prob_v_BOUE_BT = 1.0;  //probabilite pour que BOUE devienne BT lors d'un transport vertical
  double prob_v_BOUE_BOUE = (1.0-prob_v_BOUE_BT);
  double prob_h_BT_BOUE = 0.1;  //probabilite pour que BT devienne BOUE lors d'un transport horizontal
  double prob_h_BT_BT = (1.0-prob_h_BT_BOUE);

  LogPrintf("coef_v_L = %f\n",coef_v_L);
  LogPrintf("coef_h_L = %f\n",coef_h_L);
  LogPrintf("coef_h_C = %f\n",coef_h_C);

  //diffusion BOUE/EAU
#ifndef LISSAGE
  trans(  VERTICAL, EAU,  BOUE,  BOUE, EAU, lambda_D/gravite );
#endif
#ifndef ETAT_BT
  trans_ref(0,  VERTICAL,  BOUE, EAU, EAU,  BOUE, lambda_D*gravite );
#else
  trans_ref(0,  VERTICAL,  BOUE, EAU, EAU,  BOUE, lambda_D*gravite*prob_v_BOUE_BOUE );
  trans_ref(0,  VERTICAL,  BOUE, EAU, EAU,  BT, lambda_D*gravite*prob_v_BOUE_BT );
#endif
  //trans_ref(1, HORIZONTAL, EAU,  BOUE,  BOUE, EAU, lambda_D+flux1 );
  //trans_ref(2, HORIZONTAL,  BOUE, EAU, EAU,  BOUE, lambda_D+flux2 );
  if (flux_eo == flux_ns){
    trans_ref(1, HORIZONTAL, EAU,  BOUE,  BOUE, EAU, lambda_D+flux_o );
    trans_ref(2, HORIZONTAL,  BOUE, EAU, EAU,  BOUE, lambda_D+flux_e );
  }
  else{
    trans_ref(1, EST_OUEST, EAU,  BOUE,  BOUE, EAU, lambda_D+flux_o );
    trans_ref(2, EST_OUEST,  BOUE, EAU, EAU,  BOUE, lambda_D+flux_e );
    trans(NORD_SUD, EAU,  BOUE,  BOUE, EAU, lambda_D+flux_n );
    trans(NORD_SUD,  BOUE, EAU, EAU,  BOUE, lambda_D+flux_s );
  }
  //dissolution
  trans_ref(3,  VERTICAL,  TERRE, EAU,  BOUE,  BOUE, lambda_L*coef_v_L );
  trans_ref(4,  VERTICAL,  TERRE, EAU,  TERRE,  BOUE, lambda_L*coef_v_L );
  trans_ref(5,  VERTICAL, EAU,  TERRE,  BOUE,  BOUE, lambda_L*coef_v_L );
  trans_ref(6,  VERTICAL, EAU,  TERRE,  BOUE,  TERRE, lambda_L*coef_v_L );
  trans_ref(7, HORIZONTAL,  TERRE, EAU,  BOUE,  BOUE, lambda_L );
  trans_ref(8, HORIZONTAL,  TERRE, EAU,  TERRE,  BOUE, lambda_L );
  trans_ref(9, HORIZONTAL, EAU,  TERRE,  BOUE,  BOUE, lambda_L*coef_h_L );
  trans_ref(10, HORIZONTAL, EAU,  TERRE,  BOUE,  TERRE, lambda_L*coef_h_L );
  //cristallisation
  trans(  VERTICAL,  TERRE,  BOUE,  TERRE,  TERRE, lambda_C*coef_C);
  trans(  VERTICAL,  BOUE,  TERRE,  TERRE,  TERRE, lambda_C*coef_C);
#ifndef LISSAGE
  trans(  VERTICAL,  TERRE,  BOUE,  TERRE, EAU, lambda_C*2*coef_C);
  trans(  VERTICAL,  BOUE,  TERRE, EAU,  TERRE, lambda_C*2*coef_C);
#endif
  trans(HORIZONTAL,  TERRE,  BOUE,  TERRE,  TERRE, lambda_C*coef_h_C*coef_C);
  trans(HORIZONTAL,  BOUE,  TERRE,  TERRE,  TERRE, lambda_C*coef_h_C*coef_C);
#ifndef LISSAGE
  trans(HORIZONTAL,  TERRE,  BOUE,  TERRE, EAU, lambda_C*2*coef_h_C*coef_C);
  trans(HORIZONTAL,  BOUE,  TERRE, EAU,  TERRE, lambda_C*2*coef_h_C*coef_C);
#endif
/*
  trans(  VERTICAL,  BOUE,  BOUE,  BOUE,  TERRE, lambda_C ); //test cristallisation BOUE+BOUE
  trans(  VERTICAL,  BOUE,  BOUE,  EAU,  BOUE, lambda_C*2 ); //test cristallisation BOUE+BOUE
  trans(HORIZONTAL,  BOUE,  BOUE,  TERRE,  BOUE, lambda_C*coef_h_C2 ); //test cristallisation BOUE+BOUE
  trans(HORIZONTAL,  BOUE,  BOUE,  BOUE,  TERRE, lambda_C*coef_h_C2 ); //test cristallisation BOUE+BOUE
  trans(HORIZONTAL,  BOUE,  BOUE,  BOUE,  EAU, lambda_C*2*coef_h_C2 ); //test cristallisation BOUE+BOUE
  trans(HORIZONTAL,  BOUE,  BOUE,  EAU,  BOUE, lambda_C*2*coef_h_C2 ); //test cristallisation BOUE+BOUE
*/

#ifdef ETAT_BT
  trans(  VERTICAL, EAU,  BT,  BT, EAU, lambda_D/gravite );
  trans_ref(20,  VERTICAL,  BT, EAU, EAU,  BT, lambda_D*gravite );
  //trans_ref(1, HORIZONTAL, EAU,  BOUE,  BOUE, EAU, lambda_D+flux1 );
  //trans_ref(2, HORIZONTAL,  BOUE, EAU, EAU,  BOUE, lambda_D+flux2 );
  if (flux_eo == flux_ns){
    trans_ref(21, HORIZONTAL, EAU,  BT,  BT, EAU, (lambda_D+flux_o)*coef_D_BT*prob_h_BT_BT );
    trans_ref(22, HORIZONTAL,  BT, EAU, EAU,  BT, (lambda_D+flux_e)*coef_D_BT*prob_h_BT_BT );
    trans_ref(23, HORIZONTAL, EAU,  BT,  BOUE, EAU, (lambda_D+flux_o)*coef_D_BT*prob_h_BT_BOUE );
    trans_ref(24, HORIZONTAL,  BT, EAU, EAU,  BOUE, (lambda_D+flux_e)*coef_D_BT*prob_h_BT_BOUE );
  }
  else{
    trans_ref(21, EST_OUEST, EAU,  BT,  BT, EAU, (lambda_D+flux_o)*coef_D_BT*prob_h_BT_BT );
    trans_ref(22, EST_OUEST,  BT, EAU, EAU,  BT, (lambda_D+flux_e)*coef_D_BT*prob_h_BT_BT );
    trans_ref(23, EST_OUEST, EAU,  BT,  BOUE, EAU, (lambda_D+flux_o)*coef_D_BT*prob_h_BT_BOUE );
    trans_ref(24, EST_OUEST,  BT, EAU, EAU,  BOUE, (lambda_D+flux_e)*coef_D_BT*prob_h_BT_BOUE );
    trans(NORD_SUD, EAU,  BT,  BT, EAU, (lambda_D+flux_n)*coef_D_BT*prob_h_BT_BT );
    trans(NORD_SUD,  BT, EAU, EAU,  BT, (lambda_D+flux_s)*coef_D_BT*prob_h_BT_BT );
    trans(NORD_SUD, EAU,  BT,  BOUE, EAU, (lambda_D+flux_n)*coef_D_BT*prob_h_BT_BOUE );
    trans(NORD_SUD,  BT, EAU, EAU,  BOUE, (lambda_D+flux_s)*coef_D_BT*prob_h_BT_BOUE );
}

  //cristallisation
  trans(  VERTICAL,  TERRE,  BT,  TERRE,  TERRE, lambda_C*coef_C_BT);
  trans(  VERTICAL,  BT,  TERRE,  TERRE,  TERRE, lambda_C*coef_C_BT);
  trans(  VERTICAL,  TERRE,  BT,  TERRE, EAU, lambda_C*2*coef_C_BT);
  trans(  VERTICAL,  BT,  TERRE, EAU,  TERRE, lambda_C*2*coef_C_BT);
  trans(  VERTICAL,  BT,  BT,  BT,  TERRE, lambda_C*coef_C_BT*10 ); //test cristallisation BT+BT
  trans(  VERTICAL,  BT,  BT,  EAU,  BT, lambda_C*2*coef_C_BT*10 ); //test cristallisation BT+BT

  //erosion-deposition (test ... weird !)
  //trans(VERTICAL,  BT, TERRE,  BOUE,  BOUE, lambda_D*lambda_E_v*2 );
  //trans(VERTICAL,  BT, TERRE,  TERRE,  TERRE, lambda_D*lambda_E_v );
#endif

#ifdef TRANSITION_BORD
  //trans(HORIZONTAL,  BORD,  BOUE, BORD,  EAU, (!flux2)?lambda_D:0 );
  //trans(HORIZONTAL,  BOUE,  BORD, EAU,  BORD, (!flux2)?lambda_D:flux2 /*/2*/ );
  if (flux_eo == flux_ns){
    trans(HORIZONTAL,  BORD,  BOUE, BORD,  EAU, (!flux_o)?lambda_D*gravite:flux_o );
    trans(HORIZONTAL,  BOUE,  BORD, EAU,  BORD, (!flux_e)?lambda_D*gravite:flux_e /*/2*/ );
  }
  else{
    trans(EST_OUEST,  BORD,  BOUE, BORD,  EAU, (!flux_o)?lambda_D:flux_o );
    trans(EST_OUEST,  BOUE,  BORD, EAU,  BORD, (!flux_e)?lambda_D:flux_e /*/2*/ );
    trans(NORD_SUD,  BORD,  BOUE, BORD,  EAU, (!flux_n)?lambda_D:flux_n );
    trans(NORD_SUD,  BOUE,  BORD, EAU,  BORD, (!flux_s)?lambda_D:flux_s /*/2*/ );
  }
  //trans(VERTICAL,  BOUE,  BORD, EAU,  BORD, lambda_D*gravite );
  //trans(VERTICAL,  BOUE,  BORD, DUM,  BORD, lambda_D*gravite );
  //trans(VERTICAL,  BOUE,  BORD, TERRE,  BORD, lambda_C );

  trans(HORIZONTAL, TERRE, BORD, EAU, BORD, 1); //avalanches aux bords
  trans(HORIZONTAL, BORD, TERRE, BORD, EAU, 1);
#endif
  //injection IN->BOUE
  trans(  VERTICAL,  IN,  EAU,  IN, BOUE, debit );
  trans(  VERTICAL,  IN,  TERRE,  IN, BOUE, debit );  //pour ne pas boucher la source

  //gravite TERRE
  trans(  VERTICAL, TERRE, EAU, EAU, TERRE, lambda_D*gravite/*0.002*/);
  trans(  VERTICAL, TERRE, BOUE, BOUE, TERRE, lambda_D*gravite/*0.002*/);
#ifdef ETAT_BT
  trans(  VERTICAL, TERRE, BT, BT, TERRE, lambda_D*gravite/*0.002*/);
#endif
  //dilution BOUE/EAU
  //trans(  VERTICAL,  EAU,  BOUE, EAU,  EAU, 0.01 );

//#ifndef ETAT_BT
  //erosion par transport horizontal, correlations verticales
  if (flux_e == 0){
    trans_link(1, 5, 2, lambda_E_v);
    trans_link(1, 6, 2, lambda_E_v);
  }
  if (flux_o == 0){
    trans_link(2, 5, 1, lambda_E_v);
    trans_link(2, 6, 1, lambda_E_v);
  }

  //erosion par transport horizontal, correlations horizontales
  if (flux_e == 0){
    trans_link(1, 7, 2, lambda_E_h);
    trans_link(1, 8, 2, lambda_E_h);
    trans_link(1, 9, 2, lambda_E_h);
    trans_link(1, 10, 2, lambda_E_h);
  }
  if (flux_o == 0){
    trans_link(2, 7, 1, lambda_E_h);
    trans_link(2, 8, 1, lambda_E_h);
    trans_link(2, 9, 1, lambda_E_h);
    trans_link(2, 10, 1, lambda_E_h);
  }

  //erosion par transport vertical, correlations horizontales
  /*trans_link(0, 7, lambda_E_h);
  trans_link(0, 8, lambda_E_h);
  trans_link(0, 9, lambda_E_h);
  trans_link(0, 10, lambda_E_h);*/

//#else
#ifdef ETAT_BT

  //erosion par transport horizontal de cellules BT, correlations verticales
  //trans_link(21, 5, lambda_E_v*5);
  //trans_link(21, 6, lambda_E_v*5);
  //trans_link(22, 5, lambda_E_v*5);
  //trans_link(22, 6, lambda_E_v*5);
  //trans_link(23, 5, lambda_E_v*5);
  //trans_link(23, 6, lambda_E_v*5);
  //trans_link(24, 5, lambda_E_v*5);
  //trans_link(24, 6, lambda_E_v*5);

  //trans_link(23, 5, lambda_E_v);
  //trans_link(23, 6, lambda_E_v);
  //trans_link(24, 5, lambda_E_v);
  //trans_link(24, 6, lambda_E_v);

  trans_link(23, 5, 2, 0.5);
  trans_link(23, 6, 2, 0.5);
  trans_link(24, 5, 1, 0.5);
  trans_link(24, 6, 1, 0.5);

  //erosion par transport horizontal de cellules BT, correlations hoprizontales
  /*trans_link(23, 7, 0.05);
  trans_link(23, 8, 0.05);
  trans_link(23, 9, 0.05);
  trans_link(23, 10, 0.05);
  trans_link(24, 7, 0.05);
  trans_link(24, 8, 0.05);
  trans_link(24, 9, 0.05);
  trans_link(24, 10, 0.05);*/
#endif

  //erosion par transport vertical, correlations horizontales
  /*trans_link(20, 7, lambda_E_h);
  trans_link(20, 8, lambda_E_h);
  trans_link(20, 9, lambda_E_h);
  trans_link(20, 10, lambda_E_h);*/ //test
//#endif

  //gravite PIERRE/*
  trans(  VERTICAL, PIERRE, EAU, EAU, PIERRE, 10);
  trans(  VERTICAL, PIERRE, BOUE, BOUE, PIERRE, 1);
#endif

#ifdef MODELE_CHE

  //diffusion BOUE/EAU
  trans(  VERTICAL, EAU,  BOUE,  BOUE, EAU, lambda_D/gravite /* /1.5 */ );
  trans(  VERTICAL,  BOUE, EAU, EAU,  BOUE, lambda_D*gravite /* *1.5 */ );
  trans(HORIZONTAL, EAU,  BOUE,  BOUE, EAU, lambda_D+flux1 );
  trans(HORIZONTAL,  BOUE, EAU, EAU,  BOUE, lambda_D+flux2 /*+tr_h/50*/ );
  //dissolution
  trans(  VERTICAL,  TERRE, EAU,  BOUE,  BOUE, lambda_L );
  trans(  VERTICAL,  TERRE, EAU,  TERRE,  BOUE, lambda_L );
  trans(  VERTICAL, EAU,  TERRE,  BOUE,  BOUE, lambda_L );
  trans(  VERTICAL, EAU,  TERRE,  BOUE,  TERRE, lambda_L );
  trans(HORIZONTAL,  TERRE, EAU,  BOUE,  BOUE, lambda_L );
  trans(HORIZONTAL,  TERRE, EAU,  TERRE,  BOUE, lambda_L );
  trans(HORIZONTAL, EAU,  TERRE,  BOUE,  BOUE, lambda_L );
  trans(HORIZONTAL, EAU,  TERRE,  BOUE,  TERRE, lambda_L );
  //cristallisation
  trans(  VERTICAL,  TERRE,  BOUE,  TERRE,  TERRE, lambda_C );
  trans(  VERTICAL,  TERRE,  BOUE,  TERRE, EAU, lambda_C*2 );
  trans(  VERTICAL,  BOUE,  TERRE,  TERRE,  TERRE, lambda_C );
  trans(  VERTICAL,  BOUE,  TERRE, EAU,  TERRE, lambda_C*2 );
  trans(HORIZONTAL,  TERRE,  BOUE,  TERRE,  TERRE, lambda_C );
  trans(HORIZONTAL,  TERRE,  BOUE,  TERRE, EAU, lambda_C*2 );
  trans(HORIZONTAL,  BOUE,  TERRE,  TERRE,  TERRE, lambda_C );
  trans(HORIZONTAL,  BOUE,  TERRE, EAU,  TERRE, lambda_C*2 );

#ifdef DIFFUSION_TERRE_EAU
  //diffusion TERRE/EAU
  trans(  VERTICAL,  TERRE, EAU, EAU,  TERRE, lambda_D*2 );
  trans(  VERTICAL, EAU,  TERRE,  TERRE, EAU, lambda_D/2 );
  trans(HORIZONTAL, EAU,  TERRE,  TERRE, EAU, lambda_D+flux1);
  trans(HORIZONTAL,  TERRE, EAU, EAU,  TERRE, lambda_D+flux2 );
#endif

  //gravite TERRE/AIR
  trans(  VERTICAL, TERRE, AIR, AIR, TERRE, tr_v);

  //diffusion EAU/AIR
  trans(  VERTICAL, EAU,   AIR,   AIR, EAU, tr_v );
  trans(HORIZONTAL, EAU,   AIR,   AIR, EAU, diff_h+tr_h );
  trans(HORIZONTAL,   AIR, EAU, EAU,   AIR, diff_h );

  //diffusion BOUE/AIR
  trans(  VERTICAL,  BOUE,   AIR,   AIR,  BOUE, tr_v );
  trans(HORIZONTAL,  BOUE,   AIR,   AIR,  BOUE, (diff_h+tr_h)/5 );
  trans(HORIZONTAL,   AIR,  BOUE,  BOUE,   AIR, diff_h/5 );

  //injection IN->EAU
  trans(  VERTICAL,  IN,   AIR,  IN, EAU, debit /* tr_v/50 */ );/* /100000 pour la pluie*/
#ifdef TRANSITION_BORD
  //echappement EAU+BOUE+TERRE
  trans(HORIZONTAL,  EAU,  BORD,   AIR,  BORD, (diff_h+tr_h)/5 );
  trans(HORIZONTAL,  BOUE, BORD,   AIR,  BORD, (diff_h+tr_h)/10 );
#ifdef DIFFUSION_TERRE_EAU
  trans(HORIZONTAL,  TERRE,  BORD,   EAU,  BORD, (lambda_D+flux2)/2 );
#endif
#endif

  //gravite PIERRE/*
  //trans(  VERTICAL, PIERRE, AIR, AIR, PIERRE, tr_v);
  //trans(  VERTICAL, PIERRE, EAU, EAU, PIERRE, tr_v/2);
  //trans(  VERTICAL, PIERRE, BOUE, BOUE, PIERRE, tr_v/4);
#endif

#ifdef MODELE_CRI

  //cristallisation
  trans(VERTICAL, PLUS, ZERO, PLUS, PLUS, lambda_C/100000);
  trans(VERTICAL, ZERO, PLUS, PLUS, PLUS, lambda_C/100000);
  trans(HORIZONTAL, PLUS, ZERO, PLUS, PLUS, lambda_C);
  trans(HORIZONTAL, ZERO, PLUS, PLUS, PLUS, lambda_C);
  //diffusion ZERO/AIR
  trans(VERTICAL, ZERO, AIR, AIR, ZERO, lambda_D);
  trans(VERTICAL, AIR, ZERO, ZERO, AIR, lambda_D);
  trans(HORIZONTAL, ZERO, AIR, AIR, ZERO, lambda_D+flux2);
  trans(HORIZONTAL, AIR, ZERO, ZERO, AIR, lambda_D);
  //injection ZERO
  trans(VERTICAL, IN, AIR, IN, ZERO, lambda_D/100);
#ifdef TRANSITION_BORD
  //echappement ZERO
  trans(HORIZONTAL, ZERO, BORD, AIR, BORD, lambda_D+flux2);
#endif

#endif

#ifdef MODELE_DIF

  //diffusion
  trans(VERTICAL, EAU, AIR, AIR, EAU, 1000);
  trans(HORIZONTAL, EAU, AIR, AIR, EAU, 1);
  trans(HORIZONTAL, AIR, EAU, EAU, AIR, 1);

  //infiltration
  //trans(VERTICAL, EAU, TERRE, EAU, BOUE, 1000);
  trans(VERTICAL, EAU, TERRE, AIR, BOUE, lambda_D);
  trans(VERTICAL, EAU, BOUE, AIR, BOUE, lambda_D);

  //source
  trans(VERTICAL, IN, AIR, IN, EAU, debit);

#ifdef TRANSITION_BORD
  trans(HORIZONTAL, EAU, BORD, AIR, BORD, 1);
  trans(HORIZONTAL, BORD, EAU, BORD, AIR, 1);
#endif
#endif

  //LogPrintf("gravite = %f\n",gravite);

#ifdef MODELE_LIFE

  //birth
  //trans_ref(0,  ISOTROPE,  ALIVE, DEAD,  ALIVE,  ALIVE, lambda_B );
  trans_ref(0,  HORIZONTAL,  ALIVE, DEAD,  ALIVE,  ALIVE, lambda_B );
  trans_ref(1,  HORIZONTAL,  DEAD,  ALIVE, ALIVE,  ALIVE, lambda_B );
  trans_check(0, check_birth, 2, TIME_EVOL);
  trans_check(1, check_birth, 1, TIME_EVOL);
  //death
  //trans_ref(1,  ISOTROPE,  ALIVE, DEAD,  DEAD,  DEAD, lambda_D );
  //trans_ref(2,  ISOTROPE,  ALIVE, ALIVE,  ALIVE,  DEAD, lambda_D );
  trans_ref(2,  HORIZONTAL,  ALIVE, DEAD,  DEAD,  DEAD, lambda_D );
  trans_ref(3,  HORIZONTAL,  DEAD,  ALIVE, DEAD,  DEAD, lambda_D );
  trans_ref(4,  HORIZONTAL,  ALIVE, ALIVE, ALIVE, DEAD, lambda_D );
  trans_ref(5,  HORIZONTAL,  ALIVE, ALIVE, DEAD,  ALIVE,lambda_D );
  trans_check(2, check_death, 1, TIME_EVOL);
  trans_check(3, check_death, 2, TIME_EVOL);
  trans_check(4, check_death, 2, TIME_EVOL);
  trans_check(5, check_death, 1, TIME_EVOL);

#endif
}

#ifdef MODELE_DUN

int callback_bord_dun(void *data)
{
  static int Ncel_GR_debut = 0;
  int Ncel_GR;
  double *ref_intensite = (double*)data;

  Ncel_GR = Ncel[GR] + Ncel[GRJ];

  if (!Ncel_GR_debut){
    Ncel_GR_debut = Ncel_GR;
    coef_injection = *ref_intensite;
    *ref_intensite = 0.0;
  }
  else{
    //si la quantite de grain diminue, on reinjecte de la matiere
    //*ref_intensite = (Ncel_GR < Ncel_GR_debut) ? 1.0 : 0.0;
    *ref_intensite = (Ncel_GR < Ncel_GR_debut) ? coef_injection : 0.0;
    //if (*ref_intensite) LogPrintf("intensite injection = %f\n",coef_injection);
  }

  return 0;
}

int callback_injection_coef(void *data)
{
  static char start = 1;
  static double lambda_injection_debut;

  double *ref_intensite = (double*)data;

  if (start){
    lambda_injection_debut = *ref_intensite;
    start = 0;
  }
  *ref_intensite = lambda_injection_debut*coef_injection;
  //LogPrintf("intensite injection = %f\n", *ref_intensite);

  return 0;
}

#endif


//modulation intensite des transitions d'avalanches
int callback_ava_trans(void *data)
{
  double *ref_intensite = (double*)data;

  *ref_intensite = lambda_A;

  return 0;
}


