 MODULE contareaddat
#if defined key_substance && defined key_contaminant
   !!======================================================================
   !!                   ***  MODULE  contareaddat  ***
   !!              Ocean dynamics Bio Sed:  Read files *.dat 
   !!                  (chemical variables and reactions)
   !!======================================================================

#include "toolcpp.h"

   !! * Modules used
   use comsubstance
   USE parameters, ONLY : rsh,rlg,lchain,imin,imax,jmin,jmax,kmax,ljmin,ljmax,limin,limax  
   USE comvars2d, ONLY : iscreenlog        
#ifdef key_sedim                          
   USE parameters, ONLY : ksdmin,ksdmax    
#endif
   USE parameters, ONLY : ntypreac,nb_var,nflimcomax,nrecomax,nvactcomax   
   USE_MPI toolmpi, ONLY : MASTER

   IMPLICIT NONE
   PRIVATE

   !! * Accessibility
   PUBLIC conta_read_react,                & ! routine called by main.F90
          conta_read_data                   ! routine called by traceur.F90
          

   !! * Shared module variables
   LOGICAL,PUBLIC                                          :: l_sorfluxreac,l_flim3_re
   LOGICAL,PUBLIC                                          :: l_treat_re_after,l_treat_re_dc
   CHARACTER(len=lchain),PUBLIC                            :: filereaction
   REAL(KIND=rsh),PUBLIC                                   :: minpt_conta
   REAL(KIND=rsh),PUBLIC                                   :: minpt_mes_re5
!   REAL(KIND=rlg),PUBLIC                                   :: f_dtconta_sksc
   REAL(KIND=rlg),PUBLIC                                   :: dtconta_dc
!   REAL(KIND=rlg),PUBLIC                                   :: tconta_sksc
!   REAL(KIND=rlg),PUBLIC                                   :: tconta_sksc0
   INTEGER ,PUBLIC                                         :: nflimax,nreacmax
   INTEGER ,PUBLIC                                         :: nvaractmax
   INTEGER,dimension(:),allocatable,PUBLIC                 :: ityp_reac
   INTEGER,dimension(:),allocatable,PUBLIC                 :: nb_var_re
   INTEGER,dimension(ntypreac),PUBLIC                      :: nb_reac
   INTEGER,dimension(:,:),allocatable,PUBLIC               :: irkm_var_re
   INTEGER,dimension(:),allocatable,PUBLIC                 :: ityp_arh,ire_typ
   CHARACTER(len=lchain),PUBLIC                            :: fileconta1
   CHARACTER(len=lchain),PUBLIC                            :: fileconta2
   CHARACTER(len=lchain),PUBLIC                            :: fileconta3
   CHARACTER(len=lchain),PUBLIC                            :: fileconta4
   CHARACTER(len=lchain),PUBLIC                            :: fileconta5
   CHARACTER(len=lchain),dimension(:),allocatable,PUBLIC   :: name_reac
   CHARACTER(len=lchain),dimension(:,:),allocatable,PUBLIC :: name_var_re
   REAL(KIND=rsh),dimension(:),allocatable,PUBLIC          :: cin_max_re
   REAL(KIND=rsh),dimension(:),allocatable,PUBLIC          :: cf_teta_arh
   REAL(KIND=rsh),dimension(:),allocatable,PUBLIC          :: cf_xetsig_arh
   REAL(KIND=rsh),dimension(:),allocatable,PUBLIC          :: tempref_arh
   REAL(KIND=rsh),dimension(:,:),allocatable ,PUBLIC       :: cf_expovar_re
   REAL(KIND=rsh),dimension(:,:),allocatable ,PUBLIC       :: cf_stochio_re
   LOGICAL,dimension(:),allocatable ,PUBLIC                :: l_arh,l_treat_re
   LOGICAL,dimension(:),allocatable ,PUBLIC                :: l_cincte_re
   LOGICAL,dimension(:),allocatable ,PUBLIC                :: l_coefcte_re
   INTEGER,dimension(:),allocatable,PUBLIC                 :: nb_flim1_re
   INTEGER,dimension(:),allocatable,PUBLIC                 :: nb_flim2_re
   INTEGER,dimension(:),allocatable,PUBLIC                 :: nb_flim3_re
   INTEGER,dimension(:,:),allocatable,PUBLIC               :: irkm_flim1_re
   INTEGER,dimension(:,:),allocatable,PUBLIC               :: irkm_flim2_re
   REAL(KIND=rsh),dimension(:,:),allocatable,PUBLIC        :: cf_flim1_re
   REAL(KIND=rsh),dimension(:,:),allocatable,PUBLIC        :: cf_flim2_re
   REAL(KIND=rsh),dimension(:,:),allocatable,PUBLIC        :: cf_flim3_re
   LOGICAL,dimension(:),allocatable,PUBLIC                 :: l_pce_re5
   LOGICAL,dimension(:),allocatable,PUBLIC                 :: l_salt_re5
   INTEGER,dimension(:,:),allocatable ,PUBLIC              :: cf_pce_re5
   INTEGER,dimension(:,:),allocatable ,PUBLIC              :: cf_salt_re5
   INTEGER,dimension(:,:),allocatable ,PUBLIC              :: irkm_sub_re5
   INTEGER,dimension(:),allocatable,PUBLIC                 :: ndg_kds_re5
   LOGICAL,dimension(:),allocatable,PUBLIC                 :: l_kds_re5
   REAL(KIND=rsh),dimension(:,:),allocatable,PUBLIC        :: cf_distb_re5
   REAL(KIND=rsh),dimension(:),allocatable,PUBLIC          :: cf_kd0_re5
   REAL(KIND=rsh),dimension(:,:),allocatable,PUBLIC        :: cf_xkds_re5
   REAL(KIND=rsh),dimension(:),allocatable,PUBLIC          :: cf_rate_re3
   REAL(KIND=rsh),dimension(:),allocatable,PUBLIC          :: cv_watequ_re3
   REAL(KIND=rsh),dimension(:,:),allocatable,PUBLIC        :: cin_re2_tmp
 
   CHARACTER(len=lchain),dimension(:),allocatable,PUBLIC   :: name_conta
   CHARACTER(LEN=lchain),dimension(:,:),allocatable,PUBLIC :: name_spc
   INTEGER,PUBLIC                                          :: nb_conta,nb_spcmax
   INTEGER,dimension(:,:) ,allocatable,PUBLIC              :: nb_sorpt_nlin
   INTEGER,dimension(:,:) ,allocatable,PUBLIC              :: ispc_assoc_pcin
   INTEGER,dimension(:,:) ,allocatable,PUBLIC              :: ivar_spc
   INTEGER,dimension(:,:) ,allocatable,PUBLIC              :: isit_spc
   INTEGER,dimension(:),allocatable,PUBLIC                 :: ivdl_conta
   INTEGER,dimension(:),allocatable,PUBLIC                 :: nb_spc_diss
   INTEGER,dimension(:),allocatable,PUBLIC                 :: nb_spc_part_eq
   INTEGER,dimension(:),allocatable,PUBLIC                 :: nb_spc_part_ci
   INTEGER,dimension(:),allocatable,PUBLIC                 :: nb_spc_part_ne
   INTEGER,dimension(:),allocatable,PUBLIC                 :: nb_spc_part_tot
   INTEGER,dimension(:),allocatable,PUBLIC                 :: nb_spc_tot
   REAL(KIND=rsh),dimension(nb_var) ,PUBLIC                :: massmol,cf_henry
   REAL(KIND=rsh),dimension(nb_var), PUBLIC                 :: nschmidt
   REAL(KIND=rsh),dimension(nb_var) ,PUBLIC                :: cf_difair
   REAL(KIND=rsh),dimension(nb_var),PUBLIC                 :: cf_henry_adim
   REAL(KIND=rsh),dimension(nb_var),PUBLIC                 :: cf_temp_henry
   REAL(KIND=rsh),dimension(nb_var),PUBLIC                 :: cf_salt_henry
   REAL(KIND=rsh),dimension(nb_var),PUBLIC                 :: cf_difwat
   REAL(KIND=rsh),dimension(:,:,:),allocatable,PUBLIC      :: f_sorpt_nlin
   REAL(KIND=rsh),dimension(:,:,:),allocatable,PUBLIC      :: cf_kfkl_nlin
   REAL(KIND=rsh),dimension(:,:,:),allocatable,PUBLIC      :: cf_expqmax_nlin
   LOGICAL,PUBLIC                                          :: l_conta
   LOGICAL,dimension(:),allocatable,PUBLIC                 :: l_sorpt_nlin
   LOGICAL,dimension(:),allocatable,PUBLIC                 :: l_langmuir
   LOGICAL,dimension(:),allocatable,PUBLIC                 :: l_organic,l_metal

! indices pour reperage des especes impliquees dans les réactions 
   INTEGER,dimension(:,:),allocatable,PUBLIC               :: ispc_reac
   INTEGER,dimension(:),allocatable,PUBLIC                 :: icon_reac

! indices pour reperage des variables contenant plusieurs especes de contaminant
   INTEGER,dimension(:),allocatable,PUBLIC                 :: icon_var
#ifdef key_species
   INTEGER,dimension(:),allocatable,PUBLIC                 :: nb_frac_var
   LOGICAL,PUBLIC                                          :: l_gpspecies
   LOGICAL,PUBLIC                                          :: l_sorspecies
   INTEGER,dimension(:,:) ,allocatable,PUBLIC              :: ispc_var
   INTEGER,dimension(:,:) ,allocatable,PUBLIC              :: ifrac_var_reac
! fractions de chaque espece de contaminant dans les concentrations 
!                  totales de chaque variable d etat 
   REAL(KIND=rsh),dimension(:,:,:,:,:),allocatable,PUBLIC  :: frac_cwat
   REAL(KIND=rsh),dimension(:,:,:,:,:),allocatable,PUBLIC  :: frac_csed
   REAL(KIND=rsh),dimension(:,:,:,:,:),allocatable,PUBLIC  :: dc_spc_wat
   INTEGER,dimension(:,:),allocatable,PUBLIC               :: ifrac_var_spc
#else
   INTEGER,dimension(:),allocatable,PUBLIC                 :: ispc_var
#endif

! volatilisation
   LOGICAL,dimension(:) ,allocatable,PUBLIC     :: l_volatil
   LOGICAL,dimension(nb_var),PUBLIC             :: l_volkl_dif,l_volkl_sal
   LOGICAL,dimension(nb_var),PUBLIC             :: l_volkl_cor,l_volkg_cor
   INTEGER,dimension(nb_var),PUBLIC             :: ich_volkl_wind
   REAL(KIND=rsh),dimension(nb_var),PUBLIC      :: cf_expvol_liq,cf_expvol_gas
   REAL(KIND=rsh),PUBLIC                        :: massmol_h2o,massmol_co2
   REAL(KIND=rsh),PUBLIC                        :: massmol_o2,cf_difwat_co2
   REAL(KIND=rsh),PUBLIC                        :: cf_difwat_o2,cf_difair_h2o
   REAL(KIND=rsh),PUBLIC                        :: cst_gas_univ,thick_volat

! fractions at boundaries for each species
!   REAL(KIND=rsh),dimension(:,:),allocatable,PUBLIC        :: flim_n_spc
!   REAL(KIND=rsh),dimension(:,:),allocatable,PUBLIC        :: flim_e_spc
!   REAL(KIND=rsh),dimension(:,:),allocatable,PUBLIC        :: flim_s_spc
!   REAL(KIND=rsh),dimension(:,:),allocatable,PUBLIC        :: flim_w_spc
! flux reactifs entre variables
   REAL(KIND=rsh),dimension(:,:,:,:,:),allocatable, PUBLIC :: flux_reac_wat


   !! * Locales variables of module which will be deALLOCATE after  

   INTEGER                                     :: eof,nreacmax_sco
   INTEGER,DIMENSION(nb_var)                   :: nb_reac_pervar
   INTEGER,dimension(:,:),ALLOCATABLE          :: irkm_sub_r,irkf_var_re
   INTEGER,dimension(:) ,ALLOCATABLE           :: nb_reac_genco
   INTEGER,dimension(:,:),ALLOCATABLE          :: nb_flim1_co,nb_flim2_co
   INTEGER,dimension(:,:),ALLOCATABLE          :: nb_flim3_co
   INTEGER,dimension(:,:),ALLOCATABLE          :: ityp_arh_rco,nb_var_co
   INTEGER,dimension(:) ,ALLOCATABLE           :: ndg_kds_r
   INTEGER,dimension(:),ALLOCATABLE            :: ndg_kds_metal
   INTEGER,dimension(:),ALLOCATABLE            :: ich_volkl_windco
   INTEGER,dimension(:),allocatable            :: ityp_arh_pco
   REAL(KIND=rsh),dimension(:,:),allocatable   :: cin_hydrol,cin_photo
   REAL(KIND=rsh),dimension(:) ,allocatable    :: cf_kow_corg,solubil_co
   REAL(KIND=rsh),dimension(:,:),allocatable   :: cin_biodeg
   REAL(KIND=rsh),dimension(:),allocatable     :: tempref_arh_pco
   REAL(KIND=rsh),dimension(:),allocatable     :: cf_teta_arh_pco
   REAL(KIND=rsh),dimension(:),allocatable     :: cf_xetsig_arh_pco
   REAL(KIND=rsh),dimension(:,:),allocatable   :: cf_distb_spc
   REAL(KIND=rsh),dimension(:,:),allocatable   :: cf_rate_spc
   REAL(KIND=rsh),dimension(:,:),allocatable   :: cv_watequ_spc
   REAL(KIND=rsh),dimension(:),ALLOCATABLE     :: cf_kd0_metal
   REAL(KIND=rsh),dimension(:),ALLOCATABLE     :: cv_watequ_r
   REAL(KIND=rsh),dimension(:,:),allocatable   :: cf_pce_co,cf_salt_co
   REAL(KIND=rsh),dimension(:,:),allocatable   :: cf_pce_r,cf_salt_r
   REAL(KIND=rsh),dimension(:,:),ALLOCATABLE   :: cf_distb_r
   REAL(KIND=rsh),dimension(:),allocatable     :: cf_kd0_r,cf_rate_r
   REAL(KIND=rsh),dimension(:,:),allocatable   :: cf_xkds_r
   REAL(KIND=rsh),dimension(:),ALLOCATABLE     :: massmol_co,nschmidt_co
   REAL(KIND=rsh),dimension(:),ALLOCATABLE     :: cf_difair_co
   REAL(KIND=rsh),dimension(:),ALLOCATABLE     :: cf_difwat_co
   REAL(KIND=rsh),dimension(:),ALLOCATABLE     :: cf_temp_henry_co
   REAL(KIND=rsh),dimension(:),ALLOCATABLE     :: cf_salt_henry_co
   REAL(KIND=rsh),dimension(:),ALLOCATABLE     :: cf_henry_co
   REAL(KIND=rsh),dimension(:),ALLOCATABLE     :: cf_henry_adim_co
   REAL(KIND=rsh),dimension(:,:),ALLOCATABLE   :: cin_max_co
   REAL(KIND=rsh),dimension(:,:),ALLOCATABLE   :: cf_teta_arh_rco
   REAL(KIND=rsh),dimension(:,:),ALLOCATABLE   :: tempref_arh_rco
   REAL(KIND=rsh),dimension(:,:),ALLOCATABLE   :: cf_xkds_metal
   REAL(KIND=rsh),dimension(:,:),ALLOCATABLE   :: cf_xetsig_arh_rco
   REAL(KIND=rsh),dimension(:,:,:),allocatable :: cf_expovar_co
   REAL(KIND=rsh),dimension(:,:,:),allocatable :: cf_stochio_co
   REAL(KIND=rsh),dimension(:,:,:),allocatable :: cf_flim2_co
   REAL(KIND=rsh),dimension(:,:,:),allocatable :: cf_flim1_co
   REAL(KIND=rsh),dimension(:,:,:),allocatable :: cf_flim3_co
   REAL(KIND=rsh),dimension(:),ALLOCATABLE     :: cf_expvol_liq_co
   REAL(KIND=rsh),dimension(:),ALLOCATABLE     :: cf_expvol_gas_co
   CHARACTER(LEN=lchain),dimension(:,:),allocatable    :: name_sit_spc
   CHARACTER(LEN=lchain),dimension(:,:),allocatable    :: name_reac_co
   CHARACTER(len=lchain),dimension(:,:,:),allocatable  :: name_var_co
   CHARACTER(len=lchain),dimension(:,:,:),allocatable  :: name_var_flim1
   CHARACTER(len=lchain),dimension(:,:,:),allocatable  :: name_var_flim2
   CHARACTER(LEN=lchain),dimension(:),allocatable      :: fileconta
   LOGICAL,dimension(:),ALLOCATABLE            :: l_partag
   LOGICAL,dimension(:),ALLOCATABLE            :: l_arh_pco,l_pce_co
   LOGICAL,dimension(:),ALLOCATABLE            :: l_salt_co,l_kds_metal
   LOGICAL,dimension(:,:),allocatable          :: l_treat_spc
   LOGICAL,dimension(:),ALLOCATABLE            :: l_pce_r,l_salt_r
   LOGICAL,dimension(:),ALLOCATABLE            :: l_distb_spc
   LOGICAL,dimension(:),ALLOCATABLE            :: l_kds_r
   LOGICAL,dimension(:,:),ALLOCATABLE          :: l_arh_rco,l_treat_co
   LOGICAL,dimension(:,:),ALLOCATABLE          :: l_coefcte_co
   LOGICAL,dimension(:,:),ALLOCATABLE          :: l_cincte_co
   LOGICAL,dimension(:),ALLOCATABLE            :: l_volkl_difco
   LOGICAL,dimension(:),ALLOCATABLE            :: l_volkl_salco
   LOGICAL,dimension(:),ALLOCATABLE            :: l_volkl_corco
   LOGICAL,dimension(:),ALLOCATABLE            :: l_volkg_corco

 CONTAINS

   !!======================================================================

  SUBROUTINE conta_read_data(obc_cv_name_r,  &
             nb_frac_var_r,icon_var_r, init_cv_name_r,ispc_var_r,                &
             flx_atm_r,cv_rain_r, cini_wat_r,cini_sed_r,cini_air_r,l_out_subs_r, &
             tocd_r,ros_r,diam_r,ws_free_opt_r,ws_free_min_r,ws_free_max_r,     &
             ws_free_para_r,ws_hind_opt_r,ws_hind_para_r                         &
#ifdef key_species
                   ,frac_cwat_r,frac_csed_r   &
#endif
                                            )
                                            

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE conta_read_data  ***
   !&E
   !&E ** Purpose : Reading the file which describes chemical characteristics
   !&E              et species of the studied contaminant 
   !&E
   !&E ** Description :
   !&E
   !&E ** used ij-arrays : none
   !&E
   !&E ** modified variables : none
   !&E
   !&E ** Called by : traceur.F90
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2008-06  (B. Thouvenin)  Original code
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
    USE parameters, ONLY : epsilon
    USE comsubstance, ONLY : name_filesubs,name_var,nv_tot,nvp,nv_ncp,nv_sorb,  &
                             nv_adv,nv_dis,nv_state,nv_dri,nv_int

   !! * Arguments
   LOGICAL,dimension(nb_var),intent(inout)               :: l_out_subs_r
   CHARACTER(len=lchain),dimension(nb_var),intent(inout)     :: init_cv_name_r
   CHARACTER(len=lchain),dimension(nb_var),intent(inout)     :: obc_cv_name_r
   INTEGER,dimension(nb_var),intent(inout)               :: icon_var_r
   INTEGER,dimension(nb_var),intent(inout)               :: nb_frac_var_r
   REAL(KIND=rsh),dimension(nb_var),intent(inout)        :: flx_atm_r,diam_r
   REAL(KIND=rsh),dimension(nb_var),intent(inout)        :: tocd_r,ros_r
   REAL(KIND=rsh),dimension(nb_var),intent(inout)        :: cini_wat_r
   REAL(KIND=rsh),dimension(nb_var),intent(inout)        :: cini_sed_r,cv_rain_r
   REAL(KIND=rsh),dimension(nb_var),intent(inout)        :: cini_air_r
   REAL(KIND=rsh), DIMENSION(4,nb_var),intent(inout)     :: ws_free_para_r
   REAL(KIND=rsh), DIMENSION(nb_var),intent(inout)       :: ws_free_min_r,ws_free_max_r 
   REAL(KIND=rsh), DIMENSION(2,nb_var),intent(inout)     :: ws_hind_para_r 
   INTEGER       , DIMENSION(nb_var),intent(inout)       :: ws_hind_opt_r,ws_free_opt_r 
#ifdef key_species
   REAL(KIND=rsh),dimension(nb_var,nfracmax),intent(inout):: frac_cwat_r
   REAL(KIND=rsh),dimension(nb_var,nfracmax),intent(inout):: frac_csed_r
   INTEGER,dimension(nb_var,nfracmax),intent(inout)       :: ispc_var_r
#else
   INTEGER,dimension(nb_var),intent(inout)               :: ispc_var_r
#endif
   !! * Local declarations
   LOGICAL                                     :: ex,l_unitKd
   LOGICAL                                     :: l_varnew,l_KocKow,l_Kocsolub
   INTEGER                                     :: iv,nlmax,nfl2,ico,ifracv,nfl3
   INTEGER                                     :: nfunct_limit,nvtot_anc,iflim
   INTEGER                                     :: ispc,inl,ispcl,isubs,ir,nblu
   INTEGER                                     :: ityp_flim,nfl1,ispcp,iconta
   INTEGER                                     :: it,il,ind_white,nv_ncp_co,nv_sorb_co
   CHARACTER(LEN=1)                            :: num_part
   CHARACTER(LEN=5)                            :: comment
   CHARACTER(LEN=18)                           :: unitvar_r,nameread
   CHARACTER(LEN=18)                           :: name_spc_assoc
   REAL(KIND=rsh)                              :: cf_flim,flx_atm_t,frain_t
   REAL(KIND=rsh)                              :: finit_watt,finit_sedt,fspcw_tot,fspcs_tot
   REAL(KIND=rsh)                              :: xakoc,xbkoc,valread
   REAL(KIND=rsh),dimension(:,:),allocatable   :: f_init_wat_spc,f_init_sed_spc,f_rain_spc
   REAL(KIND=rsh),dimension(:,:),allocatable   :: f_atm_spc
   REAL(KIND=rsh),DIMENSION(:)      ,ALLOCATABLE :: ws_spc_free_min, ws_spc_free_max
   REAL(KIND=rsh),DIMENSION(:,:)    ,ALLOCATABLE :: ws_spc_free_para, ws_spc_hind_para
   INTEGER       ,DIMENSION(:)      ,ALLOCATABLE :: ws_spc_free_opt,ws_spc_hind_opt
   REAL(KIND=rsh),dimension(:),allocatable     :: tocd_spc,ros_spc,diam_spc
   
   !!----------------------------------------------------------------------
   !! * Executable part

   ! save into simu.log
   !-------------------
   IF_MPI (MASTER) THEN
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '**************************************************'
   WRITE(iscreenlog,*) '**************** MET&OR_read_CONTA.F90 ***************'
   WRITE(iscreenlog,*) '**************************************************'
   WRITE(iscreenlog,*)
   ENDIF_MPI 

!  ALLOCATE arrays of dimension nb_conta
   ALLOCATE(fileconta(nb_conta))
   ALLOCATE(nb_spc_diss(nb_conta))
   ALLOCATE(nb_spc_part_eq(nb_conta))
   ALLOCATE(nb_spc_part_ci(nb_conta))
   ALLOCATE(nb_spc_part_ne(nb_conta))
   ALLOCATE(nb_spc_part_tot(nb_conta))
   ALLOCATE(nb_spc_tot(nb_conta))
   ALLOCATE(name_conta(nb_conta))
   ALLOCATE(massmol_co(nb_conta))
   ALLOCATE(l_metal(nb_conta))
   ALLOCATE(l_organic(nb_conta))
   ALLOCATE(cf_kow_corg(nb_conta))
   ALLOCATE(solubil_co(nb_conta))
   ALLOCATE(cf_difwat_co(nb_conta))
   ALLOCATE(cf_difair_co(nb_conta))
   ALLOCATE(cf_henry_co(nb_conta))
   ALLOCATE(cf_henry_adim_co(nb_conta))
   ALLOCATE(cf_temp_henry_co(nb_conta))
   ALLOCATE(cf_salt_henry_co(nb_conta))
   ALLOCATE(nschmidt_co(nb_conta))
   ALLOCATE(l_arh_pco(nb_conta))
   ALLOCATE(ityp_arh_pco(nb_conta))
   ALLOCATE(cf_teta_arh_pco(nb_conta))
   ALLOCATE(tempref_arh_pco(nb_conta))
   ALLOCATE(cf_xetsig_arh_pco(nb_conta))
   ALLOCATE(l_pce_co(nb_conta))
   ALLOCATE(cf_pce_co(nb_conta,2))
   ALLOCATE(cf_salt_co(nb_conta,2))
   ALLOCATE(l_salt_co(nb_conta))
   ALLOCATE(l_kds_metal(nb_conta))
   ALLOCATE(ndg_kds_metal(nb_conta))
   ALLOCATE(cf_kd0_metal(nb_conta))
   ALLOCATE(l_volatil(nb_conta))
   ALLOCATE(l_sorpt_nlin(nb_conta))
   ALLOCATE(l_volkl_difco(nb_conta))
   ALLOCATE(l_volkl_salco(nb_conta))
   ALLOCATE(ich_volkl_windco(nb_conta))
   ALLOCATE(l_volkl_corco(nb_conta))
   ALLOCATE(cf_expvol_liq_co(nb_conta))
   ALLOCATE(l_volkg_corco(nb_conta))
   ALLOCATE(cf_expvol_gas_co(nb_conta))
   ALLOCATE(cf_xkds_metal(nb_conta,5))
   ALLOCATE(nb_reac_genco(nb_conta))
   ALLOCATE(l_partag(nb_conta))
   
! preliminary reading the fileconta to compute dimensions of species 
!   and ALLOCATE arrays
   DO ico=1,nb_conta
    eof = 0
    SELECT CASE(ico)
     CASE (1)
      fileconta(ico)=fileconta1
     CASE (2)
      fileconta(ico)=fileconta2
     CASE (3)
      fileconta(ico)=fileconta3
     CASE (4)
      fileconta(ico)=fileconta4
     CASE (5)
      fileconta(ico)=fileconta5
    END SELECT              
     
    INQUIRE(file=fileconta(ico),exist=ex)
    IF (ex) THEN
     OPEN(49,file = fileconta(ico),form='formatted')
     comment='debut'
     DO WHILE (comment /= '*****')
       READ(49,'(a)',iostat=eof) comment
     END DO
     READ(49,'(a)',iostat=eof) nameread
     READ(49,'(a)',iostat=eof) nameread 
     READ(49,*,iostat=eof)nb_spc_diss(ico)
     READ(49,*,iostat=eof)nb_spc_part_eq(ico)
     READ(49,*,iostat=eof)nb_spc_part_ci(ico)
     READ(49,*,iostat=eof)nb_spc_part_ne(ico)
     nb_spc_part_tot(ico)=nb_spc_part_eq(ico)+nb_spc_part_ci(ico)+  &
                          nb_spc_part_ne(ico)     
     nb_spc_tot(ico)=nb_spc_diss(ico)+nb_spc_part_tot(ico)
     CLOSE(49)     
    END IF     ! fin du test sur l existence du fichier
   ENDDO 
   nb_spcmax=maxval(nb_spc_tot)
   ALLOCATE(ivar_spc(nb_conta,nb_spcmax))
   ALLOCATE(name_spc(nb_conta,nb_spcmax))
   ALLOCATE(cin_hydrol(nb_conta,nb_spcmax))
   ALLOCATE(cin_photo(nb_conta,nb_spcmax))
   ALLOCATE(cin_biodeg(nb_conta,nb_spcmax))
   ALLOCATE(name_sit_spc(nb_conta,nb_spcmax))
   ALLOCATE(cf_distb_spc(nb_conta,nb_spcmax))
   ALLOCATE(isit_spc(nb_conta,nb_spcmax))
   ALLOCATE(f_init_wat_spc(nb_conta,nb_spcmax))
   ALLOCATE(f_init_sed_spc(nb_conta,nb_spcmax))
!   ALLOCATE(flim_n_spc(nb_conta,nb_spcmax))
!   ALLOCATE(flim_e_spc(nb_conta,nb_spcmax))
!   ALLOCATE(flim_s_spc(nb_conta,nb_spcmax))
!   ALLOCATE(flim_w_spc(nb_conta,nb_spcmax))
   ALLOCATE(f_atm_spc(nb_conta,nb_spcmax))
   ALLOCATE(f_rain_spc(nb_conta,nb_spcmax))
   ALLOCATE(nb_sorpt_nlin(nb_conta,nb_spcmax))
   ALLOCATE(f_sorpt_nlin(nb_conta,nb_spcmax,2))
   ALLOCATE(cf_kfkl_nlin(nb_conta,nb_spcmax,2))
   ALLOCATE(cf_expqmax_nlin(nb_conta,nb_spcmax,2))
   ALLOCATE(ispc_assoc_pcin(nb_conta,nb_spcmax))
   ALLOCATE(cf_rate_spc(nb_conta,nb_spcmax))
   ALLOCATE(cv_watequ_spc(nb_conta,nb_spcmax))
   ALLOCATE(l_treat_spc(nb_conta,nb_spcmax))
   ALLOCATE(cf_flim1_co(nb_conta,nrecomax,nflimcomax))
   ALLOCATE(cf_flim2_co(nb_conta,nrecomax,nflimcomax))
   ALLOCATE(cf_flim3_co(nb_conta,nrecomax,nflimcomax))
   ALLOCATE(name_reac_co(nb_conta,nrecomax))
   ALLOCATE(l_treat_co(nb_conta,nrecomax))
   ALLOCATE(nb_var_co(nb_conta,nrecomax))
   ALLOCATE(cin_max_co(nb_conta,nrecomax))
   ALLOCATE(l_cincte_co(nb_conta,nrecomax)) 
   ALLOCATE(l_coefcte_co(nb_conta,nrecomax)) 
   ALLOCATE(l_arh_rco(nb_conta,nrecomax)) 
   ALLOCATE(ityp_arh_rco(nb_conta,nrecomax)) 
   ALLOCATE(cf_teta_arh_rco(nb_conta,nrecomax)) 
   ALLOCATE(tempref_arh_rco(nb_conta,nrecomax)) 
   ALLOCATE(cf_xetsig_arh_rco(nb_conta,nrecomax)) 
   ALLOCATE(cf_stochio_co(nb_conta,nrecomax,nvactcomax)) 
   ALLOCATE(cf_expovar_co(nb_conta,nrecomax,nvactcomax)) 
   ALLOCATE(name_var_co(nb_conta,nrecomax,nvactcomax)) 
   ALLOCATE(nb_flim1_co(nb_conta,nrecomax))  
   ALLOCATE(nb_flim2_co(nb_conta,nrecomax))  
   ALLOCATE(nb_flim3_co(nb_conta,nrecomax))  
   ALLOCATE(name_var_flim1(nb_conta,nrecomax,nflimcomax))
   ALLOCATE(name_var_flim2(nb_conta,nrecomax,nflimcomax))
   ALLOCATE(ws_spc_free_opt(nb_spcmax)) 
   ALLOCATE(ws_spc_hind_opt(nb_spcmax)) 
   ALLOCATE(ws_spc_free_para(4,nb_spcmax)) 
   ALLOCATE(ws_spc_free_min(nb_spcmax)) 
   ALLOCATE(ws_spc_free_max(nb_spcmax)) 
   ALLOCATE(ws_spc_hind_para(2,nb_spcmax)) 
   ALLOCATE(ros_spc(nb_spcmax))
   ALLOCATE(tocd_spc(nb_spcmax))
   ALLOCATE(diam_spc(nb_spcmax))
#ifdef l_species
   ALLOCATE(ifrac_var_spc(nb_conta,nb_spcmax)
#endif
   
  DO ico=1,nb_conta
    iconta=ivdl_conta(ico)
    eof = 0
    INQUIRE(file=fileconta(ico),exist=ex)
    IF_MPI (MASTER) THEN
    WRITE(iscreenlog,*) 'fileconta = ',trim(fileconta(ico))
    ENDIF_MPI
    IF (ex) THEN
     OPEN(49,file = fileconta(ico),form='formatted')
     comment='debut'
     DO WHILE (comment /= '*****')
       READ(49,'(a)',iostat=eof) comment
     END DO
     READ(49,'(a)',iostat=eof) nameread
     ind_white=index(nameread,' ')
     name_conta(ico)=nameread(1:ind_white)
     IF(name_conta(ico)/=name_var(iconta))then
      WRITE(*,*)'*****  ERROR   *******'
      WRITE(*,*)'the contaminant ',name_var(iconta),'READing in the file :',&
                                    trim(name_filesubs)
      WRITE(*,*)'is not the same as the contaminant ',trim(name_conta(ico)), &
                        'READing in the file :',trim(fileconta(iconta))
      WRITE(*,*)&
         'CHECK to simulate your contaminant with the appropriate data file'
      WRITE(*,*)&
      ' either into these 2 files either in paraconta.txt (namconta)'
      STOP
     ENDIF
     READ(49,'(a)',iostat=eof) unitvar_r 
     ind_white=index(unitvar_r,' ')
     unitvar_r=unitvar_r(1:ind_white)
     unit_var(iconta)=unitvar_r
  ! number of species already readed   
     READ(49,*,iostat=eof)
     READ(49,*,iostat=eof)
     READ(49,*,iostat=eof)
     READ(49,*,iostat=eof)

! Reading chemical CHARACTERistics
! ================================
     READ(49,*,iostat=eof)massmol_co(ico)
     READ(49,*,iostat=eof)l_metal(ico)
     IF(l_metal(ico))then
       l_organic(ico)=.false.
     ELSE
       l_organic(ico)=.true.
       READ(49,*,iostat=eof)cf_kow_corg(ico)
       READ(49,*,iostat=eof)solubil_co(ico)
     ENDIF
     READ(49,*,iostat=eof)nschmidt_co(ico)
     READ(49,*,iostat=eof)cf_difwat_co(ico)
     READ(49,*,iostat=eof)cf_difair_co(ico)
     READ(49,*,iostat=eof)cf_henry_co(ico)
     READ(49,*,iostat=eof)cf_henry_adim_co(ico)
     READ(49,*,iostat=eof)cf_temp_henry_co(ico)
     READ(49,*,iostat=eof)cf_salt_henry_co(ico)
     READ(49,*,iostat=eof)

! Reading temperature, salinity,MES effects on sorption/desorption process
! =====================================================================
! arhenius Law
     READ(49,*,iostat=eof)l_arh_pco(ico)
     READ(49,*,iostat=eof)ityp_arh_pco(ico)
     READ(49,*,iostat=eof)cf_teta_arh_pco(ico)
     READ(49,*,iostat=eof)tempref_arh_pco(ico)
     READ(49,*,iostat=eof)cf_xetsig_arh_pco(ico)
! Particule Concentration Effect PCE
     READ(49,*,iostat=eof)l_pce_co(ico)
     READ(49,*,iostat=eof)cf_pce_co(ico,:)
     IF(.not.l_pce_co(ico))then
          cf_pce_co(ico,1)=0.
          cf_pce_co(ico,2)=1.
     ENDIF        
! effect salting out         
     READ(49,*,iostat=eof)l_salt_co(ico)
     READ(49,*,iostat=eof)cf_salt_co(ico,:)
     IF(.not.l_salt_co(ico))then
          cf_salt_co(ico,1)=0.
          cf_salt_co(ico,2)=0.
     ENDIF
     
! ONLY FOR METALS
!  Reading KD law if sorption/desorption process treated globally
!  with a function issued from MOCO or Turner 
! =====================================================================
     IF(l_metal(ico))then
        READ(49,*,iostat=eof)l_kds_metal(ico)
        IF(l_kds_metal(ico))then
            READ(49,*,iostat=eof)ndg_kds_metal(ico),cf_kd0_metal(ico)
            cf_xkds_metal(ico,:)=0.
            IF(ndg_kds_metal(ico).gt.0)then
              READ(49,*,iostat=eof)(cf_xkds_metal(ico,il), &
                                              il=1,ndg_kds_metal(ico))
            ELSE
              READ(49,*,iostat=eof)
            ENDIF
        ELSE
          READ(49,*,iostat=eof)
          READ(49,*,iostat=eof)
        ENDIF 
     ENDIF 
        
! Reading data for first dissolved specie (free dissolved) 
! ==========================================================
     ispc=1
     READ(49,*,iostat=eof)nameread
     ind_white=index(nameread,' ')
     name_spc(ico,ispc)=nameread(1:ind_white)
     READ(49,*,iostat=eof)cin_hydrol(ico,ispc)
     cin_hydrol(ico,ispc)=cin_hydrol(ico,ispc)/86400
     READ(49,*,iostat=eof)cin_photo(ico,ispc)
     cin_photo(ico,ispc)=cin_photo(ico,ispc)/86400
     READ(49,*,iostat=eof)cin_biodeg(ico,ispc)
     cin_biodeg(ico,ispc)=cin_biodeg(ico,ispc)/86400 
              
! Reading data for other dissolved variables
! ==========================================================
     READ(49,*,iostat=eof)
     DO ispc=2,nb_spc_diss(ico)
        READ(49,*,iostat=eof)nameread
        ind_white=index(nameread,' ')
        name_spc(ico,ispc)=nameread(1:ind_white)
        READ(49,*,iostat=eof)nameread
        ind_white=index(nameread,' ')
        name_sit_spc(ico,ispc)=nameread(1:ind_white)
        READ(49,*,iostat=eof)l_unitKd
        READ(49,*,iostat=eof)cf_distb_spc(ico,ispc)
        IF(.not.l_unitKd)then
!         coeffient de partage exprime  en log10(ml/g) converti en  L/G 
           cf_distb_spc(ico,ispc)=10**(cf_distb_spc(ico,ispc))/1000
        ENDIF
        IF(l_organic(ico))then
          READ(49,*,iostat=eof)l_KocKow
          READ(49,*,iostat=eof)l_Kocsolub
          READ(49,*,iostat=eof)xakoc,xbkoc
          IF(l_KocKow)THEN
           cf_distb_spc(ico,ispc)=xakoc*log10(cf_kow_corg(ico))+xbkoc
          ELSE IF(l_Kocsolub)then
           cf_distb_spc(ico,ispc)=xakoc*log10(solubil_co(ico))+xbkoc
          ENDIF
        ENDIF
        READ(49,*,iostat=eof)cin_hydrol(ico,ispc)
        cin_hydrol(ico,ispc)=cin_hydrol(ico,ispc)/86400
        READ(49,*,iostat=eof)cin_photo(ico,ispc)
        cin_photo(ico,ispc)=cin_photo(ico,ispc)/86400
        READ(49,*,iostat=eof)cin_biodeg(ico,ispc)
        cin_biodeg(ico,ispc)=cin_biodeg(ico,ispc)/86400          
! recherche du numero de la variable qui donne le substrat
! pour calcul partage
          DO iv=1,nb_var
            IF(name_var(iv)==name_sit_spc(ico,ispc))then
              isit_spc(ico,ispc)=iv
              EXIT
            ENDIF 
          ENDDO         
      ENDDO  
            
! Reading data for particulate variables which distribution is in equilibrium  
! ============================================================================
      READ(49,*,iostat=eof)
      l_sorpt_nlin(ico)=.false.
      l_partag(ico)=.false.
      DO ispc=nb_spc_diss(ico)+1,nb_spc_diss(ico)+nb_spc_part_eq(ico)
          l_partag(ico)=.true.
          READ(49,*,iostat=eof)nameread
          ind_white=index(nameread,' ')
          name_spc(ico,ispc)=nameread(1:ind_white)
          READ(49,*,iostat=eof)nameread
          ind_white=index(nameread,' ')
          name_sit_spc(ico,ispc)=nameread(1:ind_white)
! vitesses de chute de cette espece = celle du substrat 
! fraction initiale par rapport a la conc. totale dans l eau
          f_init_wat_spc(ico,ispc)=0.
! fraction initiale par rapport a la conc. totale dans le sediment
          READ(49,*,iostat=eof)f_init_sed_spc(ico,ispc)
! fractions in input at boundaries
!          READ(49,*,iostat=eof)flim_n_spc(ico,ispc),flim_e_spc(ico,ispc),  &
!                               flim_s_spc(ico,ispc),flim_w_spc(ico,ispc)
! fractions in atmosphere (flux and rain)
          READ(49,*,iostat=eof)f_atm_spc(ico,ispc)
          READ(49,*,iostat=eof)f_rain_spc(ico,ispc)
          READ(49,*,iostat=eof)l_unitKd
          READ(49,*,iostat=eof)cf_distb_spc(ico,ispc)
          IF(.not.l_unitKd)then
!           coeffient de partage exprime  log10 et en ml/g
           cf_distb_spc(ico,ispc)=10**(cf_distb_spc(ico,ispc))/1000
          ENDIF
        IF(l_organic(ico))then
          READ(49,*,iostat=eof)l_KocKow
          READ(49,*,iostat=eof)l_Kocsolub
          READ(49,*,iostat=eof)xakoc,xbkoc
          IF(l_KocKow)THEN
           cf_distb_spc(ico,ispc)=xakoc*log10(cf_kow_corg(ico))+xbkoc
          ELSE IF(l_Kocsolub)then
           cf_distb_spc(ico,ispc)=xakoc*log10(solubil_co(ico))+xbkoc
          ENDIF
        ENDIF
          READ(49,*,iostat=eof)nb_sorpt_nlin(ico,ispc)
          READ(49,*,iostat=eof)
          IF(nb_sorpt_nlin(ico,ispc) /= 0)then
           READ(49,*,iostat=eof)l_langmuir(ico)
           l_sorpt_nlin(ico)=.true.
           IF((l_langmuir(ico)).and.(nb_sorpt_nlin(ico,ispc).gt.1.   &
                 .or.nb_spc_part_eq(ico).gt.1))then
             WRITE(*,*)'ATTENTION ::  si on choisit le modele DMM '
             WRITE(*,*)'il ne peut y avoir qu UNE SEULE espece '
             WRITE(*,*)'particulaire a l equilibre '
             WRITE(*,*)' avec UN SEUL domaine non lineaire'
             STOP
           ENDIF
           IF(nb_sorpt_nlin(ico,ispc).gt.2)then
             WRITE(*,*)'ATTENTION ::  le modele actuel ne permet pas de'
             WRITE(*,*)' traiter plus de 2 domaines non lineaires'
             STOP
           ENDIF
           DO inl=1,nb_sorpt_nlin(ico,ispc)
             READ(49,*,iostat=eof)f_sorpt_nlin(ico,ispc,inl)
             READ(49,*,iostat=eof)l_unitKd
             READ(49,*,iostat=eof)cf_kfkl_nlin(ico,ispc,inl)         
             READ(49,*,iostat=eof)cf_expqmax_nlin(ico,ispc,inl)         
             IF(l_unitKd)then
               cf_kfkl_nlin(ico,ispc,inl)=10**cf_kfkl_nlin(ico,ispc,inl)/1000.
             ENDIF
           ENDDO
          ENDIF
          READ(49,*,iostat=eof)cin_hydrol(ico,ispc)
          cin_hydrol(ico,ispc)=cin_hydrol(ico,ispc)/86400
          READ(49,*,iostat=eof)cin_photo(ico,ispc)
          cin_photo(ico,ispc)=cin_photo(ico,ispc)/86400
          READ(49,*,iostat=eof)cin_biodeg(ico,ispc)
          cin_biodeg(ico,ispc)=cin_biodeg(ico,ispc)/86400          
! recherche du numero de la variable qui donne le substrat
! pour reaction de partage et egalisation des parameters sedimento
          DO iv=1,nb_var
            IF(name_var(iv)==name_sit_spc(ico,ispc))then
              isit_spc(ico,ispc)=iv
              EXIT
            ENDIF 
          ENDDO         
      ENDDO        
      READ(49,*,iostat=eof)
        
! Reading data for particulate variables with adsorption/desorption CINETIC  
! ===========================================================================
      ispcp=nb_spc_diss(ico)+nb_spc_part_eq(ico)
      DO ispc=ispcp+1,ispcp+nb_spc_part_ci(ico)
          READ(49,*,iostat=eof)nameread
          ind_white=index(nameread,' ')
          name_spc(ico,ispc)=nameread(1:ind_white)
          READ(49,*,iostat=eof)nameread
          ind_white=index(nameread,' ')
          name_spc_assoc=nameread(1:ind_white)
          DO ispcl=nb_spc_diss(ico)+1,nb_spc_diss(ico)+nb_spc_part_eq(ico)
              IF(name_spc(ico,ispcl).eq.name_spc_assoc)then
                ispc_assoc_pcin(ico,ispc)=ispcl
              ENDIF
          ENDDO
! fraction initiale par rapport a la conc. totale dans l eau
          READ(49,*,iostat=eof)f_init_wat_spc(ico,ispc)
! fraction initiale par rapport a la conc. totale dans le sediment
          READ(49,*,iostat=eof)f_init_sed_spc(ico,ispc)
! fractions in input at boundaries
!          READ(49,*,iostat=eof)flim_n_spc(ico,ispc),flim_e_spc(ico,ispc),   &
!                               flim_s_spc(ico,ispc),flim_w_spc(ico,ispc)
! fractions in atmosphere (flux and rain)
          READ(49,*,iostat=eof)f_atm_spc(ico,ispc)
          READ(49,*,iostat=eof)f_rain_spc(ico,ispc)
          READ(49,*,iostat=eof)cf_rate_spc(ico,ispc)
          cf_rate_spc(ico,ispc)=cf_rate_spc(ico,ispc)/86400
          name_sit_spc(ico,ispc)=name_sit_spc(ico,ispc_assoc_pcin(ico,ispc))
!! concentration in water at equilibrium for the first active variable
!!  ATTENTION unit M/M of particules
          READ(49,*,iostat=eof)cv_watequ_spc(ico,ispc)
          READ(49,*,iostat=eof)l_treat_spc(ico,ispc)
          cin_hydrol(ico,ispc)=cin_hydrol(ico,ispc_assoc_pcin(ico,ispc))
          cin_photo(ico,ispc)=cin_photo(ico,ispc_assoc_pcin(ico,ispc))
          cin_biodeg(ico,ispc)=cin_biodeg(ico,ispc_assoc_pcin(ico,ispc))
! recherche du numero de la variable qui donne le substrat
! egalisation des parameters sedimento
          DO iv=1,nb_var
            IF(name_var(iv)==name_sit_spc(ico,ispc))then
              isit_spc(ico,ispc)=iv
              EXIT
            ENDIF 
          ENDDO         
      ENDDO        
      READ(49,*,iostat=eof)
        
! Reading data for particulate variables without sorption/desorption 
!                            (NON ECHANGeable)
! =====================================================================
      ispcp=nb_spc_diss(ico)+nb_spc_part_eq(ico)+nb_spc_part_ci(ico)
      DO ispc=ispcp+1,ispcp+nb_spc_part_ne(ico)
          READ(49,*,iostat=eof)nameread
          ind_white=index(nameread,' ')
          name_spc(ico,ispc)=nameread(1:ind_white)
          READ(49,*,iostat=eof)nameread
          ind_white=index(nameread,' ')
          name_spc_assoc=nameread(1:ind_white)
! fraction intiale par rapport a la conc. totale (part et dissoute) dans l eau
          READ(49,*,iostat=eof)f_init_wat_spc(ico,ispc)
! fraction intiale par rapport a la conc. totale (part et dissoute) dans le sediment
          READ(49,*,iostat=eof)f_init_sed_spc(ico,ispc)
! fractions in input at boundaries
!          READ(49,*,iostat=eof)flim_n_spc(ico,ispc),flim_e_spc(ico,ispc),   &
!                               flim_s_spc(ico,ispc),flim_w_spc(ico,ispc)
! fractions in atmosphere (flux and rain)
          READ(49,*,iostat=eof)f_atm_spc(ico,ispc)
          READ(49,*,iostat=eof)f_rain_spc(ico,ispc)
          READ(49,*,iostat=eof)cin_hydrol(ico,ispc)
          cin_hydrol(ico,ispc)=cin_hydrol(ico,ispc)/86400
          READ(49,*,iostat=eof)cin_photo(ico,ispc)
          cin_photo(ico,ispc)=cin_photo(ico,ispc)/86400
          READ(49,*,iostat=eof)cin_biodeg(ico,ispc)
          cin_biodeg(ico,ispc)=cin_biodeg(ico,ispc)/86400          
! recherche du numero de la variable qui donne la particule associee
! egalisation des parameters sedimento
          DO iv=1,nb_var
            IF(name_var(iv)==name_spc_assoc)then
              isit_spc(ico,ispc)=iv
              EXIT
            ENDIF 
          ENDDO         
      ENDDO 
        
!  Reading  paremeters for volatilisation process
! ==========================================================
      READ(49,*,iostat=eof)
      READ(49,*,iostat=eof)l_volatil(ico)
! l_volkl_difco : choix de la formule pour coef de transfert phase liquide
!              par diffusion turbulente 
!              true pour Churchill ou false pour O Connor 
      READ(49,*,iostat=eof)l_volkl_difco(ico)
! l_volkl_salco : choix de la prise en compte de l effet de sel (true si oui)
!              par formule de Wen  pour coef de transfert phase liquide
      READ(49,*,iostat=eof)l_volkl_salco(ico)
! ich_volkl_windco : choix de la formule pour coef de transfert phase liquide
!              par turbulence engendree par le vent 
!              1 pour Schwartzenbach 
!              2 pour Wanninkhof 
!              3 pour Nightingale 
      READ(49,*,iostat=eof)ich_volkl_windco(ico)
! l_volkl_corco : choix de la propriete pour correction du coefficient  
!               de transfert phase liquide
!               par nombre de Schmidt (true) ou pas massse molaire (false) 
      READ(49,*,iostat=eof)l_volkl_corco(ico)
! cf_expvol_liq_co : exposant du rapport des proprietes pour phase liquide 
      READ(49,*,iostat=eof)cf_expvol_liq_co(ico)
! l_volkg_corco : choix de la propriete pour correction du coefficient 
!               de transfert phase gazeuse
!               par diffusivite (true) ou pas massse molaire (false) 
      READ(49,*,iostat=eof)l_volkg_corco(ico)
! cf_expvol_gas_co : exposant du rapport des proprietes pour phase gazeuse
      READ(49,*,iostat=eof)cf_expvol_gas_co(ico)
          
!  Reading others reactions (type generic) affected the species
! ==========================================================
     READ(49,*,iostat=eof)
     READ(49,*,iostat=eof)nb_reac_genco(ico)
     DO ir=1,nb_reac_genco(ico)
       READ(49,'(a)',iostat=eof)nameread
       ind_white=index(nameread,' ')
       name_reac_co(ico,ir)=nameread(1:ind_white)
       READ(49,*,iostat=eof)l_treat_co(ico,ir)
       READ(49,*,iostat=eof)nblu
       nb_var_co(ico,ir)=nblu
       READ(49,*,iostat=eof)cin_max_co(ico,ir)
       cin_max_co(ico,ir)=cin_max_co(ico,ir)/86400.
       READ(49,*,iostat=eof)l_cincte_co(ico,ir)
       READ(49,*,iostat=eof)l_coefcte_co(ico,ir)
       READ(49,*,iostat=eof)nfunct_limit
! law of arhenius
       READ(49,*,iostat=eof)l_arh_rco(ico,ir)
       READ(49,*,iostat=eof)ityp_arh_rco(ico,ir)
       READ(49,*,iostat=eof)cf_teta_arh_rco(ico,ir)
       READ(49,*,iostat=eof)tempref_arh_rco(ico,ir)
       READ(49,*,iostat=eof)cf_xetsig_arh_rco(ico,ir)
       READ(49,'(a)',iostat=eof)
! reactives variables 
       DO it=1,nb_var_co(ico,ir)
        READ(49,'(a)',iostat=eof) nameread ! name of reactives variables  
        ind_white=index(nameread,' ')
        name_var_co(ico,ir,it)=nameread(1:ind_white)
        READ(49,*,iostat=eof)cf_stochio_co(ico,ir,it)
        READ(49,*,iostat=eof)cf_expovar_co(ico,ir,it)
       ENDDO
       READ(49,'(a)',iostat=eof)
! Limiting functions 
      IF(nfunct_limit.gt.0) then
         nb_flim1_co(ico,ir)=0
         nb_flim2_co(ico,ir)=0
         nb_flim3_co(ico,ir)=0
         DO iflim=1,nfunct_limit
           READ(49,*,iostat=eof)ityp_flim
           READ(49,*,iostat=eof)cf_flim
           SELECT CASE (ityp_flim)
            CASE (1)
   ! Limiting functions type Monod 
              READ(49,*,iostat=eof)nameread
              nb_flim1_co(ico,ir)=nb_flim1_co(ico,ir)+1
              nfl1=nb_flim1_co(ico,ir)
              cf_flim1_co(ico,ir,nfl1)=cf_flim
              ind_white=index(nameread,' ')
              name_var_flim1(ico,ir,nfl1)=nameread(1:ind_white)
           CASE (2) 
   ! Limiting functions type Inhibitrice 
              READ(49,*,iostat=eof)nameread
              nb_flim2_co(ico,ir)=nb_flim2_co(ico,ir)+1
              nfl2=nb_flim2_co(ico,ir)
              cf_flim2_co(ico,ir,nfl2)=cf_flim
              ind_white=index(nameread,' ')
              name_var_flim2(ico,ir,nfl2)=nameread(1:ind_white)
           CASE (3) 
   ! Limiting functions by solar radiation 
              nb_flim3_co(ico,ir)=nb_flim3_co(ico,ir)+1
              nfl3=nb_flim3_co(ico,ir)
              ind_white=index(nameread,' ')
              cf_flim3_co(ico,ir,nfl3)=cf_flim
           END SELECT
          ENDDO
       ENDIF

     ENDDO
      
     CLOSE(49)
   END IF     ! fin du test sur l existence du fichier
   
       
   ! save into simu.log
   !-------------------
   IF_MPI (MASTER) THEN
   WRITE(iscreenlog,*)' '
   WRITE(iscreenlog,*)' '
   WRITE(iscreenlog,*)'Chemical CHARACTERistics of ',trim(name_var(iconta))
   WRITE(iscreenlog,*)'simulated contaminant number : ',ico
   WRITE(iscreenlog,*)' '
   WRITE(iscreenlog,*)'Molecular mass             =',massmol_co(ico)
   WRITE(iscreenlog,*)'Schmidt number in water    =',nschmidt_co(ico)
   WRITE(iscreenlog,*)'molecular diffusion in water    =',cf_difwat_co(ico)
   WRITE(iscreenlog,*)'Diffusivity in air         =',cf_difair_co(ico)
   IF(cf_henry_co(ico).ne.0)then
     WRITE(iscreenlog,*)'using  dimensional Henry Constant  =', &
                                                          cf_henry_co(ico)
   ELSE
     WRITE(iscreenlog,*)'using adimensional Henry Constant  =',&
                                                      cf_henry_adim_co(ico)
   ENDIF
   WRITE(iscreenlog,*)'B correc Henry with Staudinger f(T)  =',&
                                                     cf_temp_henry_co(ico)
   WRITE(iscreenlog,*)'correc Henry salting out (L/mol) f(Sal) =',&
                                                     cf_salt_henry_co(ico)
   WRITE(iscreenlog,*)' '
   WRITE(iscreenlog,*)' '
   WRITE(iscreenlog,*)'Speciation and processes of ',trim(name_var(iconta))
   WRITE(iscreenlog,*)' '
   WRITE(iscreenlog,*)nb_spc_diss(ico),' dissolved species :  '
   WRITE(iscreenlog,*)nb_spc_part_eq(ico),' sorbed species in linear Equ.'
   DO ispc=nb_spc_diss(ico)+1,nb_spc_diss(ico)+nb_spc_part_eq(ico)
      WRITE(iscreenlog,*)'in sediment :initial fractions/total contaminant species n°:',&
                                                   ispc,f_init_sed_spc(ico,ispc)
   ENDDO
   IF(l_sorpt_nlin(ico))then
    DO ispc=nb_spc_diss(ico)+1,nb_spc_diss(ico)+nb_spc_part_eq(ico)
     WRITE(iscreenlog,*)nb_sorpt_nlin(ico,ispc),&
                      'non linear domain for species number ',ispc
     IF(l_langmuir(ico))then
        WRITE(iscreenlog,*)' model DMM  Langmuir '
        DO inl=1, nb_sorpt_nlin(ico,ispc)
          WRITE(iscreenlog,*)' KL= ',cf_kfkl_nlin(ico,ispc,inl),&
                          ' Qmax= ' ,cf_expqmax_nlin(ico,ispc,inl)
        ENDDO
     ELSE
        WRITE(iscreenlog,*)' model DRM  Freundlich ' 
        DO inl=1, nb_sorpt_nlin(ico,ispc)
          WRITE(iscreenlog,*)' KF= ',cf_kfkl_nlin(ico,ispc,inl),&
                             ' n= ',cf_expqmax_nlin(ico,ispc,inl)
        ENDDO
     ENDIF
    ENDDO
   ELSE
     WRITE(iscreenlog,*)'0 non linear domain'
   ENDIF
   IF(nb_spc_part_ci(ico).ne.0.)then
     WRITE(iscreenlog,*)nb_spc_part_ci(ico),&
                     ' sorbed species with ad/desorption cinetic'
     ispcp=nb_spc_diss(ico)+nb_spc_part_eq(ico)
     DO ispc=ispcp+1,ispcp+nb_spc_part_ci(ico)
      WRITE(iscreenlog,*)'in water : initial fractions/total contaminant species n°:',&
                                                   ispc,f_init_wat_spc(ico,ispc)
      WRITE(iscreenlog,*)'in sediment :initial fractions/total contaminant species n°:',&
                                                   ispc,f_init_sed_spc(ico,ispc)
     ENDDO
   ENDIF
   IF(nb_spc_part_ne(ico).ne.0.)then
     WRITE(iscreenlog,*)nb_spc_part_ne(ico),' sorbed species NON exchangeable' 
     ispcp=nb_spc_diss(ico)+nb_spc_part_eq(ico)+nb_spc_part_ci(ico)
     DO ispc=ispcp+1,ispcp+nb_spc_part_ne(ico)
      WRITE(iscreenlog,*)'in water : initial fractions/total contaminant species n°:',&
                                                   ispc,f_init_wat_spc(ico,ispc)
      WRITE(iscreenlog,*)'in sediment :initial fractions/total contaminant species n°:',&
                                                   ispc,f_init_sed_spc(ico,ispc)
     ENDDO
   ENDIF
   IF(l_volatil(ico))then
     WRITE(iscreenlog,*)' '
     WRITE(iscreenlog,*)' '
     WRITE(iscreenlog,*)'Volatilisation process of ',trim(name_var(iconta))
     WRITE(iscreenlog,*)' '
     IF(l_volkl_difco(ico))then
       WRITE(iscreenlog,*)'Churchill formulation for KL by turbulent diffusion'
      ELSE
       WRITE(iscreenlog,*)'O Connor formulation for KL by turbulent diffusion'
      ENDIF
     IF(l_volkl_salco(ico))then
       WRITE(iscreenlog,*)&
         'effect of salinity (Wen formul) for KL by turbulent diffusion'
      ELSE
       WRITE(iscreenlog,*)'no effect of salinity for KL by turbulent diffusion'
      ENDIF
     IF(ich_volkl_windco(ico) == 1)then
       WRITE(iscreenlog,*)'Schwartzenbach formulation for KL induced by wind '
      ELSE IF(ich_volkl_windco(ico) == 2)then
       WRITE(iscreenlog,*)'Wanninkhof formulation for KL induced by wind '
      ELSE 
       WRITE(iscreenlog,*)'Nightingale formulation for KL induced by wind '
      ENDIF
     IF(l_volkl_corco(ico))then
       WRITE(iscreenlog,*)'for Kl : relation between molar mass '
      ELSE
       WRITE(iscreenlog,*)'for Kl : relation between Schmidt numbers '
      ENDIF
      WRITE(iscreenlog,*)'exponent nL= ',cf_expvol_liq_co(ico)
     IF(l_volkg_corco(ico))then
       WRITE(iscreenlog,*)'for Kg : relation between diffusivities '
      ELSE
       WRITE(iscreenlog,*)'for Kl : relation between Schmidt numbers '
      ENDIF
      WRITE(iscreenlog,*)'exponent ng= ',cf_expvol_gas_co(ico)
   ENDIF
   ENDIF_MPI
   
#ifdef key_species
!  PLUSIEURS ESPECES POSSIBLES / VARIABLE (KEY_SPECIES)
!!  IF l_gpspecies is .true. :
!    Grouping particulate species with same substrat for set up sorbed advected variables
!     only when they are distributed at equilibrium with dissolved variables
!!  IF l_gpspecies is .false. : each species becomes a new variable
!!=========================================================================================
! reperage des especes dans chaque variable
!    ispc_var_r(iv,ispc) : numero de l espece dans la variable iv
!    icon_var_r(iv) : numero du contaminant correspondant aux especes 
!                     inclus dans la variable iv
!     nb_frac_var_r(iv)   : nombre d especes comprises dans la variable iv
!                           (somme des fractions = 1)
!   tandis que ivar_spc(ico,ispc) = numero de la variable qui inclut l espece
!                                   ispc du contaminant ico
!   et ifrac_var_spc(ico,ispc) = numero de la fraction dans la variable 
!                    ivar_spc correspondant à l espece ispc du contaminant ico 
!!=========================================================================================
   frac_cwat_r(:,:)=0.
   nb_frac_var_r(:)=1
   nvtot_anc=nv_tot
   nv_ncp_co=0
   nv_sorb_co=0
   flx_atm_t=0.
   frain_t=0.
   finit_watt=0.
   finit_sedt=0.
   IF(nb_spc_part_tot(ico) /= 0)THEN
     ispc=nb_spc_diss(ico)+1
     IF(ispc == 1)then
       ivar_spc(ico,ispc)=iconta
       isubs=nv_tot
     ELSE
       ivar_spc(ico,ispc)=nv_tot+1 
       isubs=nv_tot+1
     ENDIF
     iv=ivar_spc(ico,ispc)
     ifracv=1
     nb_frac_var_r(iv)=ifracv
     ifrac_var_spc(ico,ispc)=ifracv
     ispc_var_r(iv,ifracv)=ispc
     icon_var_r(iv)=ico
     IF(itypv(isit_spc(ico,ispc))==5) THEN   ! variable substrat is a SORB
       itypv_r(iv)=5
       nv_sorb_co=nv_sorb_co+1
       name_varpc_assoc(iv)=name_var(isit_spc(ico,ispc))
     ELSE IF (itypv(isit_spc(ico,ispc))==4) THEN   ! variable substrat is a NoCP
       itypv_r(iv)=4
       nv_ncp_co=nv_ncp_co+1
       ws_spc_free_min(iv)=ws_free_min_r(isit_spc(ico,ispc))
       ws_spc_free_max(iv)=ws_free_max_r(isit_spc(ico,ispc))
       ws_spc_free_para(:,iv)=ws_free_para_r(:,isit_spc(ico,ispc))
       ws_spc_hind_para(:,iv)=ws_hind_para_r(:,isit_spc(ico,ispc))
       ws_spc_free_opt(iv)=ws_free_opt_r(isit_spc(ico,ispc))
       ws_spc_hind_opt(iv)=ws_hind_opt_r(isit_spc(ico,ispc))
     ENDIF
     l_out_subs_r(iv)=.true.
     name_var(iv)=trim(name_var(iconta))//'_P1'
     init_cv_name_r(iv)=init_cv_name_r(iconta)
     obc_cv_name_r(iv)='none'
     tocd_r(iv)=tocd_spc(ispc)
     ros_r(iv)=ros_spc(ispc)
     diam_r(iv)=diam_spc(ispc)
     flx_atm_r(iv)=flx_atm_r(iv)+f_atm_spc(ico,ispc)*flx_atm_r(iconta)
     cv_rain_r(iv)=cv_rain_r(iv)+f_rain_spc(ico,ispc)*cv_rain_r(iconta)
     cini_wat_r(iv)=cini_wat_r(iv)+f_init_wat_spc(ico,ispc)*cini_wat_r(iconta)
     cini_sed_r(iv)=cini_sed_r(iv)+f_init_sed_spc(ico,ispc)*cini_sed_r(iconta)
     cini_air_r(iv)=cini_air_r(iv)+f_init_wat_spc(ico,ispc)*cini_air_r(iconta)
     flx_atm_t=flx_atm_t+f_atm_spc(ico,ispc)
     frain_t=frain_t+f_rain_spc(ico,ispc)
     finit_watt=finit_watt+f_init_wat_spc(ico,ispc)
     finit_sedt=finit_sedt+f_init_sed_spc(ico,ispc)
     l_varnew=.false.
     DO ispc=nb_spc_diss(ico)+2,nb_spc_tot(ico)
       IF(ispc>nb_spc_diss(ico)+nb_spc_part_eq(ico).or..not.l_gpspecies )then
        l_varnew=.true.
       ELSE
        DO ispcp=nb_spc_diss(ico)+1,ispc-1
         IF(isit_spc(ico,ispc) == isit_spc(ico,ispcp)) then
          ivar_spc(ico,ispc)=ivar_spc(ico,ispcp)
          iv=ivar_spc(ico,ispc)
          nb_frac_var_r(iv)=nb_frac_var_r(iv)+1
          ifracv=nb_frac_var_r(iv)
          ifrac_var_spc(ico,ispc)=ifracv
          ispc_var_r(iv,ifracv)=ispc
          icon_var_r(iv)=ico
          finit_watt=finit_watt+f_init_wat_spc(ico,ispc)
          finit_sedt=finit_sedt+f_init_sed_spc(ico,ispc)
          l_varnew=.false.
          cini_wat_r(iv)=cini_wat_r(iv)+f_init_wat_spc(ico,ispc)*cini_wat_r(iconta)
          cini_sed_r(iv)=cini_sed_r(iv)+f_init_sed_spc(ico,ispc)*cini_sed_r(iconta)
          cini_air_r(iv)=cini_air_r(iv)+f_init_wat_spc(ico,ispc)*cini_air_r(iconta)
          flx_atm_r(iv)=flx_atm_r(iv)+f_atm_spc(ico,ispc)*flx_atm_r(iconta)
          cv_rain_r(iv)=cv_rain_r(iv)+f_rain_spc(ico,ispc)*cv_rain_r(iconta)
          flx_atm_t=flx_atm_t+f_atm_spc(ico,ispc)
          frain_t=frain_t+f_rain_spc(ico,ispc)
          l_varnew=.false.
          EXIT
         ELSE
          l_varnew=.true.
         ENDIF
        ENDDO
       ENDIF
       IF(l_varnew)then
        ivar_spc(ico,ispc)=isubs+1
        isubs=isubs+1
        iv=ivar_spc(ico,ispc)
        ifracv=1
        nb_frac_var_r(iv)=ifracv
        ifrac_var_spc(ico,ispc)=ifracv
        ispc_var_r(iv,ifracv)=ispc
        icon_var_r(iv)=ico
        IF(itypv(isit_spc(ico,ispc))==5) THEN   ! variable substrat is a SORB
          itypv_r(iv)=5
          nv_sorb_co=nv_sorb_co+1       
          name_varpc_assoc(iv)=name_var(isit_spc(ico,ispc))
          write(num_part,'(i1)')nv_sorb_co
          name_var(iv)=trim(name_var(iconta))//'_Ps'//num_part
        ELSE IF (itypv(isit_spc(ico,ispc))==4) THEN   ! variable substrat is a NoCP
          itypv_r(iv)=4
          nv_ncp_co=nv_ncp_co+1   
          write(num_part,'(i1)')nv_ncp_co 
          name_var(iv)=trim(name_var(iconta))//'_Pn'//num_part   
          ws_spc_free_min(iv)=ws_free_min_r(isit_spc(ico,ispc))
          ws_spc_free_max(iv)=ws_free_max_r(isit_spc(ico,ispc))
          ws_spc_free_para(:,iv)=ws_free_para_r(:,isit_spc(ico,ispc))
          ws_spc_hind_para(:,iv)=ws_hind_para_r(:,isit_spc(ico,ispc))
          ws_spc_free_opt(iv)=ws_free_opt_r(isit_spc(ico,ispc))
          ws_spc_hind_opt(iv)=ws_hind_opt_r(isit_spc(ico,ispc))
        ENDIF
        l_out_subs_r(iv)=.true.
        unit_var(iv)=trim(unit_var(iconta))
        init_cv_name_r(iv)=init_cv_name_r(iconta)
        obc_cv_name_r(iv)='none'
        tocd_r(iv)=tocd_spc(isit_spc(ico,ispc))
        ros_r(iv)=ros_spc(isit_spc(ico,ispc))
        diam_r(iv)=diam_spc(isit_spc(ico,ispc))
        flx_atm_r(iv)=flx_atm_r(ivar_spc(ico,ispc))+&
                      f_atm_spc(ico,ispc)*flx_atm_r(iconta)
        cv_rain_r(iv)=cv_rain_r(iv)+f_rain_spc(ico,ispc)*cv_rain_r(iconta)
        frac_cwat_r(iv,ifracv)=f_init_wat_spc(ico,ispc)
        frac_csed_r(iv,ifracv)=f_init_sed_spc(ico,ispc)
        cini_wat_r(iv)=cini_wat_r(iv)+f_init_wat_spc(ico,ispc)*cini_wat_r(iconta)
        cini_sed_r(iv)=cini_sed_r(iv)+f_init_sed_spc(ico,ispc)*cini_sed_r(iconta)
        cini_air_r(iv)=cini_air_r(iv)+f_init_wat_spc(ico,ispc)*cini_air_r(iconta)
        flx_atm_t=flx_atm_t+f_atm_spc(ico,ispc)
        frain_t=frain_t+f_rain_spc(ico,ispc)
        finit_watt=finiwatt_t+f_init_wat_spc(ico,ispc)
        finit_sedt=finisedt_t+f_init_sed_spc(ico,ispc)
      ENDIF
     ENDDO
   ENDIF
   IF(l_gpspecies)THEN
    DO iv=nv_tot+1,nv_tot+nv_ncp_co+nv_sorb_co
      fspcw_tot=0.
      fspcs_tot=0.
      IF(nb_frac_var_r(iv) == 1)then
        frac_cwat_r(iv,1)=1.
        frac_csed_r(iv,1)=1.
      ELSE  
       DO ispc=1,nb_frac_var_r(iv)-1
         frac_cwat_r(iv,ispc)=f_init_wat_spc(ico,ispc_var_r(iv,ispc))&
                              *cini_wat_r(iconta)/(cini_wat_r(iv)+epsilon)
         frac_csed_r(iv,ispc)=f_init_sed_spc(ico,ispc_var_r(iv,ispc))&
                              *cini_sed_r(iconta)/(cini_sed_r(iv)+epsilon)
         fspcw_tot=fspcw_tot+frac_cwat_r(iv,ispc)
         fspcs_tot=fspcs_tot+frac_csed_r(iv,ispc)
       ENDDO
       frac_cwat_r(iv,nb_frac_var_r(iv))=1.-fspcw_tot
       frac_csed_r(iv,nb_frac_var_r(iv))=1.-fspcs_tot
     ENDIF
    ENDDO
   ENDIF
   nv_ncp=nv_ncp_co+nv_ncp
   nv_sorb=nv_sorb_co+nv_sorb              
   
!!  Grouping dissolved species for set up the dissolved advected variable number "iconta" and suivantes
!!=========================================================================================
   IF(nb_spc_diss(ico).ne.0)then
     itypv_r(iconta)=6
     ivar_spc(ico,1)=iconta
     ifrac_var_spc(ico,1)=1
     icon_var_r(iconta)=ico
     init_cv_name_r(iconta)=init_cv_name_r(iconta)
     obc_cv_name_r(iconta)='none'
     flx_atm_r(iconta)=(1-flx_atm_t)*flx_atm_r(iconta)
     cv_rain_r(iconta)=(1-frain_t)*cv_rain_r(iconta)
     ifracv=1
     ispc_var_r(iconta,ifracv)=1
     f_init_wat_spc(ico,ifracv)=1-finit_watt
     cini_wat_r(iconta)=(1-finit_watt)*cini_wat_r(iconta)
     f_init_sed_spc(ico,ifracv)=1-finit_sedt
     cini_sed_r(iconta)=(1-finit_sedt)*cini_sed_r(iconta)
     IF (l_gpspecies)THEN
      nb_frac_var_r(iconta)=nb_spc_diss(ico)
      name_var(iconta)=trim(name_var(iconta))//'_D'
      frac_cwat_r(iconta,ifracv)=1.
      frac_csed_r(iconta,ifracv)=1.
      DO ispc=2,nb_spc_diss(ico)
       ifracv=ifracv+1
       ivar_spc(ico,ispc)=iconta
       ifrac_var_spc(ico,ispc)=ifracv
       ispc_var_r(iconta,ifracv)=ifracv
       frac_cwat_r(iconta,ifracv)=0.
       frac_csed_r(iconta,ifracv)=0.
       f_init_wat_spc(ico,ifracv)=0.
       f_init_sed_spc(ico,ifracv)=0.
      ENDDO
     ELSE
      nb_frac_var_r(iconta)=1
      name_var(iconta)=trim(name_spc(ico,1))
      DO ispc=2,nb_spc_diss(ico)
        isubs=isubs+1
        l_out_subs_r(isubs)=.true.
        itypv_r(isubs)=6
        ivar_spc(ico,ispc)=isubs
        ifrac_var_spc(ico,ispc)=1
        nb_frac_var_r(isubs)=1
        icon_var_r(isubs)=ico
        name_var(isubs)=trim(name_spc(ico,ispc))
        init_cv_name_r(isubs)=init_cv_name_r(iconta)
        obc_cv_name_r(isubs)=name_var(isubs)
        flx_atm_r(isubs)=0.
        cv_rain_r(isubs)=0.
        ifracv=1
        ispc_var_r(isubs,ifracv)=1
        frac_cwat_r(isubs,ifracv)=1.
        frac_csed_r(isubs,ifracv)=1.
        f_init_wat_spc(ico,ifracv)=0.
        f_init_sed_spc(ico,ifracv)=0.
        cini_wat_r(isubs)=0.
        cini_sed_r(isubs)=0.
        nv_dis=nv_dis+1
      ENDDO
    ENDIF  
   ELSE
     nv_dis=nv_dis-1
   ENDIF
    
!  UNE ESPECE / VARIABLE (PAS DE KEY_SPECIES)
! relation  especes / variable
!    icon_var_r(iv) : numero du contaminant correspondant a l espece dans la variable iv
!     nb_frac_var_r(iv)=1   : une seule espece par variable iv (somme des fractions = 1)
#else
   nb_frac_var_r(:)=1
   nvtot_anc=nv_tot
   nv_ncp_co=0
   nv_sorb_co=0
   flx_atm_t=0.
   frain_t=0.
   finit_watt=0.
   finit_sedt=0.
   IF(nb_spc_part_tot(ico) /= 0)THEN
     ispc=nb_spc_diss(ico)+1
     IF(ispc == 1)Then
       isubs=nv_tot
       iv=iconta
     ELSE
       isubs=nv_tot+1  
       iv=nv_tot+1
     ENDIF
     icon_var_r(iv)=ico
     ivar_spc(ico,ispc)=iv
     ispc_var_r(iv)=ispc
     l_out_subs_r(iv)=.true.
     name_var(iv)=trim(name_spc(ico,ispc))
     unit_var(iv)=trim(unit_var(iconta))
     init_cv_name_r(iv)=init_cv_name_r(iconta)
     obc_cv_name_r(iv)='none'
     IF(itypv(isit_spc(ico,ispc))==5) THEN   ! variable substrat is a SORB
          itypv_r(iv)=5
          nv_sorb_co=nv_sorb_co+1       
          name_varpc_assoc(iv)=name_var(isit_spc(ico,ispc))
     ELSE IF (itypv(isit_spc(ico,ispc))==4) THEN   ! variable substrat is a NoCP
          itypv_r(iv)=4
          nv_ncp_co=nv_ncp_co+1   
          ws_free_min_r(iv)=ws_free_min_r(isit_spc(ico,ispc))
          ws_free_max_r(iv)=ws_free_max_r(isit_spc(ico,ispc))
          ws_free_para_r(:,iv)=ws_free_para_r(:,isit_spc(ico,ispc))
          ws_hind_para_r(:,iv)=ws_hind_para_r(:,isit_spc(ico,ispc))
          ws_free_opt_r(iv)=ws_free_opt_r(isit_spc(ico,ispc))
          ws_hind_opt_r(iv)=ws_hind_opt_r(isit_spc(ico,ispc))
     ENDIF
     tocd_r(iv)=tocd_spc(ispc)
     ros_r(iv)=ros_spc(ispc)
     diam_r(iv)=diam_spc(ispc)
     flx_atm_r(iv)=flx_atm_r(iv)+f_atm_spc(ico,ispc)*flx_atm_r(iconta)
     cv_rain_r(iv)=cv_rain_r(iv)+f_rain_spc(ico,ispc)*cv_rain_r(iconta)
     cini_wat_r(iv)=cini_wat_r(iv)+f_init_wat_spc(ico,ispc)*cini_wat_r(iconta)
     cini_sed_r(iv)=cini_sed_r(iv)+f_init_sed_spc(ico,ispc)*cini_sed_r(iconta)
     cini_air_r(iv)=cini_air_r(iv)+f_init_wat_spc(ico,ispc)*cini_air_r(iconta)
     flx_atm_t=flx_atm_t+f_atm_spc(ico,ispc)
     frain_t=frain_t+f_rain_spc(ico,ispc)
     finit_watt=finit_watt+f_init_wat_spc(ico,ispc)
     finit_sedt=finit_sedt+f_init_sed_spc(ico,ispc)
     DO ispc=nb_spc_diss(ico)+2,nb_spc_tot(ico)
        nv_ncp_co=nv_ncp_co+1       
        iv=isubs+1
        isubs=isubs+1
        icon_var_r(iv)=ico
        ispc_var_r(iv)=ispc
        ivar_spc(ico,ispc)=iv
        itypv_r(iv)=4
        l_out_subs_r(iv)=.true.
        name_var(iv)=trim(name_spc(ico,ispc))
        unit_var(iv)=trim(unit_var(iconta))
        init_cv_name_r(iv)=init_cv_name_r(iconta)
        obc_cv_name_r(iv)='none'
        IF(itypv(isit_spc(ico,ispc))==5) THEN   ! variable substrat is a SORB
          itypv_r(iv)=5
          nv_sorb_co=nv_sorb_co+1       
          name_varpc_assoc(iv)=name_var(isit_spc(ico,ispc))
        ELSE IF (itypv(isit_spc(ico,ispc))==4) THEN   ! variable substrat is a NoCP
          itypv_r(iv)=4
          nv_ncp_co=nv_ncp_co+1   
          ws_free_min_r(iv)=ws_free_min_r(isit_spc(ico,ispc))
          ws_free_max_r(iv)=ws_free_max_r(isit_spc(ico,ispc))
          ws_free_para_r(:,iv)=ws_free_para_r(:,isit_spc(ico,ispc))
          ws_hind_para_r(:,iv)=ws_hind_para_r(:,isit_spc(ico,ispc))
          ws_free_opt_r(iv)=ws_free_opt_r(isit_spc(ico,ispc))
          ws_hind_opt_r(iv)=ws_hind_opt_r(isit_spc(ico,ispc))
        ENDIF
        tocd_r(iv)=tocd_spc(ispc)
        ros_r(iv)=ros_spc(ispc)
        diam_r(iv)=diam_spc(ispc)
        flx_atm_r(iv)=flx_atm_r(iv)+f_atm_spc(ico,ispc)*flx_atm_r(iconta)
        cv_rain_r(iv)=cv_rain_r(iv)+f_rain_spc(ico,ispc)*cv_rain_r(iconta)
        cini_wat_r(iv)=cini_wat_r(iv)+f_init_wat_spc(ico,ispc)*cini_wat_r(iconta)
        cini_sed_r(iv)=cini_sed_r(iv)+f_init_sed_spc(ico,ispc)*cini_sed_r(iconta)
        cini_air_r(iv)=cini_air_r(iv)+f_init_wat_spc(ico,ispc)*cini_air_r(iconta)
        flx_atm_t=flx_atm_t+f_atm_spc(ico,ispc)
        frain_t=frain_t+f_rain_spc(ico,ispc)
        finit_watt=finit_watt+f_init_wat_spc(ico,ispc)
        finit_sedt=finit_sedt+f_init_sed_spc(ico,ispc)
     ENDDO
   ENDIF
   nv_ncp=nv_ncp_co+nv_ncp 
   nv_sorb=nv_sorb_co+nv_sorb             
    
!!  Dissolved species for set up the dissolved advected variable number iconta and suivantes
!!=========================================================================================
   IF(nb_spc_diss(ico).ne.0)then
     itypv_r(iconta)=6
     ivar_spc(ico,1)=iconta
     icon_var_r(iconta)=ico
     ispc_var_r(iconta)=1
     init_cv_name_r(iconta)=init_cv_name_r(iconta)
     obc_cv_name_r(iconta)='none'
     flx_atm_r(iconta)=(1-flx_atm_t)*flx_atm_r(iconta)
     cv_rain_r(iconta)=(1-frain_t)*cv_rain_r(iconta)
     f_init_wat_spc(ico,1)=1-finit_watt
     cini_wat_r(iconta)=(1-finit_watt)*cini_wat_r(iconta)
     f_init_sed_spc(ico,1)=1-finit_sedt
     cini_sed_r(iconta)=(1-finit_sedt)*cini_sed_r(iconta)
     name_var(iconta)=trim(name_spc(ico,1))
     DO ispc=2,nb_spc_diss(ico)
        isubs=isubs+1
        l_out_subs_r(isubs)=.true.
        itypv_r(isubs)=6
        ispc_var_r(isubs)=ispc
        ivar_spc(ico,ispc)=isubs
        icon_var_r(isubs)=ico
        name_var(isubs)=trim(name_spc(ico,ispc))
        unit_var(isubs)=trim(unit_var(iconta))
        init_cv_name_r(isubs)=init_cv_name_r(iconta)
        obc_cv_name_r(isubs)='none'
        flx_atm_r(isubs)=0.
        cv_rain_r(isubs)=0.
        f_init_wat_spc(ico,ispc)=0.
        f_init_sed_spc(ico,ispc)=0.
        cini_wat_r(isubs)=0.
        cini_sed_r(isubs)=0.
        nv_dis=nv_dis+1
      ENDDO
   ELSE
     nv_dis=nv_dis-1
   ENDIF
#endif

  ! reactualisation des nouveaux nombres de variable
  ! (une espece dissoute et 1 ou plusieurs especes particulaires selon la diversite des substrats)
        nvp=nvpc+nv_ncp+nv_sorb
        nv_adv=nvp+nv_dis
        nv_state=nv_adv+nv_fix
        nv_tot=nv_state+nv_dri+nv_int
! fin boucle sur ico
  ENDDO
  
  END SUBROUTINE conta_read_data

   !!======================================================================

  SUBROUTINE conta_read_react

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE conta_read_react  ***
   !&E
   !&E ** Purpose : Reading parameters for biochemical reactions 
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : main.F90
   !&E
   !&E ** used ij-arrays : none
   !&E
   !&E ** modified variables : none
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2008-06  (B. Thouvenin)  Original code
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
!!   USE comvars2d, ONLY : iscreenlog        
    USE comsubstance, ONLY : name_filesubs,nv_state,nv_tot,name_var
     
   !! * Arguments

   !! * Local declarations
   LOGICAL                                     :: ex,l_tmp
   INTEGER                                     :: it,nblu,nbvaractmax,ico,ispc
   INTEGER                                     :: ir,ireac,ndgmax,irr,ll,ind_white
   CHARACTER(LEN=5)                            :: comment
   CHARACTER(LEN=10)                           :: commentlu
   CHARACTER(LEN=50)                           :: name_reacr
   CHARACTER(LEN=10)                           :: name_typ_reacr

   !!----------------------------------------------------------------------
   !! * Executable part
   
   ! save into simu.log
   !-------------------
   IF_MPI (MASTER) THEN
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '**************************************************'
   WRITE(iscreenlog,*) '**************** MET&OR_read_REACT.F90 ***************'
   WRITE(iscreenlog,*) '**************************************************'
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) 'file for definition of reactions : ',trim(filereaction)
   ENDIF_MPI

   eof=0
   nb_reac_pervar(:)=0
! reading file in order to compute the number of reactions
   INQUIRE(file= filereaction,exist=ex)
   IF(ex) then
     OPEN(49,file = filereaction,form='formatted')
     comment='debut'
     DO while (comment /= '*****')
      READ(49,'(a)',iostat=eof)comment
     ENDDO
     ireac=0
     DO while(eof==0)
       READ(49,'(a)',iostat=eof)name_reacr
       IF(eof==0) then
        READ(49,'(a)',iostat=eof)comment
        DO while (comment /= '=====' .AND. comment /= '****')
          READ(49,'(a)',iostat=eof)comment
        ENDDO
        ireac=ireac+1
       ENDIF
     ENDDO           ! fin boucle de lecture des reactions jusqu a la fin du fichier 
   ENDIF
   CLOSE(49)
   nreacmax_sco=ireac
   
 ! adding reactions from contaminant
  IF(l_conta)then
   DO ico=1,nb_conta
    DO ispc=1,nb_spc_tot(ico)
     IF(cin_hydrol(ico,ispc) /= 0.)then
      ireac=ireac+1
     ENDIF
     IF(cin_biodeg(ico,ispc) /= 0.)then
      ireac=ireac+1
     ENDIF
     IF(cin_photo(ico,ispc) /= 0.)then
      ireac=ireac+1
     ENDIF
    ENDDO
    ireac=ireac+nb_reac_genco(ico)+nb_spc_part_ci(ico)
    IF(l_volatil(ico))ireac=ireac+1
    IF(l_partag(ico))ireac=ireac+1
   ENDDO
  ENDIF
   nreacmax=ireac
   l_treat_re_dc=.false.
   l_treat_re_after=.false.
  IF(nreacmax /= 0)then
  
! allocation of arrays with nreacmax =real max number of reactions
   ALLOCATE(ityp_reac(nreacmax))
   ALLOCATE(nb_var_re(nreacmax))
   ALLOCATE(ityp_arh(nreacmax))
   ALLOCATE(ire_typ(nreacmax))
   ALLOCATE(name_reac(nreacmax))
   ALLOCATE(cin_max_re(nreacmax))
   ALLOCATE(cf_teta_arh(nreacmax))
   ALLOCATE(cf_xetsig_arh(nreacmax))
   ALLOCATE(tempref_arh(nreacmax))
   ALLOCATE(l_arh(nreacmax))
   ALLOCATE(l_treat_re(nreacmax))
   ALLOCATE(l_cincte_re(nreacmax))
   ALLOCATE(l_coefcte_re(nreacmax))
   ALLOCATE(nb_flim1_re(nreacmax))
   ALLOCATE(nb_flim2_re(nreacmax))
   ALLOCATE(nb_flim3_re(nreacmax))
   ALLOCATE(ndg_kds_r(nreacmax))
   ALLOCATE(cf_pce_r(nreacmax,2))
   ALLOCATE(cf_salt_r(nreacmax,2))
   ALLOCATE(cf_kd0_r(nreacmax))
   ALLOCATE(cf_rate_r(nreacmax))
   ALLOCATE(cv_watequ_r(nreacmax))
   ALLOCATE(cf_xkds_r(nreacmax,5))
   ALLOCATE(l_pce_r(nreacmax))
   ALLOCATE(l_kds_r(nreacmax))
   ALLOCATE(l_salt_r(nreacmax))
   ALLOCATE(l_distb_spc(nreacmax))
   ALLOCATE(icon_reac(nreacmax))
     ! temporary allocation for max nb-active variables =5
   nvaractmax=2   
   ALLOCATE(irkm_var_re(nreacmax,nvaractmax))
   ALLOCATE(irkf_var_re(nreacmax,nvaractmax))
   ALLOCATE(irkm_sub_r(nreacmax,nvaractmax))
   ALLOCATE( name_var_re(nreacmax,nvaractmax))
   ALLOCATE( cf_stochio_re(nreacmax,nvaractmax))
   ALLOCATE( cf_expovar_re(nreacmax,nvaractmax))
   ALLOCATE( cf_distb_r(nreacmax,nvaractmax))
   ALLOCATE(ispc_reac(nreacmax,nvaractmax))
#ifdef key_species
   ALLOCATE(ifrac_var_reac(nreacmax,nvaractmax))
#endif
     ! temporary allocation for max nb-limiting functions =5
   nflimax=2
   ALLOCATE(irkm_flim1_re(nreacmax,nflimax))
   ALLOCATE(irkm_flim2_re(nreacmax,nflimax))
   ALLOCATE(cf_flim1_re(nreacmax,nflimax))
   ALLOCATE(cf_flim2_re(nreacmax,nflimax))
   ALLOCATE(cf_flim3_re(nreacmax,nflimax))
   l_flim3_re=.false.
   
   eof=0
   OPEN(49,file = filereaction,form='formatted')
   comment='debut'
   DO while (comment /= '*****')
     READ(49,'(a)',iostat=eof)comment
   ENDDO
   ireac=0
   nb_reac(:)=0
   DO while(eof==0)
     READ(49,'(a)',iostat=eof)name_reacr
     IF(eof==0) then
        READ(49,'(a)',iostat=eof)commentlu
        name_typ_reacr=trim(adjustl(adjustr(commentlu)))
        ll=len_trim(name_typ_reacr)
!  name_typ_reacr longueur de chaines=10
        ireac=ireac+1
        IF(name_typ_reacr(1:ll) ==      'REAC_1')then
         ityp_reac(ireac)=1
         nb_reac(1)=nb_reac(1)+1
        ELSE IF(name_typ_reacr(1:ll) == 'GENERIC')then
         ityp_reac(ireac)=2
         nb_reac(2)=nb_reac(2)+1
        ELSE IF(name_typ_reacr(1:ll) == 'REVERSIBLE')then
         ityp_reac(ireac)=3
         nb_reac(3)=nb_reac(3)+1
        ELSE IF(name_typ_reacr(1:ll) == 'VOLATIL')then
         ityp_reac(ireac)=4
         nb_reac(4)=nb_reac(4)+1
        ELSE IF(name_typ_reacr(1:ll) == 'PARTAG_EQ')then
         ityp_reac(ireac)=5
         nb_reac(5)=nb_reac(5)+1
         l_distb_spc(nb_reac(5))=.false.
        ELSE IF(name_typ_reacr(1:ll) == 'MOCO')then
         ityp_reac(ireac)=6
         nb_reac(6)=nb_reac(6)+1
        ELSE IF(name_typ_reacr(1:ll) == 'OTHER')then
         ityp_reac(ireac)=7
         nb_reac(7)=nb_reac(7)+1
        ENDIF
        ind_white=index(name_reacr,' ')
        name_reac(ireac)=name_reacr(1:ind_white)
        READ(49,*,iostat=eof)l_treat_re(ireac)
    ! numero du contaminant et de l espece de contaminant concerne, ici=0 
        ispc_reac(ireac,:)=0
        icon_reac(ireac)=0
        
        SELECT CASE(ityp_reac(ireac))
        
         CASE (1)
          CALL conta_read_param_REAC1(ireac)
          ire_typ(ireac)=nb_reac(1)
          
         CASE (2)
          CALL conta_read_param_GEN(ireac)
          ire_typ(ireac)=nb_reac(2)
          
         CASE (3)
          CALL conta_read_param_REV(ireac)
          ire_typ(ireac)=nb_reac(3)
           
         CASE (4)
          CALL conta_read_param_VOL(ireac)
          ire_typ(ireac)=nb_reac(4)
          
         CASE (5)
          CALL conta_read_param_PART(ireac)
          ire_typ(ireac)=nb_reac(5)
          
         CASE (6)
!          CALL conta_read_param_MOCO(ireac)
          ire_typ(ireac)=nb_reac(6)
          
!         CASE (7)
!          CALL conta_read_param_OTH(ireac)
        END SELECT 
     ENDIF
   ENDDO           ! fin boucle de lecture des reactions jusqu a la fin du fichier 
11 CONTINUE
   CLOSE(49)
   
! if we simulate some contaminant defined by chemical CHARACTERistics
! addition of correspondant reactions
   IF(l_conta)then
      CALL conta_add_reac
   ENDIF

! save into simu.log
!-------------------
   IF_MPI (MASTER) THEN
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '**************************************************'
   WRITE(iscreenlog,*) '**************** REACTIONS *********************'
   WRITE(iscreenlog,*) '**************************************************'
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*)'REACTIONS DEFINED IN FILE :',filereaction
   WRITE(iscreenlog,*)'il y a  ',nb_reac(1),' reactions de type 1  (REAC1)'
   WRITE(iscreenlog,*)'il y a  ',nb_reac(2),' reactions de type 2  (GENERIC)'
   WRITE(iscreenlog,*)'il y a  ',nb_reac(3),' reactions de type 3  (REVERSIBLE)'
   WRITE(iscreenlog,*)'il y a  ',nb_reac(4),' reactions de type 4  (VOLATIL)'
   WRITE(iscreenlog,*)'il y a  ',nb_reac(5),' reactions de type 5  (PARTAG_EQ)'
   WRITE(iscreenlog,*)'il y a  ',nb_reac(6),' reactions de type 6   (MOCO)'   
   WRITE(iscreenlog,*)'il y a  ',nb_reac(7),' reactions de type 7   (OTHER)'   
   WRITE(iscreenlog,*)

   DO it=1,nv_state
    WRITE(iscreenlog,*)
    WRITE(iscreenlog,*)' variable :',trim(name_var(it)),'is affected by ', nb_reac_pervar(it),' biogeochemical reactions '
   ENDDO
   ENDIF_MPI

  nbvaractmax=0
  ndgmax=0
  nbvaractmax=maxval(nb_var_re)
  ndgmax=maxval(ndg_kds_r(1:nb_reac(5)))
  IF(nb_reac(5).gt.0)then
    ALLOCATE(l_pce_re5(nb_reac(5)))
    ALLOCATE(l_salt_re5(nb_reac(5)))
    ALLOCATE(cf_pce_re5(nb_reac(5),2))
    ALLOCATE(cf_salt_re5(nb_reac(5),2))
    ALLOCATE(cf_distb_re5(nb_reac(5),nbvaractmax))
    ALLOCATE(irkm_sub_re5(nb_reac(5),nbvaractmax))
    ALLOCATE(l_kds_re5(nb_reac(5)))
    l_pce_re5(:)=l_pce_r(1:nb_reac(5))
    l_salt_re5(:)=l_salt_r(1:nb_reac(5))
    cf_pce_re5(:,:)=cf_pce_r(1:nb_reac(5),:)
    cf_salt_re5(:,:)=cf_salt_r(1:nb_reac(5),:)
    irkm_sub_re5(:,:)=irkm_sub_r(1:nb_reac(5),1:nbvaractmax)
    cf_distb_re5(:,:)=cf_distb_r(1:nb_reac(5),1:nbvaractmax)
    l_kds_re5(:)=l_kds_r(1:nb_reac(5))
    ALLOCATE(ndg_kds_re5(nb_reac(5)))
    ALLOCATE(cf_kd0_re5(nb_reac(5)))
    ALLOCATE(cf_xkds_re5(nb_reac(5),ndgmax))
    ndg_kds_re5(:)=ndg_kds_r(1:nb_reac(5))
    cf_kd0_re5(:)=cf_kd0_r(1:nb_reac(5))
    cf_xkds_re5(:,:)=cf_xkds_r(1:nb_reac(5),1:ndgmax)
  ENDIF
  
  IF(nb_reac(2).gt.0)then
    ALLOCATE(cin_re2_tmp(limin:limax,ljmin:ljmax))   
  ENDIF

  IF(nb_reac(3).gt.0)then
    ALLOCATE(cf_rate_re3(nb_reac(3)))
    ALLOCATE(cv_watequ_re3(nb_reac(3)))
    cf_rate_re3(:)=cf_rate_r(1:nb_reac(3))
    cv_watequ_re3(:)=cv_watequ_r(1:nb_reac(3))
  ENDIF
  
! sortie des flux reactifs entre variables
   IF(l_sorfluxreac)then
    ALLOCATE(flux_reac_wat(nreacmax,nbvaractmax,kmax,limin:limax,ljmin:ljmax))   
   ENDIF
    
! desallocation of areas which have been used only into this module
   if(l_conta)call conta_deALLOCATE_data
  ELSE
   IF_MPI (MASTER) THEN
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '**************************************************'
   WRITE(iscreenlog,*) '**************** NO REACTIONS *********************'
   WRITE(iscreenlog,*) '**************************************************'
   WRITE(iscreenlog,*) ' '
   ENDIF_MPI
  ENDIF
  DO ireac=1,nreacmax
! definition de l_treat_re_dc = true si il y a au moins une reaction traitee dans l equation
!   et de l_treat_re_after = true si il y a au moins une reaction traitee en time splitting
    IF(l_treat_re(ireac))then
            l_treat_re_dc=.true.
    ELSE
            l_treat_re_after=.true.
    ENDIF 
  ENDDO                  

  END SUBROUTINE conta_read_react
   !!======================================================================
      
  SUBROUTINE conta_read_param_REAC1(ir)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE conta_read_param_REAC1 ***
   !&E
   !&E ** Purpose : Reading parameters for biochemical reactions type REAC1
   !&E
   !&E ** Description :
   !&E
   !&E ** used ij-arrays : none
   !&E
   !&E ** modified variables : none
   !&E
   !&E ** Called by : conta_read_react.F90
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2008-06  (B. Thouvenin)  Original code
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
    USE comsubstance, ONLY : name_filesubs,irk_mod,name_var,nv_tot

   !! * Arguments
   INTEGER, intent(in)                         :: ir

   !! * Local declarations
   LOGICAL                                     :: ex,l_var_sim,l_unitKd
   INTEGER                                     :: ire,it,isubs,nfl1,ind_white
   INTEGER                                     :: iflim,nfunct_limit,  &
                                                  ityp_flim,nfl2,nfl3
   CHARACTER(LEN=5)                            :: comment
   CHARACTER(LEN=20)                           :: commentlu,namcop,       &
                                                    name_var_flim
   CHARACTER(LEN=10)                           :: name_typ_reacr
   REAL(KIND=rsh)                                   :: cf_flim_re
     
   !!----------------------------------------------------------------------
   !! * Executable part

!   =======================================================
!   Reaction  type 1
!   cinetic simple order 1 (degradation : dC/dt=µC)
!   une seule variable active
!   =======================================================
    nb_var_re(ir)=1
    IF(eof==0) then
      READ(49,'(a)',iostat=eof) commentlu     ! nom de la variable active
      ind_white=index(commentlu,' ')
      name_var_re(ir,1)=commentlu(1:ind_white)
      READ(49,*,iostat=eof)cin_max_re(ir)
      cin_max_re(ir)=cin_max_re(ir)/86400.
      READ(49,*,iostat=eof)cf_stochio_re(ir,1)
      READ(49,*,iostat=eof)l_arh(ir)
      READ(49,*,iostat=eof)ityp_arh(ir)
      READ(49,*,iostat=eof)cf_teta_arh(ir)
      READ(49,*,iostat=eof)tempref_arh(ir)
      READ(49,*,iostat=eof)cf_xetsig_arh(ir)
      DO isubs=1,nb_var
       IF(name_var_re(ir,1) == name_var(isubs)) then
         irkm_var_re(ir,1)=irk_mod(isubs)
         irkf_var_re(ir,1)=isubs
         nb_reac_pervar(isubs)=nb_reac_pervar(isubs)+1
       ENDIF
      ENDDO
      IF(irkm_var_re(ir,1).eq.0)then
        WRITE(*,*)'********ERROR ***********'
        WRITE(*,*)' The active variable of the reaction ',ir,' (',trim(name_reac(ir)),')'
        WRITE(*,*)'is not simulated (probably not the same name) !!'
        WRITE(*,*)'Simulated  variables  are :'
        DO isubs=1,nv_tot
          WRITE(*,*)name_var(isubs)
        ENDDO
        WRITE(*,*)' Modify her name :',name_var_re(ir,1)
        WRITE(*,*)' in the file of reaction parameters'
        STOP
      ENDIF
      READ(49,'(a)',iostat=eof)
    ELSE
      WRITE(*,*)'End of the file : ',filereaction
      WRITE(*,*)'in the reaction number : ',ir 
      STOP
    ENDIF
       
   
! save into simu.log
!-------------------
   IF_MPI (MASTER) THEN
   WRITE(iscreenlog,*)'REACTION : ',name_reac(ir),' TYPE :',ityp_reac(ir)
   WRITE(iscreenlog,*)'maxi cinetic ',cin_max_re(ir)
   IF(l_treat_re(ir)) then
      WRITE(iscreenlog,*)'integration of sink & sources due to this reaction '
      WRITE(iscreenlog,*)'in transport equations '
   ELSE
      WRITE(iscreenlog,*)'sink & sources due to this reaction computed '
      WRITE(iscreenlog,*)' after (operator splitting)'
   ENDIF
   WRITE(iscreenlog,*)'maxi cinetic ',cin_max_re(ir)
   WRITE(iscreenlog,*)'active variable ',name_var_re(ir,1)
   WRITE(iscreenlog,*)'     stochiometric coefficient  ',cf_stochio_re(ir,1)
   WRITE(iscreenlog,*)'     model s order number of this variable : ',irkm_var_re(ir,1)
   WRITE(iscreenlog,*)'     data file s order number of this variable  : ', irkf_var_re(ir,1)  
   ENDIF_MPI

  END SUBROUTINE conta_read_param_REAC1
  !!======================================================================

  SUBROUTINE conta_read_param_GEN(ir)
 

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE conta_read_param_GEN ***
   !&E
   !&E ** Purpose : Reading parameters for biochemical reactions type GENERIC
   !&E
   !&E ** Description :
   !&E
   !&E ** used ij-arrays : none
   !&E
   !&E ** modified variables : none
   !&E
   !&E ** Called by : conta_read_react.F90
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2008-06  (B. Thouvenin)  Original code
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
    USE comsubstance, ONLY : name_filesubs,irk_mod, name_var,nv_tot

   !! * Arguments
   INTEGER, intent(in)                           :: ir

   !! * Local declarations
   LOGICAL                                     :: l_var_sim
   INTEGER                                     :: nvaractmax_new, nflimax_new
   INTEGER                                     :: ire,it,isubs,nfl1,nblu
   INTEGER                                     :: ind_white
   INTEGER                                     :: iflim,nfunct_limit,  &
                                                  ityp_flim,nfl2,nfl3
   CHARACTER(LEN=5)                            :: comment
   CHARACTER(LEN=20)                           :: commentlu,namcop,name_var_flim
   CHARACTER(LEN=10)                           :: name_typ_reacr
   REAL(KIND=rsh)                              :: cf_flim_re
   INTEGER,dimension(:,:),allocatable               :: irkm_var_tmp,irkf_var_tmp
   CHARACTER(len=lchain),dimension(:,:),allocatable :: name_var_tmp
   REAL(KIND=rsh),dimension(:,:),allocatable        :: cf_expovar_tmp
   REAL(KIND=rsh),dimension(:,:),allocatable        :: cf_stochio_tmp
   INTEGER,dimension(:,:),allocatable               :: irkm_flim1_tmp
   INTEGER,dimension(:,:),allocatable               :: irkm_flim2_tmp
   REAL(KIND=rsh),dimension(:,:),allocatable        :: cf_flim1_tmp
   REAL(KIND=rsh),dimension(:,:),allocatable        :: cf_flim2_tmp
   REAL(KIND=rsh),dimension(:,:),allocatable        :: cf_flim3_tmp
     
   !!----------------------------------------------------------------------
   !! * Executable part
   
!   =======================================================
!   Reaction of type 2
!   generic : dC(iv)/dt=cf_stochio * maxi cinetic * limiting functions * C(iv1)^n * C(iv2)^m
!   =======================================================

! reALLOCATE arrays if nvaractmax too small
             IF(l_conta .and. nb_var_re(ir).gt.nvaractmax) then
               ALLOCATE(name_var_tmp(nreacmax,nvaractmax))
               ALLOCATE(irkm_var_tmp(nreacmax,nvaractmax))
               ALLOCATE(irkf_var_tmp(nreacmax,nvaractmax))
               ALLOCATE(cf_stochio_tmp(nreacmax,nvaractmax))
               ALLOCATE(cf_expovar_tmp(nreacmax,nvaractmax))
               name_var_tmp(:,:)=name_var_re(:,:)        
               irkm_var_tmp(:,:)=irkm_var_re(:,:) 
               irkf_var_tmp(:,:)=irkf_var_re(:,:) 
               cf_stochio_tmp(:,:)=cf_stochio_re(:,:) 
               cf_expovar_tmp(:,:)=cf_expovar_re(:,:) 
               deALLOCATE(name_var_re)
               deALLOCATE(irkm_var_re)
               deALLOCATE(irkf_var_re)
               deALLOCATE(cf_stochio_re)
               deALLOCATE(cf_expovar_re)
               nvaractmax_new=nb_var_re(ir)
               ALLOCATE(name_var_re(nreacmax,nvaractmax_new))
               ALLOCATE(irkm_var_re(nreacmax,nvaractmax_new))
               ALLOCATE(irkf_var_re(nreacmax,nvaractmax_new))
               ALLOCATE(cf_stochio_re(nreacmax,nvaractmax_new))
               ALLOCATE(cf_expovar_re(nreacmax,nvaractmax_new))
               name_var_re(:,1:nvaractmax)=name_var_tmp(:,:)        
               irkm_var_re(:,1:nvaractmax)=irkm_var_tmp(:,:)        
               irkf_var_re(:,1:nvaractmax)=irkf_var_tmp(:,:)        
               cf_stochio_re(:,1:nvaractmax)=cf_stochio_tmp(:,:)        
               cf_expovar_re(:,1:nvaractmax)=cf_expovar_tmp(:,:)        
               nvaractmax=nvaractmax_new                      
               deALLOCATE(name_var_tmp)
               deALLOCATE(irkm_var_tmp)
               deALLOCATE(irkf_var_tmp)
               deALLOCATE(cf_stochio_tmp)
               deALLOCATE(cf_expovar_tmp)
             ENDIF

    IF(eof==0) then
     READ(49,*,iostat=eof)nblu
     nb_var_re(ir)=nblu
     READ(49,*,iostat=eof)cin_max_re(ir)
     cin_max_re(ir)=cin_max_re(ir)/86400.     
     READ(49,*,iostat=eof)l_cincte_re(ir)
     READ(49,*,iostat=eof)l_coefcte_re(ir)
     READ(49,*,iostat=eof)nfunct_limit
! law of arhenius
      READ(49,*,iostat=eof)l_arh(ir)
      READ(49,*,iostat=eof)ityp_arh(ir)
      READ(49,*,iostat=eof)cf_teta_arh(ir)
      READ(49,*,iostat=eof)tempref_arh(ir)
      READ(49,*,iostat=eof)cf_xetsig_arh(ir)
      READ(49,'(a)',iostat=eof)
! reactives variables 
      DO it=1,nb_var_re(ir)
        READ(49,'(a)',iostat=eof) commentlu     ! name of reactives variables  (affected by the reaction)
        ind_white=index(commentlu,' ')
        name_var_re(ir,it)=commentlu(1:ind_white)
        l_var_sim=.false.
        DO isubs=1,nb_var
         IF(name_var_re(ir,it) == name_var(isubs)) then
           irkm_var_re(ir,it)=irk_mod(isubs)
           irkf_var_re(ir,it)=isubs
           nb_reac_pervar(isubs)=nb_reac_pervar(isubs)+1
           l_var_sim=.true.
           EXIT
         ENDIF
        ENDDO
        IF(.not.l_var_sim) then
!          si l_var_sim=.false.
          WRITE(*,*)'********ERROR ***********'
          WRITE(*,*)' The reactive variable of the reaction ',ir,' (',trim(name_reac(ir)),')'
          WRITE(*,*)'is not simulated (probably not the same name) !!'
          WRITE(*,*)'Simulated  variables  are :'
          DO isubs=1,nv_tot
            WRITE(*,*)name_var(isubs)
          ENDDO
          WRITE(*,*)' Modify her name :',name_var_re(ir,it)
          WRITE(*,*)' in the file of reaction parameters by comparison with' ,trim(name_filesubs)
          STOP
        ENDIF
        READ(49,*,iostat=eof)cf_stochio_re(ir,it)
        READ(49,*,iostat=eof)cf_expovar_re(ir,it)
      ENDDO
      READ(49,'(a)',iostat=eof)
! Limiting functions 
      nb_flim1_re(ir)=0
      nb_flim2_re(ir)=0
      nb_flim3_re(ir)=0
      IF(nfunct_limit.gt.0) then
         DO iflim=1,nfunct_limit
           READ(49,*,iostat=eof)ityp_flim
           READ(49,*,iostat=eof)cf_flim_re
           SELECT CASE (ityp_flim)
            CASE (1)
   ! Limiting functions type Monod 
             READ(49,*,iostat=eof)commentlu
             ind_white=index(commentlu,' ')
             name_var_flim=commentlu(1:ind_white)
             nb_flim1_re(ir)=nb_flim1_re(ir)+1
             nfl1=nb_flim1_re(ir)
       ! reALLOCATE arrays if nflimax too small
             IF(nb_flim1_re(ir).gt.nflimax) then
               ALLOCATE(irkm_flim1_tmp(nreacmax,nflimax))
               ALLOCATE(cf_flim1_tmp(nreacmax,nflimax))
               irkm_flim1_tmp(:,:)=irkm_flim1_re(:,:)        
               cf_flim1_tmp(:,:)=cf_flim1_re(:,:) 
               deALLOCATE(irkm_flim1_re)       
               deALLOCATE(cf_flim1_re)
               nflimax_new=nb_flim1_re(ir)
               ALLOCATE(irkm_flim1_re(nreacmax,nflimax_new))     
               ALLOCATE(cf_flim1_re(nreacmax,nflimax_new))
               irkm_flim1_re(:,1:nflimax)=irkm_flim1_tmp(:,:)        
               cf_flim1_re(:,1:nflimax)=cf_flim1_tmp(:,:)
               nflimax=nflimax_new                      
               deALLOCATE(irkm_flim1_tmp)
               deALLOCATE(cf_flim1_tmp)
             ENDIF
             cf_flim1_re(ir,nfl1)=cf_flim_re
             l_var_sim=.false.
             DO isubs=1,nb_var
               IF(name_var_flim == name_var(isubs)) then
                  irkm_flim1_re(ir,nfl1)=irk_mod(isubs)
                  l_var_sim=.true.
               ELSE IF(name_var_flim == 'sal'.or.name_var_flim == 'salinite')then
                  irkm_flim1_re(ir,nfl1)=-1
                  l_var_sim=.true.
               ENDIF
             ENDDO
           CASE (2) 
   ! Limiting functions type Inhibitrice 
             READ(49,*,iostat=eof)commentlu
             ind_white=index(commentlu,' ')
             name_var_flim=commentlu(1:ind_white)
             nb_flim2_re(ir)=nb_flim2_re(ir)+1
             nfl2=nb_flim2_re(ir)
       ! reALLOCATE arrays if nflimax too small
             IF(nb_flim2_re(ir).gt.nflimax) then
               ALLOCATE(irkm_flim2_tmp(nreacmax,nflimax))
               ALLOCATE(cf_flim2_tmp(nreacmax,nflimax))
               irkm_flim2_tmp(:,:)=irkm_flim2_re(:,:)        
               cf_flim2_tmp(:,:)=cf_flim2_re(:,:) 
               deALLOCATE(irkm_flim2_re)       
               deALLOCATE(cf_flim2_re)
               nflimax_new=nb_flim2_re(ir)
               ALLOCATE(irkm_flim2_re(nreacmax,nflimax_new))     
               ALLOCATE(cf_flim2_re(nreacmax,nflimax_new))
               irkm_flim2_re(:,1:nflimax)=irkm_flim2_tmp(:,:)        
               cf_flim2_re(:,1:nflimax)=cf_flim2_tmp(:,:)
               nflimax=nflimax_new                      
               deALLOCATE(irkm_flim2_tmp)
               deALLOCATE(cf_flim2_tmp)
             ENDIF
             cf_flim2_re(ir,nfl2)=cf_flim_re
             l_var_sim=.false.
             DO isubs=1,nb_var
               IF(name_var_flim == name_var(isubs)) then
                  irkm_flim2_re(ir,nfl2)=irk_mod(isubs)
                  l_var_sim=.true.
               ELSE IF(name_var_flim == 'sal'.or.name_var_flim == 'salinite')then
                  irkm_flim2_re(ir,nfl2)=-1
                  l_var_sim=.true.
               ENDIF
             ENDDO
           CASE (3) 
   ! Limiting functions by solar radiation 
             nb_flim3_re(ir)=nb_flim3_re(ir)+1
             nfl3=nb_flim3_re(ir)
             l_flim3_re=.true.
       ! reALLOCATE arrays if nflimax too small
             IF(nb_flim3_re(ir).gt.nflimax) then
               ALLOCATE(cf_flim3_tmp(nreacmax,nflimax))
               cf_flim3_tmp(:,:)=cf_flim3_re(:,:) 
               deALLOCATE(cf_flim3_re)
               nflimax_new=nb_flim3_re(ir)
               ALLOCATE(cf_flim3_re(nreacmax,nflimax_new))
               cf_flim3_re(:,1:nflimax)=cf_flim3_tmp(:,:)
               nflimax=nflimax_new                      
               deALLOCATE(cf_flim3_tmp)
             ENDIF
             cf_flim3_re(ir,nfl3)=cf_flim_re
           END SELECT
          ENDDO
          IF(.not.l_var_sim.and.ityp_flim.ne.3) then
!l_var_sim=.false.
             WRITE(*,*)'********ERROR ***********'
             WRITE(*,*)' The limiting variable ',name_var_flim
             WRITE(*,*)' of the reaction ',ir,' (',trim(name_reac(ir)),')'
             WRITE(*,*)'does not existed (probably not the same name) !!'
             WRITE(*,*)' Modify her name in the file of reaction parameters'
             WRITE(*,*)'  by comparison with ',trim(name_filesubs)
             STOP
          ENDIF
        ENDIF
      READ(49,'(a)',iostat=eof)
    ELSE
      WRITE(*,*)'End of the file : ',filereaction
      WRITE(*,*)'in the reaction number : ',ir 
      STOP
    ENDIF
         
! save into simu.log
!-------------------
   IF_MPI (MASTER) THEN
   WRITE(iscreenlog,*)'REACTION : ',name_reac(ir),' TYPE :',ityp_reac(ir)
   WRITE(iscreenlog,*)'maxi cinetic ',cin_max_re(ir)
   IF(l_treat_re(ir)) then
      WRITE(iscreenlog,*)'integration of sink & sources due to this reaction '
      WRITE(iscreenlog,*)'in transport equations '
   ELSE
      WRITE(iscreenlog,*)'sink & sources due to this reaction computed '
      WRITE(iscreenlog,*)' after (operator splitting)'
   ENDIF
   IF(.not.l_cincte_re(ir)) then
       WRITE(iscreenlog,*)'cinetic varies and has to be programmed in the subroutine conta_cinreac'
   ELSE
       WRITE(iscreenlog,*)'constant cinetic'
   ENDIF
   IF(.not.l_coefcte_re(ir)) then
       WRITE(iscreenlog,*)'stochiometrics coefficients varyng'
       WRITE(iscreenlog,*)'  (has to be programmed in the subroutine conta_cinreac)'
   ELSE
       WRITE(iscreenlog,*)'constants stochiometrics coefficients '
   ENDIF
   IF(l_arh(ir)) then
       WRITE(iscreenlog,*)'type of arhenius law',ityp_arh(ir)
       WRITE(iscreenlog,*)'teta ', cf_teta_arh(ir)
       WRITE(iscreenlog,*)'reference temperature  ', tempref_arh(ir)
       WRITE(iscreenlog,*)'xeta or sigma ', cf_xetsig_arh(ir)
   ENDIF          
   DO it=1,nb_var_re(ir)
       WRITE(iscreenlog,*)'for the reactive variable n° ',it, name_var_re(ir,it)
       WRITE(iscreenlog,*)'     stochiometric coefficient  ',cf_stochio_re(ir,it)
       WRITE(iscreenlog,*)'     exponent of this variable ',cf_expovar_re(ir,it)
       WRITE(iscreenlog,*)'     model s order number of this variable : ',irkm_var_re(ir,it)
       WRITE(iscreenlog,*)'     data file s order number of this variable : ', irkf_var_re(ir,it)
   ENDDO
   WRITE(iscreenlog,*)' number of limiting functions :',nb_flim1_re(ir)
   WRITE(iscreenlog,*)' number of inhibit functions  :',nb_flim2_re(ir)
   WRITE(iscreenlog,*)' number of radiation limiting functions  :',nb_flim3_re(ir)
   IF(nb_flim1_re(ir).gt.0)then
       DO iflim=1,nb_flim1_re(ir)
        WRITE(iscreenlog,*)'order number of limiting variable =',irkm_flim1_re(ir,iflim)
        WRITE(iscreenlog,*)' Michaelis constant=',cf_flim1_re(ir,iflim)
       ENDDO
   ENDIF
   IF(nb_flim2_re(ir).gt.0)then
       DO iflim=1,nb_flim2_re(ir)
        WRITE(iscreenlog,*)'order number of inhibit variable =',irkm_flim2_re(ir,iflim)
        WRITE(iscreenlog,*)'Michaelis constant=',cf_flim2_re(ir,iflim)
       ENDDO
   ENDIF
   IF(nb_flim3_re(ir).gt.0)then
       DO iflim=1,nb_flim3_re(ir)
        WRITE(iscreenlog,*)' constant a [exp(-a.extinct.z)]=',cf_flim3_re(ir,iflim)
       ENDDO
   ENDIF
   ENDIF_MPI

  END SUBROUTINE conta_read_param_GEN
  !!======================================================================

  SUBROUTINE conta_read_param_REV(ir)
 

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE conta_read_param_REV ***
   !&E
   !&E ** Purpose : Reading parameters for biochemical reactions type REVERSIBLE
   !&E
   !&E ** Description :
   !&E
   !&E ** used ij-arrays : none
   !&E
   !&E ** modified variables : none
   !&E
   !&E ** Called by : conta_read_react.F90
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2008-06  (B. Thouvenin)  Original code
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
    USE comsubstance, ONLY : name_filesubs,irk_mod,name_var,nv_tot

   !! * Arguments
   INTEGER, intent(in)                           :: ir

   !! * Local declarations
   INTEGER                                     :: ir3,it,isubs,nblu,ind_white
   CHARACTER(LEN=5)                            :: comment
   CHARACTER(LEN=20)                           :: commentlu
   LOGICAL                                     :: l_var_sim
     
   !!----------------------------------------------------------------------
   !! * Executable part
   
!   =======================================================
!   Reaction of type 3
!   reversible : dC(iv1)/dt=-  cte_vitess  * (C(iv1) -Ceq(iv1))
!                dC(iv2)/dt=+  cte_vitess  * (C(iv1) -Ceq(iv1))
!   =======================================================


    IF(eof==0) then
     ir3=nb_reac(3)
     nb_var_re(ir)=2
     READ(49,*,iostat=eof)cf_rate_r(ir3)
     cf_rate_r(ir3)=cf_rate_r(ir3)/86400.
! ATTENTION: unit of cv_watequ_re3(ir3)  : mass/mass of total particules
!            equilibrium concentration for the first active variable
     READ(49,*,iostat=eof)cv_watequ_r(ir3)
! law of arhenius
      READ(49,*,iostat=eof)l_arh(ir)
      READ(49,*,iostat=eof)ityp_arh(ir)
      READ(49,*,iostat=eof)cf_teta_arh(ir)
      READ(49,*,iostat=eof)tempref_arh(ir)
      READ(49,*,iostat=eof)cf_xetsig_arh(ir)
      READ(49,'(a)',iostat=eof)
! reactives variables 
      DO it=1,nb_var_re(ir)
        READ(49,'(a)',iostat=eof) commentlu     ! name of reactives variables  (affected by the reaction)
        ind_white=index(commentlu,' ')
        name_var_re(ir,it)=commentlu(1:ind_white)
        l_var_sim=.false.
        DO isubs=1,nb_var
         IF(name_var_re(ir,it) == name_var(isubs)) then
           irkm_var_re(ir,it)=irk_mod(isubs)
           irkf_var_re(ir,it)=isubs
           nb_reac_pervar(isubs)=nb_reac_pervar(isubs)+1
           l_var_sim=.true.
           EXIT
         ENDIF
        ENDDO
        IF(.not.l_var_sim) then
!          si l_var_sim=.false.
          WRITE(*,*)'********ERROR ***********'
          WRITE(*,*)' The reactive variable of the reaction ',ir,' (',trim(name_reac(ir)),')'
          WRITE(*,*)'is not simulated (probably not the same name) !!'
          WRITE(*,*)'Simulated  variables  are :'
          DO isubs=1,nv_tot
            WRITE(*,*)name_var(isubs)
          ENDDO
          WRITE(*,*)' Modify her name :',name_var_re(ir,it)
          WRITE(*,*)' in the file of reaction parameters by comparison with' ,trim(name_filesubs)
          STOP
        ENDIF
      ENDDO
      READ(49,'(a)',iostat=eof)
    ELSE
      WRITE(*,*)'End of the file : ',filereaction
      WRITE(*,*)'in the reaction number : ',ir 
      STOP
    ENDIF
         
! save into simu.log
!-------------------
   IF_MPI (MASTER) THEN
   WRITE(iscreenlog,*)'REACTION : ',name_reac(ir),' TYPE :',ityp_reac(ir)
   WRITE(iscreenlog,*)'rate constant ',cf_rate_r(ir3)
   WRITE(iscreenlog,*)'equilibrium concentration of the first active variable ',cv_watequ_r(ir3)
   IF(l_treat_re(ir)) then
      WRITE(iscreenlog,*)'integration of sink & sources due to this reaction '
      WRITE(iscreenlog,*)'in transport equations '
   ELSE
      WRITE(iscreenlog,*)'sink & sources due to this reaction computed '
      WRITE(iscreenlog,*)' after (operator splitting)'
   ENDIF
   DO it=1,nb_var_re(ir)
       WRITE(iscreenlog,*)'for the reactive variable n° ',it, name_var_re(ir,it)
       WRITE(iscreenlog,*)'     data file s order number of this variable : ', irkf_var_re(ir,it)
   ENDDO
   ENDIF_MPI

  END SUBROUTINE conta_read_param_REV
   !!======================================================================
      
  SUBROUTINE conta_read_param_VOL(ir)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE conta_read_param_VOL ***
   !&E
   !&E ** Purpose : Reading parameters for biochemical reactions type VOLATIL
   !&E
   !&E ** Description :
   !&E
   !&E ** used ij-arrays : none
   !&E
   !&E ** modified variables : none
   !&E
   !&E ** Called by : conta_read_react.F90
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2008-09  (B. Thouvenin)  Original code
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
    USE comsubstance, ONLY : name_filesubs,irk_mod,name_var,nv_tot

   !! * Arguments
   INTEGER, intent(in)                         :: ir

   !! * Local declarations
   LOGICAL                                     :: ex
   INTEGER                                     :: ire,it,isubs,ind_white,iv
   CHARACTER(LEN=5)                            :: comment
   CHARACTER(LEN=20)                           :: commentlu
      
   !!----------------------------------------------------------------------
   !! * Executable part

!   =======================================================
!   Reaction  type VOLATIL
!   volatilisation of one variable from water to air,
!   dependant of its atmospheric concentration and some parameters
!   choice between several formulations
!   =======================================================
    nb_var_re(ir)=1
    IF(eof==0) then
      READ(49,'(a)',iostat=eof) commentlu     ! nom de la variable active
      ind_white=index(commentlu,' ')
      name_var_re(ir,1)=commentlu(1:ind_white)
      DO isubs=1,nb_var
       IF(name_var_re(ir,1) == name_var(isubs)) then
         irkm_var_re(ir,1)=irk_mod(isubs)
         irkf_var_re(ir,1)=isubs
         nb_reac_pervar(isubs)=nb_reac_pervar(isubs)+1
       ENDIF
      ENDDO
      IF(irkm_var_re(ir,1).eq.0)then
        WRITE(*,*)'********ERROR ***********'
        WRITE(*,*)' The active variable of the reaction ',ir,' (',trim(name_reac(ir)),')'
        WRITE(*,*)'is not simulated (probably not the same name) !!'
        WRITE(*,*)'Simulated  variables  are :'
        DO isubs=1,nv_tot
          WRITE(*,*)name_var(isubs)
        ENDDO
        WRITE(*,*)' Modify her name :',name_var_re(ir,1)
        WRITE(*,*)' in the file of reaction parameters'
        STOP
      ENDIF
      iv=irkm_var_re(ir,1)
      READ(49,*,iostat=eof)massmol(iv)
      READ(49,*,iostat=eof)nschmidt(iv)
      READ(49,*,iostat=eof)cf_difwat(iv)
      READ(49,*,iostat=eof)cf_difair(iv)
      READ(49,*,iostat=eof)cf_henry(iv)
      READ(49,*,iostat=eof)cf_henry_adim(iv)
      READ(49,*,iostat=eof)cf_temp_henry(iv)
      READ(49,*,iostat=eof)cf_salt_henry(iv)
! l_volkl_dif : choix de la formule pour coef de transfert phase liquide
!              par diffusion turbulente 
!              true pour Churchill ou false pour O Connor 
      READ(49,*,iostat=eof)l_volkl_dif(iv)
! l_volkl_sal : choix de la prise en compte de l effet de sel (true si oui)
!              par formule de Wen  pour coef de transfert phase liquide
      READ(49,*,iostat=eof)l_volkl_sal(iv)
! ich_volkl_wind : choix de la formule pour coef de transfert phase liquide
!              par turbulence engendree par le vent 
!              1 pour Schwartzenbach 
!              2 pour Wanninkhof 
!              3 pour Nightingale 
      READ(49,*,iostat=eof)ich_volkl_wind(iv)
! l_volkl_cor : choix de la propriete pour correction du coefficient  
!               de transfert phase liquide
!               par nombre de Schmidt (true) ou pas massse molaire (false) 
      READ(49,*,iostat=eof)l_volkl_cor(iv)
! cf_expvol_liq : exposant du rapport des proprietes pour phase liquide 
      READ(49,*,iostat=eof)cf_expvol_liq(iv)
! l_volkg_cor : choix de la propriete pour correction du coefficient 
!               de transfert phase gazeuse
!               par diffusivite (true) ou pas massse molaire (false) 
      READ(49,*,iostat=eof)l_volkg_cor(iv)
! cf_expvol_gas : exposant du rapport des proprietes pour phase gazeuse
      READ(49,*,iostat=eof)cf_expvol_gas(iv)
      READ(49,'(a)',iostat=eof)

    ELSE
      WRITE(*,*)'End of the file : ',filereaction
      WRITE(*,*)'in the reaction number : ',ir 
      STOP
    ENDIF
       
   
! save into simu.log
!-------------------
   IF_MPI (MASTER) THEN
   WRITE(iscreenlog,*)'REACTION : ',name_reac(ir),' TYPE :',ityp_reac(ir)
      WRITE(iscreenlog,*)'integration of sink & sources due to this reaction '
      WRITE(iscreenlog,*)'in transport equations '
   WRITE(iscreenlog,*)'active variable ',name_var_re(ir,1)
   WRITE(iscreenlog,*)'     model s order number of this variable : ',irkm_var_re(ir,1)
   WRITE(iscreenlog,*)'     data file s order number of this variable  : ', irkf_var_re(ir,1)  
   WRITE(iscreenlog,*)'Molecular mass             =',massmol(iv)
   WRITE(iscreenlog,*)'Schmidt number in water    =',nschmidt(iv)
   WRITE(iscreenlog,*)'molecular diffusion in water    =',cf_difwat(iv)
   WRITE(iscreenlog,*)'Diffusivity in air       =',cf_difair(iv)
   IF(cf_henry(iv).ne.0)then
     WRITE(iscreenlog,*)'used  dimensional Henry Constant  =',cf_henry(iv)
   ELSE
     WRITE(iscreenlog,*)'used adimensional Henry Constant  =',cf_henry_adim(iv)
   ENDIF
   WRITE(iscreenlog,*)'B correc Henry with Staudinger f(T)  =',cf_temp_henry(iv)
   WRITE(iscreenlog,*)'correc Henry salting out (L/mol) f(Sal) =',cf_salt_henry(iv)
     IF(l_volkl_dif(iv))then
       WRITE(iscreenlog,*)'Churchill formulation for KL by turbulent diffusion'
      ELSE
       WRITE(iscreenlog,*)'O Connor formulation for KL by turbulent diffusion'
      ENDIF
     IF(l_volkl_sal(iv))then
       WRITE(iscreenlog,*)'effect of salinity (Wen formul) for KL by turbulent diffusion'
      ELSE
       WRITE(iscreenlog,*)'no effect of salinity for KL by turbulent diffusion'
      ENDIF
     IF(ich_volkl_wind(iv) == 1)then
       WRITE(iscreenlog,*)'Schwartzenbach formulation for KL induced by wind '
      ELSE IF(ich_volkl_wind(iv) == 2)then
       WRITE(iscreenlog,*)'Wanninkhof formulation for KL induced by wind '
      ELSE 
       WRITE(iscreenlog,*)'Nightingale formulation for KL induced by wind '
      ENDIF
     IF(l_volkl_cor(iv))then
       WRITE(iscreenlog,*)'for Kl : relation between molar mass '
      ELSE
       WRITE(iscreenlog,*)'for Kl : relation between Schmidt numbers '
      ENDIF
      WRITE(iscreenlog,*)'exponent nL= ',cf_expvol_liq(iv)
     IF(l_volkg_cor(iv))then
       WRITE(iscreenlog,*)'for Kg : relation between diffusivities '
      ELSE
       WRITE(iscreenlog,*)'for Kl : relation between Schmidt numbers '
      ENDIF
      WRITE(iscreenlog,*)'exponent ng= ',cf_expvol_gas(iv)
      WRITE(iscreenlog,*)
   ENDIF_MPI

  END SUBROUTINE conta_read_param_VOL
  !!======================================================================

  SUBROUTINE conta_read_param_PART(ir)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE conta_read_param_PART ***
   !&E
   !&E ** Purpose : Reading parameters for biochemical reactions type PART_EQ
   !&E
   !&E ** Description :
   !&E
   !&E ** used ij-arrays : none
   !&E
   !&E ** modified variables : none
   !&E
   !&E ** Called by : conta_read_react.F90
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2008-06  (B. Thouvenin)  Original code
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
    USE comsubstance, ONLY : name_filesubs,name_var,irk_mod
                             
   !! * Arguments
   INTEGER, intent(in)                         :: ir

   !! * Local declarations
   LOGICAL                                     :: l_var_sim,l_unitKd
   INTEGER                                     :: ire,it,isubs,nfl1,ir5,i,ind_white
   INTEGER                                     :: iflim,nfunct_limit,il,nblu,  &
                                                  ityp_flim,nfl2,nfl3,nvaractmax_new
   CHARACTER(LEN=5)                            :: comment
   CHARACTER(LEN=20)                           :: commentlu,namcop,       &
                                                    name_var_flim
   CHARACTER(LEN=lchain)                       :: name_sub
   CHARACTER(LEN=10)                           :: name_typ_reacr
   REAL(KIND=rsh)                              :: cf_flim_re
   INTEGER,dimension(:,:),allocatable          :: irkm_var_tmp,irkf_var_tmp
   INTEGER,dimension(:,:),allocatable          :: irkm_sub_tmp
   CHARACTER(len=lchain),dimension(:,:),allocatable :: name_var_tmp
   REAL(KIND=rsh),dimension(:,:),allocatable        :: cf_distb_tmp
     
   !!----------------------------------------------------------------------
   !! * Executable part
   
!   =======================================================
!   Reaction de type 5
!   partage a l equilibre
!   =======================================================

    IF(eof==0) then
     READ(49,*,iostat=eof)nblu
     nb_var_re(ir)=nblu
     ir5=nb_reac(5)
     
! reALLOCATE arrays if nvaractmax too small
             IF(nb_var_re(ir).gt.nvaractmax) then
               ALLOCATE(name_var_tmp(nreacmax,nvaractmax))
               ALLOCATE(irkm_var_tmp(nreacmax,nvaractmax))
               ALLOCATE(irkf_var_tmp(nreacmax,nvaractmax))
               ALLOCATE(irkm_sub_tmp(nreacmax,nvaractmax))
               ALLOCATE(cf_distb_tmp(nreacmax,nvaractmax))
               name_var_tmp(:,:)=name_var_re(:,:)        
               irkm_var_tmp(:,:)=irkm_var_re(:,:) 
               irkf_var_tmp(:,:)=irkf_var_re(:,:) 
               irkm_sub_tmp(:,:)=irkm_sub_r(:,:) 
               cf_distb_tmp(:,:)=cf_distb_r(:,:) 
               deALLOCATE(name_var_re)
               deALLOCATE(irkm_var_re)
               deALLOCATE(irkf_var_re)
               deALLOCATE(irkm_sub_r)
               deALLOCATE(cf_distb_r)
               nvaractmax_new=nb_var_re(ir)
               ALLOCATE(name_var_re(nreacmax,nvaractmax_new))
               ALLOCATE(irkm_var_re(nreacmax,nvaractmax_new))
               ALLOCATE(irkf_var_re(nreacmax,nvaractmax_new))
               ALLOCATE(irkm_sub_r(nreacmax,nvaractmax_new))
               ALLOCATE(cf_distb_r(nreacmax,nvaractmax_new))
               name_var_re(:,1:nvaractmax)=name_var_tmp(:,:)        
               irkm_var_re(:,1:nvaractmax)=irkm_var_tmp(:,:)        
               irkf_var_re(:,1:nvaractmax)=irkf_var_tmp(:,:)        
               irkm_sub_r(:,1:nvaractmax)=irkm_sub_tmp(:,:)        
               cf_distb_r(:,1:nvaractmax)=cf_distb_tmp(:,:)        
               nvaractmax=nvaractmax_new                      
               deALLOCATE(name_var_tmp)
               deALLOCATE(irkm_var_tmp)
               deALLOCATE(irkf_var_tmp)
               deALLOCATE(irkm_sub_tmp)
               deALLOCATE(cf_distb_tmp)
             ENDIF
     
! arhenius Law
     READ(49,*,iostat=eof)l_arh(ir)
     READ(49,*,iostat=eof)ityp_arh(ir)
     READ(49,*,iostat=eof)cf_teta_arh(ir)
     READ(49,*,iostat=eof)tempref_arh(ir)
     READ(49,*,iostat=eof)cf_xetsig_arh(ir)
! Particule Concentration Effect PCE
     READ(49,*,iostat=eof)l_pce_r(ir5)
     READ(49,*,iostat=eof)cf_pce_r(ir5,:)
     IF(.not.l_pce_r(ir5))then
          cf_pce_r(ir5,1)=0.
          cf_pce_r(ir5,2)=1.
     ENDIF
        
! effect salting out         
     READ(49,*,iostat=eof)l_salt_r(ir5)
     READ(49,*,iostat=eof)cf_salt_r(ir5,:)
        IF(.not.l_salt_r(ir5))then
          cf_salt_r(ir5,1)=0.
          cf_salt_r(ir5,2)=0.
        ENDIF

! KD variable with salinity (complexation with chlorure (MOCO) or Turner s law..
      READ(49,*,iostat=eof)l_kds_r(ir5)
      IF(l_kds_r(ir5))then
        READ(49,*,iostat=eof)ndg_kds_r(ir5),cf_kd0_r(ir5)
        cf_xkds_r(ir5,:)=0.
        IF(ndg_kds_r(ir5).gt.0)then
          READ(49,*,iostat=eof)(cf_xkds_r(ir5,il),il=1,ndg_kds_r(ir5))
        ELSE
          READ(49,*,iostat=eof)
        ENDIF
      ELSE
        ndg_kds_r(ir5)=1
        READ(49,*,iostat=eof)
        READ(49,*,iostat=eof)
      ENDIF 

! FREE Dissolved variable  
      it=1
      READ(49,'(a)',iostat=eof) commentlu     ! nom de la variables dissoute (affectee par la reaction)
      ind_white=index(commentlu,' ')
      name_var_re(ir,it)=commentlu(1:ind_white)
      l_var_sim=.false.
      DO isubs=1,nb_var
         IF(name_var_re(ir,it) == name_var(isubs)) then
           irkm_var_re(ir,it)=irk_mod(isubs)
           irkf_var_re(ir,it)=isubs
           nb_reac_pervar(isubs)=nb_reac_pervar(isubs)+1
           l_var_sim=.true.
           EXIT
         ENDIF
      ENDDO
      IF(.not.l_var_sim) then
!          si l_var_sim=.false.
          WRITE(*,*)'********ERROR ***********'
          WRITE(*,*)' The reactive variable of the reaction ',ir,' (',trim(name_reac(ir)),')'
          WRITE(*,*)'does not exist (probably not the same name) !!'
          WRITE(*,*)' Modify her name :',name_var_re(ir,it)
          WRITE(*,*)' in the file of reaction parameters by comparison with ',trim(name_filesubs)
           STOP
       ENDIF
       READ(49,'(a)',iostat=eof) 
       
! OTHERS variables in equilibrium , bound with a substrate
      DO it=2,nb_var_re(ir)
        READ(49,'(a)',iostat=eof) commentlu     ! nom des variables reactives (affectees par la reaction)
        ind_white=index(commentlu,' ')
        name_var_re(ir,it)=commentlu(1:ind_white)
        l_var_sim=.false.
        DO isubs=1,nb_var
         IF(name_var_re(ir,it) == name_var(isubs)) then
           irkm_var_re(ir,it)=irk_mod(isubs)
           irkf_var_re(ir,it)=isubs
           nb_reac_pervar(isubs)=nb_reac_pervar(isubs)+1
           l_var_sim=.true.
           EXIT
         ENDIF
        ENDDO
        IF(.not.l_var_sim) then
!          si l_var_sim=.false.
          WRITE(*,*)'********ERROR ***********'
          WRITE(*,*)' The reactive variable of the reaction ',ir,' (',trim(name_reac(ir)),')'
          WRITE(*,*)'does not exist (probably not the same name) !!'
          WRITE(*,*)' Modify her name :',name_var_re(ir,it)
          WRITE(*,*)' in the file of reaction parameters by comparison with ',trim(name_filesubs)
           STOP
        ENDIF
        READ(49,*,iostat=eof)commentlu
        ind_white=index(commentlu,' ')
        name_sub=commentlu(1:ind_white)
! if substrat=SMtot the ratio SMtot/SM will be equal to 1, no normalization
        l_var_sim=.false.
        DO isubs=1,nb_var
          IF(name_sub == name_var(isubs)) then
             irkm_sub_r(ir5,it)=irk_mod(isubs)
             l_var_sim=.true.
          ELSE IF(name_sub == 'SMtot') then
             l_var_sim=.true. 
             irkm_sub_r(ir5,it)=-1
! mis à -1 pour reperer que c est SMtot
          ENDIF
        ENDDO
        IF(.not.l_var_sim) then
!          si l_var_sim=.false.
            WRITE(*,*)'********ERROR ***********'
            WRITE(*,*)' The substrate ',name_sub
            WRITE(*,*)'of the variable ',trim(name_var_re(ir,it))
            WRITE(*,*)'in the reaction ',ir,' (',trim(name_reac(ir)),')'
            WRITE(*,*)'does not exist (probably not the same name) !!'
            WRITE(*,*)' Modify her name in the file of reaction parameters'
            WRITE(*,*)'  by comparison with ',trim(name_filesubs)
            STOP
        ENDIF
        READ(49,*,iostat=eof)l_unitKd
        READ(49,*,iostat=eof)cf_distb_r(ir5,it)
        IF(.not.l_unitKd)then
!         coeffient de partage exprime  en log10(ml/g) converti en  L/G 
           cf_distb_r(ir5,it)=10**(cf_distb_r(ir5,it))/1000
        ENDIF
        
     ENDDO
      READ(49,'(a)',iostat=eof)
    ELSE
      WRITE(*,*)'End of the file : ',filereaction
      WRITE(*,*)'in the reaction number : ',ir 
      WRITE(*,*)ir5,' reaction type 5' 
      STOP
    ENDIF
       
! save into simu.log
!-------------------
   IF_MPI (MASTER) THEN
   WRITE(iscreenlog,*)'REACTION : ',name_reac(ir),' TYPE :',ityp_reac(ir)
   l_treat_re(ir)=.false.
   WRITE(iscreenlog,*)'sink & sources due to this reaction computed '
   WRITE(iscreenlog,*)' after (operator splitting)'
   IF(l_arh(ir)) then
       WRITE(iscreenlog,*)'type of arhenius law',ityp_arh(ir)
       WRITE(iscreenlog,*)'teta ', cf_teta_arh(ir)
       WRITE(iscreenlog,*)'reference temperature  ', tempref_arh(ir)
       WRITE(iscreenlog,*)'xeta or sigma ', cf_xetsig_arh(ir)
   ENDIF          
   WRITE(iscreenlog,*)'the distribution treated by this reaction affect :',nb_var_re(ir),' variables'
   WRITE(iscreenlog,*)'the free dissolved reactive variable  is :', name_var_re(ir,1)
   WRITE(iscreenlog,*)'her model s order number is : ',irkm_var_re(ir,1)
   WRITE(iscreenlog,*)'her file s order number is : ',irkf_var_re(ir,1)
   DO it=2,nb_var_re(ir)
       WRITE(iscreenlog,*)'the reactive variable :', name_var_re(ir,it)
       WRITE(iscreenlog,*)'with model s order number  : ',irkm_var_re(ir,it)
       WRITE(iscreenlog,*)'and file s order number  :  ',irkf_var_re(ir,it)
       WRITE(iscreenlog,*)'associated substrate (model s order number : ',irkm_sub_r(ir5,it)
       WRITE(iscreenlog,*)'normalised partition coefficient  : ',cf_distb_r(ir5,it)
   ENDDO
   WRITE(iscreenlog,*)
   ENDIF_MPI


  END SUBROUTINE conta_read_param_PART
  !!======================================================================
      
  SUBROUTINE conta_add_reac

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE conta_add_reac  ***
   !&E
   !&E ** Purpose : addition of biochemical reactions definied in the file
   !&E              of chemical CHARACTERistics and contaminant speciation  
   !&E
   !&E ** Description :
   !&E
   !&E ** used ij-arrays : none
   !&E
   !&E ** modified variables : none
   !&E
   !&E ** Called by : conta_read_react.F90
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2008-06  (B. Thouvenin)  Original code
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
    USE comsubstance, ONLY : name_filesubs
   !! * Arguments

   !! * Local declarations
   INTEGER                                     :: iv,ispc,ire,ico,nb_reac_add,it,isubs,iflim,ivdl
   INTEGER                                     ::   ivs,ir,ireac,ifracv,nvaractmax_new,nflimax_new
   LOGICAL                                     :: l_var_sim
   CHARACTER(LEN=1)                            :: numcon
   INTEGER,dimension(:,:),allocatable          :: irkm_var_tmp,irkf_var_tmp
   INTEGER,dimension(:,:),allocatable          :: irkm_sub_tmp,ispc_reac_tmp
   INTEGER,dimension(:,:),allocatable          :: ifrac_var_tmp
   INTEGER,dimension(:,:),allocatable          :: irkm_flim1_tmp,irkm_flim2_tmp
   CHARACTER(len=lchain),dimension(:,:),allocatable :: name_var_tmp
   REAL(KIND=rsh),dimension(:,:),allocatable   :: cf_distb_tmp,cf_stochio_tmp
   REAL(KIND=rsh),dimension(:,:),allocatable   :: cf_expovar_tmp,cf_flim1_tmp
   REAL(KIND=rsh),dimension(:,:),allocatable   :: cf_flim2_tmp,cf_flim3_tmp
     
   !!----------------------------------------------------------------------
   !! * Executable part
   ire=nreacmax_sco
   nb_reac_add=nreacmax_sco
   
   ! save into simu.log
   !-------------------
   IF_MPI (MASTER) THEN
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '**************************************************'
   WRITE(iscreenlog,*) '**************** MET&OR_add_react.F90 ***************'
   WRITE(iscreenlog,*) '**************************************************'
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) 'ADDITION OF REACTIONS DEFINED FOR CONTAMINANTS'
   ENDIF_MPI

! Boucle si il y a plusieurs contaminants simules
  DO ico=1,nb_conta
   IF_MPI (MASTER) THEN
     WRITE(iscreenlog,*) 'IN FILE : ',trim(fileconta(ico))
   ENDIF_MPI
    
! addition of reactions type 1 (hydrolyse, biodegradation, photolyse simple)
! ===========================================================================
   DO ispc=1,nb_spc_tot(ico)
    IF(cin_hydrol(ico,ispc) /= 0.)then
     ire=ire+1
     name_reac(ire)='HYDROLYSE_'//name_spc(ico,ispc)
     ityp_reac(ire)=1
     nb_reac(1)=nb_reac(1)+1
     ire_typ(ire)=nb_reac(1)
     cin_max_re(ire)=cin_hydrol(ico,ispc)
     l_arh(ire)=.false.
     l_treat_re(ire)=.true.
     nb_var_re(ire)=1
     cf_stochio_re(ire,1)=-1.
     ispc_reac(ire,1)=ispc
     icon_reac(ire)=ico
#ifdef key_species
     DO iv=1,nv_adv
       IF (icon_reac(ire) == icon_var(iv)) then
        DO ifracv=1,nb_frac_var(iv) 
         IF(ispc_reac(ire,1) == ispc_var(iv,ifracv))then
            ifrac_var_reac(ire,1)=ifracv
            nb_reac_pervar(irk_fil(iv))=nb_reac_pervar(irk_fil(iv))+1
            name_var_re(ire,1)=name_var(irk_fil(iv))
            irkm_var_re(ire,1)=iv
         ENDIF
        ENDDO
       ENDIF
     ENDDO
#else
     DO iv=1,nv_adv
       IF (icon_reac(ire) == icon_var(iv) .and. ispc_reac(ire,1) == ispc_var(iv))then
            nb_reac_pervar(irk_fil(iv))=nb_reac_pervar(irk_fil(iv))+1
            name_var_re(ire,1)=name_var(irk_fil(iv))
            irkm_var_re(ire,1)=iv
       ENDIF
     ENDDO
#endif
    ENDIF
    IF(cin_biodeg(ico,ispc) /= 0.)then
     ire=ire+1
     name_reac(ire)='BIODEGRADATION_'//name_spc(ico,ispc)
     ityp_reac(ire)=1
     nb_reac(1)=nb_reac(1)+1
     ire_typ(ire)=nb_reac(1)
     cin_max_re(ire)=cin_biodeg(ico,ispc)
     l_arh(ire)=.false.
     l_treat_re(ire)=.true.
     nb_var_re(ire)=1
     cf_stochio_re(ire,1)=-1.
     ispc_reac(ire,1)=ispc
     icon_reac(ire)=ico
#ifdef key_species
     DO iv=1,nv_adv
       IF (icon_reac(ire) == icon_var(iv)) then
        DO ifracv=1,nb_frac_var(iv) 
         IF(ispc_reac(ire,1) == ispc_var(iv,ifracv))then
            ifrac_var_reac(ire,1)=ifracv
            nb_reac_pervar(irk_fil(iv))=nb_reac_pervar(irk_fil(iv))+1
            name_var_re(ire,1)=name_var(irk_fil(iv))
            irkm_var_re(ire,1)=iv
         ENDIF
        ENDDO
       ENDIF
     ENDDO
#else
     DO iv=1,nv_adv
       IF (icon_reac(ire) == icon_var(iv) .and. ispc_reac(ire,1) == ispc_var(iv))then
            nb_reac_pervar(irk_fil(iv))=nb_reac_pervar(irk_fil(iv))+1
            name_var_re(ire,1)=name_var(irk_fil(iv))
            irkm_var_re(ire,1)=iv
       ENDIF
     ENDDO
#endif
    ENDIF
    IF(cin_photo(ico,ispc) /= 0.)then
     ire=ire+1
     name_reac(ire)='PHOTOLYSE_'//name_spc(ico,ispc)
     ityp_reac(ire)=1
     nb_reac(1)=nb_reac(1)+1
     ire_typ(ire)=nb_reac(1)
     cin_max_re(ire)=cin_photo(ico,ispc)
     l_arh(ire)=.false.
     l_treat_re(ire)=.true.
     nb_var_re(ire)=1
     cf_stochio_re(ire,1)=-1.
     ispc_reac(ire,1)=ispc
     icon_reac(ire)=ico
#ifdef key_species
     DO iv=1,nv_adv
       IF (icon_reac(ire) == icon_var(iv)) then
        DO ifracv=1,nb_frac_var(iv) 
         IF(ispc_reac(ire,1) == ispc_var(iv,ifracv))then
             ifrac_var_reac(ire,1)=ifracv
             nb_reac_pervar(irk_fil(iv))=nb_reac_pervar(irk_fil(iv))+1
             name_var_re(ire,1)=name_var(irk_fil(iv))
             irkm_var_re(ire,1)=iv
         ENDIF
        ENDDO
       ENDIF
     ENDDO
#else
     DO iv=1,nv_adv
       IF (icon_reac(ire) == icon_var(iv) .and. ispc_reac(ire,1) == ispc_var(iv))then
            nb_reac_pervar(irk_fil(iv))=nb_reac_pervar(irk_fil(iv))+1
            name_var_re(ire,1)=name_var(irk_fil(iv))
            irkm_var_re(ire,1)=iv
       ENDIF
     ENDDO
#endif
    ENDIF
   ENDDO

  nb_reac_add=ire-nb_reac_add
  
! save into simu.log
!-------------------
  DO ir=nb_reac(1)-nb_reac_add+1,nb_reac(1)
   IF_MPI (MASTER) THEN
   WRITE(iscreenlog,*)
   WRITE(iscreenlog,*)'ADDITIONAL REACTION FOR CONTAMINANT '
   WRITE(iscreenlog,*)'REACTION : ',trim(name_reac(ir)),' TYPE :',ityp_reac(ir)
   WRITE(iscreenlog,*)'maxi cinetic ',cin_max_re(ir)
   IF(l_treat_re(ir)) then
      WRITE(iscreenlog,*)'integration of sink & sources due to this reaction '
      WRITE(iscreenlog,*)'in transport equations '
   ELSE
      WRITE(iscreenlog,*)'sink & sources due to this reaction computed '
      WRITE(iscreenlog,*)' after (operator splitting)'
   ENDIF
   WRITE(iscreenlog,*)'maxi cinetic ',cin_max_re(ir)
   WRITE(iscreenlog,*)'active variable ',trim(name_var_re(ir,1))
   WRITE(iscreenlog,*)'only one active specie  number  ',ispc_reac(ire,1)
   WRITE(iscreenlog,*)'     stochiometric coefficient  ',cf_stochio_re(ir,1)
   WRITE(iscreenlog,*)'     model s order number of this variable : ',irkm_var_re(ir,1)
   WRITE(iscreenlog,*)'     data file s order number of this variable  : ', irkf_var_re(ir,1)  
   ENDIF_MPI
  ENDDO
     
! addition of reactions type 2 (GENERIC)
! ===========================================================================
  DO ireac=1,nb_reac_genco(ico)
    ire=ire+1
    ityp_reac(ire)=2
    nb_reac(2)=nb_reac(2)+1
    ire_typ(ire)=nb_reac(2)
    name_reac(ire)=name_reac_co(ico,ireac)
    l_treat_re(ire)=l_treat_co(ico,ireac)
    nb_var_re(ire)=nb_var_co(ico,ireac)
    cin_max_re(ire)= cin_max_co(ico,ireac)
    l_cincte_re(ire)=l_cincte_co(ico,ireac)
    l_coefcte_re(ire)=l_coefcte_co(ico,ireac)
    l_arh(ire)=l_arh_rco(ico,ireac)
    ityp_arh(ire)=ityp_arh_rco(ico,ireac)
    cf_teta_arh(ire)=cf_teta_arh_rco(ico,ireac)
    tempref_arh(ire)=tempref_arh_rco(ico,ireac)
    cf_xetsig_arh(ire)=cf_xetsig_arh_rco(ico,ireac)
! reALLOCATE arrays if nvaractmax too small
    IF(nb_var_re(ire).gt.nvaractmax) then
               ALLOCATE(name_var_tmp(nreacmax,nvaractmax))
               ALLOCATE(irkm_var_tmp(nreacmax,nvaractmax))
               ALLOCATE(irkf_var_tmp(nreacmax,nvaractmax))
               ALLOCATE(cf_stochio_tmp(nreacmax,nvaractmax))
               ALLOCATE(cf_expovar_tmp(nreacmax,nvaractmax))
               ALLOCATE(ispc_reac_tmp(nreacmax,nvaractmax))
               name_var_tmp(:,:)=name_var_re(:,:)        
               irkm_var_tmp(:,:)=irkm_var_re(:,:) 
               irkf_var_tmp(:,:)=irkf_var_re(:,:) 
               cf_stochio_tmp(:,:)=cf_stochio_re(:,:) 
               cf_expovar_tmp(:,:)=cf_expovar_re(:,:) 
               ispc_reac_tmp(:,:)=ispc_reac(:,:) 
               deALLOCATE(name_var_re)
               deALLOCATE(irkm_var_re)
               deALLOCATE(irkf_var_re)
               deALLOCATE(cf_stochio_re)
               deALLOCATE(cf_expovar_re)
               deALLOCATE(ispc_reac)
               nvaractmax_new=nb_var_re(ire)
               ALLOCATE(name_var_re(nreacmax,nvaractmax_new))
               ALLOCATE(irkm_var_re(nreacmax,nvaractmax_new))
               ALLOCATE(irkf_var_re(nreacmax,nvaractmax_new))
               ALLOCATE(cf_stochio_re(nreacmax,nvaractmax_new))
               ALLOCATE(cf_expovar_re(nreacmax,nvaractmax_new))
               ALLOCATE(ispc_reac(nreacmax,nvaractmax_new))
               name_var_re(:,1:nvaractmax)=name_var_tmp(:,:)        
               irkm_var_re(:,1:nvaractmax)=irkm_var_tmp(:,:)        
               irkf_var_re(:,1:nvaractmax)=irkf_var_tmp(:,:)        
               cf_stochio_re(:,1:nvaractmax)=cf_stochio_tmp(:,:)        
               cf_expovar_re(:,1:nvaractmax)=cf_expovar_tmp(:,:)        
               ispc_reac(:,1:nvaractmax)=ispc_reac_tmp(:,:)        
               deALLOCATE(name_var_tmp)
               deALLOCATE(irkm_var_tmp)
               deALLOCATE(irkf_var_tmp)
               deALLOCATE(cf_stochio_tmp)
               deALLOCATE(cf_expovar_tmp)
               deALLOCATE(ispc_reac_tmp)
#ifdef key_species
               ALLOCATE(ifrac_var_tmp(nreacmax,nvaractmax))
               ifrac_var_tmp(:,:)=ifrac_var_reac(:,:) 
               deALLOCATE(ifrac_var_reac)
               ALLOCATE(ifrac_var_reac(nreacmax,nvaractmax_new))
               ifrac_var_reac(:,1:nvaractmax)=ifrac_var_tmp(:,:)        
               deALLOCATE(ifrac_var_tmp)
#endif
               nvaractmax=nvaractmax_new                      
    ENDIF
    DO it=1,nb_var_re(ire)
      name_var_re(ire,it)=name_var_co(ico,ireac,it)
      l_var_sim=.false.
      DO isubs=1,nb_var
        IF(name_var_re(ire,it) == name_var(isubs)) then
           irkm_var_re(ire,it)=irk_mod(isubs)
           irkf_var_re(ire,it)=isubs
           nb_reac_pervar(isubs)=nb_reac_pervar(isubs)+1
           l_var_sim=.true.
           ispc_reac(ire,it)=0.
           icon_reac(ire)=ico
           DO iv=1,nv_adv
              IF (icon_reac(ire) == icon_var(iv)) then
#ifdef key_species
                DO ifracv=1,nb_frac_var(iv) 
                  IF(ispc_reac(ire,it) == ispc_var(iv,ifracv))then
                     ifrac_var_reac(ire,it)=ifracv
                     nb_reac_pervar(irk_fil(iv))=nb_reac_pervar(irk_fil(iv))+1
                  ENDIF
                ENDDO
#else
                  IF(ispc_reac(ire,it) == ispc_var(iv))then
                     nb_reac_pervar(irk_fil(iv))=nb_reac_pervar(irk_fil(iv))+1
                  ENDIF
#endif
              ENDIF
           ENDDO
           EXIT
        ENDIF
      ENDDO
      DO ispc=1,nb_spc_tot(ico)
        IF(name_var_re(ire,it) == name_spc(ico,ispc)) then
           irkm_var_re(ire,it)=irk_mod(ivar_spc(ico,ispc))
           irkf_var_re(ire,it)=ivar_spc(ico,ispc)
           l_var_sim=.true.
           ispc_reac(ire,it)=ispc
           icon_reac(ire)=ico
           DO iv=1,nv_adv
              IF (icon_reac(ire) == icon_var(iv)) then
#ifdef key_species
                DO ifracv=1,nb_frac_var(iv) 
                  IF(ispc_reac(ire,it) == ispc_var(iv,ifracv))then
                      ifrac_var_reac(ire,it)=ifracv
                     nb_reac_pervar(irk_fil(iv))=nb_reac_pervar(irk_fil(iv))+1
                  ENDIF
                ENDDO
#else
                  IF(ispc_reac(ire,it) == ispc_var(iv))then
                     nb_reac_pervar(irk_fil(iv))=nb_reac_pervar(irk_fil(iv))+1
                  ENDIF
#endif
              ENDIF
           ENDDO
           EXIT
        ENDIF
      ENDDO
      IF(.not.l_var_sim) then
!          si l_var_sim=.false.
          WRITE(*,*)'********ERROR ***********'
          WRITE(*,*)' The reactive variable of the reaction ',ire,' (',trim(name_reac(ire)),')'
          WRITE(*,*)'is not simulated (probably not the same name) !!'
          WRITE(*,*)'Simulated  variables  are :'
          DO iv=1,nv_tot
            WRITE(*,*)trim(name_var(iv))
          ENDDO
          WRITE(*,*)' Modify her name :',trim(name_var_re(ire,it))
          WRITE(*,*)' in the file of contaminant CHARACTERistics  by comparison with' ,trim(name_filesubs)
          STOP
      ENDIF
      cf_stochio_re(ire,it)=cf_stochio_co(ico,ireac,it)
      cf_expovar_re(ire,it)=cf_expovar_co(ico,ireac,it)
    ENDDO
    nb_flim1_re(ire)=nb_flim1_co(ico,ireac)
    nb_flim2_re(ire)=nb_flim2_co(ico,ireac)
    nb_flim3_re(ire)=nb_flim3_co(ico,ireac)
  ! reALLOCATE arrays if nflimax too small
             IF(nb_flim1_re(ire).gt.nflimax) then
               ALLOCATE(irkm_flim1_tmp(nreacmax,nflimax))
               ALLOCATE(cf_flim1_tmp(nreacmax,nflimax))
               irkm_flim1_tmp(:,:)=irkm_flim1_re(:,:)        
               cf_flim1_tmp(:,:)=cf_flim1_re(:,:) 
               deALLOCATE(irkm_flim1_re)       
               deALLOCATE(cf_flim1_re)
               nflimax_new=nb_flim1_re(ire)
               ALLOCATE(irkm_flim1_re(nreacmax,nflimax_new))     
               ALLOCATE(cf_flim1_re(nreacmax,nflimax_new))
               irkm_flim1_re(:,1:nflimax)=irkm_flim1_tmp(:,:)        
               cf_flim1_re(:,1:nflimax)=cf_flim1_tmp(:,:)
               nflimax=nflimax_new                      
               deALLOCATE(irkm_flim1_tmp)
               deALLOCATE(cf_flim1_tmp)
             ENDIF
             IF(nb_flim2_re(ire).gt.nflimax) then
               ALLOCATE(irkm_flim2_tmp(nreacmax,nflimax))
               ALLOCATE(cf_flim2_tmp(nreacmax,nflimax))
               irkm_flim2_tmp(:,:)=irkm_flim2_re(:,:)        
               cf_flim2_tmp(:,:)=cf_flim2_re(:,:) 
               deALLOCATE(irkm_flim2_re)       
               deALLOCATE(cf_flim2_re)
               nflimax_new=nb_flim2_re(ire)
               ALLOCATE(irkm_flim2_re(nreacmax,nflimax_new))     
               ALLOCATE(cf_flim2_re(nreacmax,nflimax_new))
               irkm_flim2_re(:,1:nflimax)=irkm_flim2_tmp(:,:)        
               cf_flim2_re(:,1:nflimax)=cf_flim2_tmp(:,:)
               nflimax=nflimax_new                      
               deALLOCATE(irkm_flim2_tmp)
               deALLOCATE(cf_flim2_tmp)
             ENDIF
             IF(nb_flim3_re(ire).gt.nflimax) then
               ALLOCATE(cf_flim3_tmp(nreacmax,nflimax))
               cf_flim3_tmp(:,:)=cf_flim3_re(:,:) 
               deALLOCATE(cf_flim3_re)
               nflimax_new=nb_flim3_re(ire)
               ALLOCATE(cf_flim3_re(nreacmax,nflimax_new))
               cf_flim3_re(:,1:nflimax)=cf_flim3_tmp(:,:)
               nflimax=nflimax_new                      
               deALLOCATE(cf_flim3_tmp)
             ENDIF
    DO iflim=1,nb_flim1_re(ire)
     cf_flim1_re(ire,iflim)=cf_flim1_co(ico,ireac,iflim)
     DO isubs=1,nb_var
        IF(name_var_flim1(ico,ireac,iflim) == name_var(isubs)) then          
            irkm_flim1_re(ire,iflim)=irk_mod(isubs)
        ELSE IF(name_var_flim1(ico,ireac,iflim) == 'sal'   &
             .or.name_var_flim1(ico,ireac,iflim) == 'salinite')then
            irkm_flim1_re(ire,iflim)=-1
        ENDIF
     ENDDO
    ENDDO
    DO iflim=1,nb_flim2_re(ire)
     cf_flim2_re(ire,iflim)=cf_flim2_co(ico,ireac,iflim)
     DO isubs=1,nb_var
        IF(name_var_flim2(ico,ireac,iflim) == name_var(isubs)) then          
            irkm_flim2_re(ire,iflim)=irk_mod(isubs)
        ELSE IF(name_var_flim2(ico,ireac,iflim) == 'sal'   &
             .or.name_var_flim2(ico,ireac,iflim) == 'salinite')then
            irkm_flim2_re(ire,iflim)=-1
        ENDIF
     ENDDO
    ENDDO
    DO iflim=1,nb_flim3_re(ire)
     cf_flim3_re(ire,iflim)=cf_flim3_co(ico,ireac,iflim)
     l_flim3_re=.true.
    ENDDO
! save into simu.log
!-------------------
      IF_MPI (MASTER) THEN
      WRITE(iscreenlog,*)
      WRITE(iscreenlog,*)'ADDITIONAL REACTION FOR CONTAMINANT '
      WRITE(iscreenlog,*)'REACTION : ',trim(name_reac(ire)),' TYPE :',ityp_reac(ire)
      WRITE(iscreenlog,*)'maxi cinetic ',cin_max_re(ire)
      IF(l_treat_re(ire)) then
        WRITE(iscreenlog,*)'integration of sink & sources due to this reaction '
        WRITE(iscreenlog,*)'in transport equations '
      ELSE
        WRITE(iscreenlog,*)'sink & sources due to this reaction computed '
        WRITE(iscreenlog,*)' after (operator splitting)'
      ENDIF
      IF(.not.l_cincte_re(ire)) then
       WRITE(iscreenlog,*)'cinetic varies and has to be programmed in the subroutine conta_cinreac'
      ELSE
       WRITE(iscreenlog,*)'constant cinetic'
      ENDIF
      IF(.not.l_coefcte_re(ire)) then
       WRITE(iscreenlog,*)'stochiometrics coefficients varyng'
       WRITE(iscreenlog,*)'  (has to be programmed in the subroutine conta_cinreac)'
      ELSE
       WRITE(iscreenlog,*)'constants stochiometrics coefficients '
      ENDIF
      IF(l_arh(ire)) then
       WRITE(iscreenlog,*)'type of arhenius law',ityp_arh(ire)
       WRITE(iscreenlog,*)'teta ', cf_teta_arh(ire)
       WRITE(iscreenlog,*)'reference temperature  ', tempref_arh(ire)
       WRITE(iscreenlog,*)'xeta or sigma ', cf_xetsig_arh(ire)
      ENDIF          
      DO it=1,nb_var_re(ire)
       WRITE(iscreenlog,*)'for the reactive variable n° ',it, trim(name_var_re(ire,it))
       if(ispc_reac(ire,it) > 0)WRITE(iscreenlog,*)' one active specie  number  ',ispc_reac(ire,it)
       WRITE(iscreenlog,*)'     stochiometric coefficient  ',cf_stochio_re(ire,it)
       WRITE(iscreenlog,*)'     exponent of this variable ',cf_expovar_re(ire,it)
       WRITE(iscreenlog,*)'     model s order number of this variable : ',irkm_var_re(ire,it)
       WRITE(iscreenlog,*)'     data file s order number of this variable : ', irkf_var_re(ire,it)
      ENDDO
      WRITE(iscreenlog,*)' number of limiting functions :',nb_flim1_re(ire)
      WRITE(iscreenlog,*)' number of inhibit functions  :',nb_flim2_re(ire)
      WRITE(iscreenlog,*)' number of radiation limiting functions  :',nb_flim3_re(ire)
      IF(nb_flim1_re(ire).gt.0)then
       DO iflim=1,nb_flim1_re(ire)
        WRITE(iscreenlog,*)'order number of limiting variable =',irkm_flim1_re(ire,iflim)
        WRITE(iscreenlog,*)' Michaelis constant=',cf_flim1_re(ire,iflim)
       ENDDO
      ENDIF
      IF(nb_flim2_re(ire).gt.0)then
       DO iflim=1,nb_flim2_re(ire)
        WRITE(iscreenlog,*)'order number of inhibit variable =',irkm_flim2_re(ire,iflim)
        WRITE(iscreenlog,*)'Michaelis constant=',cf_flim2_re(ire,iflim)
       ENDDO
      ENDIF
      IF(nb_flim3_re(ire).gt.0)then
       DO iflim=1,nb_flim3_re(ire)
        WRITE(iscreenlog,*)' constant a [exp(-a.extinct.z)]=',cf_flim3_re(ire,iflim)
       ENDDO
      ENDIF
      ENDIF_MPI
    ENDDO
  nb_reac_add=nb_reac_add+nb_reac_genco(ico)

! addition of reactions type 3 (REVERS)
! ===========================================================================
  it=nb_spc_diss(ico)+nb_spc_part_eq(ico)
  DO ispc=it+1,it+nb_spc_part_ci(ico)
   ir=nb_reac(3)+1
   nb_reac(3)=nb_reac(3)+1
   ire=ire+1
   nb_var_re(ire)=2
   cf_rate_r(ir)=cf_rate_spc(ico,ispc)
   cv_watequ_r(ir)=cv_watequ_spc(ico,ispc)
   l_arh(ire)=.false.
   l_treat_re(ire)=l_treat_spc(ico,ispc)
   ityp_reac(ire)=3
   ire_typ(ire)=nb_reac(3)
   irkm_var_re(ire,1)=irk_mod(ivar_spc(ico,ispc))
   irkf_var_re(ire,1)=ivar_spc(ico,ispc)
   nb_reac_pervar(ivar_spc(ico,ispc))=nb_reac_pervar(ivar_spc(ico,ispc))+1
   name_var_re(ire,1)=name_var(ivar_spc(ico,ispc))
   iv=ivar_spc(ico,ispc_assoc_pcin(ico,ispc))
   irkm_var_re(ire,2)=irk_mod(iv)
   irkf_var_re(ire,2)=iv
   nb_reac_pervar(iv)=nb_reac_pervar(iv)+1
   name_var_re(ire,2)=name_var(iv) 
   nb_reac_add=nb_reac_add+1
   name_reac(ire)='REVERSIBLE_'//TRIM(name_var_re(ire,1))//'_'//trim(name_var_re(ire,2))
   
! save into simu.log
!-------------------
   IF_MPI (MASTER) THEN
    WRITE(iscreenlog,*)
    WRITE(iscreenlog,*)'ADDITIONAL REACTION FOR CONTAMINANT '
    WRITE(iscreenlog,*)'REACTION : ',trim(name_reac(ire)),' TYPE :',ityp_reac(ire)
    WRITE(iscreenlog,*)'between : ',trim(name_var_re(ire,1))
    WRITE(iscreenlog,*)'    and : ',trim(name_var_re(ire,2))
    IF(l_treat_re(ire)) then
        WRITE(iscreenlog,*)'integration of sink & sources due to this reaction '
        WRITE(iscreenlog,*)'in transport equations '
    ELSE
        WRITE(iscreenlog,*)'sink & sources due to this reaction computed '
        WRITE(iscreenlog,*)' after (operator splitting)'
    ENDIF
    WRITE(iscreenlog,*)'  rate constant ',cf_rate_r(ir)
    WRITE(iscreenlog,*)'  equilibrium concentration ',cv_watequ_r(ir)
   ENDIF_MPI
  ENDDO
  
! addition of reactions type 4 (VOLATIL)
! ===========================================================================
   IF(l_volatil(ico))then
! always here volatilisation of the first species (free species)
     ire=ire+1
     nb_reac(4)=nb_reac(4)+1
     nb_var_re(ire)=1
     l_treat_re(ire)=.true.
     name_var_re(ire,1)=name_spc(ico,1)
     irkm_var_re(ire,1)=irk_mod(ivar_spc(ico,1))
     iv=irkm_var_re(ire,1)
     ityp_reac(ire)=4
     massmol(iv)=massmol_co(ico)
     nschmidt(iv)=nschmidt_co(ico)
     cf_difwat(iv)=cf_difwat_co(ico)
     cf_difair(iv)=cf_difair_co(ico)
     cf_henry(iv)=cf_henry_co(ico)
     cf_henry_adim(iv)=cf_henry_adim_co(ico)
     cf_temp_henry(iv)=cf_temp_henry_co(ico)
     cf_salt_henry(iv)=cf_salt_henry_co(ico)
     l_volkl_dif(iv)=l_volkl_difco(ico)
     l_volkl_sal(iv)=l_volkl_salco(ico)
     ich_volkl_wind(iv)=ich_volkl_windco(ico)
     l_volkl_cor(iv)=l_volkl_corco(ico)
     cf_expvol_liq(iv)=cf_expvol_liq_co(ico)
     l_volkg_cor(iv)=l_volkg_corco(ico)
     cf_expvol_gas(iv)=cf_expvol_gas_co(ico)    
     nb_reac_pervar(irk_fil(iv))=nb_reac_pervar(irk_fil(iv))+1
     nb_reac_add=nb_reac_add+1
     name_reac(ire)='VOLAILISATION_'//trim(name_var_re(ire,1))
     
! save into simu.log
!-------------------
   IF_MPI (MASTER) THEN
      WRITE(iscreenlog,*)
      WRITE(iscreenlog,*)'ADDITIONAL REACTION FOR CONTAMINANT '
      WRITE(iscreenlog,*)'REACTION : ',trim(name_reac(ire)),' TYPE :',ityp_reac(ire)
      WRITE(iscreenlog,*)'volatilisation of : ',trim(name_var_re(ire,1))
      WRITE(iscreenlog,*)'integration of sink & sources due to this reaction '
      WRITE(iscreenlog,*)'in transport equations '
   ENDIF_MPI
   ENDIF
     
! addition of reactions type 5 (PART_EQ)
! ===========================================================================
  IF(l_partag(ico))then
    ire=ire+1
    write(numcon,'(i1)')ico
    name_reac(ire)='PARTEq_'//numcon
    nb_var_re(ire)=nb_spc_part_eq(ico)+nb_spc_diss(ico)
    l_treat_re_after=.true.
! reALLOCATE arrays if nvaractmax too small
             IF(nb_var_re(ire).gt.nvaractmax) then
               ALLOCATE(name_var_tmp(nreacmax,nvaractmax))
               ALLOCATE(irkm_var_tmp(nreacmax,nvaractmax))
               ALLOCATE(irkf_var_tmp(nreacmax,nvaractmax))
               ALLOCATE(irkm_sub_tmp(nreacmax,nvaractmax))
               ALLOCATE(cf_distb_tmp(nreacmax,nvaractmax))
               ALLOCATE(ispc_reac_tmp(nreacmax,nvaractmax))
               name_var_tmp(:,:)=name_var_re(:,:)        
               irkm_var_tmp(:,:)=irkm_var_re(:,:) 
               irkf_var_tmp(:,:)=irkf_var_re(:,:) 
               irkm_sub_tmp(:,:)=irkm_sub_r(:,:) 
               cf_distb_tmp(:,:)=cf_distb_r(:,:) 
               ispc_reac_tmp(:,:)=ispc_reac(:,:) 
               deALLOCATE(name_var_re)
               deALLOCATE(irkm_var_re)
               deALLOCATE(irkf_var_re)
               deALLOCATE(irkm_sub_r)
               deALLOCATE(cf_distb_r)
               deALLOCATE(ispc_reac)
               nvaractmax_new=nb_var_re(ire)
               ALLOCATE(name_var_re(nreacmax,nvaractmax_new))
               ALLOCATE(irkm_var_re(nreacmax,nvaractmax_new))
               ALLOCATE(irkf_var_re(nreacmax,nvaractmax_new))
               ALLOCATE(irkm_sub_r(nreacmax,nvaractmax_new))
               ALLOCATE(cf_distb_r(nreacmax,nvaractmax_new))
               ALLOCATE(ispc_reac(nreacmax,nvaractmax_new))
               name_var_re(:,1:nvaractmax)=name_var_tmp(:,:)        
               irkm_var_re(:,1:nvaractmax)=irkm_var_tmp(:,:)        
               irkf_var_re(:,1:nvaractmax)=irkf_var_tmp(:,:)        
               irkm_sub_r(:,1:nvaractmax)=irkm_sub_tmp(:,:)        
               cf_distb_r(:,1:nvaractmax)=cf_distb_tmp(:,:)        
               ispc_reac(:,1:nvaractmax)=ispc_reac_tmp(:,:)        
               deALLOCATE(name_var_tmp)
               deALLOCATE(irkm_var_tmp)
               deALLOCATE(irkf_var_tmp)
               deALLOCATE(irkm_sub_tmp)
               deALLOCATE(cf_distb_tmp)
               deALLOCATE(ispc_reac_tmp)
#ifdef key_species
               ALLOCATE(ifrac_var_tmp(nreacmax,nvaractmax))
               ifrac_var_tmp(:,:)=ifrac_var_reac(:,:) 
               deALLOCATE(ifrac_var_reac)
               ALLOCATE(ifrac_var_reac(nreacmax,nvaractmax_new))
               ifrac_var_reac(:,1:nvaractmax)=ifrac_var_tmp(:,:)        
               deALLOCATE(ifracfracact_tmp)
#endif
               nvaractmax=nvaractmax_new                      
             ENDIF
    icon_reac(ire)=ico
    l_arh(ire)=l_arh_pco(ico)
    l_treat_re(ire)=.false.
    ityp_reac(ire)=5
    ityp_arh(ire)=ityp_arh_pco(ico)
    nb_reac(5)=nb_reac(5)+1
    ire_typ(ire)=nb_reac(5)
    cf_teta_arh(ire)=cf_teta_arh_pco(ico)
    tempref_arh(ire)=tempref_arh_pco(ico)
    cf_xetsig_arh(ire)=cf_xetsig_arh_pco(ico)
    l_pce_r(ire_typ(ire))=l_pce_co(ico)
    l_salt_r(ire_typ(ire))=l_salt_co(ico)
    cf_salt_r(ire_typ(ire),:)=cf_salt_co(ico,:)
    cf_pce_r(ire_typ(ire),:)=cf_pce_co(ico,:)
    l_kds_r(ire_typ(ire))= l_kds_metal(ico) 
    IF(l_metal(ico) .and. l_kds_metal(ico))then       
       ndg_kds_r(ire_typ(ire))= ndg_kds_metal(ico)        
       cf_kd0_r(ire_typ(ire))= cf_kd0_metal(ico)        
       cf_xkds_r(ire_typ(ire),1:ndg_kds_metal(ico))= cf_xkds_metal(ico,1:ndg_kds_metal(ico))
    ELSE
       ndg_kds_r(ire_typ(ire))= 1        
       cf_kd0_r(ire_typ(ire))= 0.       
       cf_xkds_r(ire_typ(ire),1)= 0.     
    ENDIF       
    l_distb_spc(ire)=.true.
    ivdl=ivdl_conta(ico)
    name_var_re(ire,1)=name_var(ivdl)
    irkm_var_re(ire,1)=irk_mod(ivdl)
    irkf_var_re(ire,1)=ivdl
    nb_reac_pervar(ivdl)=nb_reac_pervar(ivdl)+1
    it=1
#ifdef key_species
    ifrac_var_reac(ire,it)=1
    DO ispc=2,nb_spc_diss(ico)+nb_spc_part_eq(ico)
       it=it+1 
       irkm_var_re(ire,it)=irk_mod(ivar_spc(ico,ispc))
       irkf_var_re(ire,it)=ivar_spc(ico,ispc)
       ispc_reac(ire,it)=ispc
       cf_distb_r(ire_typ(ire),it)=cf_distb_spc(ico,ispc)
       irkm_sub_r(ire_typ(ire),it)=irk_mod(isit_spc(ico,ispc))
       DO iv=1,nv_adv
         IF (icon_reac(ire) == icon_var(iv)) then
            DO ifracv=1,nb_frac_var(iv) 
               IF(ispc_reac(ire,it) == ispc_var(iv,ifracv))then
                  ifrac_var_reac(ire,it)=ifracv
                  nb_reac_pervar(irk_fil(iv))=nb_reac_pervar(irk_fil(iv))+1
               ENDIF
            ENDDO
         ENDIF
       ENDDO
    ENDDO
#else
    DO ispc=2,nb_spc_diss(ico)+nb_spc_part_eq(ico)
       it=it+1 
       name_var_re(ire,it)=name_var(ivar_spc(ico,ispc))
       irkm_var_re(ire,it)=irk_mod(ivar_spc(ico,ispc))
       irkf_var_re(ire,it)=ivar_spc(ico,ispc)
       ispc_reac(ire,it)=ispc
       cf_distb_r(ire_typ(ire),it)=cf_distb_spc(ico,ispc)
       irkm_sub_r(ire_typ(ire),it)=irk_mod(isit_spc(ico,ispc))
       DO iv=1,nv_adv
         IF (icon_reac(ire) == icon_var(iv) .and. ispc_reac(ire,it) == ispc_var(iv))then
                  nb_reac_pervar(irk_fil(iv))=nb_reac_pervar(irk_fil(iv))+1
         ENDIF
       ENDDO
    ENDDO
#endif
    nb_reac_add=nb_reac_add+1
    
! save into simu.log
!-------------------
   IF_MPI (MASTER) THEN
    WRITE(iscreenlog,*)
    WRITE(iscreenlog,*)'ADDITIONAL REACTION FOR CONTAMINANT '
    WRITE(iscreenlog,*)'REACTION : ',trim(name_reac(ire)),' TYPE :',ityp_reac(ire)
    WRITE(iscreenlog,*)'sink & sources due to this reaction computed '
    WRITE(iscreenlog,*)' after (operator splitting)'
    IF(l_arh(ire)) then
       WRITE(iscreenlog,*)'type of arhenius law',ityp_arh(ire)
       WRITE(iscreenlog,*)'teta ', cf_teta_arh(ire)
       WRITE(iscreenlog,*)'reference temperature  ', tempref_arh(ire)
       WRITE(iscreenlog,*)'xeta or sigma ', cf_xetsig_arh(ire)
    ENDIF 
    DO it=1,nb_var_re(ire)
      iv=irkm_var_re(ire,it)
      write(iscreenlog,*)'variable active :',trim(name_var(irk_fil(iv)))
      write(iscreenlog,*)'distribution coefficient :',cf_distb_r(ire_typ(ire),it)
    ENDDO 
   ENDIF_MPI 
  ENDIF          
  
! fin boucle sur ico
 ENDDO
   
  END SUBROUTINE conta_add_reac
   !!======================================================================
      
  SUBROUTINE conta_deALLOCATE_data

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE conta_deALLOCATE_data  ***
   !&E
   !&E ** Purpose : deallocation of arrays used in contareaddat and unused after  
   !&E
   !&E ** Description :
   !&E
   !&E ** used ij-arrays : none
   !&E
   !&E ** modified variables : none
   !&E
   !&E ** Called by : conta_read_react.F90
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2008-06  (B. Thouvenin)  Original code
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   !! * Arguments

   !! * Local declarations

      deALLOCATE(irkm_sub_r)
      deALLOCATE(irkf_var_re)
      deALLOCATE(nb_reac_genco)
      deALLOCATE(nb_flim1_co)
      deALLOCATE(nb_flim2_co)
      deALLOCATE(nb_flim3_co)
      deALLOCATE(ityp_arh_rco)
      deALLOCATE(nb_var_co)
      deALLOCATE(ndg_kds_r)
      deALLOCATE(ndg_kds_metal)
      deALLOCATE(ich_volkl_windco)
      deALLOCATE( ityp_arh_pco)
      deALLOCATE(cin_hydrol)
      deALLOCATE(cin_photo)
      deALLOCATE(cf_kow_corg)
      deALLOCATE(solubil_co)
      deALLOCATE(cin_biodeg)
      deALLOCATE(cf_teta_arh_pco)
      deALLOCATE(tempref_arh_pco)
      deALLOCATE( cf_xetsig_arh_pco)
      deALLOCATE(cf_distb_spc)
      deALLOCATE(cf_rate_spc)
      deALLOCATE(cv_watequ_spc)
      deALLOCATE(cf_kd0_metal)
      deALLOCATE(cv_watequ_r)
      deALLOCATE(cf_pce_co)
      deALLOCATE(cf_salt_co)
      deALLOCATE(cf_pce_r)
      deALLOCATE(cf_salt_r)
      deALLOCATE(cf_distb_r)
      deALLOCATE(cf_kd0_r)
      deALLOCATE(cf_rate_r)
      deALLOCATE(cf_xkds_r)
      deALLOCATE(massmol_co)
      deALLOCATE(nschmidt_co)
      deALLOCATE(cf_difair_co)
      deALLOCATE(cf_difwat_co)
      deALLOCATE(cf_temp_henry_co)
      deALLOCATE(cf_salt_henry_co)
      deALLOCATE(cf_henry_co)
      deALLOCATE(cf_henry_adim_co)
      deALLOCATE(cin_max_co)
      deALLOCATE(cf_teta_arh_rco)
      deALLOCATE(tempref_arh_rco)
      deALLOCATE(cf_xkds_metal)
      deALLOCATE(cf_xetsig_arh_rco)
      deALLOCATE(cf_stochio_co)
      deALLOCATE(cf_expovar_co)
      deALLOCATE(cf_flim1_co)
      deALLOCATE(cf_flim2_co)
      deALLOCATE(cf_flim3_co)
      deALLOCATE(cf_expvol_liq_co)
      deALLOCATE(cf_expvol_gas_co)
      deALLOCATE(name_sit_spc)
      deALLOCATE(name_reac_co)
      deALLOCATE(name_var_co)
      deALLOCATE(name_var_flim1)
      deALLOCATE(name_var_flim2)
      deALLOCATE(fileconta)
      deALLOCATE(l_partag)
      deALLOCATE(l_arh_pco)
      deALLOCATE(l_pce_co)
      deALLOCATE(l_salt_co)
      deALLOCATE(l_kds_metal)
      deALLOCATE(l_treat_spc)
      deALLOCATE(l_pce_r)
      deALLOCATE(l_salt_r)
      deALLOCATE(l_distb_spc)
      deALLOCATE(l_kds_r)
      deALLOCATE(l_arh_rco)
      deALLOCATE(l_treat_co)
      deALLOCATE(l_cincte_co)
      deALLOCATE(l_coefcte_co)
      deALLOCATE(l_volkl_difco)
      deALLOCATE(l_volkl_salco)
      deALLOCATE(l_volkl_corco)
      deALLOCATE(l_volkg_corco)

  END SUBROUTINE conta_deALLOCATE_data
   !!======================================================================
#endif
 END MODULE contareaddat

