/* ReSCAL - Stochastic simulation engine
 *
 * Copyright (C) 2011-2013
 *
 * Author: Olivier Rozier <rozier@ipgp.fr>
 *
 * Code based on dissol program,
 * by Eduardo Sepulveda <edo@espci.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.
 */



#include <stdio.h>
#include <stdlib.h>
#include <memory.h>
#include <math.h>

#include "defs.h"
#include "macros.h"
#include "param.h"
#include "format.h"
#include "cells.h"
#include "space.h"
#include "surface.h"
#include "doublets.h"
#ifdef PARALLEL
#include "synchro.h"
#endif
#include "transitions.h"
#include "simul.h"
#include "trace.h"
#ifdef LGCA
#include "lgca.h"
#endif
#include "callbacks.h"
#include "interface.h"
#include "view.h"

extern unsigned char opt_info, opt_nv;
extern int H, L, P, HL, HLP;       // les dimensions de la terre
extern int LN, LS, LNS, HLN;    //couloir est-ouest (limite nord, limite sud, largeur nord-sud, ...)
extern Cell  *TE;	           // la 'terre'
extern int *db_pos[];                   // les tableaux contenant la position des doublets actifs
extern Doublet t_doub[MAX_DB];               // la description des doublets
extern int Ncel[];                       // nombre de cellules par type
extern int Ndb[];                        // nombre doublets actifs par type
extern int nb_trans_db, nb_trans_cel, nb_trans; // nombre de transitions de 2 cellules, d'1 cellule et au total
extern TransitionDb t_trans[];            // la description des transitions de doublets
extern TransitionCel t_trans_cel[];           // la description des transitions de cellules
extern LienTransDb   t_lien[];        // la description des liens de correlation entre transitions de 2 cellules
extern int nb_liens;             // nombre de liens entre transitions
extern const char *classes_db[];        // noms des classes de doublet
extern const char *etats[];             // les noms des types de cellules
#ifdef INFO_TRANS
extern int           cpt_lien[];      // compteur des correlations effectives
#endif
extern int           direction[];             // conversion classe de transition -> direction
#ifdef LGCA
extern int col_iter;    //nombre de cycles de collisions
extern int CLEO;
extern char *nom_fic_mvt;
#ifdef CGV
extern float grdvc_max; //seuil max pour le controle de l'erosion en fonction du gradient de vitesse
extern float grdvc_min; //seuil min pour le controle de l'erosion en fonction du gradient de vitesse
#ifdef CELL_COLOR
extern float grdvc_max_col; //seuil max pour le controle de l'erosion des grains colores en fonction du gradient de vitesse
extern float grdvc_min_col; //seuil min pour le controle de l'erosion des grains colores en fonction du gradient de vitesse
#endif
#endif
#endif
extern float lambda_A;
extern float lambda_A_unstable;
#ifdef CELL_COLOR
extern float color_coef;         // coefficient associe a la callback check_color
#endif
#ifdef PARALLEL
extern int mode_par;    //mode parallele
#endif

#ifdef REFDB_PTR
extern RefDoublets *RefDB;       // references des cellules de la terre vers les doublets actifs
#else
extern RefDoublets_Type *RefDB_Type;       // references des cellules de la terre vers les doublets actifs
extern RefDoublets_Ind *RefDB_Ind;       // references des cellules de la terre vers les doublets actifs
#endif

#ifdef CGV
extern float cgv_coef;
#endif

extern int vdir_mode; //display mode of the current orientation
extern unsigned char reorient_flag;

unsigned long int iter=0, md_iter=0;  // nombre d'iterations (unites, milliards d'unites)
double csp_time=0.0;                  // time variable for the cellular space
//double temps=0.0;                  // temps reel simule
#ifdef TIME_SCALE
double real_time=0.0;       // real time with time scale
#endif
char end_of_simul=0;
int cpt_trans_blocked[MAX_TRANSITIONS_DB];  //nombre de transitions blockees
double PII, PII_time;     //probabilités totales
double PI_ij[MAX_TRANSITIONS];                  // les PI_ij
double cum[MAX_TRANSITIONS];     //probabilites cumulees
char *boundary_str = NULL;
#ifdef CYCLAGE_HOR
int boundary = BC_PERIODIC;  //boundary conditions
int pbc_mode = 1; //periodic boundary conditions
#else
int boundary = BC_REINJECTION;  //boundary conditions
int pbc_mode = 0; //periodic boundary conditions
#endif
int ava_mode = AVA_SYNC;  //mode d'avalanches
char *ava_mode_str = NULL;
int ava_trans = 0;        //transitions d'avalanches
int ava_norm = 0;       //avalanches avec normales
int ava_h_lim = 1; //10;	// hauteur limite avant avalanche
int ava_nb_cel_max = 1;  // nb de cellules qui tombent simultanement
float ava_delay = 0.0; //delay between avalanches
float ava_duration = 1.0; //default value for DUN model
//float ava_delay_long = 0.0; //8000.0;
float ava_angle = 0.0; //angle of avalanches (degrees)
float ava_angle_col = 0.0; //angle of avalanches (degrees) for colored grains
float ava_angle_stable = 0.0;
float ava_angle_unstable = 0.0;
float ava_angle_stable_col = 0.0;
float ava_angle_unstable_col = 0.0;
float lambda_A_stable = 0.0;
int ava_anisotropy = 0;
//#endif
unsigned char simul_dump_flag = 0;
unsigned char csphpp_flag = 0;
float dump_delay_png = 0.0;
float dump_delay_csp = 0.0;
float stop_delay_t0 = 0.0;

#ifdef SURRECTIONS
float sur_delay = 0.0;  //10000.0;              // delai entre deux surrections successives (0 = pas de surrection)
int sur_mode = 0;                // type de surrection (uniforme ou locale)
char *sur_mode_str=NULL;
#endif

int rot_mode = ROT_CYCL;  //default rotation mode
#ifdef LGCA
int rot_mvt_nb_cycles = -1;  //nombre de cycles de stabilisation apres rotation
#endif
#ifdef ROTATIONS
float rot_angle = 0.0;
float rot_delay = 0.0;
char *rot_mode_str=NULL;
#endif

#ifdef LGCA
int init_mvt_nb_cycles = -1;  //nombre de cycles de stabilisation du flux au depart
double lgca_delay = 0.0; //delai entre deux cycles de gaz sur reseau
int lgca_reset = 0;
int lgca_ready = 0;
#endif
//#ifdef WIND_DATA
char *wind_data_filename = NULL;
//#endif
#ifdef CGV
char *qsat_data_filename = NULL; //"real_data/PDF.data";
#endif
#ifdef TIME_SCALE
char *phys_prop_filename = NULL; //"real_data/desert_earth.prop";
float time_scale = 1.0; //t_0
#endif

int simul_check(int , int , int );
void simul_lien_trans(int , int , int );
#ifdef LGCA
void stabilisation_collisions(int );
#endif
void simul_dump();

void params_simul()
{
  //csp_time
  parametre("Time", "initial physical time - optional", &csp_time, PARAM_DOUBLE);

#ifdef CYCLAGE_HOR
#ifdef MODELE_DUN
  parametre("Boundary", "boundary conditions (PERIODIC|OPEN|REINJECTION), PERIODIC by default", &boundary_str, PARAM_STRING);
#else
  parametre("Boundary", "boundary conditions (PERIODIC|OPEN), PERIODIC by default", &boundary_str, PARAM_STRING);
#endif
#elif defined(MODELE_DUN)
  parametre("Boundary", "boundary conditions (OPEN|REINJECTION), REINJECTION by default", &boundary_str, PARAM_STRING);
#endif

#ifdef AVALANCHES
  parametre("Ava_mode", "mode of avalanches (NONE|SYNC|TRANS|PROPAG, SYNC by default)", &ava_mode_str, PARAM_STRING);
  parametre("Ava_h_lim", "height limit for avalanches (1 by default)", &ava_h_lim, PARAM_INT);
  parametre("Ava_nb_cel", "number of cells falling together in SYNC mode (1 by default) - optional", &ava_nb_cel_max, PARAM_INT);
  parametre("Ava_delay", "delay between avalanches - optional", &ava_delay, PARAM_FLOAT);
  parametre("Ava_angle", "angle of repose (degrees) - optional", &ava_angle, PARAM_FLOAT);
  parametre("Ava_duration", "duration of avalanches - optional", &ava_duration, PARAM_FLOAT);
  parametre("Ava_angle_unstable", "unstable angle of repose (degrees) - optional", &ava_angle_unstable, PARAM_FLOAT);
  parametre("Ava_anisotropy", "anisotropic conditions in TRANS mode (0|1) - optional", &ava_anisotropy, PARAM_INT);
#ifdef CELL_COLOR
  parametre("Ava_angle_col", "angle of repose for white grains (degrees) - optional", &ava_angle_col, PARAM_FLOAT);
  parametre("Ava_angle_unstable_col", "unstable angle of repose for white grains (degrees) - optional", &ava_angle_unstable_col, PARAM_FLOAT);
#endif
#endif

#ifdef CELL_COLOR
  parametre("Color_coef", "control coefficient of colored grains - optional", &color_coef, PARAM_FLOAT);
#endif

#ifdef SURRECTIONS
  parametre("Sur_mode", "surrection mode (UNIFORM|CONE|TECTO)", &sur_mode_str, PARAM_STRING);
  parametre("Sur_delay", "surrection delay", &sur_delay, PARAM_FLOAT);
#endif

#ifdef ROTATIONS
  parametre("Rot_mode", "rotation mode (DISK|CYCLE) - optional", &rot_mode_str, PARAM_STRING);
  parametre("Rot_angle", "rotation angle - optional", &rot_angle, PARAM_FLOAT);
  parametre("Rot_delay", "rotation delay - optional", &rot_delay, PARAM_FLOAT);
#ifdef LGCA
  parametre("Rot_ncycl", "number of flow cycles without transitions after rotation - optional", &rot_mvt_nb_cycles, PARAM_INT);
#endif
#endif

#ifdef LGCA
  parametre("Init_ncycl", "initial number of flow cycles without transitions - optional", &init_mvt_nb_cycles, PARAM_INT);
  parametre("Lgca_delay", "physical time delay between flow cycles - optional", &lgca_delay, PARAM_DOUBLE);
  parametre("Lgca_reset", "reset of the flow after rotation (0|1) - optional", &lgca_reset, PARAM_INT);
#endif

#ifdef CGV
//  parametre("Grdvc", "critical shear stress for relationship between erosion rate and shear stress - optional", &grdvc, PARAM_FLOAT);
  parametre("Tau_min", "min shear stress value in the linear relationship between erosion rate and shear stress - optional", &grdvc_min, PARAM_FLOAT);
  parametre("Tau_max", "max shear stress value in the linear relationship between erosion rate and shear stress - optional", &grdvc_max, PARAM_FLOAT);
  parametre_cache("Grdvc_min", "Grdvc_min parameter is obsolete - Please use Tau_min parameter instead !", &grdvc_min, PARAM_FLOAT);
  parametre_cache("Grdvc_max", "Grdvc_max parameter is obsolete - Please use Tau_max parameter instead !", &grdvc_max, PARAM_FLOAT);
#ifdef CELL_COLOR
  parametre("Tau_min_col", "min shear stress value for white grains in the linear relationship between erosion rate and shear stress - optional", &grdvc_min_col, PARAM_FLOAT);
  parametre("Tau_max_col", "max shear stress value for white grains in the linear relationship between erosion rate and shear stress - optional", &grdvc_max_col, PARAM_FLOAT);
  parametre_cache("Grdvc_min_col", "Grdvc_min_col parameter is obsolete - Please use Tau_min_col parameter instead !", &grdvc_min_col, PARAM_FLOAT);
  parametre_cache("Grdvc_max_col", "Grdvc_max_col parameter is obsolete - Please use Tau_max_col parameter instead !", &grdvc_max_col, PARAM_FLOAT);
#endif
  parametre("Qsat_file", "real PDF data between Qsat and tau1 - optional", &qsat_data_filename, PARAM_STRING);
#endif
#ifdef WIND_DATA
  parametre("Wind_file", "filename for wind data - optional", &wind_data_filename, PARAM_STRING);
#endif
#ifdef TIME_SCALE
  parametre("Phys_prop_file", "filename for physical parameters - optional", &phys_prop_filename, PARAM_STRING);
#endif
}

void simul_parse()
{
  if (boundary_str){
    if (!strcmp(boundary_str, "PERIODIC")){
#ifdef CYCLAGE_HOR
      boundary = BC_PERIODIC;
#else
      ErrPrintf("Incorrect value for Boundary : %s, please define CYCLAGE_HOR in defs.h and recompile\n", boundary_str);
      exit(-2);
#endif
    }
    else if (!strcmp(boundary_str, "OPEN")){
      boundary = BC_OPEN;
    }
    else if (!strcmp(boundary_str, "REINJECTION")){
#ifdef MODELE_DUN
      boundary = BC_REINJECTION;
#else
      ErrPrintf("Incorrect value for Boundary : %s, not compatible with model %s\n", boundary_str, MOD_NAME);
      exit(-2);
#endif
    }
    else{
      ErrPrintf("Incorrect value for Boundary : %s\n", boundary_str);
      exit(-2);
    }
    pbc_mode = (boundary == BC_PERIODIC);
#ifdef PARALLEL
    if (mode_par) pbc_mode = 0;
#endif
  }
  if (ava_mode_str){
    if (!strcmp(ava_mode_str, "NONE")){
      ava_mode = AVA_NONE;
    }
    else if (!strcmp(ava_mode_str, "SYNC")){
      ava_mode = AVA_SYNC;
    }
    else if (!strcmp(ava_mode_str, "TRANS")){
      ava_mode = AVA_TRANS;
    }
    else if (!strcmp(ava_mode_str, "PROPAG")){
      ava_mode = AVA_PROPAG;
    }
    else{
      ErrPrintf("Incorrect value for Ava_mode : %s\n", ava_mode_str);
      exit(-2);
    }
  }
#ifdef SURRECTIONS
  if (sur_mode_str){
    if (!strcmp(sur_mode_str, "UNIFORM")){
      sur_mode = SUR_MODE_UNIFORM;
    }
    else if (!strcmp(sur_mode_str, "CONE")){
      sur_mode = SUR_MODE_CONE;
    }
    else if (!strcmp(sur_mode_str, "TECTO")){
      sur_mode = SUR_MODE_TECTO;
    }
  }
#endif
#ifdef ROTATIONS
  if (rot_mode_str){
    if (!strcmp(rot_mode_str, "DISK")){
      rot_mode = ROT_DISK;
    }
    else if (!strcmp(rot_mode_str, "CYCLE")){
      rot_mode = ROT_CYCL;
    }
  }
#endif
}

void init_simul()
{
  //params_simul();
  //params_modele();
  simul_parse();
  ava_trans = (ava_mode == AVA_TRANS) || (ava_mode == AVA_PROPAG);
  ava_norm = (ava_angle != 0);
  simul_dump_flag = (dump_delay_png>0) || (dump_delay_csp>0);
}


// computation of the probability distribution
inline void compute_prob_dist()
{
  int i;
  double PI_ij_time, ccc;
  TransitionDb *ptr;
  TransitionCel *ptr_cel;

  //LogPrintf("calcul_prob: start\n");

  // doublet transitions
  ptr = t_trans;
  for(i=0; i<nb_trans_db; i++, ptr++){
    //regulation of a transition
    if (ptr->regul) ptr->regul(&ptr->intensite);
    PI_ij[i] = (double) Ndb[ptr->depart]*ptr->intensite;
#ifdef CGV
    //if (t_trans[i].check) PI_ij[i] *= cgv_coef;
#endif
    //if (PI_ij[i]) LogPrintf("tr = %d   rate = %f   time = %f\n", i, ptr->intensite, csp_time); //debug
  }

  // cell transitions
  ptr_cel = t_trans_cel;
  for(; i<nb_trans; i++,ptr_cel++){
    if (ptr_cel->regul) ptr_cel->regul(&ptr_cel->intensite); //regulation de la transition
    PI_ij[i] = (double) Ncel[ptr_cel->depart] * ptr_cel->intensite;
  }

  // calcul PII : overall probability for any transition to occur in a unit of time
  // calcul PII_time : the same only for the transitions with time evolution (and time correction if transitions are blocked)
  for(i=0, PII=0.0, ccc=0.0, PII_time = 0.0; i<nb_trans; i++){
    PII  += PI_ij[i];
    if (i<nb_trans_db){
      if (t_trans[i].time_mode == TIME_EVOL)
        PI_ij_time = PI_ij[i];
      else if (t_trans[i].time_mode == TIME_CORR){
#ifdef CGV
        PI_ij_time = PI_ij[i]*cgv_coef;
#else
        ErrPrintf("calcule_prob: no coefficient for time correction\n");
        exit(-1);
#endif
      }
      else
        PI_ij_time = 0;
    }
    else{
      PI_ij_time = PI_ij[i];
    }
/*#ifdef CGV
    //if ((i<nb_trans_db) && (t_trans[i].check != NULL)) PII_time  += PI_ij[i]*cgv_coef; else PII_time  += PI_ij[i];
    if ((i<nb_trans_db) && (t_trans[i].time_mode == TIME_CORR)) PII_time  += PI_ij[i]*cgv_coef; else PII_time  += PI_ij[i];
#endif*/
    cum[i] = PI_ij[i] + ccc;
    ccc    = cum[i];
    PII_time += PI_ij_time;
  }
/*#ifndef CGV
  PII_time = PII;
#endif*/

  //LogPrintf("PII = %f\n", PII);
  //return PII;
}


void simul_time()
{
  double aleat, t;

  //computation of csp and real time evolution
  //using Poisson formula
  if (PII_time){
    aleat = drand48();
    if (!aleat) {ErrPrintf("ERROR: simul_time - aleat = 0 !\n"); exit(-1); }
    t = -log(aleat)/PII_time;
    csp_time += t;
#ifdef TIME_SCALE
    real_time += t*time_scale;
#endif
  }
  else{
    LogPrintf("simul_time: PII_time = 0\n");
  }
}


inline int simul_trans()
{
  int i, elt, tr;
  double ccc, aleat;
  int dir;
  int cpt_blk=0;
  int max_blk=HLP;
  char flag_evol_time=0;
  char flag_trans_done=0;

  while ((!flag_trans_done) && (cpt_blk < max_blk)){

    //we choose randomly which transition will occur
    tr = 0;
    aleat = drand48();
    ccc = PII * aleat;
    while (ccc >= cum[tr]) tr++;

    //LogPrintf("tr=%d\n",tr);//debug

    if (tr > nb_trans){
      ErrPrintf("ERROR: array overflow in t_trans (%d)\n", tr);
      exit(-1);
    }

    if (tr < nb_trans_db){
      /// DOUBLET TRANSITION
      //we choose randomly which doublet will undergo the transition
      int db_depart = t_trans[tr].depart;
      elt = (int) floor( drand48() * Ndb[db_depart]);

      int ix = db_pos[db_depart][elt];
  /*#ifdef TRACE
      if (ii>10) trace_test(ix);
  #endif*/


      //direction (class) of the transition
      if (t_trans[tr].classe == HORIZONTAL){
  #ifdef REFDB_PTR
        dir = (RefDB[ix][EST] == &db_pos[db_depart][elt]) ? EST : SUD;
  #else
        dir = ((RefDB_Type[ix][EST] == db_depart) && (RefDB_Ind[ix][EST] == elt)) ? EST : SUD;
  #endif
      }
      else
        dir = direction[t_trans[tr].classe];

      //control of the transition
      if (t_trans[tr].nb_chk && !simul_check(tr, ix, dir)){
        //blocked transition
        cpt_blk++;
        continue;
      }

      do_trans_db(tr, ix, dir);
      flag_trans_done = 1;
      iter++; //LogPrintf("iter = %ld\n", iter);
      if (iter > 1000000000L){
        iter -= 1000000000L;
        md_iter++;
      }

  #ifdef LIENS_TRANSITIONS
      //linked transitions
      if (t_trans[tr].lien) simul_lien_trans(tr, ix, dir);
  #endif

      flag_evol_time = (t_trans[tr].time_mode != TIME_NO_EVOL);
    }

    else{
      /// CELL TRANSITION
      static int imax = 0;
      int ix;
      //LogPrintf("tr = %d   aleat = %f   ccc = %f   PII = %f\n", tr, aleat, ccc, PII);

      tr -= nb_trans_db;
      i=0;
      //we look for a cell of the desired type
      do{
        ix = (int) floor( drand48() * HLP);
        i++;
      }
      while (TE[ix].celltype != t_trans_cel[tr].depart);

      if (i>imax){
        imax = i;
        if (imax>100){
          LogPrintf ("WARNING: cell transition was difficult (%d tries)\n", imax);
        }
      }

      do_trans_cel(tr, ix);
      iter++; //LogPrintf("iter = %ld\n", iter);
      flag_trans_done = 1;
      flag_evol_time = 1;
    }
  }

  //LogPrintf("cpt_blk = %d\n", cpt_blk);

  if (!flag_trans_done){
    LogPrintf("simul_trans: too many blocked transitions (%d)\n", cpt_blk);
    cpt_blk = 0;
  }

  //time evolution
  if (flag_evol_time) simul_time();

  return flag_trans_done;
}


#ifdef TIME_SCALE
extern float lambda_C;
extern float lambda_T;

//#ifdef CGV
void simul_time_scale(float wind_value)
{
  FILE *fp;
  static int nb_rows;
  static float *tau1_values; //=grdvc_min
  static float *erosion_rates;
  static float *qsat_ratios;
  static float kappa; //Karman constant
  static float z_0; //roughness length
  static float z; //height of observation (velovity gradient)
  static float rho_air; //density of air
  static float rho_sand; //density of sand
  static float g; //gravity constant
  static float ut; //threshold velocity (erosion)
  static float d; //diameter of grain
  static float l_0; //length scale
  static char start=1;
  int i;

  if (!phys_prop_filename) return;

  if (start){
    //read physical parameters file
    LogPrintf("simul_time_scale: read physical parameters file %s\n", phys_prop_filename);

    if( (fp = fopen(phys_prop_filename,"r")) == NULL ) {
      ErrPrintf("ERROR: cannot open file %s\n",phys_prop_filename);
      exit(-4);
    }

    if (fscanf(fp,"kappa=%f\n",&kappa) != 1){
      ErrPrintf("ERROR: cannot read parameter kappa in %s\n",phys_prop_filename);
      exit(-4);
    }
    LogPrintf("kappa=%f\n", kappa);

    if (fscanf(fp,"z_0=%f\n",&z_0) != 1){
      ErrPrintf("ERROR: cannot read parameter z_0 in %s\n",phys_prop_filename);
      exit(-4);
    }
    LogPrintf("z_0=%f\n", z_0);

    if (fscanf(fp,"z=%f\n",&z) != 1){
      ErrPrintf("ERROR: cannot read parameter z in %s\n",phys_prop_filename);
      exit(-4);
    }
    LogPrintf("z=%f\n", z);

    if (fscanf(fp,"rhoair=%f\n",&rho_air) != 1){
      ErrPrintf("ERROR: cannot read parameter rho_air in %s\n",phys_prop_filename);
      exit(-4);
    }
    LogPrintf("rhoair=%f\n", rho_air);

    if (fscanf(fp,"rhosand=%f\n",&rho_sand) != 1){
      ErrPrintf("ERROR: cannot read parameter rho_sand in %s\n",phys_prop_filename);
      exit(-4);
    }
    LogPrintf("rhosand=%f\n", rho_sand);

    if (fscanf(fp,"g=%f\n",&g) != 1){
      ErrPrintf("ERROR: cannot read parameter g in %s\n",phys_prop_filename);
      exit(-4);
    }
    LogPrintf("g=%f\n", g);


    if (fscanf(fp,"d=%f\n",&d) != 1){
      ErrPrintf("ERROR: cannot read parameter d in %s\n",phys_prop_filename);
      exit(-4);
    }
    LogPrintf("d=%f\n", d);


    fclose(fp);

    ut = 0.1*sqrt((rho_sand/rho_air)*g*d);
    LogPrintf("ut=%f\n", ut);

    l_0=5*(rho_sand/rho_air)*d/4.;
    LogPrintf("l_0=%f\n", l_0);

    //read qsat data file
    LogPrintf("simul_time_scale: read qsat data file %s\n", qsat_data_filename);

    if( (fp = fopen(qsat_data_filename,"r")) == NULL ) {
      ErrPrintf("ERROR: cannot open file %s\n",qsat_data_filename);
      exit(-4);
    }

    //number of points
    if (fscanf(fp,"%d",&nb_rows) != 1){
      ErrPrintf("ERROR: cannot read the number of points in %s\n", qsat_data_filename);
      exit(-4);
    };
    LogPrintf("number of rows: %d\n", nb_rows);

    //read qsat data
    AllocMemory(tau1_values, float, nb_rows);
    AllocMemory(erosion_rates, float, nb_rows);
    AllocMemory(qsat_ratios, float, nb_rows);
    for(i=0; i<nb_rows; i++){
      if (fscanf(fp,"%f %f %f",&tau1_values[i], &erosion_rates[i], &qsat_ratios[i]) != 3){
        ErrPrintf("ERROR: cannot read qsat data in %s (row %d)\n",qsat_data_filename, i);
        exit(-4);
      };
      //LogPrintf("%d\t%f\t%f\t%f\n", i, tau1_values[i], erosion_rates[i], qsat_ratios[i]); //debug
    }

    fclose(fp);
  }

  //computation of the time_scale and tau1
  float ustar; //friction velovity
  float Qsat_ratio; //ratio of sand fluxes
  //float tau1; //threshold of the shear stress
  float qs; //sand flux in the nature
  float erosion_rate; //erosion rate
  float transport_rate; //transport rate
  float deposition_rate; //deposition rate
  float qsat; //sand flux in the model
  //float time_0; //time scale

  if (wind_value >= 0){
    //ustar= wind_value*kappa/log(z/z_0);
    ustar= wind_value*kappa/log(z/z_0);
    Qsat_ratio=1-(ut/ustar)*(ut/ustar);
    for (i=0; (i<nb_rows) && (Qsat_ratio<qsat_ratios[i]); i++);
    if (i<nb_rows){
      //grdvc_min = (tau1_values[i-1]+tau1_values[i])/2; //here we use PDF.data to get tau1 according to the value Qsat_ratio[i]
      //linear interpolation
      float t = (Qsat_ratio - qsat_ratios[i])/(qsat_ratios[i-1] - qsat_ratios[i]);
      grdvc_min = tau1_values[i-1]*t + tau1_values[i]*(1-t); //here we use PDF.data to get tau1 according to the value Qsat_ratio[i]
      erosion_rate = erosion_rates[i];
    }
    else{
      grdvc_min = tau1_values[nb_rows-1];
      erosion_rate = erosion_rates[nb_rows-1];
    }
    if (grdvc_min == 0) grdvc_min = 1.0; //value 0 is not allowed
    grdvc_max = grdvc_min + 1000.0;
    LogPrintf("grdvc_min = %f\n", grdvc_min);
    //LogPrintf("time_scale = %f\n", time_scale);
  }
  else{
    //if no wind data, then we compute ustar and erosion rate from tau1
    if (grdvc_min == 0) grdvc_min = 1.0; //value 0 is not allowed
    for (i=0; (i<nb_rows) && (grdvc_min>tau1_values[i]); i++);
    Qsat_ratio = qsat_ratios[i];
    erosion_rate = erosion_rates[i];
    if (Qsat_ratio < 1){
      ustar = ut/sqrtf(1 - Qsat_ratio);
    }
    else{
      ustar = 0; //no erosion
    }
  }
  LogPrintf("ustar=%f\n",ustar);
  if (ustar>ut){
    qs = 22.0*(rho_air/rho_sand)*sqrtf(d/g)*(ustar*ustar-ut*ut);  //qsat in the natural
    transport_rate = lambda_T;
    deposition_rate = lambda_C;
    qsat = transport_rate*erosion_rate/deposition_rate; //in model
    time_scale = qsat*(l_0*l_0)/qs; //t_0
  }
  else{
    //no erosion in this case
    //grdvc_min = 1000;
    //grdvc_max = grdvc_min + 1000.0;
    time_scale = 0; //not possible to set the time scale
    if (start) ErrPrintf("Warning: no time scale\n");
  }
  LogPrintf("time_scale = %g\n", time_scale);
  start = 0;
}
#endif

void simul_correl(int tr1, int ix, int ix2, int dir)
{
  static double cum2[MAX_LIENS];
  static int ln[MAX_LIENS], ind[MAX_LIENS];
  double ccc, aleat;
  int i, db, db2, n;
  int trans2_db_depart;
  //determination du doublet de depart
  db = type_doublet(ix, dir);
  db2 = type_doublet(ix2, dir);

  //construction du tableau des transitions correlees
  for(i=0, ccc=0.0, n=0; i<nb_liens; i++){
    if (t_lien[i].trans1 == tr1){
      trans2_db_depart = t_trans[t_lien[i].trans2].depart;
      if (t_lien[i].cel == 1){
        if (trans2_db_depart == db){
          ccc += t_lien[i].intensite;
          cum2[n] = ccc;
          ln[n] = i;
          ind[n] = ix;
          n++;
        }
      }
      else{
        if (trans2_db_depart == db2){
          ccc += t_lien[i].intensite;
          cum2[n] = ccc;
          ln[n] = i;
          ind[n] = ix2;
          n++;
        }
      }
    }
  }
  if (n > 0){ //test tableau non-vide ?
    //choix d'une transition selon les probabilites cumulees
    if (ccc > 1.0){
      ErrPrintf("ERROR : cumulated probabilites = %f > 1.0\n",ccc);
      exit(-1);
    }
    aleat = drand48();
    i=0;
    while ((i < n) && (aleat > cum2[i])) i++;
    if (i < n){
      //LogPrintf("Lien %d -> %d\n", tr1, ltr[i]);
      int tr2 = t_lien[ln[i]].trans2;
      /*if (t_trans[tr2].nb_chk){
        //controle de la transition
        if (!simul_check(tr2, ind[i], dir)) return; //transition bloquee
      }*/
      do_trans_db(tr2, ind[i], dir);
      iter++;
#ifdef INFO_TRANS
      cpt_lien[ln[i]]++;
#endif
      if (t_trans[tr2].lien) simul_lien_trans(tr2, ind[i], dir);
    }
  }
}


void simul_lien_trans(int tr, int ix, int dir)
{
  int ix2;
  int ix_ouest, ix2_ouest;
  int ix_haut, ix2_haut;
  int ix_nord, ix2_nord;

  if (t_trans[tr].classe == VERTICAL){
    // PREMIERE TRANSITION VERTICALE
    ix2 = ix+L;

    // liens horizontaux, doublet est
    simul_correl(tr, ix, ix2, EST);

    // liens horizontaux, doublet ouest
    ix_ouest = ix-1;
    ix2_ouest = ix2-1;
#ifdef CYCLAGE_HOR
    if (pbc_mode && (TE[ix_ouest].celltype == BORD)){
      ix_ouest += (L-2);
      ix2_ouest += (L-2);
    }
#endif
    simul_correl(tr, ix_ouest, ix2_ouest, EST);

    // liens horizontaux, doublet sud
    simul_correl(tr, ix, ix2, SUD);

    // liens horizontaux, doublet nord
    ix_nord = ix-HL;
    ix2_nord = ix2-HL;
#ifdef CYCLAGE_HOR
    if (pbc_mode && (TE[ix_nord].celltype == BORD)){
      ix_nord += (P-2)*HL;
      ix2_nord += (P-2)*HL;
    }
#endif
    simul_correl(tr, ix_nord, ix2_nord, SUD);
  }
  else {
    // PREMIERE TRANSITION HORIZONTALE
    int flag_eo = (dir == EST/*EST_OUEST*/);
    ix2 = flag_eo ? ix+1 : ix+HL;
#ifdef CYCLAGE_HOR
    if (pbc_mode && (TE[ix2].celltype == BORD)){
      ix2 -= flag_eo ? (L-2) : (P-2)*HL;
    }
#endif

    // liens verticaux, doublet bas
    simul_correl(tr, ix, ix2, BAS);

    // liens verticaux, doublet haut
    ix_haut = ix-L;
    ix2_haut = ix2-L;
    simul_correl(tr, ix_haut, ix2_haut, BAS);

    if (!flag_eo){
      // liens horizontaux, doublet est
      simul_correl(tr, ix, ix2, EST);

      // liens horizontaux, doublet ouest
      ix_ouest = ix-1;
      ix2_ouest = ix2-1;
#ifdef CYCLAGE_HOR
      if (pbc_mode && (TE[ix_ouest].celltype == BORD)){
        ix_ouest += (L-2);
        ix2_ouest += (L-2);
      }
#endif
      simul_correl(tr, ix_ouest, ix2_ouest, EST);
    }
    else{
      // liens horizontaux, doublet sud
      simul_correl(tr, ix, ix2, SUD);

      // liens horizontaux, doublet nord
      ix_nord = ix-HL;
      ix2_nord = ix2-HL;
#ifdef CYCLAGE_HOR
      if (pbc_mode && (TE[ix_nord].celltype == BORD)){
        ix_nord += (P-2)*HL;
        ix2_nord += (P-2)*HL;
      }
#endif
      simul_correl(tr, ix_nord, ix2_nord, SUD);
    }
  }
}


int simul_check(int tr, int ix, int dir)
{
  int i;
  int chk = 1;
  DataCheck *pdc;

  pdc = t_trans[tr].checks;

  for(i=0;  (i<t_trans[tr].nb_chk) && chk; i++, pdc++){
    pdc->dir = dir;

    //Appel de la callback de controle (si elle est definie)
    if (pdc->cel == 1){ // test sur la premiere cellule
      chk = pdc->func(ix, (void*)pdc);
    }
    else{ // test sur la deuxieme cellule
      int ix2;
      if (dir == EST){
        ix2 = ix+1;
#ifdef CYCLAGE_HOR
        if (pbc_mode && (TE[ix2].celltype == BORD)) ix2 -= (L-2);
#endif
      }
      else if (dir == SUD){ //NORD_SUD
        ix2 = ix+HL;
#ifdef CYCLAGE_HOR
        if (pbc_mode && (TE[ix2].celltype == BORD)) ix2 -= (P-2)*HL;
#endif
      }
      else { //BAS
        ix2 = ix+L;
      }
      chk = pdc->func(ix2, (void*)pdc);
    }
  }

  if (!chk){
    cpt_trans_blocked[tr]++;
#ifdef TRACE_TRANS
    //trace_trans_blocked(tr, ix);
#endif
  }

  return chk;
}


#ifdef SURRECTIONS
void simul_sur()
{
  //static int seuil = HLP*10;  //10000000L;
  static float time_threshold = 0.0;

  if (sur_delay){
    if (!time_threshold){
      time_threshold = sur_delay*ceil(csp_time/sur_delay);
      LogPrintf("delai de surrection = %f\n", sur_delay);
    }
    if (csp_time >= time_threshold){
      //regul_niveau();
      //if (calcule_hmin() > H/10) surrection_locale();
      //int nsur;
      //for(nsur=0; nsur<3; nsur++){
      //if (calcule_hmin() > 3) surrection();/*(type_sur == SUR_UNI) ? surrection() : surrection_locale();*/
      if (calcule_alti_max() < H-3) surrection();/*(type_sur == SUR_UNI) ? surrection() : surrection_locale();*/
      else { LogPrintf("surrection bloquee, sommet atteint !\n");  }
      //}
      time_threshold += sur_delay;
    }
  }
}
#endif

#ifdef ROTATIONS
//extern double lambda_E;
extern int Ncel[];           // nombre de cellules par type
extern float coef_injection; //regulation de l'injection de grain

void simul_rot()
{
  static char start = 1;
  static float time_threshold = 0.0;
  static int cpt_rot = 0;
  static short lg_seq = 5; //2
  static short seq_rot[] = {0, 1, 2, 3, 4};
  static short seq_del[] = {1, 1, 1, 1, 1};
  static float seq_inj[] = {1.0, 1.0};
  static int total_seq_rot = 0;
  static int seq_flag = 0; //flag for using sequences of angles/delays/...
  int nseq;
  //static double lambda_E_debut = 0.0;

  //if (!rot_delay || !rot_angle) return;
  if (!rot_delay) return;

  if (start){
    //seq_flag = (lg_seq > 1);
    if (!seq_flag) time_threshold = rot_delay*ceil(csp_time/rot_delay);
    LogPrintf("mode rotation = %s\n", (rot_mode == ROT_DISK) ? "table tournante" : "cyclage");
    LogPrintf("delai de rotation (par defaut) = %f\n", rot_delay);
    if (seq_flag) LogPrintf("sequence de periode %d\n", lg_seq);
#ifdef LGCA
    LogPrintf("nombre de cycles de stabilisation apres rotation = %d\n", rot_mvt_nb_cycles);
    if (rot_mvt_nb_cycles < VSTEP_TIME) ErrPrintf("WARNING: stabilization after rotation is too short ( Rot_ncycl = %d < VSTEP_TIME = %d )\n", rot_mvt_nb_cycles, VSTEP_TIME);
#endif
    //conditions_bords();
    start = 0;
  }
  if (csp_time >= time_threshold){
    //translation
    /*int alea_dx = (drand48()-0.5)*L/5;
    int alea_dz = (drand48()-0.5)*L/5;
    translation(alea_dx, alea_dz);*/

    if (seq_flag){
      //sequences de rotations
      //nseq = (short)lg_seq*drand48(); //random direction in sequence
      nseq = cpt_rot % lg_seq;
      LogPrintf("nseq=%d   seq_rot=%d   seq_del=%d\n", nseq, seq_rot[nseq], seq_del[nseq]);
      rotation(seq_rot[nseq]*rot_angle, rot_mode);
      total_seq_rot += seq_rot[nseq];
      //coef_injection = seq_inj[nseq];
    }
    else{
      //rotation d'angle rot_angle
      //short alea = 1 + (short)4.0*drand48();
      //short alea = (drand48()<0.5) ? 1 : 4;
      //rotation(alea*rot_angle, rot_mode);
      if (1) //(rot_angle)
        rotation(rot_angle, rot_mode);
      else
        rotation(360*drand48(), rot_mode);
    }
    cpt_rot++;

#if 0
//#ifdef CGV
    //intensité du vent
    //grdvc_min = drand48()*500.0;
    grdvc_min = seq_gvc[nseq];
    grdvc_max = grdvc_min + 1000.0;
    LogPrintf("variation d'intensite du vent : grdvc_min = %f   grdvc_max = %f   csp_time = %f\n", grdvc_min, grdvc_max, csp_time);
#endif

    //conditions_bords();
    //float delai_rot_alea = delai_rot*1.0*seq_del[nseq];
    //float delai_rot_alea = rot_delay*drand48();
    //float delai_rot_alea = delai_rot*0.5*(1.0+drand48());
    //time_threshold += delai_rot_alea;
    if (seq_flag){
      time_threshold += rot_delay*seq_del[nseq];
    }
    else
      time_threshold += rot_delay;
    //time_threshold = 10000000L;

    //LogPrintf("cpt_rot = %d   delai_rot = %f   total_seq_rot = %d   coef_injection = %f\n", cpt_rot, rot_delay, total_seq_rot, coef_injection);
    LogPrintf("cpt_rot = %d\n", cpt_rot);

#if defined(GR) || defined(GRJ)
    LogPrintf("nombre de grains : %d\n", Ncel[GR] + Ncel[GRJ]);
#ifdef AVALANCHES
    if (ava_mode == AVA_SYNC){
      avalanches(GR, ava_h_lim, ava_nb_cel_max, ALTI_MODE_BAS);
    }
#endif
#endif

#ifdef LGCA
    stabilisation_collisions(rot_mvt_nb_cycles);
#endif
  }
}
#endif

#ifdef AVALANCHES
//synchronized avalanches (SYNC mode)
void simul_ava_sync()
{
  static float time_threshold = 0.0;
  static unsigned long iter_0 = 0, md_iter_0;
  static unsigned long inter_iter_ava = 20000; //200000; //1000000; //nb min de transitions entre 2 avalanches
  static float ava_delay_sync = 0.0;
  static char start = 1;

  //if (!lambda_A) return;

  if (start){
#ifdef MODELE_AVA
    //ava_delay = 0.1;
    inter_iter_ava = 1;
#endif
#ifdef DUMP_SIGMA
    inter_iter_ava = 1;
#endif
    //inter_iter_ava *= (float)(LNS/100.0);
    //ava_delay_sync = 1.0 / lambda_A;
    ava_delay_sync = ava_delay;

    time_threshold = ava_delay_sync*ceil(csp_time/ava_delay_sync);//*drand48();
    iter_0 = iter;
    md_iter_0 = md_iter;

    LogPrintf("seuil temps avalanches = %f\n", time_threshold);
    LogPrintf("delai avalanches = %f\n", ava_delay_sync);
    LogPrintf("nb min iterations avant avalanches = %lu\n", inter_iter_ava);
    LogPrintf("ava_h_lim = %d\n", ava_h_lim);
    LogPrintf("ava_nb_cel_max = %d\n", ava_nb_cel_max);
    start = 0;
  }
//#if !defined(DUMP_SIGMA) && !defined(DUMP_AUTOCORREL)
#if 0
  if ((!csp_time) || ((csp_time >= time_threshold) && ((iter-iter_0 >= inter_iter_ava) || (md_iter-md_iter_0 > 0)))){
#else
  if ((!csp_time) || (csp_time >= time_threshold)){
#endif

#ifdef DUMP_SIGMA
//    dump_sigma_alti();
#endif
    avalanches(ALTI, ava_h_lim, ava_nb_cel_max, ALTI_MODE_BAS);
    if (ava_angle){
      avalanches_norm(ALTI, ava_nb_cel_max, ALTI_MODE_BAS);
    }

#ifdef DUMP_SIGMA
    dump_sigma_alti();
#endif
#ifdef DUMP_AUTOCORREL
    dump_autocorrel();
#endif
    //ava_delay_sync = 1.0 / lambda_A;
    time_threshold = csp_time + ava_delay_sync;
    iter_0 = iter;
    md_iter_0 = md_iter;
  }
}


// dynamical regulation of the rate and angles of avalanches
void simul_ava_dynamics()
{
  static float time_threshold = 0.0;
  static char start = 1;

#ifdef MODELE_AVA
// valeurs 2D
//  static float ava_delay_short = 500.0;
//  static float ava_delay_long = 4000.0;
// valeurs 3D
//  static float ava_delay_short = 1000.0;
//  static float ava_delay_long = 8000.0;
#else
#ifndef LGCA
// valeurs 2D
//  static float ava_delay_short = 4000.0;
//  static float ava_delay_long = 10000.0;
// valeurs 3D
//  static float ava_delay_short = 1000.0;
//  static float ava_delay_long = 4000.0;
#else
// valeurs 2D
//  static float ava_delay_short = 6000.0;
//  static float ava_delay_long = 12000.0;
// valeurs 3D
//  static float ava_delay_short = 1500.0;
//  static float ava_delay_long = 3000.0;
#endif
#endif

  if (ava_mode == AVA_PROPAG){
    if (start){
      time_threshold = 1 + ava_delay*ceil(csp_time/ava_delay);
      ava_angle_stable = ava_angle;
      ava_angle_stable_col = ava_angle_col;
      lambda_A_stable = lambda_A;
      //we set the unstable values for the transitions, if not null
      if (ava_angle_unstable) ava_angle = ava_angle_unstable;
#ifdef CELL_COLOR
      if (ava_angle_unstable_col) ava_angle_col = ava_angle_unstable_col;
#endif
      if (lambda_A_unstable) lambda_A = lambda_A_unstable;
      start = 0;
    }
    return;
  }

  if (start){
    time_threshold = ava_duration*ceil(csp_time/ava_duration);//*drand48();
    ava_angle_stable = ava_angle;
    lambda_A_stable = lambda_A;
    //ava_angle_unstable = atan(1.5*tan(ava_angle*PI/180.0))*180/PI; //pente * 1.5
    ava_angle_stable_col = ava_angle_col;
    //ava_angle_unstable_col = atan(1.2*tan(ava_angle_col*PI/180.0))*180/PI;
    LogPrintf("simul_ava_dynamics():\n");
    //LogPrintf("seuil temps avalanches = %f\n", time_threshold);
    LogPrintf("duration of avalanches = %f\n", ava_duration);
    LogPrintf("delay between avalanches = %f\n", ava_delay);
    LogPrintf("ava_angle = %f (stable)| %f (unstable)\n", ava_angle_stable, ava_angle_unstable);
#ifdef CELL_COLOR
    LogPrintf("ava_angle_col = %f (stable)| %f (unstable)\n", ava_angle_stable_col, ava_angle_unstable_col);
#endif
    //LogPrintf("nb min iterations avant avalanches = %lu\n", inter_iter_ava);
    start = 0;
  }

  if ((!csp_time) || (csp_time >= time_threshold)){
    if (ava_angle == ava_angle_stable){
      ava_angle = ava_angle_unstable;
      //lambda_A /= 3;
      lambda_A = lambda_A_unstable;
      time_threshold = csp_time + ava_delay;
    }
    else{
      ava_angle = ava_angle_stable;
      //lambda_A *= 3;
      lambda_A = lambda_A_stable;
      time_threshold = csp_time + ava_duration;
    }
    //LogPrintf("ava_angle = %f\n",ava_angle);
    //LogPrintf("lambda_A = %f\n",lambda_A);
#ifdef CELL_COLOR
    if (ava_angle_col == ava_angle_stable_col){
      ava_angle_col = ava_angle_unstable_col;
    }
    else{
      ava_angle_col = ava_angle_stable_col;
    }
    //LogPrintf("ava_angle_col = %f\n",ava_angle_col);
#endif //CELL_COLOR
  }
}
#endif //AVALANCHES

#ifdef CGV
//extern float lambda_E; //erosion rate
extern float grdvc_max; //seuil max pour le controle de l'erosion en fonction du gradient de vitesse
extern float grdvc_min; //seuil min pour le controle de l'erosion en fonction du gradient de vitesse

void simul_intensite_vent()
{
  static float delai_vent_max = 200.0; //1000.0; //10.0;
  static float time_threshold = 0.0;
  //static float pente_erosion_grdv = 0.0;
  //static float lambda_E_debut = 0.0;
  static char start = 1;

  if (start){
    time_threshold = delai_vent_max*ceil(csp_time/delai_vent_max); //delai_vent_max*0.5;
    //lambda_E_debut = lambda_E;
    //pente_erosion_grdv = lambda_E/(grdvc_max-grdvc_min);
    //LogPrintf("parametres d'intensite du vent au debut : lambda_E = %f   grdvc_min = %f   grdvc_max = %f\n", lambda_E, grdvc_min, grdvc_max);
    //LogPrintf("pente_erosion_grdv = %f\n", pente_erosion_grdv);
    LogPrintf("parametres d'intensite du vent au debut : grdvc_min = %f   grdvc_max = %f\n", grdvc_min, grdvc_max);
    LogPrintf("delai max de variation du vent = %f\n", delai_vent_max);
    start = 0;
  }
  if (csp_time >= time_threshold){
    //lambda_E = 2.0*lambda_E_debut*drand48();
    //lambda_E = lambda_E_debut/(1.1-drand48());
    //grdvc_min = grdvc_max - lambda_E/pente_erosion_grdv;
    grdvc_min = drand48()*500.0;
    grdvc_max = grdvc_min + 1000.0;
    float delai_vent_alea = delai_vent_max*drand48();
    //if (delai_vent_max>=500){ //affichage de l'intensite seulement en frequence basse
      //LogPrintf("variation d'intensite du vent : lambda_E = %f   grdvc_min = %f   csp_time = %f\n", lambda_E, grdvc_min, csp_time);
      LogPrintf("variation d'intensite du vent : grdvc_min = %f   grdvc_max = %f   csp_time = %f\n", grdvc_min, grdvc_max, csp_time);
      LogPrintf("delai de variation du vent = %f\n", delai_vent_alea);

    //}
    time_threshold += delai_vent_alea;
  }
}
#endif //CGV

#ifdef WIND_DATA
// load wind data from an ascii file
// format:
// nb_wind
// delay angle value
// ...
void simul_wind_data()
{
  static float *wind_delays = NULL;
  static float *wind_angles = NULL;
  static float *wind_values = NULL;
  static int nb_wind=0;
  static int time_unit=0;
  static int cpt_wind=0;
  static float time_threshold = 0.0;
  static float angle_dir = 0.0;
  float angle = 0.0;
  //float alpha = 0.0;
  float seq_duration = 0.0;
  char time_unit_str[128];
  static char start = 1;

  if (start){
    LogPrintf("simul_wind_data: read file %s\n", wind_data_filename);

    //open wind data file
    FILE *fp;
    if( (fp = fopen(wind_data_filename,"r")) == NULL ) {
      ErrPrintf("ERROR: cannot open file %s\n",wind_data_filename);
      exit(-4);
    }

    //number of wind regimes
    if (fscanf(fp,"%d,%s",&nb_wind,time_unit_str) != 2){
      ErrPrintf("ERROR: cannot read the number of wind regimes in %s\n",wind_data_filename);
      exit(-4);
    };
    LogPrintf("number of winds: %d\n",nb_wind);
    //LogPrintf("%s\n",time_unit_str);
    if (!strcmp(time_unit_str,"t0")){
      time_unit = UNIT_T0;
      LogPrintf("time unit is t0\n");
    } else if (!strcmp(time_unit_str,"rt")){
      time_unit = UNIT_RT;
      LogPrintf("time unit is physical\n");
#ifndef TIME_SCALE
      ErrPrintf("ERROR: no time scale, not possible to compute wind data from file %s\n",wind_data_filename);
      exit(-4);
#endif
    } else {
      ErrPrintf("ERROR: cannot read time unit in %s\n",wind_data_filename);
      exit(-4);
    }

    //read wind data
    int i;
    AllocMemory(wind_delays, float, nb_wind);
    AllocMemory(wind_angles, float, nb_wind);
    AllocMemory(wind_values, float, nb_wind);
    for(i=0; i<nb_wind; i++){
      if (fscanf(fp,"%f %f %f",&wind_delays[i], &wind_angles[i], &wind_values[i]) != 3){
        ErrPrintf("ERROR: cannot read wind data in %s\n",wind_data_filename);
        exit(-4);
      };
      //duration of the wind sequence in t0 unit
#ifdef TIME_SCALE
      if (time_unit == UNIT_RT){
        simul_time_scale(wind_values[i]);
        seq_duration += wind_delays[i]/time_scale;
      }
#endif
      if (time_unit == UNIT_T0){
        seq_duration += wind_delays[i];
      }
      if (i<5) LogPrintf("%d\t%f\t%f\t%f\n", i, wind_delays[i], wind_angles[i], wind_values[i]);
    }
    if (i>5) LogPrintf("...\n");
    LogPrintf("duration of the wind sequence = %f t0\n", seq_duration);

    fclose(fp);

    //process wind data
#ifdef TIME_SCALE
    if (time_unit == UNIT_RT){
      //real time
      //initial wind and time_threshold
      //time_threshold = wind_delays[cpt_wind]/time_scale;
      time_threshold = seq_duration*floor(csp_time/seq_duration);
      while (time_threshold <= csp_time){
        simul_time_scale(wind_values[cpt_wind]);
        time_threshold += wind_delays[cpt_wind++]/time_scale;
      }
      cpt_wind--;
      LogPrintf("cpt_wind = %d   time_threshold = %f\n", cpt_wind, time_threshold);
      //first wind orientation (north = 0 degree)
      angle = wind_angles[cpt_wind] + 90;
    }
#endif

    if (time_unit == UNIT_T0){
      //virtual time unit t0
      //initial wind and time_threshold
      time_threshold = seq_duration*floor(csp_time/seq_duration);
      while (time_threshold <= csp_time){
        time_threshold += wind_delays[cpt_wind++];
      }
      cpt_wind--;
      LogPrintf("cpt_wind = %d   time_threshold = %f\n", cpt_wind, time_threshold);
      //first wind orientation (west = 0 degree)
      angle = wind_angles[cpt_wind];
#ifdef CGV
      //first value of min shear stress
      grdvc_min = wind_values[cpt_wind];
#ifdef TIME_SCALE
      //time scale
      simul_time_scale(-1);
#endif
      LogPrintf("grdvc_min = %f\n", grdvc_min);
#endif
    }

    LogPrintf("simul_wind_data: direction = %f, rotation = %f, time = %e\n", wind_angles[cpt_wind], angle, csp_time);
    rotation(-angle, rot_mode);
    angle_dir = wind_angles[cpt_wind]; //first wind direction
    //LogPrintf("simul_wind_data: first direction = %f, time_threshold = %f\n", angle_dir,time_threshold);

    start = 0;
  }

  if (csp_time >= time_threshold){
    //we change the wind counter
    cpt_wind++;
    if (cpt_wind == nb_wind) cpt_wind = 0;
    //wind rotation
    angle = wind_angles[cpt_wind] - angle_dir;
    //angle = wind_angles[cpt_wind] ;
    LogPrintf("simul_wind_data: direction = %f, rotation = %f, time = %e\n", wind_angles[cpt_wind], angle, csp_time);
    rotation(-angle, rot_mode);
    //new wind direction
    angle_dir = wind_angles[cpt_wind];
#ifdef TIME_SCALE
    if (time_unit == UNIT_RT){
      //time scale
      simul_time_scale(wind_values[cpt_wind]);
      //real time
      time_threshold += wind_delays[cpt_wind]/time_scale;
      LogPrintf("real_time = %e\n", real_time);
    }
#endif

    if (time_unit == UNIT_T0){
      //virtual time unit t0
      time_threshold += wind_delays[cpt_wind];
#ifdef CGV
      //next value of min shear stress
      grdvc_min = wind_values[cpt_wind];
#ifdef TIME_SCALE
      //time scale
      simul_time_scale(-1);
      LogPrintf("real_time = %e\n", real_time);
#endif
      LogPrintf("grdvc_min = %f\n", grdvc_min);
#endif
    }
#ifdef LGCA
    if (angle) stabilisation_collisions(rot_mvt_nb_cycles);
#endif
  }
}
#endif //WIND_DATA


#ifdef LGCA
extern unsigned char opt_vss;
extern float *grdv;

void stabilisation_collisions(int nb_cyc)
{
  int i;
  char str_status[100];

//#ifdef SF_REAL
#ifdef CGV
  calcule_alti(ALTI, ALTI_MODE_BAS);
  calcule_normales();
#endif

  if (lgca_reset) init_mvt();

  LogPrintf("stabilisation des collisions : %d cycles\n", nb_cyc);
  push_status("flow stabilization (0%)",0);
  lgca_ready = 0;

  for (i=0; i<nb_cyc; i++) {
    //int j;
    //for (j=0; j<100000/L; j++) calcule_alti(GR, ALTI_MODE_BAS); //pour ralentir ...
    do_collisions();
    do_propagations();
#ifdef CGV
    if (opt_vss && grdv) calcule_grad_vel();
#endif
    if (!opt_nv){
      sprintf(str_status, "flow stabilization (%d%%)", (100*i)/nb_cyc);
      pop_status(0);
      push_status(str_status,0);
    }
    do_thread_sched();
  }
  lgca_ready = 1;
  LogPrintf("fin stabilisation des collisions\n");
  pop_status(0);
}
#endif

#ifdef ALTI
extern unsigned char opt_ls;

void simul_norm_lum()
{
  static float seuil_calcule_normales = 0.0;
  static float delai_calcule_normales = 1.0;

  if ((opt_ls) && (csp_time >= seuil_calcule_normales)){
    //recalcul des normales pour le rendu avec eclairage
    calcule_normales();
    seuil_calcule_normales += delai_calcule_normales;
  }
}
#endif

int simul_csp()
{
  int ii;
  int res;

#ifdef ALTI
  calcule_alti(ALTI, ALTI_MODE_BAS);
  simul_norm_lum();
#endif

  if (simul_dump_flag){
    lock_csp(0); //a rotation can occur here
    simul_dump();
    unlock_csp(0);
  }

#ifdef WIND_DATA
  if (wind_data_filename){
    lock_csp(0); //a rotation can occur here
    simul_wind_data();
    unlock_csp(0);
  }
#endif

#ifdef TIME_SCALE
  // time scale computation
  if (!wind_data_filename) simul_time_scale(-1);
#endif

#ifdef LGCA
  double lgca_time_threshold = csp_time;
  if (!lgca_delay){
#ifdef MODELE_DUN
    lgca_delay = 1.0/NB_MVT_EO; //1.0;
    //lgca_delay = 1.5/NB_MVT_EO; //1.0;
#else
    //we compute lgca delay from transition rates
    double lambda_max=0;
    for(i=0; i<nb_trans_db; i++){
      if (t_trans[i].intensite > lambda_max) lambda_max = t_trans[i].intensite;
    }
    lgca_delay = 1.0 / lambda_max;
#endif
  }
  LogPrintf("lgca_delay = %f\n", lgca_delay);
#endif //LGCA

  /// CSP ENGINE START
#ifdef DUMP_SIGNATURE
  if (opt_info) dump_signature(0);
#endif

#ifdef DUMP_SIGMA
//  dump_sigma_alti();
#endif

  lock_csp(0);

#ifdef AVALANCHES
  if (ava_mode == AVA_SYNC) simul_ava_sync();

  if (ava_delay) simul_ava_dynamics();
#endif //AVALANCHES

#ifdef LGCA
//#if 0
//stabilisation des collisions
  if (init_mvt_nb_cycles<0) init_mvt_nb_cycles = 2*CLEO; //default value
  else if (init_mvt_nb_cycles<VSTEP_TIME) init_mvt_nb_cycles = VSTEP_TIME; //min value
  if (nom_fic_mvt) init_mvt_nb_cycles = VSTEP_TIME;
  stabilisation_collisions(init_mvt_nb_cycles);
  if (rot_mvt_nb_cycles<0) rot_mvt_nb_cycles = lgca_reset ? 2*CLEO : 3*VSTEP_TIME;
#ifdef CGV
  //init_cgv();
  if (opt_info) dump_cgv();
#endif
#endif //LGCA

  ii = 0;
  ResetMemory(cpt_trans_blocked, int, MAX_TRANSITIONS_DB);

  /// MAIN LOOP
  while(!end_of_simul)
  {

    /*if (!opt_nv)*/
    do_thread_sched();

    //for (i=0; i<100000/L; i++) calcule_alti(GR, ALTI_MODE_BAS); //pour ralentir ...

#ifdef LGCA
    if ((PII == 0) || (csp_time >= lgca_time_threshold)){
      /// LATTICE GAZ
#ifdef ALTI
      //calcule_alti(GR, ALTI_MODE_BAS);
      calcule_normales();
#endif
#ifdef CGV
      calcule_grad_vel();
      compute_coef_cgv();
#endif
      //LogPrintf("debut collisions...");
      do_collisions();
      //LogPrintf("propagations...");
      do_propagations();
      //LogPrintf("fin collisions...\n");

      lgca_time_threshold += lgca_delay;
      while (lgca_time_threshold < csp_time) lgca_time_threshold += lgca_delay;
    }
#elif defined(ALTI)
    simul_norm_lum();
#endif

    //we compute the probability distribution
    compute_prob_dist();

#ifndef LGCA
    //we check for the presence of active doublets
    if (PII == 0){
      LogPrintf("no active doublets\n");
      end_of_simul = 1;
      break;
    }
#else
    //if no active doublets, then we only compute the lattice gaz
    if (PII == 0) continue;
#endif

    //stochastic transition
    res = simul_trans();
    if (!res){
#ifdef LGCA
      LogPrintf("simul_csp: time evolution forced (%f)\n", csp_time);
      if (csp_time < lgca_time_threshold){
        csp_time = lgca_time_threshold;
      }
      else{
        //lgca speedup
        while (lgca_time_threshold < csp_time) lgca_time_threshold += lgca_delay;
      }
#endif
#ifdef MODELE_LIFE
      LogPrintf("simul_csp: stopping simulation ...\n");
      end_of_simul = 1;
#endif
    }
    if (stop_delay_t0 && (csp_time > stop_delay_t0)){
      LogPrintf("time to quit\n");
      push_status("stopping ...", 0);
      sleep(1);
      exit(0);
    }

    if (simul_dump_flag) simul_dump();

#ifdef SURRECTIONS
    simul_sur();
#endif

#ifdef ROTATIONS
    simul_rot();
#endif

#ifdef AVALANCHES
    if (ava_mode == AVA_SYNC) simul_ava_sync();

    if (ava_delay) simul_ava_dynamics();
#endif //AVALANCHES

#ifdef CGV
    //simul_intensite_vent();
#endif
#ifdef WIND_DATA
    if (wind_data_filename) simul_wind_data();
#endif

    ii++;

#ifdef DUMP_SIGNATURE
    if (opt_info &&  !(ii & 0x000fffff)){ //toutes les 1048576 transitions
    	//empreinte terrestre
      dump_signature(ii);
    }
#endif

#if defined(CELL_COLOR) && defined(MODELE_AVA)
   /*if (!(ii & 0x00007fff)){ //toutes les 65536 transitions
     init_color();
   }*/
#endif

#ifdef TRACE3D_CEL
    //if (!(ii & 0x000001fff)) trace3d_cel();  //toutes les 8192 transitions
    //if (!(ii & 0x00000000)) trace3d_cel();  //toutes les 256 transitions
    trace3d_cel();  //a chaque transition
#endif

#ifdef DUMP_SIGMA
    //if (!(ii & 0x000000fff)) dump_sigma_alti(); //toutes les 4096 transitions
    //if (!(ii & 0x0000000ff)) dump_sigma_alti(); //toutes les 256 transitions
#endif

#ifdef DUMP_AUTOCORREL
    if (!(ii & 0x00000ffff)) dump_autocorrel();
#endif

#ifdef PARALLEL
    //if (mode_par && !(ii & 0xf)) {
    if (mode_par && !(ii & 0xff)) {
      //if (synchro_terre()) {
        //LogPrintf("pause synchro\n");
        while (synchro_terre());
        //LogPrintf("synchro ok\n");
      //}
    }
#endif

  }
  //END OF MAIN LOOP

  unlock_csp(0);

  LogPrintf("end of simulation\n");

#ifdef LISSAGE
  lissage_terre();
#endif

  //return (PII != 0);
  return (end_of_simul);
}


void simul_dump()
{
  static float dump_delay=0.0;
  static double time_threshold = 0.0;
  static double time_threshold_dpng = 0.0;
  static double time_threshold_dcsp = 0.0;
  static int cpt_dump = 0;
  static char name[100];
  static char str[100];
  static double t0=0, t1=0;
  static char start=1;

  if (start){
#ifdef PARALLEL
    extern int proc_id;
    sprintf(str,"_%d",proc_id);
#else
    *str = 0;
#endif
    //time_threshold = csp_time;
    //set dump_delay = minimal non-zero delay
    if (dump_delay_png > 0){
      //time_threshold_dpng = csp_time;
      time_threshold_dpng = dump_delay_png*ceil(csp_time/dump_delay_png); //time_threshold_dpng is mulitple of dump_delay_png
      dump_delay = dump_delay_png;
      LogPrintf("dump_delay_png = %f (t0)\n", dump_delay_png);
    }
    if (dump_delay_csp > 0){
      //time_threshold_dcsp = csp_time;
      time_threshold_dcsp = dump_delay_csp*ceil(csp_time/dump_delay_csp); //time_threshold_dcsp is mulitple of dump_delay_csp
      if (!dump_delay || (dump_delay_csp < dump_delay)) dump_delay = dump_delay_csp;
      LogPrintf("dump_delay_csp = %f (t0)\n", dump_delay_csp);
    }
    time_threshold = Min(time_threshold_dpng, time_threshold_dcsp);
    LogPrintf("time_threshold = %f\n", time_threshold);
    LogPrintf("dump_delay = %f (t0)\n", dump_delay);
    start = 0;
  }

  if ((csp_time >= time_threshold) || end_of_simul){
    //no more than 10 files per sec., for security
    elapsed(&t1);

    if (!t0 || (t1-t0 >= 0.1)){
      lock_display(0);

#ifdef REORIENT_AUTO
      reorient_flag = 1;
      rotation(0, rot_mode | ROT_REORIENT_TEMP);
      //display wind direction
      vdir_mode = VDIR_WIND;
#endif

      if (dump_delay_png && (time_threshold == time_threshold_dpng)){
        //dump PNG file
#if !defined(USE_LIBPNG) && !defined(USE_GD)
        gcallback_flash(0); //update image buffer
#endif
        sprintf(name,"%s%05d_t0%s.png", MOD_NAME, cpt_dump, str);
        dump_image(name, "png");
        while (time_threshold_dpng <= csp_time) time_threshold_dpng += dump_delay_png;
      }

      if (dump_delay_csp && (time_threshold == time_threshold_dcsp)){
        //dump CSP file
        //LogPrintf("dump CSP, cpt_dump = %d, csp_time = %f t0\n", cpt_dump, csp_time);
        //sprintf(name,"%s%05d_t0.csp", MOD_NAME, cpt_dump);
        //dump_terre(name, DUMP_CSP);
        dump_terre(DUMP_CSP, cpt_dump, UNIT_T0);
#ifdef ALTI
        //sprintf(name,"ALTI%05d.data", cpt_dump);
        //dump_surface(name);
        dump_surface(cpt_dump, UNIT_T0);
#endif
#if defined(LGCA) //&& !defined(ROTATIONS)
        if (csphpp_flag){
          dump_mvt(cpt_dump, UNIT_T0);
          dump_grad_vel(cpt_dump, UNIT_T0);
        }
#endif
        while (time_threshold_dcsp <= csp_time) time_threshold_dcsp += dump_delay_csp;
      }

#ifdef REORIENT_AUTO
      rotation(0, rot_mode | ROT_REORIENT_UNDO);
      vdir_mode = VDIR_NONE;
      reorient_flag = 0;
#endif

      unlock_display(0);

      t0 = t1;

      //update cpt
      cpt_dump++;
      cpt_dump = cpt_dump % 100000; //no more than 100000 files, for security
    }
    else{
      if (time_threshold == time_threshold_dpng) time_threshold_dpng += dump_delay_png;
      if (time_threshold == time_threshold_dcsp) time_threshold_dcsp += dump_delay_csp;
      //LogPrintf("simul_dump: no dump)n");
    }

    //next time threshold
    while (time_threshold <= csp_time) time_threshold += dump_delay;
    if (dump_delay_png && (time_threshold_dpng < time_threshold)) time_threshold = time_threshold_dpng;
    if (dump_delay_csp && (time_threshold_dcsp < time_threshold)) time_threshold = time_threshold_dcsp;
    //LogPrintf("simul_dump: time_threshold = %f\n",time_threshold);
  }
}


void dump_time()
{
  static int cpt = 0;
  static unsigned long int md_iter_0 = 0, iter_0 = 0;
  static double csp_time_0 = 0;
  static double real_time_0 = 0;
  long int delta_iter, delta_md_iter;
	FILE *fp;

  fp = fopen("TIME.log","a");
  if (!cpt){
    fprintf(fp,"      nb trans.        delta trans.    time            delta time       ");
#ifdef LGCA
    fprintf(fp,"lgca cyc.    ");
#endif
    csp_time_0 = csp_time;
#ifdef TIME_SCALE
    fprintf(fp,"real time      delta real time");
    real_time_0 = real_time;
#endif
  }
  delta_iter = iter - iter_0;
  delta_md_iter = md_iter - md_iter_0;
  if ((delta_md_iter > 0) && (delta_iter < 0))
  {
    delta_iter += 1000000000;
    delta_md_iter--;
  }

  if (delta_md_iter)
    fprintf(fp,"\n%04d: %04lu%09lu %03lu%09lu       %e    %e     ", cpt++, md_iter, iter, delta_md_iter, delta_iter, csp_time, csp_time - csp_time_0);
  else
    fprintf(fp,"\n%04d: %04lu%09lu    %09lu       %e    %e     ", cpt++, md_iter, iter, delta_iter, csp_time, csp_time - csp_time_0);
#ifdef LGCA
  fprintf(fp,"%09d    ", col_iter);
#endif
#ifdef TIME_SCALE
  fprintf(fp,"%e   %e", real_time, real_time - real_time_0);
#endif
  fclose(fp);

  md_iter_0 = md_iter;
  iter_0 = iter;
  csp_time_0 = csp_time;
#ifdef TIME_SCALE
  real_time_0 = real_time;
#endif
}
