diff --git a/exploration/GR3D_Rdescription/GR3Dmetapopulation.R b/exploration/GR3D_Rdescription/GR3Dmetapopulation.R
new file mode 100644
index 0000000000000000000000000000000000000000..cd1dbe14251d668fe6a760618649a1ac750ab9d8
--- /dev/null
+++ b/exploration/GR3D_Rdescription/GR3Dmetapopulation.R
@@ -0,0 +1,101 @@
+
+#=============================================================================================================
+# metapopulation identification
+# =============================================================================================================
+computeAutochtonousRate = function(data){
+  data %>% filter(migrationBasin == originBasin) %>% 
+    select(-originBasin) %>% 
+    inner_join(data %>% group_by(year, migrationBasin) %>% 
+                 summarise(totalRun = sum(effective), .groups = 'drop'), by = c('migrationBasin', 'year')) %>% 
+    mutate(autochtonousRate = effective / totalRun) %>% 
+    select(year, migrationBasin, autochtonousRate)
+}
+
+computeHomingRate = function(data){
+  data %>% filter(migrationBasin == originBasin) %>% 
+    select(-migrationBasin) %>% 
+    inner_join(data %>% group_by(year, originBasin) %>% 
+                 summarise(production = sum(effective), .groups = 'drop'), by = c('originBasin', 'year')) %>% 
+    mutate(homingRate = effective / production) %>% 
+    select(year, originBasin, homingRate)
+}
+
+# comparison between autochtonous rate and homing rate
+# computeAutochtonousRate(exchangesData) %>% 
+#   inner_join(computeHomingRate(data = exchangesData), by =  c('migrationBasin' = 'originBasin', 'year'))
+
+
+
+identyMetapopulation = function(exchangesData, threshold = .95, verbose = FALSE) {
+  # exchangesData {year, migrationBasin, originBasin)}
+  
+  exchangesDataUpdated = exchangesData
+  
+  # initialise the connection table between basin and metapopulation
+  metapopulation = exchangesDataUpdated %>% 
+    distinct(year, basin = migrationBasin) %>% 
+    mutate(metapop = basin)
+  
+  iteration = 1
+  
+  # find the basin with the minimum autochtonous rate for the first time,
+  basinWithMinAutochtonousRate = computeAutochtonousRate(data = exchangesDataUpdated) %>%
+    group_by(year) %>%
+    slice(which.min(autochtonousRate))
+  
+  # loop while autochtonousRate is still < 0.95
+  
+  while (basinWithMinAutochtonousRate$autochtonousRate < threshold ) {
+    #  while (iteration <= 4 ) {
+    if (verbose) cat(iteration, ": ", basinWithMinAutochtonousRate$autochtonousRate, '\n' )
+    
+    # basinWithMinAutochtonousRate will be merged with origin basin sending the maximum number of fish in this basin
+    metapopsToBeMerged <-   exchangesDataUpdated %>% 
+      inner_join(basinWithMinAutochtonousRate %>% select(-autochtonousRate), 
+                 by = c('year', 'migrationBasin' )) %>%
+      filter(migrationBasin != originBasin) %>% # to avoid self merging
+      group_by(year) %>% 
+      slice(which.max(effective)) %>% 
+      ungroup() %>% 
+      select(year, migrationBasin, originBasin)
+    
+    # update the connection table between basin and metapopulation by merging metapopulation
+    for (i in 1:nrow(metapopsToBeMerged)) {
+      mergedName = paste0("M", iteration, '_',i)
+      metapopsToBeMerged_i <- metapopsToBeMerged %>% slice(i) 
+      if (verbose) cat("\t", metapopsToBeMerged_i$migrationBasin, ' + ', metapopsToBeMerged_i$originBasin, ' = ', mergedName, "\n")
+      
+      metapopulation <- metapopulation %>% 
+        mutate(metapop = if_else(year == metapopsToBeMerged_i$year &
+                                   metapop %in%  (metapopsToBeMerged_i %>% select(ends_with("Basin")) %>% unlist(use.names = FALSE)), 
+                                 mergedName,
+                                 metapop)) 
+    }
+    # sum by migration and origin basins according to updated metapopulations
+    # computed with initial exchangesData
+    exchangesDataUpdated <- exchangesData %>% 
+      #sum up on migration basins
+      inner_join(metapopulation, 
+                 by = c('year', 'migrationBasin' = 'basin')) %>% 
+      group_by(year, metapop, originBasin) %>% 
+      summarise(effective = sum(effective), .groups = 'drop') %>% 
+      rename(migrationBasin = metapop) %>% 
+      # sum on origine basins
+      inner_join(metapopulation, 
+                 by = c('year', 'originBasin' = 'basin')) %>% 
+      group_by(year, metapop, migrationBasin) %>% 
+      summarise(effective = sum(effective), .groups = 'drop') %>% 
+      rename(originBasin = metapop)
+    
+    # minimum of autochtonous rate
+    basinWithMinAutochtonousRate = computeAutochtonousRate(data = exchangesDataUpdated) %>%
+      group_by(year) %>%
+      slice(which.min(autochtonousRate))
+    
+    iteration = iteration + 1
+  }
+  return(list(metapopulation = metapopulation, exchangesData = exchangesDataUpdated))
+}
+
+
+
diff --git a/exploration/GR3D_Rdescription/thermalRanges.R b/exploration/GR3D_Rdescription/thermalRanges.R
new file mode 100644
index 0000000000000000000000000000000000000000..cf5d758f1293a4138f5a282ba81361ce7fbab178
--- /dev/null
+++ b/exploration/GR3D_Rdescription/thermalRanges.R
@@ -0,0 +1,17 @@
+library(tidyverse)
+
+#TODO extract from xml
+
+thermalRange = tibble(process = "Grow", Tmin = 1.7, Topt = 4.5, Tmax = 27.9, type = "Rosso") %>% 
+  bind_rows(tibble(process = "SurvivalSpawnerInRiv", Tmin = 5.4, Topt = 16.7, Tmax = 27.5, type = "Rosso") ) %>% 
+  bind_rows(tibble(process = "Reproduction", Tmin = 5.1, Topt = 13.5, Tmax = 24.5, type = "Rosso") ) %>% 
+  bind_rows(tibble(process = "SurviveAfterReproduction", Tmin = -Inf, Topt = 19.58 - log(19)/0.58, Tmax = 19.58 + log(19)/0.58, type = "logit") )
+
+
+thermalRange %>% mutate(numero = row_number(),
+                        process = factor(process, 
+                                         levels = c("Grow","SurvivalSpawnerInRiv", "SurviveAfterReproduction", "Reproduction" ))) %>% 
+  ggplot() +
+  geom_segment(aes(x = Tmin, y = process, xend = Tmax, yend = process)) +
+  geom_point(aes(x = Topt, y = process)) 
+
diff --git a/exploration/NEA_sensitivity_analysis/extractNumeroProcess.R b/exploration/NEA_sensitivity_analysis/extractNumeroProcess.R
index 6d6e57d9e7678f890b6bbc00a1c1384ad50c6485..80bd310c0ba9bfc3ffd044c6c2ada6c667002b25 100644
--- a/exploration/NEA_sensitivity_analysis/extractNumeroProcess.R
+++ b/exploration/NEA_sensitivity_analysis/extractNumeroProcess.R
@@ -1,4 +1,5 @@
 library(tidyverse)
+library(XML)
 #source("../GR3D_Rdescription/GR3D_NEA_XML_parameters.R")
 
 setwd("~/Documents/workspace/GR3D/exploration/NEA_sensitivity_analysis")
diff --git a/exploration/NEA_sensitivity_analysis/metapopulationIdentification.R b/exploration/NEA_sensitivity_analysis/metapopulationIdentification.R
index 61ddeb400ff91bf41b8badc9f4859fe62288aaee..03b1c22efb62721f7255a9082822393dade46e3e 100644
--- a/exploration/NEA_sensitivity_analysis/metapopulationIdentification.R
+++ b/exploration/NEA_sensitivity_analysis/metapopulationIdentification.R
@@ -1,102 +1,11 @@
-library(tidyverse)
-
-computeAutochtonousRate = function(data){
-  data %>% filter(migrationBasin == originBasin) %>% 
-    select(-originBasin) %>% 
-    inner_join(data %>% group_by(year, migrationBasin) %>% 
-                 summarise(totalRun = sum(effective), .groups = 'drop'), by = c('migrationBasin', 'year')) %>% 
-    mutate(autochtonousRate = effective / totalRun) %>% 
-    select(year, migrationBasin, autochtonousRate)
-}
-
-computeHomingRate = function(data){
-  data %>% filter(migrationBasin == originBasin) %>% 
-    select(-migrationBasin) %>% 
-    inner_join(data %>% group_by(year, originBasin) %>% 
-                 summarise(production = sum(effective), .groups = 'drop'), by = c('originBasin', 'year')) %>% 
-    mutate(homingRate = effective / production) %>% 
-    select(year, originBasin, homingRate)
-}
-
-# comparison between autochtonous rate and homing rate
-# computeAutochtonousRate(exchangesData) %>% 
-#   inner_join(computeHomingRate(data = exchangesData), by =  c('migrationBasin' = 'originBasin', 'year'))
-
-
-
-identyMetapopulation = function(exchangesData, threshold = .95, verbose = FALSE) {
-  exchangesDataUpdated = exchangesData
-  
-  # initialise the connection table between basin and metapopulation
-  metapopulation = exchangesDataUpdated %>% 
-    distinct(year, basin = migrationBasin) %>% 
-    mutate(metapop = basin)
-  
-  iteration = 1
-  
-  # find the basin with the minimum autochtonous rate for the first time,
-  basinWithMinAutochtonousRate = computeAutochtonousRate(data = exchangesDataUpdated) %>%
-    group_by(year) %>%
-    slice(which.min(autochtonousRate))
-  
-  # loop while autochtonousRate is still < 0.95
-  
-  while (basinWithMinAutochtonousRate$autochtonousRate < threshold ) {
-    #  while (iteration <= 4 ) {
-    if (verbose) cat(iteration, ": ", basinWithMinAutochtonousRate$autochtonousRate, '\n' )
-    
-    # basinWithMinAutochtonousRate will be merged with origin basin sending the maximum number of fish in this basin
-    metapopsToBeMerged <-   exchangesDataUpdated %>% 
-      inner_join(basinWithMinAutochtonousRate %>% select(-autochtonousRate), 
-                 by = c('year', 'migrationBasin' )) %>%
-      filter(migrationBasin != originBasin) %>% # to avoid self merging
-      group_by(year) %>% 
-      slice(which.max(effective)) %>% 
-      ungroup() %>% 
-      select(year, migrationBasin, originBasin)
-    
-    # update the connection table between basin and metapopulation by merging metapopulation
-    for (i in 1:nrow(metapopsToBeMerged)) {
-      mergedName = paste0("M", iteration, '_',i)
-      metapopsToBeMerged_i <- metapopsToBeMerged %>% slice(i) 
-      if (verbose) cat("\t", metapopsToBeMerged_i$migrationBasin, ' + ', metapopsToBeMerged_i$originBasin, ' = ', mergedName, "\n")
-      
-      metapopulation <- metapopulation %>% 
-        mutate(metapop = if_else(year == metapopsToBeMerged_i$year &
-                                   metapop %in%  (metapopsToBeMerged_i %>% select(ends_with("Basin")) %>% unlist(use.names = FALSE)), 
-                                 mergedName,
-                                 metapop)) 
-    }
-    # sum by migration and origin basins according to updated metapopulations
-    # computed with initial exchangesData
-    exchangesDataUpdated <- exchangesData %>% 
-      #sum up on migration basins
-      inner_join(metapopulation, 
-                 by = c('year', 'migrationBasin' = 'basin')) %>% 
-      group_by(year, metapop, originBasin) %>% 
-      summarise(effective = sum(effective), .groups = 'drop') %>% 
-      rename(migrationBasin = metapop) %>% 
-      # sum on origine basins
-      inner_join(metapopulation, 
-                 by = c('year', 'originBasin' = 'basin')) %>% 
-      group_by(year, metapop, migrationBasin) %>% 
-      summarise(effective = sum(effective), .groups = 'drop') %>% 
-      rename(originBasin = metapop)
-    
-    # minimum of autochtonous rate
-    basinWithMinAutochtonousRate = computeAutochtonousRate(data = exchangesDataUpdated) %>%
-      group_by(year) %>%
-      slice(which.min(autochtonousRate))
-    
-    iteration = iteration + 1
-  }
-  return(list(metapopulation = metapopulation, exchangesData = exchangesDataUpdated))
-}
+require(tidyverse)
 
+source(GR3Dmetapopulation)
 # upload data =====
 # longer table of effective of different origin basins in each migration basin 
 exchangesData <- read_csv("../../data/output/northeastamerica/effectiveFluxes_1-observed.csv")
 
+
 exchangesData %>% group_by(year, basin = originBasin) %>% 
   summarise(production = sum(effective), .groups = 'drop') %>% 
   inner_join(exchangesData %>% group_by(year, basin = migrationBasin) %>% 
diff --git a/exploration/NEA_sensitivity_analysis/tryFromJava.R b/exploration/NEA_sensitivity_analysis/tryFromJava.R
index a47b94dedc2b72a1aa5feae142ea4c9596f71862..1a4a7b1648f75ccf275c6f5c987624763ffc8b0d 100644
--- a/exploration/NEA_sensitivity_analysis/tryFromJava.R
+++ b/exploration/NEA_sensitivity_analysis/tryFromJava.R
@@ -3,15 +3,16 @@ library(rJava) # see J4R
 #library(yaml)
 library(tidyverse)
 library(jdx)
+library(tictoc)
 
-rm(list = ls())
 
+rm(list = ls())
 
 ## to have the same working directory as GR3D
 #setwd("~/Documents/workspace/GR3D/exploration/NEA_sensitivity_analysis")
 setwd("../..")
 
-
+source("exploration/GR3D_Rdescription/GR3Dmetapopulation.R")
 # ====================================== LOCAL FUNCTIONS =====================================
 convertFromSALtoR = function(java_out){
   require(jdx)
@@ -22,8 +23,8 @@ convertFromSALtoR = function(java_out){
   R_out = lapply(seq(2, length(java_out)), FUN = function(i){convertToR(java_out[[i]])}) %>% 
     unlist() %>% 
     matrix(nrow = length(java_out) - 1, byrow = TRUE, dimnames = list(NULL, headers) ) %>%
-    data.frame() %>% 
-    suppressMessages(type_convert()) 
+    as_tibble() %>% 
+    type_convert()
   
   return(R_out)
 }
@@ -41,7 +42,7 @@ exploreEasyRun = function(){
   # 1. northernmost populated basin_id   
   # 2. centroid latitude
   out_range = .jcall("analysis.EasyRun","[[Ljava/lang/String;","getValuesFromAquanismGroupProcess","processes.processesEachStep.10.exportToR")
-  range = convertFromSALtoR(out_range)
+  range = suppressMessages(convertFromSALtoR(out_range))
   
   # 3. log-likelihood ----
   #       the proportion of reproductions leading to a recruitment > 50 over the 1920-1950 period was used 
@@ -55,66 +56,79 @@ exploreEasyRun = function(){
   #     then in R: annual slope of the linear regression
   #     and average over the period
   # 5. mean age of female and age
-  #     annual mean age in each basin, 
-  #     then in R: average per year,
-  #     and average over the period (1920 and 1950)
+  #     annual mean age over the period (1920 and 1950) in each basin, 
+  #     and average over basin
   out_spawners_features = .jcall("analysis.EasyRun","[[Ljava/lang/String;","getValuesFromAquanismGroupProcess","processes.processesEachStep.9.exportToR")
-  spawners_features = convertFromSALtoR(out_spawners_features)
+  spawners_features =suppressMessages (convertFromSALtoR(out_spawners_features))
   
-  ## ICI coder pente de la droite moyenne des ages
+  slopePrimiparous =  lm(pct_primiparous~latitude, data = spawners_features)$coeff[2]
+  meanAgeFemale = mean(spawners_features$mean_age_female, na.rm = TRUE)
+  meanAgeMale = mean(spawners_features$mean_age_male, na.rm = TRUE)
   
   # ------------------------------------------------------------ #
   # 5. number of metapopulation based on the exchange matrix ----
   #     mean annual exchange matrix between 1920 and 1950
   out_exchanges = .jcall("analysis.EasyRun","[[Ljava/lang/String;","getValuesFromAquanismGroupProcess","processes.processesEachStep.16.exportToR")
-  exchanges = convertFromSALtoR(out_exchanges)
+  exchanges = suppressMessages(convertFromSALtoR(out_exchanges))
+  
+  # keep only basin with a totalRun  (sum over migrationBasin) >0
+  selectedBasin <- exchanges %>% group_by(basin = migrationBasin) %>% 
+    summarise(totalRun = sum(effective), .groups = 'drop') %>% 
+    filter(totalRun > 0) %>% select(basin)
   
-  ## ICI  lancer metapopulationIdentification
+  metapopulations <- exchanges %>% inner_join(selectedBasin, by=c('migrationBasin' = 'basin')) %>%
+    inner_join(selectedBasin, by=c('natalBasin' = 'basin')) %>% 
+    mutate(year =0) %>% rename(originBasin = natalBasin) %>% 
+    identyMetapopulation(., threshold = .95, verbose = FALSE)
   
+  nbMetapopulation = metapopulations$metapopulation %>% 
+    distinct(metapop) %>% 
+    summarise(n = n()) %>% 
+    unlist(use.names = FALSE)
   
-  return(list(northern_basin_id = range$northern_basin_id, 
-              effective_centroid_latitude = range$effective_centroid_latitude,
-              logL = logL,
-              spawners_features = spawners_features,
-              exchanges = exchanges))
+  return(tibble(northern_basin_id = range$northern_basin_id, 
+                effective_centroid_latitude = range$effective_centroid_latitude,
+                logL = logL,
+                slopePrimiparous = slopePrimiparous,
+                meanAgeFemale = meanAgeFemale,
+                meanAgeMale = meanAgeMale, 
+                nbMetapopulation = nbMetapopulation))
 }
 
 
 
 # GR3D arguments -----------------------------------------
-jarfile = "target/GR3D-3.2-SNAPSHOT.jar"
-simDuration = 150
-simBegin = 1
-timeStepDuration = 1
-
-seed = 1
-#'-q',
-arguments = c('-simDuration', simDuration, '-simBegin',simBegin,
-              '-timeStepDuration', timeStepDuration, 
-              '-RNGStatusIndex', format(seed,scientific = FALSE),
-              '-groups',"data/input/northeastamerica/fishRIOBasin_Sapidissima_Rjava.xml",
-              '-env',"data/input/northeastamerica/RIOBNneaBasins_Rjava.xml",
-              '-observers',"data/input/northeastamerica/RIO_obs_empty.xml")
-
-
-
+#' jarfile = "target/GR3D-3.2-SNAPSHOT.jar"
+#' simDuration = 150
+#' simBegin = 1
+#' timeStepDuration = 1
+#' 
+#' seed = 1
+#' #'-q',
+#' arguments = c('-simDuration', simDuration, '-simBegin',simBegin,
+#'               '-timeStepDuration', timeStepDuration, 
+#'               '-RNGStatusIndex', format(seed,scientific = FALSE),
+#'               '-groups',"data/input/northeastamerica/fishRIOBasin_Sapidissima_Rjava.xml",
+#'               '-env',"data/input/northeastamerica/RIOBNneaBasins_Rjava.xml",
+#'               '-observers',"data/input/northeastamerica/RIO_obs_empty.xml")
 
 
 # -----------------------------------------------------
 #  initializes the Java Virtual Machine -------------------------------------
+jarfile = "target/GR3D-3.2-SNAPSHOT.jar"
 .jinit(NULL, force.init = TRUE)
 .jinit(classpath = jarfile, force.init = TRUE)
 .jclassPath()
 
 #  modification of xlm parameters -----------------------
-parametersNamesANG = c("processes.processesEachStep.11.Soffset") 
-thetasANG = c(-1)
-
-parametersNamesENV = c("simulationName")
-thetasENV = c("Allee")
+# parametersNamesANG = c("processes.processesEachStep.11.Soffset") 
+# thetasANG = c(-1)
+# 
+# parametersNamesENV = c("simulationName")
+# thetasENV = c("Allee")
 
 # run the GR3D with updated parameter values in processes --------------------------------------------
-runEasyRun(arguments, parametersNamesANG,thetasANG, parametersNamesENV, thetasENV )
+# runEasyRun(arguments, parametersNamesANG,thetasANG, parametersNamesENV, thetasENV )
 
 # .jcall("analysis.EasyRun", "V", "runSimulation", arguments, 
 #        .jarray(parametersNamesANG), .jarray(thetasANG),
@@ -122,42 +136,60 @@ runEasyRun(arguments, parametersNamesANG,thetasANG, parametersNamesENV, thetasEN
 
 
 # call GR3D method to get model outputs ----------------------------------------------
-res1 = exploreEasyRun()
-res1$spawners_features
 
 
-#
-listeSimulation <- list(list(replicat = 1, 
-                             parametersNamesANG = c("processes.processesEachStep.11.Soffset"), 
-                             thetasANG = c(-1),
-                             parametersNamesENV = c("simulationName"), 
-                             thetasENV = c("Allee")),
-                        list(replicat = 1, 
-                             parametersNamesANG = c("processes.processesEachStep.11.Soffset"), 
-                             thetasANG = c(0),
-                             parametersNamesENV = c("simulationName"), 
-                             thetasENV = c("no Allee")))
-
-for (i in 1:length(listeSimulation)) {
-  arguments = c('-simDuration', simDuration, '-simBegin',simBegin,
-                '-timeStepDuration', timeStepDuration, 
-                '-RNGStatusIndex', format(listeSimulation[[i]]$replicat,scientific = FALSE), 
+runSimulation = function(simulationParameter) {
+  arguments = c('-simDuration', 600, '-simBegin',2,
+                '-timeStepDuration', 1, 
+                '-RNGStatusIndex', format(simulationParameter$replicat,scientific = FALSE), 
                 '-groups',"data/input/northeastamerica/fishRIOBasin_Sapidissima_Rjava.xml",
                 '-env',"data/input/northeastamerica/RIOBNneaBasins_Rjava.xml",
                 '-observers',"data/input/northeastamerica/RIO_obs_empty.xml")
   
-  parametersNamesANG = listeSimulation[[i]]$parametersNamesANG
-  thetasANG = listeSimulation[[i]]$thetasANG
-  parametersNamesENV = listeSimulation[[i]]$parametersNamesENV
-  thetasENV = listeSimulation[[i]]$thetasENV
-  
   runEasyRun(arguments, 
-             parametersNamesANG,
-             thetasANG, 
-             parametersNamesENV , 
-             thetasENV )
-  res1 = exploreEasyRun()
-  cat(res1$effective_centroid_latitude, "\n")
+             parametersNamesANG = simulationParameter$parametersNamesANG,
+             thetasANG = simulationParameter$thetasANG, 
+             parametersNamesENV = simulationParameter$parametersNamesENV, 
+             thetasENV = simulationParameter$thetasENV )
+  
+  #tic()
+  result = exploreEasyRun()
+  #toc()
+  return(result)
 }
 
-toto = listeSimulation[[i]]$parametersNamesANG
+
+#============================
+listeSimulation <- list(list(replicat = 1, 
+                             parametersNamesANG = c("processes.processesEachStep.11.Soffset",
+                                                    "processes.processesEachStep.6.pHomingAfterEquil"), 
+                             thetasANG = c(-1, 0.97),
+                             parametersNamesENV = c("simulationName"), 
+                             thetasENV = c("Allee_highHoming")),
+                        
+                        list(replicat = 1, 
+                             parametersNamesANG = c("processes.processesEachStep.11.Soffset",
+                                                    "processes.processesEachStep.6.pHomingAfterEquil"), 
+                             thetasANG = c(-1, 0.75),
+                             parametersNamesENV = c("simulationName"), 
+                             thetasENV = c("Allee_lowHoming")),
+                        
+                        list(replicat = 1, 
+                             parametersNamesANG = c("processes.processesEachStep.11.Soffset",
+                                                    "processes.processesEachStep.6.pHomingAfterEquil"), 
+                             thetasANG = c(0, 0.97),
+                             parametersNamesENV = c("simulationName"), 
+                             thetasENV = c("noAllee_highHoming")),
+                        
+                        list(replicat = 1, 
+                             parametersNamesANG = c("processes.processesEachStep.11.Soffset",
+                                                    "processes.processesEachStep.6.pHomingAfterEquil"), 
+                             thetasANG = c(0, 0.75),
+                             parametersNamesENV = c("simulationName"), 
+                             thetasENV = c("noAllee_lowHoming"))
+                        )
+tic()
+resSimulation = lapply(listeSimulation, runSimulation)
+toc()
+
+save(listeSimulation, resSimulation, file = 'simulation.Rdata')
diff --git a/src/main/java/analysis/AnalyseFishDistribution.java b/src/main/java/analysis/AnalyseFishDistribution.java
index a1839c52104b209d7d68795514e6f1c1cdf4af3b..dc037c1f0420f50c3449dc537b03ce9e3401ca60 100644
--- a/src/main/java/analysis/AnalyseFishDistribution.java
+++ b/src/main/java/analysis/AnalyseFishDistribution.java
@@ -254,40 +254,34 @@ public class AnalyseFishDistribution extends AquaNismsGroupProcess<DiadromousFis
 		// initialise the distribution range
 		double southernEdge = northernLimit;
 		double northernEdge = southernLimit;
-		int southernId = -1;
-		int northernId = -1;
-		String southernName = "";
-		String northernName = "";
-		// String southernName = "";
-		// String northernName = "";
+
+		RiverBasin northernBasin = null;
+		RiverBasin southernBasin = null;
+
 		// initialise for centroid
 		TreeMapForCentile latitudeEffective = new TreeMapForCentile();
 		TreeMapForCentile latitudePresence = new TreeMapForCentile();
 
-		String[][] export = new String[2][6];
-
 		for (RiverBasin riverBasin : riverBasins) {
 			double meanLastRecruitment = riverBasin.getLastRecruitments().getMean();
 			if (meanLastRecruitment >= minimumRecruitsForPopulatedBasin) {
 				// the river basin is considered populated
 
-				// NOTE : recruiit number is calulated after thermal tolerance impact
+				// CHECK : recruiit number is calulated after thermal tolerance impact
 				// (which is intreaged in the stock-recruitment relationship)
 
 				// southern edge
 				if (riverBasin.getLatitude() < southernEdge) {
 					// the basin is the new southern edge
 					southernEdge = riverBasin.getLatitude();
-					southernId = riverBasin.getBasin_id();
-					southernName = riverBasin.getName();
+					southernBasin = riverBasin;
 				}
 
 				// northern edge
 				if (riverBasin.getLatitude() > northernEdge) {
 					// the basin is the new northern edge
 					northernEdge = riverBasin.getLatitude();
-					northernId = riverBasin.getBasin_id();
-					northernName = riverBasin.getName();
+					northernBasin = riverBasin;
 				}
 
 				// for distribution centroide computation
@@ -297,29 +291,41 @@ public class AnalyseFishDistribution extends AquaNismsGroupProcess<DiadromousFis
 			}
 		}
 
-		// if the universe is empty (none of basins is populated)
-		if (southernEdge == northernLimit & northernEdge == southernLimit) {
-			southernEdge = (northernLimit + southernLimit) / 2.;
-			northernEdge = southernEdge;
-		}
-
 		// distribution centroids computation
 		double basinCentroid = latitudePresence.calculateMedian();
 		double effectiveCentroid = latitudeEffective.calculateMedian();
 
+		// export
+		String[][] export = new String[2][8];
 		export[0][0] = "southern_basin_id";
 		export[0][1] = "southern_basin_name";
-		export[0][2] = "northern_basin_id";
-		export[0][3] = "northern_basin_name";
-		export[0][4] = "basin_centroid_latitude";
-		export[0][5] = "effective_centroid_latitude";
-
-		export[1][0] = String.valueOf(southernId);
-		export[1][1] = southernName;
-		export[1][2] = String.valueOf(northernId);
-		export[1][3] = northernName;
-		export[1][4] = String.valueOf(basinCentroid);
-		export[1][5] = String.valueOf(effectiveCentroid);
+		export[0][2] = "southern_latitude";
+		export[0][3] = "northern_basin_id";
+		export[0][4] = "northern_basin_name";
+		export[0][5] = "northern_latitude";
+		export[0][6] = "basin_centroid_latitude";
+		export[0][7] = "effective_centroid_latitude";
+
+		if (southernBasin != null & northernBasin != null) {
+			export[1][0] = String.valueOf(southernBasin.getBasin_id());
+			export[1][1] = southernBasin.getName();
+			export[1][2] = String.valueOf(southernBasin.getLatitude());
+			export[1][3] = String.valueOf(northernBasin.getBasin_id());
+			export[1][4] = northernBasin.getName();
+			export[1][5] = String.valueOf(northernBasin.getLatitude());
+			export[1][6] = String.valueOf(basinCentroid);
+			export[1][7] = String.valueOf(effectiveCentroid);
+		} else {
+			export[1][0] = "NA";
+			export[1][1] = "NA";
+			export[1][2] = "NA";
+			export[1][3] = "NA";
+			export[1][4] = "NA";
+			export[1][5] = "NA";
+			export[1][6] = "NA";
+			export[1][7] = "NA";
+		}
+
 		return export;
 	}
 
diff --git a/src/main/java/analysis/AnalyseSpawnerFeatures.java b/src/main/java/analysis/AnalyseSpawnerFeatures.java
index 79faa13861041368eda95d03c731c40caf842ebe..50c7bc17f66b4dc7f7865a91ac1423253476e234 100644
--- a/src/main/java/analysis/AnalyseSpawnerFeatures.java
+++ b/src/main/java/analysis/AnalyseSpawnerFeatures.java
@@ -328,7 +328,7 @@ public class AnalyseSpawnerFeatures extends AquaNismsGroupProcess<DiadromousFish
 		result[0][1] = "basin_name";
 		result[0][2] = "latitude";
 		result[0][3] = "mean_age_female"; // primiparous and multiparous
-		result[0][4] = "mean_age _male"; // primiparous and multiparous
+		result[0][4] = "mean_age_male"; // primiparous and multiparous
 		result[0][5] = "pct_primiparous";
 
 		int i = 1;
diff --git a/src/main/java/analysis/EasyRun.java b/src/main/java/analysis/EasyRun.java
index d528cd319699cc746bb77e276cf80ef4cd2e636b..e4bbbf13fa9b89581e6acbc42fed70baeae6f3f7 100644
--- a/src/main/java/analysis/EasyRun.java
+++ b/src/main/java/analysis/EasyRun.java
@@ -55,7 +55,6 @@ public class EasyRun {
 	public static void runSimulation(String[] batchArgs, String[] paramNameaANG, double[] paramValuesANG, String[] paramNameaENV,
 			String[] paramValuesENV) throws Exception {
 
-		// TODO parse string to other type if needed
 		try {
 			pilot = new Pilot();
 			BatchRunner runner = new BatchRunner(pilot);
@@ -129,21 +128,29 @@ public class EasyRun {
 		// double[] parameterValues = new double[] { 10, 2, 0.7 };
 
 		// =====================================================
-		System.out.println("WITHOUT ALLEE EFFECT");
-		String[] parameterNamesANG = new String[] { "processes.processesEachStep.11.Soffset" };
-		double[] parameterValuesANG = new double[] { 0. };
+		// System.out.println("WITHOUT ALLEE EFFECT");
+		// String[] parameterNamesANG = new String[] { "processes.processesEachStep.11.Soffset" };
+		// double[] parameterValuesANG = new double[] { 0. };
+		// // runSimulation(batchArguments, parameterNamesANG, parameterValuesANG);
+		// String[] parameterNamesENV = new String[] { "simulationName" };
+		// String[] parameterValuesENV = new String[] { "noAllee" };
+
+		System.out.println("reference");
+		String[] parameterNamesANG = new String[0];
+		double[] parameterValuesANG = new double[0];
 		// runSimulation(batchArguments, parameterNamesANG, parameterValuesANG);
-		String[] parameterNamesENV = new String[] { "simulationName" };
-		String[] parameterValuesENV = new String[] { "noAllee" };
+		String[] parameterNamesENV = new String[0];
+		String[] parameterValuesENV = new String[0];
+
 		runSimulation(batchArguments, parameterNamesANG, parameterValuesANG, parameterNamesENV, parameterValuesENV);
 
-		System.out.println("\n== AnalyseSpawnerFeatures ==");
-		String[][] spawnerRunResults = getValuesFromAquanismGroupProcess("processes.processesEachStep.9.exportToR");
-		for (String[] record : spawnerRunResults) {
-			for (String value : record)
-				System.out.print(value + "\t");
-			System.out.println();
-		}
+		// System.out.println("\n== AnalyseSpawnerFeatures ==");
+		// String[][] spawnerRunResults = getValuesFromAquanismGroupProcess("processes.processesEachStep.9.exportToR");
+		// for (String[] record : spawnerRunResults) {
+		// for (String value : record)
+		// System.out.print(value + "\t");
+		// System.out.println();
+		// }
 
 		System.out.println("\n== AnalyseFishDistribution ==");
 		String[][] distributionResults = getValuesFromAquanismGroupProcess("processes.processesEachStep.10.exportToR");
@@ -156,13 +163,14 @@ public class EasyRun {
 		System.out.println("\n== AnalyseLikelihoodOfPresence ==");
 		System.out.println(getSingleValueFromAquanismGroupProcess("processes.processesAtEnd.0.getLogLikelihood"));
 
-		System.out.println("\n== ExportExchangeSBetweenCatchments ==");
-		String[][] exchangeSBetweenCatchments = getValuesFromAquanismGroupProcess("processes.processesEachStep.16.exportToR");
-		for (String[] record : exchangeSBetweenCatchments) {
-			for (String value : record)
-				System.out.print(value + "\t");
-			System.out.println();
-		}
+		// System.out.println("\n== ExportExchangeSBetweenCatchments ==");
+		// String[][] exchangeSBetweenCatchments =
+		// getValuesFromAquanismGroupProcess("processes.processesEachStep.16.exportToR");
+		// for (String[] record : exchangeSBetweenCatchments) {
+		// for (String value : record)
+		// System.out.print(value + "\t");
+		// System.out.println();
+		// }
 
 	}
 }
diff --git a/src/main/java/analysis/ExportExchangesBetweenCatchments.java b/src/main/java/analysis/ExportExchangesBetweenCatchments.java
index 88e894f00e53fec00102d0709510b98497138dd0..67efe43a1e8dac913381060e0ed90d32b7ab9749 100644
--- a/src/main/java/analysis/ExportExchangesBetweenCatchments.java
+++ b/src/main/java/analysis/ExportExchangesBetweenCatchments.java
@@ -241,22 +241,17 @@ public class ExportExchangesBetweenCatchments extends AquaNismsGroupProcess<Diad
 	public String[][] exportToR() {
 
 		int i, j;
-		String[][] export = new String[riverBasins.length + 1][riverBasins.length + 1];
+		String[][] export = new String[riverBasins.length * riverBasins.length + 1][3];
 
 		// headers
 		export[0][0] = "migrationBasin";
-		j = 1;
-		for (RiverBasin natalBasin : riverBasins) {
-			export[0][j] = natalBasin.getName();
-			j++;
-		}
-
+		export[0][1] = "natalBasin";
+		export[0][2] = "effective";
 		i = 1;
 		for (RiverBasin migrationBasin : riverBasins) {
-			export[i][0] = migrationBasin.getName();
-			j = 1;
-			for (RiverBasin natalBasin : riverBasins) {
 
+			for (RiverBasin natalBasin : riverBasins) {
+				// compute the mean over years
 				Map<Long, Double> year_effective = exchanges.get(migrationBasin).get(natalBasin);
 				int n = 0;
 				double sum = 0;
@@ -266,10 +261,12 @@ public class ExportExchangesBetweenCatchments extends AquaNismsGroupProcess<Diad
 						sum += effective;
 					}
 				}
-				export[i][j] = String.valueOf(Math.round(10. * sum / n) / 10);
-				j++;
+				// fill in export array
+				export[i][0] = migrationBasin.getName();
+				export[i][1] = natalBasin.getName();
+				export[i][2] = String.valueOf(Math.round(10. * sum / n) / 10.);
+				i++;
 			}
-			i++;
 		}
 
 		return export;
diff --git a/src/main/java/environment/BasinNetworkObserverPresence.java b/src/main/java/environment/BasinNetworkObserverPresence.java
index 1246c13c568b2a760e1f942ca7aecca0b545990c..4b27f34f7c6af45f96c730d599a8d6d734269111 100644
--- a/src/main/java/environment/BasinNetworkObserverPresence.java
+++ b/src/main/java/environment/BasinNetworkObserverPresence.java
@@ -37,8 +37,7 @@ import fr.cemagref.simaqualife.kernel.util.TransientParameters;
 import fr.cemagref.simaqualife.pilot.Pilot;
 
 @SuppressWarnings("serial")
-public class BasinNetworkObserverPresence extends ObserverListener
-		implements Configurable, Drawable, MouseMotionListener {
+public class BasinNetworkObserverPresence extends ObserverListener implements Configurable, Drawable, MouseMotionListener {
 
 	private String title;
 
@@ -46,6 +45,8 @@ public class BasinNetworkObserverPresence extends ObserverListener
 
 	private String period = "obs_1751_1850";
 
+	private Boolean fillInshoreBasin = false;
+
 	@Description(name = "Color scale", tooltip = "")
 	public ColorScaleEnum colorScaleEnum = ColorScaleEnum.BluesScale;
 
@@ -73,11 +74,13 @@ public class BasinNetworkObserverPresence extends ObserverListener
 	// basin under the mouse
 	private transient Basin basinUnderMouse;
 
-
 	@TransientParameters.InitTransientParameters
 	public void init(Pilot pilot) {
 		this.pilot = pilot;
 
+		if (fillInshoreBasin == null)
+			fillInshoreBasin = false;
+
 		if (this.colorScaleEnum == null) {
 			this.colorScaleEnum = ColorScaleEnum.RedsScale;
 
@@ -113,13 +116,13 @@ public class BasinNetworkObserverPresence extends ObserverListener
 						}
 					}
 				}
-				// reader.close();
-				// scanner.close();
+				reader.close();
+				scanner.close();
 			} catch (Exception e) {
 				e.printStackTrace();
 			}
-
 		}
+
 		// the Jpanal that holds all the components to be displayed
 		display = new JPanel(new BorderLayout());
 
@@ -142,7 +145,6 @@ public class BasinNetworkObserverPresence extends ObserverListener
 
 	}
 
-
 	/*
 	 * @Override public void valueChanged(ObservablesHandler arg0, Object arg1, long arg2) {
 	 * 
@@ -152,6 +154,7 @@ public class BasinNetworkObserverPresence extends ObserverListener
 	 * label.setText(txt); display.repaint(); }
 	 */
 
+
 	@Override
 	public void mouseDragged(MouseEvent arg0) {
 		// nothing to do
@@ -327,12 +330,21 @@ public class BasinNetworkObserverPresence extends ObserverListener
 				 * else g.setColor(Color.LIGHT_GRAY);
 				 */
 
-				if (presence > 0)
-					g.setColor(colorScaleEnum.getScale().getColor(presence));
-				else if (presence == 0)
-					g.setColor(Color.WHITE);
-				else
-					g.setColor(Color.GRAY);
+				if (basin instanceof RiverBasin) {
+					if (presence > 0)
+						g.setColor(colorScaleEnum.getScale().getColor(presence));
+					else if (presence == 0)
+						g.setColor(Color.WHITE);
+					else
+						g.setColor(Color.GRAY);
+				} else if (basin instanceof InshoreBasin & fillInshoreBasin) {
+					if (presence > 0)
+						g.setColor(colorScaleEnum.getScale().getColor(presence));
+					else if (presence == 0)
+						g.setColor(Color.WHITE);
+					else
+						g.setColor(Color.GRAY);
+				}
 
 				g2d.fill(displayShape);
 			}
@@ -341,10 +353,9 @@ public class BasinNetworkObserverPresence extends ObserverListener
 
 	public enum ColorScaleEnum {
 
-		RedsScale(new RedsScale()), BluesScale(new BluesScale()), BicolorScale(new BicolorScale()), GraysScale(
-				new GraysScale());
-		private ColorScale scale;
+		RedsScale(new RedsScale()), BluesScale(new BluesScale()), BicolorScale(new BicolorScale()), GraysScale(new GraysScale());
 
+		private ColorScale scale;
 
 		ColorScaleEnum(ColorScale colorScale) {
 			scale = colorScale;
diff --git a/src/main/java/environment/RIOBasinNetworkObserverPresence.java b/src/main/java/environment/RIOBasinNetworkObserverPresence.java
index f30f81cba30d1c80005eac5badf329f6a91bb210..a51b111c72a4b4f8ce74e7747e96a56696b8e158 100644
--- a/src/main/java/environment/RIOBasinNetworkObserverPresence.java
+++ b/src/main/java/environment/RIOBasinNetworkObserverPresence.java
@@ -45,6 +45,8 @@ public class RIOBasinNetworkObserverPresence extends ObserverListener implements
 
 	private String period = "obs_1751_1850";
 
+	private Boolean fillInshoreBasin = false;
+
 	@Description(name = "Color scale", tooltip = "")
 	public ColorScaleEnum colorScaleEnum = ColorScaleEnum.BluesScale;
 
@@ -72,11 +74,13 @@ public class RIOBasinNetworkObserverPresence extends ObserverListener implements
 	// basin under the mouse
 	private transient Basin basinUnderMouse;
 
-
 	@TransientParameters.InitTransientParameters
 	public void init(Pilot pilot) {
 		this.pilot = pilot;
 
+		if (fillInshoreBasin == null)
+			fillInshoreBasin = false;
+
 		if (this.colorScaleEnum == null)
 			this.colorScaleEnum = ColorScaleEnum.RedsScale;
 
@@ -139,7 +143,6 @@ public class RIOBasinNetworkObserverPresence extends ObserverListener implements
 		bn = (RIOBasinNetworkWithContinent) pilot.getAquaticWorld().getEnvironment();
 	}
 
-
 	/*
 	 * @Override public void valueChanged(ObservablesHandler arg0, Object arg1, long arg2) {
 	 * 
@@ -149,6 +152,7 @@ public class RIOBasinNetworkObserverPresence extends ObserverListener implements
 	 * label.setText(txt); display.repaint(); }
 	 */
 
+
 	@Override
 	public void mouseDragged(MouseEvent arg0) {
 		// nothing to do
@@ -318,12 +322,22 @@ public class RIOBasinNetworkObserverPresence extends ObserverListener implements
 					else
 						System.out.println(basinName + " is not present in the presence file");
 
-					if (presence > 0)
-						g.setColor(colorScaleEnum.getScale().getColor(presence));
-					else if (presence == 0)
-						g.setColor(Color.WHITE);
-					else
-						g.setColor(Color.GRAY);
+					if (basin instanceof RiverBasin) {
+						if (presence > 0)
+							g.setColor(colorScaleEnum.getScale().getColor(presence));
+						else if (presence == 0)
+							g.setColor(Color.WHITE);
+						else
+							g.setColor(Color.GRAY);
+					} else if (basin instanceof InshoreBasin & fillInshoreBasin) {
+						if (presence > 0)
+							g.setColor(colorScaleEnum.getScale().getColor(presence));
+						else if (presence == 0)
+							g.setColor(Color.WHITE);
+						else
+							g.setColor(Color.GRAY);
+					} else
+						g.setColor(Color.BLUE);
 
 					g2d.fill(displayShape);
 				}
@@ -334,8 +348,8 @@ public class RIOBasinNetworkObserverPresence extends ObserverListener implements
 	public enum ColorScaleEnum {
 
 		RedsScale(new RedsScale()), BluesScale(new BluesScale()), BicolorScale(new BicolorScale()), GraysScale(new GraysScale());
-		private ColorScale scale;
 
+		private ColorScale scale;
 
 		ColorScaleEnum(ColorScale colorScale) {
 			scale = colorScale;
diff --git a/src/main/java/species/ReproduceWithDiagnose.java b/src/main/java/species/ReproduceWithDiagnose.java
index e45c074cfc423f068750531bf6244b48ef6009ee..0eb64bb656f71299a894b0ba7c05dfbf174dff72 100644
--- a/src/main/java/species/ReproduceWithDiagnose.java
+++ b/src/main/java/species/ReproduceWithDiagnose.java
@@ -25,7 +25,6 @@ import environment.Time.Season;
 import fr.cemagref.simaqualife.kernel.processes.AquaNismsGroupProcess;
 import fr.cemagref.simaqualife.kernel.util.TransientParameters.InitTransientParameters;
 import fr.cemagref.simaqualife.pilot.Pilot;
-import miscellaneous.BinomialForSuperIndividualGen;
 import miscellaneous.Miscellaneous;
 import miscellaneous.Trio;
 import species.DiadromousFish.Gender;
@@ -83,12 +82,12 @@ public class ReproduceWithDiagnose extends AquaNismsGroupProcess<DiadromousFish,
 
 	private transient NormalGen genNormal;
 
-	/**
-	 * the random numbers generator for binomial draws
-	 * 
-	 * @unit --
-	 */
-	private transient BinomialForSuperIndividualGen aleaGen;
+	// /**
+	// * the random numbers generator for binomial draws
+	// *
+	// * @unit --
+	// */
+	// private transient BinomialForSuperIndividualGen aleaGen;
 
 	private transient MortalityFunction mortalityFunction;
 
@@ -111,7 +110,7 @@ public class ReproduceWithDiagnose extends AquaNismsGroupProcess<DiadromousFish,
 		super.initTransientParameters(pilot);
 		genNormal = new NormalGen(pilot.getRandomStream(), new NormalDist(0., 1.));
 
-		aleaGen = new BinomialForSuperIndividualGen(pilot.getRandomStream());
+		// aleaGen = new BinomialForSuperIndividualGen(pilot.getRandomStream());
 
 		mortalityFunction = new MortalityFunction();
 		stockRecruitmentRelationship = new StockRecruitmentRelationship();
@@ -205,11 +204,11 @@ public class ReproduceWithDiagnose extends AquaNismsGroupProcess<DiadromousFish,
 					amountPerSuperIndividual = alpha / maxNumberOfSuperIndividualPerReproduction;
 
 					// Compute the Allee effect parameters S95 and S50
-					if (Soffset >= 0.) {// Allee effect independant of catchment size (including
+					if (Soffset >= 0.) {// Allee effect independant of catchment size (including no Allee Effect)
 						S95 = Soffset;
 						S50 = Soffset;
 					} else {
-						S95 = eta * riverBasin.getAccessibleSurface(); // corresponds to S* in the rougier publication
+						S95 = eta * riverBasin.getAccessibleSurface(); // corresponds to S* in the rougier etal 2015
 						S50 = S95 / ratioS95_S50;
 					}
 					// initilisation of the stock recruitment relationship
@@ -535,7 +534,7 @@ class StockRecruitmentRelationship implements UnivariateFunction {
 		// BH Stock-Recruitment relationship with logistic depensation
 		double meanNumberOfRecruit = 0.;
 		double effectiveStock = getEffectiveStock(stock);
-		if (stock > 0)
+		if (effectiveStock > 0)
 			meanNumberOfRecruit = Math.round(alpha * effectiveStock) / (beta + effectiveStock);
 		return meanNumberOfRecruit;
 	}