Index: LMDZ6/trunk/libf/dyn3d/check_isotopes.F90
===================================================================
--- LMDZ6/trunk/libf/dyn3d/check_isotopes.F90	(revision 4319)
+++ LMDZ6/trunk/libf/dyn3d/check_isotopes.F90	(revision 4325)
@@ -2,5 +2,5 @@
    USE strings_mod, ONLY: maxlen, msg, strIdx, strStack, int2str, real2str
    USE infotrac,    ONLY: nqtot, niso, nphas, isotope, isoCheck, iqIsoPha, isoSelect, &
-                          ntiso, iH2O, nzone, tracers, isoName,  itZonIso, tnat
+                          ntiso, iH2O, nzone, tracers, isoName,  itZonIso, getKey
    IMPLICIT NONE
    include "dimensions.h"
@@ -10,5 +10,6 @@
    CHARACTER(LEN=maxlen) :: modname, msg1, nm(2)
    INTEGER :: ixt, ipha, k, i, iq, iiso, izon, ieau, iqeau, iqpar
-   INTEGER, ALLOCATABLE :: ix(:)
+   INTEGER, ALLOCATABLE ::   ix(:)
+   REAL,    ALLOCATABLE :: tnat(:)
    REAL    :: xtractot, xiiso, deltaD, q1, q2
    REAL, PARAMETER :: borne     = 1e19,  &
@@ -33,4 +34,5 @@
       iso_O17 = strIdx(isoName,'H2[17]O')
       iso_HTO = strIdx(isoName,'H[3]HO')
+      IF(getKey('tnat', tnat)) CALL abort_gcm(modname, 'missing isotopic parameter', 1)
       first = .FALSE.
    END IF
Index: LMDZ6/trunk/libf/dyn3d/dynetat0.F90
===================================================================
--- LMDZ6/trunk/libf/dyn3d/dynetat0.F90	(revision 4319)
+++ LMDZ6/trunk/libf/dyn3d/dynetat0.F90	(revision 4325)
@@ -6,9 +6,9 @@
 ! Purpose: Initial state reading.
 !-------------------------------------------------------------------------------
-  USE infotrac,    ONLY: nqtot, tracers, niso, iqIsoPha, tnat, alpha_ideal, iH2O
+  USE infotrac,    ONLY: nqtot, tracers, niso, iqIsoPha, iH2O, isoName
   USE strings_mod, ONLY: maxlen, msg, strStack, real2str, int2str
   USE netcdf,      ONLY: NF90_OPEN,  NF90_NOWRITE, NF90_INQ_VARID, &
                          NF90_CLOSE, NF90_GET_VAR, NF90_NoErr
-  USE readTracFiles_mod, ONLY: new2oldH2O, newHNO3, oldHNO3
+  USE readTracFiles_mod, ONLY: new2oldH2O, newHNO3, oldHNO3, getKey
   USE control_mod, ONLY: planet_type
   USE assert_eq_m, ONLY: assert_eq
@@ -41,5 +41,5 @@
   INTEGER, PARAMETER :: length=100
   INTEGER :: iq, fID, vID, idecal, iqParent, iName, iZone, iPhase
-  REAL    :: time, tab_cntrl(length)               !--- RUN PARAMS TABLE
+  REAL    :: time, tnat, alpha_ideal, tab_cntrl(length)    !--- RUN PARAMS TABLE
   LOGICAL :: lSkip, ll
 !-------------------------------------------------------------------------------
@@ -155,6 +155,8 @@
       iqParent = tracers(iq)%iqParent
       IF(tracers(iq)%iso_iZone == 0) THEN
+         IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) &
+            CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1)
          CALL msg('Tracer <'//TRIM(var)//'> is missing => initialized with a simplified Rayleigh distillation law.', modname)
-         q(:,:,:,iq) = q(:,:,:,iqParent)*tnat(iName)*(q(:,:,:,iqParent)/30.e-3)**(alpha_ideal(iName)-1.)
+         q(:,:,:,iq) = q(:,:,:,iqParent)*tnat*(q(:,:,:,iqParent)/30.e-3)**(alpha_ideal-1.)
       ELSE
          CALL msg('Tracer <'//TRIM(var)//'> is missing => initialized to its parent isotope concentration.', modname)
Index: LMDZ6/trunk/libf/dyn3d/gcm.F90
===================================================================
--- LMDZ6/trunk/libf/dyn3d/gcm.F90	(revision 4319)
+++ LMDZ6/trunk/libf/dyn3d/gcm.F90	(revision 4325)
@@ -20,5 +20,5 @@
 
   USE filtreg_mod
-  USE infotrac
+  USE infotrac, ONLY: nqtot, init_infotrac
   USE control_mod
   USE mod_const_mpi, ONLY: COMM_LMDZ
@@ -205,5 +205,5 @@
   !  Choix du nombre de traceurs et du schema pour l'advection
   !  dans fichier traceur.def, par default ou via INCA
-  call infotrac_init
+  call init_infotrac
 
   ! Allocation de la tableau q : champs advectes   
Index: LMDZ6/trunk/libf/dyn3d/iniacademic.F90
===================================================================
--- LMDZ6/trunk/libf/dyn3d/iniacademic.F90	(revision 4319)
+++ LMDZ6/trunk/libf/dyn3d/iniacademic.F90	(revision 4325)
@@ -5,5 +5,5 @@
 
   USE filtreg_mod, ONLY: inifilr
-  USE infotrac,    ONLY: nqtot, niso, tnat, alpha_ideal, iqIsoPha, tracers
+  USE infotrac,    ONLY: nqtot, niso, iqIsoPha, tracers, getKey, isoName
   USE control_mod, ONLY: day_step,planet_type
   use exner_hyb_m, only: exner_hyb
@@ -73,5 +73,5 @@
   integer idum
 
-  REAL zdtvr
+  REAL zdtvr, tnat, alpha_ideal
   
   character(len=*),parameter :: modname="iniacademic"
@@ -286,5 +286,7 @@
               iqParent = tracers(iq)%iqParent
               IF(tracers(iq)%iso_iZone == 0) THEN
-                 q(:,:,iq) = q(:,:,iqParent)*tnat(iName)*(q(:,:,iqParent)/30.e-3)**(alpha_ideal(iName)-1.)
+                 IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) &
+                    CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1)
+                 q(:,:,iq) = q(:,:,iqParent)*tnat*(q(:,:,iqParent)/30.e-3)**(alpha_ideal-1.)
               ELSE
                  q(:,:,iq) = q(:,:,iqIsoPha(iName,iPhase))
Index: LMDZ6/trunk/libf/dyn3d/vlsplt.F
===================================================================
--- LMDZ6/trunk/libf/dyn3d/vlsplt.F	(revision 4319)
+++ LMDZ6/trunk/libf/dyn3d/vlsplt.F	(revision 4325)
@@ -437,5 +437,5 @@
         enddo
       enddo
-      do ifils=1,tracers(iq)%nqChilds
+      do ifils=1,tracers(iq)%nqChildren
         iq2=tracers(iq)%iqDescen(ifils)
         call vlx(Ratio,pente_max,masseq,u_mq,iq2)
@@ -969,5 +969,5 @@
 ! CRisi: appel récursif de l'advection sur les fils.
 ! Il faut faire ça avant d'avoir mis à jour q et masse
-      !write(*,*) 'vlsplt 942: iq,nqChilds(iq)=',iq,nqChilds(iq)
+      !write(*,*) 'vlsplt 942: iq,nqChildren(iq)=',iq,nqChildren(iq)
       do ifils=1,tracers(iq)%nqDescen
         iq2=tracers(iq)%iqDescen(ifils)
@@ -987,5 +987,5 @@
       enddo
         
-      do ifils=1,tracers(iq)%nqChilds
+      do ifils=1,tracers(iq)%nqChildren
         iq2=tracers(iq)%iqDescen(ifils)
         call vlz(Ratio,pente_max,masseq,wq,iq2)
Index: LMDZ6/trunk/libf/dyn3d/vlspltqs.F
===================================================================
--- LMDZ6/trunk/libf/dyn3d/vlspltqs.F	(revision 4319)
+++ LMDZ6/trunk/libf/dyn3d/vlspltqs.F	(revision 4325)
@@ -479,5 +479,6 @@
 ! CRisi: appel récursif de l'advection sur les fils.
 ! Il faut faire ça avant d'avoir mis à jour q et masse
-      !write(*,*) 'vlspltqs 326: iq,nqChilds(iq)=',iq,tracers(iq)%nqChilds
+      !write(*,*) 'vlspltqs 326: iq,nqChildren(iq)=',iq,
+!     &                 tracers(iq)%nqChildren
       
       do ifils=1,tracers(iq)%nqDescen
@@ -491,5 +492,5 @@
         enddo
       enddo
-      do ifils=1,tracers(iq)%nqChilds
+      do ifils=1,tracers(iq)%nqChildren
         iq2=tracers(iq)%iqDescen(ifils)
         call vlx(Ratio,pente_max,masseq,u_mq,iq2)
@@ -786,5 +787,6 @@
 ! CRisi: appel récursif de l'advection sur les fils.
 ! Il faut faire ça avant d'avoir mis à jour q et masse
-      !write(*,*) 'vlyqs 689: iq,nqChilds(iq)=',iq,tracers(iq)%nqChilds
+      !write(*,*) 'vlyqs 689: iq,nqChildren(iq)=',iq,
+!     &              tracers(iq)%nqChildren
    
       do ifils=1,tracers(iq)%nqDescen
@@ -797,5 +799,5 @@
         enddo
       enddo
-      do ifils=1,tracers(iq)%nqChilds
+      do ifils=1,tracers(iq)%nqChildren
         iq2=tracers(iq)%iqDescen(ifils)
         !write(*,*) 'vlyqs 783: appel rec de vly, iq2=',iq2
Index: LMDZ6/trunk/libf/dyn3d_common/infotrac.F90
===================================================================
--- LMDZ6/trunk/libf/dyn3d_common/infotrac.F90	(revision 4319)
+++ LMDZ6/trunk/libf/dyn3d_common/infotrac.F90	(revision 4325)
@@ -3,9 +3,8 @@
 MODULE infotrac
 
-   USE       strings_mod, ONLY: msg, find, strIdx,  strFind, strParse, dispTable, int2str,  reduceExpr, &
-                          cat, fmsg, test, strTail, strHead, strStack, strReduce, bool2str, maxlen, testFile
-   USE readTracFiles_mod, ONLY: trac_type, readTracersFiles, indexUpdate, addPhase, getKey, maxTableWidth, keys_type, &
-                                isot_type, setGeneration,   initIsotopes, delPhase, getKey_init, ancestor, tran0
-                                
+   USE       strings_mod, ONLY: msg, fmsg, maxlen, cat, dispTable, int2str, bool2str, strStack, strParse
+   USE readTracFiles_mod, ONLY: trac_type, readTracersFiles, tracers, setGeneration, itZonIso, nbIso, tran0, delPhase, &
+                        getKey, isot_type, readIsotopesFile, isotope, maxTableWidth, iqIsoPha, ntiso, ixIso, addPhase, &
+                   indexUpdate, isoSelect, isoPhas, isoZone, isoName, isoKeys, iH2O, isoCheck, nphas, nzone, niso
    IMPLICIT NONE
 
@@ -13,5 +12,5 @@
 
    !=== FOR TRACERS:
-   PUBLIC :: infotrac_init                                 !--- Initialization of the tracers
+   PUBLIC :: init_infotrac                                 !--- Initialization of the tracers
    PUBLIC :: tracers, type_trac, types_trac                !--- Full tracers database, tracers type keyword
    PUBLIC :: nqtot,   nbtr,   nqo,   nqCO2,   nqtottr      !--- Main dimensions
@@ -19,8 +18,8 @@
 
    !=== FOR ISOTOPES: General
-   PUBLIC :: isotopes,  nbIso                              !--- Derived type, full isotopes families database + nb of families
+   PUBLIC :: isot_type, nbIso                              !--- Derived type, full isotopes families database + nb of families
    PUBLIC :: isoSelect, ixIso                              !--- Isotopes family selection tool + selected family index
    !=== FOR ISOTOPES: Specific to water
-   PUBLIC :: iH2O, tnat, alpha_ideal                       !--- H2O isotopes index, natural abundance, fractionning coeff.
+   PUBLIC :: iH2O                                          !--- H2O isotopes class index
    PUBLIC :: min_qParent, min_qMass, min_ratio             !--- Min. values for various isotopic quantities
    !=== FOR ISOTOPES: Depending on the selected isotopes family
@@ -33,6 +32,4 @@
    !=== FOR BOTH TRACERS AND ISOTOPES
    PUBLIC :: getKey                                        !--- Get a key from "tracers" or "isotope"
-
-   INTERFACE isoSelect; MODULE PROCEDURE isoSelectByIndex, isoSelectByName; END INTERFACE isoSelect
 
 !=== CONVENTIONS FOR TRACERS NUMBERS:
@@ -77,5 +74,5 @@
 !  | iqDescen    | Indexes of the childs       (all generations)        | iqfils      | 1:nqtot                |
 !  | nqDescen    | Number of the descendants   (all generations)        | nqdesc      | 1:nqtot                |
-!  | nqChilds    | Number of childs            (1st generation only)    | nqfils      | 1:nqtot                |
+!  | nqChildren  | Number of childs            (1st generation only)    | nqfils      | 1:nqtot                |
 !  | iso_iGroup  | Isotopes group index in isotopes(:)                  | /           | 1:nbIso                |
 !  | iso_iName   | Isotope  name  index in isotopes(iso_iGroup)%trac(:) | iso_indnum  | 1:niso                 |
@@ -102,39 +99,19 @@
 
    !=== DIMENSIONS OF THE TRACERS TABLES AND OTHER SCALAR VARIABLES
-   INTEGER,                 SAVE :: nqtot,  &                   !--- Tracers nb in dynamics (incl. higher moments + H2O)
-                                    nbtr,   &                   !--- Tracers nb in physics  (excl. higher moments + H2O)
-                                    nqo,    &                   !--- Number of water phases
-                                    nbIso,  &                   !--- Number of available isotopes family
-                                    nqtottr, &                  !--- Number of tracers passed to phytrac (TO BE DELETED ?)
-                                    nqCO2                       !--- Number of tracers of CO2  (ThL)
-   CHARACTER(LEN=maxlen),   SAVE :: type_trac                   !--- Keyword for tracers type(s)
-   CHARACTER(LEN=maxlen),   SAVE, ALLOCATABLE :: types_trac(:)  !--- Keyword for tracers type(s), parsed version
-
-   !=== DERIVED TYPES EMBEDDING MOST INFORMATIONS ABOUT TRACERS AND ISOTOPES FAMILIES
-   TYPE(trac_type), TARGET, SAVE, ALLOCATABLE ::  tracers(:)    !=== TRACERS DESCRIPTORS VECTOR
-   TYPE(isot_type), TARGET, SAVE, ALLOCATABLE :: isotopes(:)    !=== ISOTOPES PARAMETERS VECTOR
-
-   !=== ALIASES FOR CURRENTLY SELECTED ISOTOPES FAMILY OF VARIABLES EMBEDDED IN THE VECTOR "isotopes"
-   TYPE(isot_type),         SAVE, POINTER :: isotope            !--- CURRENTLY SELECTED ISOTOPES FAMILY DESCRIPTOR
-   INTEGER,                 SAVE          :: ixIso, iH2O        !--- Index of the selected isotopes family and H2O family
-   LOGICAL,                 SAVE          :: isoCheck           !--- Flag to trigger the checking routines
-   TYPE(keys_type),         SAVE, POINTER :: isoKeys(:)         !--- ONE SET OF KEYS FOR EACH ISOTOPE (LISTED IN isoName)
-   CHARACTER(LEN=maxlen),   SAVE, POINTER :: isoName(:),   &    !--- ISOTOPES NAMES FOR THE CURRENTLY SELECTED FAMILY
-                                             isoZone(:),   &    !--- TAGGING ZONES  FOR THE CURRENTLY SELECTED FAMILY
-                                             isoPhas            !--- USED PHASES    FOR THE CURRENTLY SELECTED FAMILY
-   INTEGER,                 SAVE          ::  niso, nzone, &    !--- NUMBER OF ISOTOPES, TAGGING ZONES AND PHASES
-                                             nphas, ntiso       !--- NUMBER OF PHASES AND ISOTOPES + ISOTOPIC TAGGING TRACERS
-   INTEGER,                 SAVE, POINTER ::itZonIso(:,:), &    !--- INDEX IN "isoTrac" AS f(tagging zone idx,  isotope idx)
-                                            iqIsoPha(:,:)       !--- INDEX IN "qx"      AS f(isotopic tracer idx, phase idx)
-
-   !=== VARIABLES FOR ISOTOPES INITIALIZATION AND FOR INCA
-   REAL,                SAVE, ALLOCATABLE ::     tnat(:), &     !--- Natural relative abundance of water isotope        (niso)
-                                          alpha_ideal(:)        !--- Ideal fractionning coefficient (for initial state) (niso)
-   INTEGER,             SAVE, ALLOCATABLE :: conv_flg(:), &     !--- Convection     activation ; needed for INCA        (nbtr)
-                                              pbl_flg(:)        !--- Boundary layer activation ; needed for INCA        (nbtr)
+   INTEGER,               SAVE :: nqtot,  &                     !--- Tracers nb in dynamics (incl. higher moments + H2O)
+                                  nbtr,   &                     !--- Tracers nb in physics  (excl. higher moments + H2O)
+                                  nqo,    &                     !--- Number of water phases
+                                  nqtottr, &                    !--- Number of tracers passed to phytrac (TO BE DELETED ?)
+                                  nqCO2                         !--- Number of tracers of CO2  (ThL)
+   CHARACTER(LEN=maxlen), SAVE :: type_trac                     !--- Keyword for tracers type(s)
+   CHARACTER(LEN=maxlen), SAVE, ALLOCATABLE :: types_trac(:)    !--- Keyword for tracers type(s), parsed version
+
+   !=== VARIABLES FOR INCA
+   INTEGER,               SAVE, ALLOCATABLE :: conv_flg(:), &   !--- Convection     activation ; needed for INCA        (nbtr)
+                                                pbl_flg(:)      !--- Boundary layer activation ; needed for INCA        (nbtr)
 
 CONTAINS
 
-SUBROUTINE infotrac_init
+SUBROUTINE init_infotrac
    USE control_mod, ONLY: planet_type, config_inca
 #ifdef REPROBUS
@@ -180,13 +157,12 @@
    INTEGER :: nqtrue                                                 !--- Tracers nb from tracer.def (no higher order moments)
    INTEGER :: iad                                                    !--- Advection scheme number
-   INTEGER :: ic, ip, np, iq, jq, it, nt, im, nm, ix, iz, nz, k      !--- Indexes and temporary variables
+   INTEGER :: ic, iq, jq, it, nt, im, nm, iz, k                      !--- Indexes and temporary variables
    LOGICAL :: lerr, ll, lRepr
    CHARACTER(LEN=1) :: p
    TYPE(trac_type), ALLOCATABLE, TARGET :: ttr(:)
    TYPE(trac_type), POINTER             :: t1, t(:)
-   TYPE(isot_type), POINTER             :: iso
    INTEGER :: ierr
 
-   CHARACTER(LEN=*), PARAMETER :: modname="infotrac_init"
+   CHARACTER(LEN=*), PARAMETER :: modname="init_infotrac"
 !------------------------------------------------------------------------------------------------------------------------------
 ! Initialization :
@@ -249,5 +225,5 @@
 !==============================================================================================================================
    lRepr = ANY(types_trac(:) == 'repr')
-   IF(readTracersFiles(type_trac, fType, tracers, lRepr)) CALL abort_gcm(modname, 'problem with tracers file(s)',1)
+   IF(readTracersFiles(type_trac, fType, lRepr)) CALL abort_gcm(modname, 'problem with tracers file(s)',1)
    !---------------------------------------------------------------------------------------------------------------------------
    IF(fType == 0) CALL abort_gcm(modname, 'Missing tracers file: "traceur.def", "tracer.def" or "tracer_<keyword>.def file.',1)
@@ -297,6 +273,4 @@
    !---------------------------------------------------------------------------------------------------------------------------
 
-   CALL getKey_init(tracers)
-
    !--- Transfert the number of tracers to Reprobus
 #ifdef REPROBUS
@@ -377,5 +351,5 @@
    CALL MOVE_ALLOC(FROM=ttr, TO=tracers)
 
-   !--- SET FIELDS %iqParent, %nqChilds, %iGeneration, %iqDescen, %nqDescen
+   !--- SET FIELDS %iqParent, %nqChildren, %iGeneration, %iqDescen, %nqDescen
    CALL indexUpdate(tracers)
 
@@ -401,36 +375,7 @@
    END DO
 
-   niso = 0; nzone=0; nphas=nqo; ntiso = 0; isoCheck=.FALSE.
-   IF(initIsotopes(tracers, isotopes)) CALL abort_gcm(modname, 'Problem when reading isotopes parameters', 1)
-   nbIso = SIZE(isotopes)
-   nqtottr = nqtot - COUNT(delPhase(tracers%gen0Name) == 'H2O' .AND. tracers%component == 'lmdz')
-   IF(nbIso/=0) THEN                        !--- ISOTOPES FOUND
-
-      !=== READ PHYSICAL PARAMETERS FROM "isotopes_params.def" FILE SPECIFIC TO WATER ISOTOPES
-      !    DONE HERE, AND NOT ONLY IN "infotrac_phy", BECAUSE SOME PHYSICAL PARAMS ARE NEEDED FOR RESTARTS (tnat, alpha_ideal)
-      CALL getKey_init(tracers, isotopes)
-      IF(isoSelect('H2O')) RETURN           !--- Select water isotopes ; finished if no water isotopes
-      iH2O = ixIso                          !--- Keep track of water family index
-      IF(getKey('tnat' , tnat,        isoName(1:niso))) CALL abort_gcm(modname, 'can''t read "tnat"', 1)
-      IF(getKey('alpha', alpha_ideal, isoName(1:niso))) CALL abort_gcm(modname, 'can''t read "alpha_ideal"', 1)
-
-      !=== MAKE SURE THE MEMBERS OF AN ISOTOPES FAMILY ARE PRESENT IN THE SAME PHASES
-      DO ix = 1, nbIso
-         iso => isotopes(ix)
-         !--- Check whether each isotope and tagging isotopic tracer is present in the same number of phases
-         DO it = 1, iso%ntiso
-            np = SUM([(COUNT(tracers(:)%name == addPhase(iso%trac(it), iso%phase(ip:ip))), ip=1, iso%nphas)])
-            IF(np == iso%nphas) CYCLE
-            WRITE(msg1,'("Found ",i0," phases for ",a," instead of ",i0)')np, TRIM(iso%trac(it)), iso%nphas
-            CALL abort_gcm(modname, msg1, 1)
-         END DO
-         DO it = 1, iso%niso
-            nz = SUM([(COUNT(iso%trac == TRIM(iso%trac(it))//'_'//iso%zone(iz)), iz=1, iso%nzone)])
-            IF(nz == iso%nzone) CYCLE
-            WRITE(msg1,'("Found ",i0," tagging zones for ",a," instead of ",i0)')nz, TRIM(iso%trac(it)), iso%nzone
-            CALL abort_gcm(modname, msg1, 1)
-         END DO
-      END DO
-   END IF
+   !=== READ PHYSICAL PARAMETERS FOR ISOTOPES ; DONE HERE BECAUSE dynetat0 AND iniacademic NEED "tnat" AND "alpha_ideal"
+   niso = 0; nzone = 0; nphas = nqo; ntiso = 0; isoCheck = .FALSE.
+   IF(readIsotopesFile()) CALL abort_gcm(modname, 'Problem when reading isotopes parameters', 1)
 
    !--- Convection / boundary layer activation for all tracers
@@ -439,31 +384,34 @@
 
    !--- Note: nqtottr can differ from nbtr when nmom/=0
+   nqtottr = nqtot - COUNT(delPhase(tracers%gen0Name) == 'H2O' .AND. tracers%component == 'lmdz')
    IF(COUNT(tracers%iso_iName == 0) - COUNT(delPhase(tracers%name) == 'H2O' .AND. tracers%component == 'lmdz') /= nqtottr) &
-      CALL abort_gcm('infotrac_init', 'pb dans le calcul de nqtottr', 1)
+      CALL abort_gcm(modname, 'pb dans le calcul de nqtottr', 1)
 
    !=== DISPLAY THE RESULTS
-   CALL msg('nqo    = '//TRIM(int2str(nqo)),    modname)
-   CALL msg('nbtr   = '//TRIM(int2str(nbtr)),   modname)
-   CALL msg('nqtrue = '//TRIM(int2str(nqtrue)), modname)
-   CALL msg('nqtot  = '//TRIM(int2str(nqtot)),  modname)
-   CALL msg('niso   = '//TRIM(int2str(niso)),   modname)
-   CALL msg('ntiso  = '//TRIM(int2str(ntiso)),  modname)
+   IF(prt_level > 1) THEN
+      CALL msg('nqo    = '//TRIM(int2str(nqo)),    modname)
+      CALL msg('nbtr   = '//TRIM(int2str(nbtr)),   modname)
+      CALL msg('nqtrue = '//TRIM(int2str(nqtrue)), modname)
+      CALL msg('nqtot  = '//TRIM(int2str(nqtot)),  modname)
+      CALL msg('niso   = '//TRIM(int2str(niso)),   modname)
+      CALL msg('ntiso  = '//TRIM(int2str(ntiso)),  modname)
 #ifdef INCA
-   CALL msg('nqCO2  = '//TRIM(int2str(nqCO2)),  modname)
-   CALL msg('nqINCA = '//TRIM(int2str(nqINCA)), modname)
-#endif
+      CALL msg('nqCO2  = '//TRIM(int2str(nqCO2)),  modname)
+      CALL msg('nqINCA = '//TRIM(int2str(nqINCA)), modname)
+#endif
+   END IF
    t => tracers
    CALL msg('Information stored in infotrac :', modname)
    IF(dispTable('isssssssssiiiiiiiii', &
-      ['iq    ', 'name  ', 'lName ', 'gen0N ', 'parent', 'type  ', 'phase ', 'compon', 'isAdv ', 'isPhy ', &
+      ['iq    ', 'name  ', 'lName ', 'gen0N ', 'parent', 'type  ', 'phase ', 'compon', 'isPhy ', 'isAdv ', &
        'iadv  ', 'iGen  ', 'iqPar ', 'nqDes ', 'nqChld', 'iGroup', 'iName ', 'iZone ', 'iPhase'],          &
-      cat(t%name, t%longName, t%gen0Name, t%parent, t%type, t%phase, t%component, bool2str(t%isAdvected),  &
-                                                                                  bool2str(t%isInPhysics)),&
-      cat([(iq, iq=1, nqtot)], t%iadv, t%iGeneration, t%iqParent, t%nqDescen, t%nqChilds, t%iso_iGroup,    &
+      cat(t%name, t%longName, t%gen0Name, t%parent, t%type, t%phase, t%component, bool2str(t%isInPhysics), &
+                                                                                  bool2str(t%isAdvected)), &
+      cat([(iq, iq=1, nqtot)], t%iadv, t%iGeneration, t%iqParent, t%nqDescen, t%nqChildren, t%iso_iGroup,  &
                   t%iso_iName, t%iso_iZone, t%iso_iPhase), nColMax=maxTableWidth, nHead=2, sub=modname))   &
       CALL abort_gcm(modname, "problem with the tracers table content", 1)
    IF(niso > 0) THEN
       CALL msg('Where, for isotopes family "'//TRIM(isotope%parent)//'":', modname)
-      CALL msg('  isoKeys = '//strStack(isoKeys%name), modname)
+      CALL msg('  isoKeys%name = '//strStack(isoKeys%name), modname)
       CALL msg('  isoName = '//strStack(isoName),      modname)
       CALL msg('  isoZone = '//strStack(isoZone),      modname)
@@ -474,50 +422,5 @@
    CALL msg('end', modname)
 
-END SUBROUTINE infotrac_init
-
-
-!==============================================================================================================================
-!=== THE ROUTINE isoSelect IS USED TO SWITCH FROM AN ISOTOPE FAMILY TO ANOTHER: ISOTOPES DEPENDENT PARAMETERS ARE UPDATED
-!     Single generic "isoSelect" routine, using the predefined index of the parent (fast version) or its name (first call).
-!==============================================================================================================================
-LOGICAL FUNCTION isoSelectByName(iName, lVerbose) RESULT(lerr)
-   IMPLICIT NONE
-   CHARACTER(LEN=*),  INTENT(IN) :: iName
-   LOGICAL, OPTIONAL, INTENT(IN) :: lVerbose
-   INTEGER :: iIso
-   LOGICAL :: lV
-   lV = .FALSE.; IF(PRESENT(lVerbose)) lV = lVerbose
-   iIso = strIdx(isotopes(:)%parent, iName)
-   lerr = iIso == 0
-   IF(lerr) THEN
-      niso = 0; ntiso = 0; nzone=0; nphas=nqo; isoCheck=.FALSE.
-      CALL msg('no isotope family named "'//TRIM(iName)//'"', ll=lV)
-      RETURN
-   END IF
-   lerr = isoSelectByIndex(iIso, lV)
-END FUNCTION isoSelectByName
-!==============================================================================================================================
-LOGICAL FUNCTION isoSelectByIndex(iIso, lVerbose) RESULT(lerr)
-   IMPLICIT NONE
-   INTEGER,           INTENT(IN) :: iIso
-   LOGICAL, OPTIONAL, INTENT(IN) :: lVerbose
-   LOGICAL :: lv
-   lv = .FALSE.; IF(PRESENT(lVerbose)) lv = lVerbose
-   lerr = .FALSE.
-   IF(iIso == ixIso) RETURN                                          !--- Nothing to do if the index is already OK
-   lerr = iIso<=0 .OR. iIso>nbIso
-   CALL msg('Inconsistent isotopes family index '//TRIM(int2str(iIso))//': should be > 0 and <= '//TRIM(int2str(nbIso))//'"',&
-            ll=lerr .AND. lV)
-   IF(lerr) RETURN
-   ixIso = iIso                                                      !--- Update currently selected family index
-   isotope  => isotopes(ixIso)                                       !--- Select corresponding component
-   isoKeys  => isotope%keys;     niso     = isotope%niso
-   isoName  => isotope%trac;     ntiso    = isotope%ntiso
-   isoZone  => isotope%zone;     nzone    = isotope%nzone
-   isoPhas  => isotope%phase;    nphas    = isotope%nphas
-   itZonIso => isotope%itZonIso; isoCheck = isotope%check
-   iqIsoPha => isotope%iqIsoPha
-END FUNCTION isoSelectByIndex
-!==============================================================================================================================
+END SUBROUTINE init_infotrac
 
 END MODULE infotrac
Index: LMDZ6/trunk/libf/dyn3d_common/iso_verif_dyn.F
===================================================================
--- LMDZ6/trunk/libf/dyn3d_common/iso_verif_dyn.F	(revision 4319)
+++ LMDZ6/trunk/libf/dyn3d_common/iso_verif_dyn.F	(revision 4325)
@@ -64,5 +64,5 @@
         function iso_verif_aberrant_nostop
      :           (x,iso,q,err_msg)
-        USE infotrac, ONLY: tnat
+        USE infotrac, ONLY: isoName, getKey
         implicit none
         
@@ -74,5 +74,5 @@
         ! locals
         real qmin,deltaD
-        real deltaDmax,deltaDmin
+        real deltaDmax,deltaDmin,tnat
         parameter (qmin=1e-11)
         parameter (deltaDmax=200.0,deltaDmin=-999.9)
@@ -85,5 +85,10 @@
         ! verifier que HDO est raisonable
          if (q.gt.qmin) then
-             deltaD=(x/q/tnat(iso)-1)*1000
+             IF(getKey('tnat', tnat, isoName(iso))) THEN
+                  err_msg = 'Missing isotopic parameter "tnat"'
+                  iso_verif_aberrant_nostop=1
+                  RETURN
+             END IF
+             deltaD=(x/q/tnat-1)*1000
              if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then
                   write(*,*) 'erreur detectee par iso_verif_aberrant:'
Index: LMDZ6/trunk/libf/dyn3dmem/check_isotopes_loc.F90
===================================================================
--- LMDZ6/trunk/libf/dyn3dmem/check_isotopes_loc.F90	(revision 4319)
+++ LMDZ6/trunk/libf/dyn3dmem/check_isotopes_loc.F90	(revision 4325)
@@ -3,5 +3,5 @@
    USE strings_mod, ONLY: maxlen, msg, strIdx, strStack, int2str, real2str
    USE infotrac,    ONLY: nqtot, niso, nphas, isotope, isoCheck, iqIsoPha, isoSelect, &
-                          ntiso, iH2O, nzone, tracers, isoName,  itZonIso, tnat
+                          ntiso, iH2O, nzone, tracers, isoName,  itZonIso, getKey
    IMPLICIT NONE
    include "dimensions.h"
@@ -11,5 +11,6 @@
    CHARACTER(LEN=maxlen) :: modname, msg1, nm(2)
    INTEGER :: ixt, ipha, k, i, iq, iiso, izon, ieau, iqeau, iqpar
-   INTEGER, ALLOCATABLE :: ix(:)
+   INTEGER, ALLOCATABLE ::   ix(:)
+   REAL,    ALLOCATABLE :: tnat(:)               !--- OpenMP shared variable
    REAL    :: xtractot, xiiso, deltaD, q1, q2
    REAL, PARAMETER :: borne     = 1e19,  &
@@ -36,4 +37,5 @@
       iso_O17 = strIdx(isoName,'H2[17]O')
       iso_HTO = strIdx(isoName,'H[3]HO')
+      IF(getKey('tnat', tnat)) CALL abort_gcm(modname, 'missing isotopic parameter', 1)
 !$OMP END MASTER
 !$OMP BARRIER
Index: LMDZ6/trunk/libf/dyn3dmem/dynetat0_loc.F90
===================================================================
--- LMDZ6/trunk/libf/dyn3dmem/dynetat0_loc.F90	(revision 4319)
+++ LMDZ6/trunk/libf/dyn3dmem/dynetat0_loc.F90	(revision 4325)
@@ -7,9 +7,9 @@
 !-------------------------------------------------------------------------------
   USE parallel_lmdz
-  USE infotrac,    ONLY: nqtot, tracers, niso, iqIsoPha, tnat, alpha_ideal, iH2O
+  USE infotrac,    ONLY: nqtot, tracers, niso, iqIsoPha, iH2O, isoName
   USE strings_mod, ONLY: maxlen, msg, strStack, real2str, int2str
   USE netcdf,      ONLY: NF90_OPEN,  NF90_NOWRITE, NF90_INQUIRE_DIMENSION, NF90_INQ_VARID, &
                          NF90_CLOSE, NF90_GET_VAR, NF90_INQUIRE_VARIABLE,  NF90_NoErr
-  USE readTracFiles_mod, ONLY: new2oldH2O, newHNO3, oldHNO3
+  USE readTracFiles_mod, ONLY: new2oldH2O, newHNO3, oldHNO3, getKey
   USE control_mod, ONLY: planet_type
   USE assert_eq_m, ONLY: assert_eq
@@ -42,5 +42,5 @@
   INTEGER, PARAMETER :: length=100
   INTEGER :: iq, fID, vID, idecal, ierr, iqParent, iName, iZone, iPhase, ix
-  REAL    :: time, tab_cntrl(length)               !--- RUN PARAMS TABLE
+  REAL    :: time, tnat, alpha_ideal, tab_cntrl(length)    !--- RUN PARAMS TABLE
   REAL,             ALLOCATABLE :: vcov_glo(:,:),masse_glo(:,:),   ps_glo(:)
   REAL,             ALLOCATABLE :: ucov_glo(:,:),    q_glo(:,:), phis_glo(:)
@@ -179,6 +179,8 @@
       iqParent = tracers(iq)%iqParent
       IF(tracers(iq)%iso_iZone == 0) THEN
+         IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) &
+            CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1)
          CALL msg('Tracer <'//TRIM(var)//'> is missing => initialized with a simplified Rayleigh distillation law.', modname)
-         q(ijb_u:ije_u,:,iq)= q(ijb_u:ije_u,:,iqParent)*tnat(iName)*(q(ijb_u:ije_u,:,iqParent)/30.e-3)**(alpha_ideal(iName)-1.)
+         q(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqParent)*tnat*(q(ijb_u:ije_u,:,iqParent)/30.e-3)**(alpha_ideal-1.)
       ELSE
          CALL msg('Tracer <'//TRIM(var)//'> is missing => initialized to its parent isotope concentration.', modname)
Index: LMDZ6/trunk/libf/dyn3dmem/gcm.F90
===================================================================
--- LMDZ6/trunk/libf/dyn3dmem/gcm.F90	(revision 4319)
+++ LMDZ6/trunk/libf/dyn3dmem/gcm.F90	(revision 4325)
@@ -11,5 +11,5 @@
   USE mod_const_mpi, ONLY: init_const_mpi
   USE parallel_lmdz
-  USE infotrac, ONLY: nqtot, infotrac_init
+  USE infotrac, ONLY: nqtot, init_infotrac
 !#ifdef CPP_PHYS
 !  USE mod_interface_dyn_phys, ONLY: init_interface_dyn_phys
@@ -205,5 +205,5 @@
   !  Choix du nombre de traceurs et du schema pour l'advection
   !  dans fichier traceur.def, par default ou via INCA
-  call infotrac_init
+  call init_infotrac
 
   ! Allocation de la tableau q : champs advectes   
Index: LMDZ6/trunk/libf/dyn3dmem/iniacademic_loc.F90
===================================================================
--- LMDZ6/trunk/libf/dyn3dmem/iniacademic_loc.F90	(revision 4319)
+++ LMDZ6/trunk/libf/dyn3dmem/iniacademic_loc.F90	(revision 4325)
@@ -5,5 +5,5 @@
 
   USE filtreg_mod, ONLY: inifilr
-  USE infotrac,    ONLY: nqtot, niso, tnat, alpha_ideal, iqIsoPha, tracers
+  USE infotrac,    ONLY: nqtot, niso, iqIsoPha, tracers, getKey, isoName
   USE control_mod, ONLY: day_step,planet_type
   use exner_hyb_m, only: exner_hyb
@@ -77,5 +77,5 @@
   integer idum
 
-  REAL zdtvr
+  REAL zdtvr, tnat, alpha_ideal
   
   character(len=*),parameter :: modname="iniacademic"
@@ -290,6 +290,7 @@
               iqParent = tracers(iq)%iqParent
               IF(tracers(iq)%iso_iZone == 0) THEN
-                 q(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqParent)*tnat(iName) &
-                                     *(q(ijb_u:ije_u,:,iqParent)/30.e-3)**(alpha_ideal(iName)-1.)
+                 IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) &
+                    CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1)
+                 q(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqParent)*tnat*(q(ijb_u:ije_u,:,iqParent)/30.e-3)**(alpha_ideal-1.)
               ELSE
                  q(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqIsoPha(iName,iPhase))
Index: LMDZ6/trunk/libf/dyn3dmem/vlsplt_loc.F
===================================================================
--- LMDZ6/trunk/libf/dyn3dmem/vlsplt_loc.F	(revision 4319)
+++ LMDZ6/trunk/libf/dyn3dmem/vlsplt_loc.F	(revision 4325)
@@ -351,5 +351,5 @@
 c$OMP END DO NOWAIT
       enddo !do ifils=1,tracers(iq)%nqDescen
-      do ifils=1,tracers(iq)%nqChilds
+      do ifils=1,tracers(iq)%nqChildren
         iq2=tracers(iq)%iqDescen(ifils)
         call vlx_loc(Ratio,pente_max,masse,u_mq,ijb_x,ije_x,iq2)
@@ -726,5 +726,5 @@
 ! CRisi: appel récursif de l'advection sur les fils.
 ! Il faut faire ça avant d'avoir mis à jour q et masse
-      !write(*,*) 'vly 689: iq,nqChilds(iq)=',iq,tracers(iq)%nqChilds
+!     write(*,*)'vly 689: iq,nqChildren(iq)=',iq,tracers(iq)%nqChildren
 
       ijb=ij_begin-2*iip1
@@ -761,5 +761,5 @@
       enddo
 
-      do ifils=1,tracers(iq)%nqChilds
+      do ifils=1,tracers(iq)%nqChildren
         iq2=tracers(iq)%iqDescen(ifils)
         call vly_loc(Ratio,pente_max,masse,qbyv,iq2)
@@ -1148,5 +1148,5 @@
 ! CRisi: appel récursif de l'advection sur les fils.
 ! Il faut faire ça avant d'avoir mis à jour q et masse
-      !write(*,*) 'vlsplt 942: iq,nqChilds(iq)=',iq,tracers(iq)%nqChilds
+!     write(*,*)'vlsplt 942: iq,nqChildren(iq)=',iq,tracers(iq)%nqChildren
       do ifils=1,tracers(iq)%nqDescen
         iq2=tracers(iq)%iqDescen(ifils)
@@ -1169,5 +1169,5 @@
 c$OMP BARRIER
 
-      do ifils=1,tracers(iq)%nqChilds
+      do ifils=1,tracers(iq)%nqChildren
         iq2=tracers(iq)%iqDescen(ifils)
         call vlz_loc(Ratio,pente_max,masse,w,ijb_x,ije_x,iq2)
Index: LMDZ6/trunk/libf/dyn3dmem/vlspltqs_loc.F
===================================================================
--- LMDZ6/trunk/libf/dyn3dmem/vlspltqs_loc.F	(revision 4319)
+++ LMDZ6/trunk/libf/dyn3dmem/vlspltqs_loc.F	(revision 4325)
@@ -337,6 +337,6 @@
 ! CRisi: appel recursif de l'advection sur les fils.
 ! Il faut faire ca avant d'avoir mis a jour q et masse
-      !write(*,*) 'vlspltqs 336: iq,ijb_x,nqChilds(iq)=',
-!     &     iq,ijb_x,tracers(iq)%nqChilds
+      !write(*,*) 'vlspltqs 336: iq,ijb_x,nqChildren(iq)=',
+!     &     iq,ijb_x,tracers(iq)%nqChildren
 
       do ifils=1,tracers(iq)%nqDescen
@@ -356,5 +356,5 @@
 c$OMP END DO NOWAIT
       enddo
-      do ifils=1,tracers(iq)%nqChilds
+      do ifils=1,tracers(iq)%nqChildren
         iq2=tracers(iq)%iqDescen(ifils)
         !write(*,*) 'vlxqs 349: on appelle vlx pour iq2=',iq2
@@ -729,5 +729,6 @@
 ! CRisi: appel recursif de l'advection sur les fils.
 ! Il faut faire ca avant d'avoir mis a jour q et masse
-      !write(*,*) 'vlyqs 689: iq,nqChilds(iq)=',iq,tracers(iq)%nqChilds
+!     write(*,*)'vlyqs 689: iq,nqChildren(iq)=',iq,
+!    &             tracers(iq)%nqChildren
       
       ijb=ij_begin-2*iip1
@@ -767,5 +768,5 @@
 c$OMP END DO NOWAIT
       enddo
-      do ifils=1,tracers(iq)%nqChilds
+      do ifils=1,tracers(iq)%nqChildren
         iq2=tracers(iq)%iqDescen(ifils)
         !write(lunout,*) 'vly: appel recursiv vly iq2=',iq2
Index: LMDZ6/trunk/libf/dynphy_lonlat/phylmd/ce0l.F90
===================================================================
--- LMDZ6/trunk/libf/dynphy_lonlat/phylmd/ce0l.F90	(revision 4319)
+++ LMDZ6/trunk/libf/dynphy_lonlat/phylmd/ce0l.F90	(revision 4325)
@@ -23,5 +23,5 @@
   USE netcdf,         ONLY: NF90_OPEN, NF90_NOWRITE, NF90_CLOSE, NF90_NOERR,    &
          NF90_INQUIRE_DIMENSION, NF90_INQ_DIMID, NF90_INQ_VARID, NF90_GET_VAR
-  USE infotrac,       ONLY: type_trac, infotrac_init
+  USE infotrac,       ONLY: type_trac, init_infotrac
   USE dimphy,         ONLY: klon
   USE test_disvert_m, ONLY: test_disvert
@@ -132,5 +132,5 @@
 
 !--- Tracers initializations
-  CALL infotrac_init()
+  CALL init_infotrac()
 
   CALL inifilr()
Index: LMDZ6/trunk/libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90
===================================================================
--- LMDZ6/trunk/libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90	(revision 4319)
+++ LMDZ6/trunk/libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90	(revision 4325)
@@ -16,5 +16,5 @@
   USE mod_phys_lmdz_para, ONLY: klon_omp ! number of columns (on local omp grid)
   USE vertical_layers_mod, ONLY : init_vertical_layers
-  USE infotrac, ONLY: nbtr,nqCO2,tracers,isotopes,type_trac,types_trac,conv_flg,pbl_flg,nqtottr
+  USE infotrac, ONLY: nbtr, type_trac, types_trac
 #ifdef CPP_StratAer
   USE infotrac_phy, ONLY: nbtr_bin, nbtr_sulgas, id_OCS_strat, &
@@ -137,5 +137,5 @@
 
   ! Initialize tracer names, numbers, etc. for physics
-  CALL init_infotrac_phy(type_trac, tracers, isotopes, nqtottr, nqCO2, pbl_flg, conv_flg)
+  CALL init_infotrac_phy
 
   ! Initializations for Reprobus
Index: LMDZ6/trunk/libf/misc/readTracFiles_mod.f90
===================================================================
--- LMDZ6/trunk/libf/misc/readTracFiles_mod.f90	(revision 4319)
+++ LMDZ6/trunk/libf/misc/readTracFiles_mod.f90	(revision 4325)
@@ -1,7 +1,6 @@
 MODULE readTracFiles_mod
 
-  USE strings_mod,    ONLY: msg, find, get_in, str2int, dispTable, testFile, strReduce,  strFind, strStack, strHead, &
+  USE strings_mod,    ONLY: msg, find, get_in, str2int, dispTable, testFile, strReduce,  strFind, strStack, strHead,  &
        test, removeComment, cat, fmsg, maxlen, int2str, checkList, strParse, strReplace, strTail, strCount, strIdx, reduceExpr
-  USE trac_types_mod, ONLY: trac_type, isot_type, keys_type
 
   IMPLICIT NONE
@@ -9,37 +8,101 @@
   PRIVATE
 
-  PUBLIC :: maxlen                                                   !--- PARAMETER FOR CASUAL STRING LENGTH
-  PUBLIC :: trac_type, readTracersFiles, setGeneration, indexUpdate  !--- TRACERS  DESCRIPTION ASSOCIATED TOOLS
-  PUBLIC :: keys_type, getKey, fGetKey,  setDirectKeys, getKey_init  !--- TOOLS TO GET/SET KEYS FROM/TO tracers & isotopes
-
-  PUBLIC :: addPhase, getiPhase,  old_phases, phases_sep, nphases, & !--- FUNCTIONS RELATED TO THE PHASES
-            delPhase, getPhase, known_phases, phases_names           !--- + ASSOCIATED VARIABLES
-
-  PUBLIC :: oldH2OIso, newH2OIso, old2newH2O, new2oldH2O             !--- H2O ISOTOPES BACKWARD COMPATIBILITY (OLD traceur.def)
-  PUBLIC :: oldHNO3,   newHNO3                                       !--- HNO3 REPRO   BACKWARD COMPATIBILITY (OLD start.nc)
-
-  PUBLIC :: tran0, idxAncestor, ancestor                             !--- GENERATION 0 TRACER + TOOLS FOR GENERATIONS
+  PUBLIC :: maxlen                                              !--- PARAMETER FOR CASUAL STRING LENGTH
+  PUBLIC :: tracers                                             !--- TRACERS  DESCRIPTION DATABASE
+  PUBLIC :: trac_type, setGeneration, indexUpdate               !--- TRACERS  DESCRIPTION ASSOCIATED TOOLS
+  PUBLIC :: testTracersFiles, readTracersFiles                  !--- TRACERS FILES READING ROUTINES
+  PUBLIC :: getKey, fGetKey, setDirectKeys                      !--- TOOLS TO GET/SET KEYS FROM/TO  tracers & isotopes
+  PUBLIC :: getKeysDBase,    setKeysDBase                       !--- TOOLS TO GET/SET THE DATABASE (tracers & isotopes)
+
+  PUBLIC :: addPhase, getiPhase,  old_phases, phases_sep, &     !--- FUNCTIONS RELATED TO THE PHASES
+   nphases, delPhase, getPhase, known_phases, phases_names      !--- + ASSOCIATED VARIABLES
+
+  PUBLIC :: oldH2OIso, newH2OIso, old2newH2O, new2oldH2O        !--- H2O ISOTOPES BACKWARD COMPATIBILITY (OLD traceur.def)
+  PUBLIC :: oldHNO3,   newHNO3                                  !--- HNO3 REPRO   BACKWARD COMPATIBILITY (OLD start.nc)
+
+  PUBLIC :: tran0, idxAncestor, ancestor                        !--- GENERATION 0 TRACER + TOOLS FOR GENERATIONS
 
   !=== FOR ISOTOPES: GENERAL
-  PUBLIC :: isot_type, readIsotopesFile, initIsotopes                !--- ISOTOPES DESCRIPTION TYPE + READING ROUTINE
+  PUBLIC :: isot_type, readIsotopesFile, isoSelect              !--- ISOTOPES DESCRIPTION TYPE + READING ROUTINE
+  PUBLIC :: ixIso, nbIso                                        !--- INDEX OF SELECTED ISOTOPES CLASS + NUMBER OF CLASSES
+
+  !=== FOR ISOTOPES: H2O FAMILY ONLY
+  PUBLIC :: iH2O
+
+  !=== FOR ISOTOPES: DEPENDING ON THE SELECTED ISOTOPES CLASS
+  PUBLIC :: isotope, isoKeys                                    !--- SELECTED ISOTOPES DATABASE + ASSOCIATED KEYS
+  PUBLIC :: isoName, isoZone, isoPhas                           !--- ISOTOPES AND TAGGING ZONES NAMES AND PHASES
+  PUBLIC :: niso,    nzone,   nphas,   ntiso                    !---  " " NUMBERS + ISOTOPES AND TAGGING TRACERS NUMBERS
+  PUBLIC :: itZonIso                                            !--- Idx IN isoName(1:niso) = f(tagging idx, isotope idx)
+  PUBLIC :: iqIsoPha                                            !--- Idx IN qx(1:nqtot)     = f(isotope idx,   phase idx)
+  PUBLIC :: isoCheck                                            !--- FLAG TO RUN ISOTOPES CHECKING ROUTINES
 
   PUBLIC :: maxTableWidth
 !------------------------------------------------------------------------------------------------------------------------------
-  TYPE :: dataBase_type                                              !=== TYPE FOR TRACERS SECTION
-    CHARACTER(LEN=maxlen)  :: name                                   !--- Section name
-    TYPE(trac_type), ALLOCATABLE :: trac(:)                          !--- Tracers descriptors
+  TYPE :: keys_type                                        !=== TYPE FOR A SET OF KEYS ASSOCIATED TO AN ELEMENT
+    CHARACTER(LEN=maxlen)              :: name             !--- Tracer name
+    CHARACTER(LEN=maxlen), ALLOCATABLE :: key(:)           !--- Keys string list
+    CHARACTER(LEN=maxlen), ALLOCATABLE :: val(:)           !--- Corresponding values string list
+  END TYPE keys_type
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE :: trac_type                                        !=== TYPE FOR A SINGLE TRACER NAMED "name"
+    CHARACTER(LEN=maxlen) :: name        = ''              !--- Name of the tracer
+    CHARACTER(LEN=maxlen) :: gen0Name    = ''              !--- First generation ancestor name
+    CHARACTER(LEN=maxlen) :: parent      = ''              !--- Parent name
+    CHARACTER(LEN=maxlen) :: longName    = ''              !--- Long name (with advection scheme suffix)
+    CHARACTER(LEN=maxlen) :: type        = 'tracer'        !--- Type  (so far: 'tracer' / 'tag')
+    CHARACTER(LEN=maxlen) :: phase       = 'g'             !--- Phase ('g'as / 'l'iquid / 's'olid)
+    CHARACTER(LEN=maxlen) :: component   = ''              !--- Coma-separated list of components (Ex: lmdz,inca)
+    INTEGER               :: iadv        = 10              !--- Advection scheme used
+    INTEGER               :: iGeneration = -1              !--- Generation number (>=0)
+    LOGICAL               :: isAdvected  = .FALSE.         !--- "true" tracers: iadv > 0.   COUNT(isAdvected )=nqtrue
+    LOGICAL               :: isInPhysics = .TRUE.          !--- "true" tracers: in tr_seri. COUNT(isInPhysics)=nqtottr
+    INTEGER               :: iqParent    = 0               !--- Parent index
+    INTEGER,  ALLOCATABLE :: iqDescen(:)                   !--- Descendants index (in growing generation order)
+    INTEGER               :: nqDescen    = 0               !--- Number of descendants (all generations)
+    INTEGER               :: nqChildren  = 0               !--- Number of children  (first generation)
+    INTEGER               :: iso_iGroup  = 0               !--- Isotopes group index in isotopes(:)
+    INTEGER               :: iso_iName   = 0               !--- Isotope  name  index in isotopes(iso_iGroup)%trac(:)
+    INTEGER               :: iso_iZone   = 0               !--- Isotope  zone  index in isotopes(iso_iGroup)%zone(:)
+    INTEGER               :: iso_iPhase  = 0               !--- Isotope  phase index in isotopes(iso_iGroup)%phase
+    TYPE(keys_type)       :: keys                          !--- <key>=<val> pairs vector
+  END TYPE trac_type
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE :: isot_type                                        !=== TYPE FOR AN ISOTOPES FAMILY DESCENDING ON TRACER "parent"
+    CHARACTER(LEN=maxlen)              :: parent           !--- Isotopes family name (parent tracer name ; ex: H2O)
+    LOGICAL                            :: check=.FALSE.    !--- Triggering of the checking routines
+    TYPE(keys_type),       ALLOCATABLE :: keys(:)          !--- Isotopes keys/values pairs list     (length: niso)
+    CHARACTER(LEN=maxlen), ALLOCATABLE :: trac(:)          !--- Isotopes + tagging tracers list     (length: ntiso)
+    CHARACTER(LEN=maxlen), ALLOCATABLE :: zone(:)          !--- Geographic tagging zones names list (length: nzone)
+    CHARACTER(LEN=maxlen)              :: phase = 'g'      !--- Phases list: [g][l][s]              (length: nphas)
+    INTEGER                            :: niso  = 0        !--- Number of isotopes, excluding tagging tracers
+    INTEGER                            :: nzone = 0        !--- Number of geographic tagging zones
+    INTEGER                            :: ntiso = 0        !--- Number of isotopes, including tagging tracers
+    INTEGER                            :: nphas = 0        !--- Number phases
+    INTEGER,               ALLOCATABLE :: iqIsoPha(:,:)    !--- Idx in "tracers(1:nqtot)" = f(name(1:ntiso)),phas)
+                                                           !---        "iqIsoPha" former name: "iqiso"
+    INTEGER,               ALLOCATABLE :: itZonIso(:,:)    !--- Idx in "trac(1:ntiso)" = f(zone, name(1:niso))
+                                                           !---        "itZonIso" former name: "index_trac"
+  END TYPE isot_type
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE :: dataBase_type                                         !=== TYPE FOR TRACERS SECTION
+    CHARACTER(LEN=maxlen)  :: name                              !--- Section name
+    TYPE(trac_type), ALLOCATABLE :: trac(:)                     !--- Tracers descriptors
   END TYPE dataBase_type
 !------------------------------------------------------------------------------------------------------------------------------
   INTERFACE getKey
     MODULE PROCEDURE getKeyByName_s1,  getKeyByName_i1,  getKeyByName_r1, &
-                     getKeyByName_sm,  getKeyByName_im,  getKeyByName_rm
+                     getKeyByName_sm,  getKeyByName_im,  getKeyByName_rm, &
+                     getKeyByName_s1m, getKeyByName_i1m, getKeyByName_r1m
   END INTERFACE getKey
 !------------------------------------------------------------------------------------------------------------------------------
-  INTERFACE fGetKey;       MODULE PROCEDURE fgetKeyByIndex_s1, fgetKeyByName_s1; END INTERFACE fGetKey
+  INTERFACE    isoSelect;  MODULE PROCEDURE  isoSelectByIndex,  isoSelectByName; END INTERFACE isoSelect
   INTERFACE  old2newH2O;   MODULE PROCEDURE  old2newH2O_1,  old2newH2O_m;        END INTERFACE old2newH2O
   INTERFACE  new2oldH2O;   MODULE PROCEDURE  new2oldH2O_1,  new2oldH2O_m;        END INTERFACE new2oldH2O
+  INTERFACE fGetKey;       MODULE PROCEDURE fgetKeyIdx_s1, fgetKeyNam_s1, fgetKey_sm;        END INTERFACE fGetKey
   INTERFACE tracersSubset; MODULE PROCEDURE trSubset_Indx, trSubset_Name, trSubset_gen0Name; END INTERFACE tracersSubset
-  INTERFACE idxAncestor;   MODULE PROCEDURE idxAncestor_1, idxAncestor_m;        END INTERFACE idxAncestor
-  INTERFACE    ancestor;   MODULE PROCEDURE    ancestor_1,    ancestor_m;        END INTERFACE    ancestor
+  INTERFACE idxAncestor;   MODULE PROCEDURE idxAncestor_1, idxAncestor_m, idxAncestor_mt;    END INTERFACE idxAncestor
+  INTERFACE    ancestor;   MODULE PROCEDURE    ancestor_1,    ancestor_m,    ancestor_mt;    END INTERFACE    ancestor
+  INTERFACE      addKey;   MODULE PROCEDURE      addKey_1,      addKey_m,     addKey_mm;     END INTERFACE addKey
   INTERFACE    addPhase;   MODULE PROCEDURE   addPhase_s1,   addPhase_sm,   addPhase_i1,   addPhase_im; END INTERFACE addPhase
 !------------------------------------------------------------------------------------------------------------------------------
@@ -54,8 +117,9 @@
   INTEGER, PARAMETER :: nphases = LEN_TRIM(known_phases)        !--- Number of phases
   CHARACTER(LEN=maxlen), SAVE      :: phases_names(nphases) &   !--- Known phases names
-                                = ['gaseous', 'liquid ', 'solid  ', 'cloud  ']
+                                             = ['gaseous', 'liquid ', 'solid  ', 'cloud  ']
   CHARACTER(LEN=1), SAVE :: phases_sep  =  '_'                  !--- Phase separator
   LOGICAL,          SAVE :: tracs_merge = .TRUE.                !--- Merge/stack tracers lists
   LOGICAL,          SAVE :: lSortByGen  = .TRUE.                !--- Sort by growing generation
+  CHARACTER(LEN=maxlen), SAVE :: isoFile = 'isotopes_params.def'!--- Name of the isotopes parameters file
 
   !--- CORRESPONDANCE BETWEEN OLD AND NEW WATER NAMES
@@ -67,9 +131,23 @@
   CHARACTER(LEN=maxlen), SAVE ::   newHNO3(2) = ['HNO3   ', 'HNO3tot']
 
-  !=== LOCAL COPIES OF THE TRACERS AND ISOTOPES DESCRIPTORS, USED BY getKey (INITIALIZED IN getKey_init)
+  !=== TRACERS AND ISOTOPES DESCRIPTORS, USED BY getKey
   TYPE(trac_type), ALLOCATABLE, TARGET, SAVE ::  tracers(:)
   TYPE(isot_type), ALLOCATABLE, TARGET, SAVE :: isotopes(:)
 
-  INTEGER,    PARAMETER :: maxTableWidth = 192                       !--- Maximum width of a table displayed with "dispTable"
+  !=== ALIASES OF VARIABLES FROM SELECTED ISOTOPES FAMILY EMBEDDED IN "isotope" (isotopes(ixIso))
+  TYPE(isot_type),         SAVE, POINTER :: isotope             !--- CURRENTLY SELECTED ISOTOPES FAMILY DESCRIPTOR
+  INTEGER,                 SAVE          :: ixIso, iH2O         !--- Index of the selected isotopes family and H2O family
+  INTEGER,                 SAVE          :: nbIso               !--- Number of isotopes classes
+  LOGICAL,                 SAVE          :: isoCheck            !--- Flag to trigger the checking routines
+  TYPE(keys_type),         SAVE, POINTER :: isoKeys(:)          !--- ONE SET OF KEYS FOR EACH ISOTOPE (LISTED IN isoName)
+  CHARACTER(LEN=maxlen),   SAVE, POINTER :: isoName(:),   &     !--- ISOTOPES NAMES FOR THE CURRENTLY SELECTED FAMILY
+                                            isoZone(:),   &     !--- TAGGING ZONES  FOR THE CURRENTLY SELECTED FAMILY
+                                            isoPhas             !--- USED PHASES    FOR THE CURRENTLY SELECTED FAMILY
+  INTEGER,                 SAVE          ::  niso, nzone, &     !--- NUMBER OF ISOTOPES, TAGGING ZONES AND PHASES
+                                            nphas, ntiso        !--- NUMBER OF PHASES AND ISOTOPES + ISOTOPIC TAGGING TRACERS
+  INTEGER,                 SAVE, POINTER ::itZonIso(:,:), &     !--- INDEX IN "isoTrac" AS f(tagging zone idx,  isotope idx)
+                                           iqIsoPha(:,:)        !--- INDEX IN "qx"      AS f(isotopic tracer idx, phase idx)
+
+  INTEGER,    PARAMETER :: maxTableWidth = 192                  !--- Maximum width of a table displayed with "dispTable"
   CHARACTER(LEN=maxlen) :: modname
 
@@ -100,15 +178,15 @@
 !     * If you need to convert a %key/%val pair into a direct-access key, add the corresponding line in "setDirectKeys".
 !==============================================================================================================================
-LOGICAL FUNCTION readTracersFiles(type_trac, fTyp, tracs, lRepr) RESULT(lerr)
-!------------------------------------------------------------------------------------------------------------------------------
-  CHARACTER(LEN=*),             INTENT(IN)  :: type_trac              !--- List of components used
-  INTEGER,         OPTIONAL,    INTENT(OUT) :: fTyp                   !--- Type of input file found
-  TYPE(trac_type), ALLOCATABLE, INTENT(OUT) :: tracs(:)
-  LOGICAL,         OPTIONAL,    INTENT(IN)  :: lRepr
+LOGICAL FUNCTION readTracersFiles(type_trac, fTyp, lRepr) RESULT(lerr)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=*),  INTENT(IN)  :: type_trac                        !--- List of components used
+  INTEGER, OPTIONAL, INTENT(OUT) :: fTyp                             !--- Type of input file found
+  LOGICAL, OPTIONAL, INTENT(IN)  :: lRepr                            !--- Activate the HNNO3 exceptions for REPROBUS
   CHARACTER(LEN=maxlen),  ALLOCATABLE :: s(:), sections(:), trac_files(:)
   CHARACTER(LEN=maxlen) :: str, fname, mesg, tname, pname, cname
-  INTEGER               :: is, nsec, ierr, it, ntrac, ns, ip, ix, fType
+  INTEGER               :: nsec, ierr, it, ntrac, ns, ip, ix, fType
   LOGICAL, ALLOCATABLE  :: ll(:), lGen3(:)
   LOGICAL :: lRep
+  TYPE(keys_type), POINTER :: k
 !------------------------------------------------------------------------------------------------------------------------------
   lerr = .FALSE.
@@ -117,29 +195,8 @@
   lRep=0; IF(PRESENT(lRepr)) lRep = lRepr
 
-  !--- Required sections + corresponding files names (new style single section case)
-  IF(test(strParse(type_trac, '|', sections), lerr)) RETURN          !--- Parse "type_trac" list
-
-  nsec = SIZE(sections, DIM=1)
-  ALLOCATE(trac_files(nsec)); DO is=1, nsec; trac_files(is) = 'tracer_'//TRIM(sections(is))//'.def'; END DO
-
-  !--- LOOK AT AVAILABLE FILES
-  ll = .NOT.testFile(trac_files)
-  fType = 0
-  IF(.NOT.testFile('traceur.def') .AND. nsec==1) fType = 1           !--- OLD STYLE FILE
-  IF(.NOT.testFile('tracer.def'))                fType = 2           !--- NEW STYLE ; SINGLE  FILE, SEVERAL SECTIONS
-  IF(ALL(ll))                                    fType = 3           !--- NEW STYLE ; SEVERAL FILES, SINGLE SECTION USED
+  !--- Required sections + corresponding files names (new style single section case) for tests
+  IF(test(testTracersFiles(modname, type_trac, fType, trac_files, sections), lerr)) RETURN
   IF(PRESENT(fTyp)) fTyp = fType
-  IF(ANY(ll) .AND. fType/=3) THEN                                    !--- MISSING FILES
-    IF(test(checkList(trac_files, .NOT.ll, 'Failed reading tracers description', 'files', 'missing'), lerr)) RETURN
-  END IF
-
-  !--- CHECK WHETHER type_trac AND FILE TYPE ARE COMPATIBLE
-  IF(test(fmsg('No multiple sections for the old format "traceur.def"', ll = SIZE(sections)>1 .AND. fType==1), lerr)) RETURN
-
-  !--- TELLS WHAT WAS IS ABOUT TO BE USED
-  IF (fmsg('No adequate tracers description file(s) found ; default values will be used',          modname, fType==0)) RETURN
-  CALL msg('Trying to read old-style tracers description file "traceur.def"',                      modname, fType==1)
-  CALL msg('Trying to read the new style multi-sections tracers description file "tracer.def"',    modname, fType==2)
-  CALL msg('Trying to read the new style single section tracers description files "tracer_*.def"', modname, fType==3)
+  nsec = SIZE(sections)
 
   !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -156,5 +213,6 @@
 
       !--- READ THE REMAINING LINES: <hadv> <vadv> <tracer> [<transporting fluid>]
-      ALLOCATE(tracs(ntrac))
+      IF(ALLOCATED(tracers)) DEALLOCATE(tracers)
+      ALLOCATE(tracers(ntrac))
       DO it=1,ntrac                                                  !=== READ RAW DATA: loop on the line/tracer number
         READ(90,'(a)',IOSTAT=ierr) str
@@ -164,4 +222,5 @@
         CALL msg('This file is for air tracers only',           modname, ns == 3 .AND. it == 1)
         CALL msg('This files specifies the transporting fluid', modname, ns == 4 .AND. it == 1)
+        k => tracers(it)%keys
 
         !=== NAME OF THE TRACER
@@ -169,11 +228,13 @@
         ix = strIdx(oldHNO3, s(3))
         IF(ix /= 0 .AND. lRep) tname = newHNO3(ix)                   !--- Exception for HNO3 (REPROBUS ONLY)
-        tracs(it)%name = tname                                       !--- Set %name
-        tracs(it)%keys%name = tname                                  !--- Copy tracers names in keys components
+        tracers(it)%name = tname                                     !--- Set %name
+        CALL addKey('name', tname, k)                                !--- Set the name of the tracer
+        tracers(it)%keys%name = tname                                !--- Copy tracers names in keys components
 
         !=== NAME OF THE COMPONENT
         cname = type_trac                                            !--- Name of the model component
         IF(ANY([(addPhase('H2O', ip), ip = 1, nphases)] == tname)) cname = 'lmdz'
-        tracs(it)%component = cname                                  !--- Set %component
+        tracers(it)%component = cname                                !--- Set %component
+        CALL addKey('component', cname, k)                           !--- Set the name of the model component
 
         !=== NAME OF THE PARENT
@@ -184,18 +245,20 @@
           IF(ix /= 0 .AND. lRep) pname = newHNO3(ix)                 !--- Exception for HNO3 (REPROBUS ONLY)
         END IF
-        tracs(it)%parent = pname                                     !--- Set %parent
+        tracers(it)%parent = pname                                   !--- Set %parent
+        CALL addKey('parent', pname, k)
 
         !=== PHASE AND ADVECTION SCHEMES NUMBERS
-        tracs(it)%phase = known_phases(ip:ip)                        !--- Set %phase:  tracer phase (default: "g"azeous)
-        tracs(it)%keys%key = ['hadv', 'vadv']                        !--- Set %keys%key
-        tracs(it)%keys%val = s(1:2)                                  !--- Set %keys%val
+        tracers(it)%phase = known_phases(ip:ip)                      !--- Set %phase:  tracer phase (default: "g"azeous)
+        CALL addKey('phase', known_phases(ip:ip), k)                 !--- Set the phase  of the tracer (default: "g"azeous)
+        CALL addKey('hadv', s(1),  k)                                !--- Set the horizontal advection schemes number
+        CALL addKey('vadv', s(2),  k)                                !--- Set the vertical   advection schemes number
       END DO
       CLOSE(90)
-      CALL setGeneration(tracs)                                      !--- Set %iGeneration and %gen0Name
-      WHERE(tracs%iGeneration == 2) tracs%type = 'tag'               !--- Set %type:        'tracer' or 'tag'
-      IF(test(checkTracers(tracs, fname, fname), lerr)) RETURN       !--- Detect orphans and check phases
-      IF(test(checkUnique (tracs, fname, fname), lerr)) RETURN       !--- Detect repeated tracers
-      CALL sortTracers  (tracs)                                      !--- Sort the tracers
-      tracs(:)%keys%name = tracs(:)%name                             !--- Copy tracers names in keys components
+      IF(test(setGeneration(tracers), lerr)) RETURN                  !--- Set %iGeneration and %gen0Name
+      WHERE(tracers%iGeneration == 2) tracers(:)%type = 'tag'        !--- Set %type:        'tracer' or 'tag'
+      CALL addKey('type', tracers(:)%type, tracers(:)%keys)          !--- Set the type of tracers
+      IF(test(checkTracers(tracers, fname, fname), lerr)) RETURN     !--- Detect orphans and check phases
+      IF(test(checkUnique (tracers, fname, fname), lerr)) RETURN     !--- Detect repeated tracers
+      CALL sortTracers    (tracers)                                  !--- Sort the tracers
     !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     CASE(2); IF(test(feedDBase(["tracer.def"], [type_trac], modname), lerr)) RETURN !=== SINGLE   FILE, MULTIPLE SECTIONS
@@ -209,14 +272,51 @@
 
   IF(nsec  == 1) THEN; 
-    tracs = dBase(1)%trac
+    tracers = dBase(1)%trac
   ELSE IF(tracs_merge) THEN
     CALL msg('The multiple required sections will be MERGED.',    modname)
-    IF(test(mergeTracers(dBase, tracs), lerr)) RETURN
+    IF(test(mergeTracers(dBase, tracers), lerr)) RETURN
   ELSE
     CALL msg('The multiple required sections will be CUMULATED.', modname)
-    IF(test(cumulTracers(dBase, tracs), lerr)) RETURN
+    IF(test(cumulTracers(dBase, tracers), lerr)) RETURN
   END IF
-  CALL setDirectKeys(tracs)                                          !--- Set %iqParent, %iqDescen, %nqDescen, %nqChilds
+  CALL setDirectKeys(tracers)                                        !--- Set %iqParent, %iqDescen, %nqDescen, %nqChildren
 END FUNCTION readTracersFiles
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+LOGICAL FUNCTION testTracersFiles(modname, type_trac, fType, tracf, sects) RESULT(lerr)
+  CHARACTER(LEN=*),                             INTENT(IN)  :: modname, type_trac
+  INTEGER,                                      INTENT(OUT) :: fType
+  CHARACTER(LEN=maxlen), ALLOCATABLE, OPTIONAL, INTENT(OUT) :: tracf(:), sects(:)
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: trac_files(:), sections(:)
+  LOGICAL, ALLOCATABLE :: ll(:)
+  INTEGER :: is, nsec
+  lerr = .FALSE.
+
+  !--- PARSE "type_trac" LIST AND DETERMINE THE TRACERS FILES NAMES (FOR CASE 3: MULTIPLE FILES, SINNGLE SECTION PER FILE)
+  IF(test(strParse(type_trac, '|', sections,  n=nsec), lerr)) RETURN !--- Parse "type_trac" list
+  IF(PRESENT(sects)) sects = sections
+  ALLOCATE(trac_files(nsec)); DO is=1, nsec; trac_files(is) = 'tracer_'//TRIM(sections(is))//'.def'; END DO
+  IF(PRESENT(tracf)) tracf = trac_files
+
+  nsec = SIZE(trac_files, DIM=1)
+  ll = .NOT.testFile(trac_files)
+  fType = 0
+  IF(.NOT.testFile('traceur.def') .AND. nsec==1) fType = 1           !--- OLD STYLE FILE
+  IF(.NOT.testFile('tracer.def'))                fType = 2           !--- NEW STYLE ; SINGLE  FILE, SEVERAL SECTIONS
+  IF(ALL(ll))                                    fType = 3           !--- NEW STYLE ; SEVERAL FILES, SINGLE SECTION USED
+  IF(ANY(ll) .AND. fType/=3) THEN                                    !--- MISSING FILES
+    IF(test(checkList(trac_files, .NOT.ll, 'Failed reading tracers description', 'files', 'missing'), lerr)) RETURN
+  END IF
+
+  !--- CHECK WHETHER type_trac AND FILE TYPE ARE COMPATIBLE
+  IF(test(fmsg('No multiple sections for the old format "traceur.def"', ll = nsec>1 .AND. fType==1), lerr)) RETURN
+
+  !--- TELLS WHAT WAS IS ABOUT TO BE USED
+  CALL msg('Trying to read old-style tracers description file "traceur.def"',                      modname, fType==1)
+  CALL msg('Trying to read the new style multi-sections tracers description file "tracer.def"',    modname, fType==2)
+  CALL msg('Trying to read the new style single section tracers description files "tracer_*.def"', modname, fType==3)
+END FUNCTION testTracersFiles
 !==============================================================================================================================
 
@@ -253,5 +353,5 @@
     lerr = ANY([(dispTraSection('RAW CONTENT OF SECTION "'//TRIM(snm)//'"', snm, modname), idb=1, SIZE(dBase))])
     IF(test(expandSection(dBase(idb)%trac, snm, fnm), lerr)) RETURN  !--- EXPAND NAMES ;  set %parent, %type, %component
-    CALL setGeneration   (dBase(idb)%trac)                           !---                 set %iGeneration,   %genOName
+    IF(test(setGeneration(dBase(idb)%trac),           lerr)) RETURN  !---                 set %iGeneration,   %genOName
     IF(test(checkTracers (dBase(idb)%trac, snm, fnm), lerr)) RETURN  !--- CHECK ORPHANS AND PHASES
     IF(test(checkUnique  (dBase(idb)%trac, snm, fnm), lerr)) RETURN  !--- CHECK TRACERS UNIQUENESS
@@ -359,6 +459,6 @@
   ky => t(jd)%keys
   DO k = 1, SIZE(ky%key)                                             !--- Loop on the keys of the tracer named "defName"
-!   CALL addKey_m(ky%key(k), ky%val(k), t(:)%keys)                   !--- Add key to all the tracers (no overwriting)
-    DO it = 1, SIZE(t); CALL addKey_1(ky%key(k), ky%val(k), t(it)%keys); END DO
+!   CALL addKey_m(ky%key(k), ky%val(k), t(:)%keys, .FALSE.)           !--- Add key to all the tracers (no overwriting)
+    DO it = 1, SIZE(t); CALL addKey_1(ky%key(k), ky%val(k), t(it)%keys, .FALSE.); END DO
   END DO
   tt = [t(1:jd-1),t(jd+1:SIZE(t))]; CALL MOVE_ALLOC(FROM=tt, TO=t)   !--- Remove the virtual tracer named "defName"
@@ -409,6 +509,6 @@
   TYPE(trac_type),       ALLOCATABLE :: ttr(:)
   CHARACTER(LEN=maxlen), ALLOCATABLE :: ta(:), pa(:)
-  CHARACTER(LEN=maxlen) :: msg1, modname
-  INTEGER :: it, nt, iq, nq, itr, ntr, ipr, npr, i
+  CHARACTER(LEN=maxlen) :: msg1, modname, tname, cname , pname
+  INTEGER :: it, nt, iq, nq, jq, itr, ntr, ipr, npr, i
   LOGICAL :: ll
   modname = 'expandSection'
@@ -423,4 +523,5 @@
     tr(it)%type      = fgetKey(it, 'type'  , tr(:)%keys, 'tracer')
     tr(it)%component = sname
+    CALL addKey('component', sname, tr(:)%keys)
 
     !--- Determine the number of tracers and parents ; coherence checking
@@ -438,5 +539,4 @@
   END DO
   !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-  CALL delKey(['parent','type  '], tr)
 
   ALLOCATE(ttr(nq))
@@ -449,16 +549,14 @@
     DO ipr=1,npr                                                     !--- Loop on parents list elts
       DO itr=1,ntr                                                   !--- Loop on tracers list elts
-        i = iq+itr-1+(ipr-1)*ntr
-        ttr(i)%name   = TRIM(ta(itr))
-        ttr(i)%parent = TRIM(pa(ipr))
-        ttr(i)%keys%name = ta(itr)
-        ttr(i)%keys%key  = tr(it)%keys%key
-        ttr(i)%keys%val  = tr(it)%keys%val
-!        ttr(i)%keys   = keys_type(ta(itr), tr(it)%keys%key, tr(it)%keys%val)
+        ttr(iq)%keys%key  = tr(it)%keys%key
+        ttr(iq)%keys%val  = tr(it)%keys%val
+        ttr(iq)%keys%name = ta(itr)
+        ttr(iq)%name      = TRIM(ta(itr));    CALL addKey('name',      ta(itr),          ttr(iq)%keys)
+        ttr(iq)%parent    = TRIM(pa(ipr));    CALL addKey('parent',    pa(ipr),          ttr(iq)%keys)
+        ttr(iq)%type      = tr(it)%type;      CALL addKey('type',      tr(it)%type,      ttr(iq)%keys)
+        ttr(iq)%component = tr(it)%component; CALL addKey('component', tr(it)%component, ttr(iq)%keys)
+        iq = iq+1
       END DO
     END DO
-    ttr(iq:iq+ntr*npr-1)%type      = tr(it)%type                     !--- Duplicating type
-    ttr(iq:iq+ntr*npr-1)%component = tr(it)%component                !--- Duplicating type
-    iq = iq + ntr*npr
   !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   END DO
@@ -471,5 +569,5 @@
 
 !==============================================================================================================================
-SUBROUTINE setGeneration(tr)
+LOGICAL FUNCTION setGeneration(tr) RESULT(lerr)
 !------------------------------------------------------------------------------------------------------------------------------
 ! Purpose: Determine, for each tracer of "tr(:)":
@@ -479,23 +577,27 @@
   TYPE(trac_type),     INTENT(INOUT) :: tr(:)                        !--- Tracer derived type vector
   INTEGER                            :: iq, nq, ig
-  LOGICAL,               ALLOCATABLE :: lg(:)
-  CHARACTER(LEN=maxlen), ALLOCATABLE :: prn(:)
-!------------------------------------------------------------------------------------------------------------------------------
-  tr(:)%iGeneration = -1                                             !--- error if -1
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: parent(:), prn(:)
+  CHARACTER(LEN=maxlen) :: gen0(SIZE(tr))
+  INTEGER               :: iGen(SIZE(tr))
+  LOGICAL               ::   lg(SIZE(tr))
+!------------------------------------------------------------------------------------------------------------------------------
+  iGen(:) = -1                                                       !--- error if -1
   nq = SIZE(tr, DIM=1)                                               !--- Number of tracers lines
-  lg = tr(:)%parent == tran0                                         !--- Flag for generation 0 tracers
-  WHERE(lg) tr(:)%iGeneration = 0                                    !--- Generation 0 tracers
+  IF(test(fmsg('missing "parent" attribute', 'setGeneration', getKey('parent', parent, ky=tr(:)%keys)), lerr)) RETURN
+  WHERE(parent == tran0) iGen(:) = 0
 
   !=== Determine generation for each tracer
   ig=-1; prn = [tran0]
   DO                                                                 !--- Update current generation flag
-    IF(ig/=-1) prn = PACK( tr(:)%name, MASK=tr(:)%iGeneration == ig)
-    lg(:) = [(ANY(prn(:) == tr(iq)%parent), iq=1, nq)]               !--- Current generation tracers flag
+    IF(ig/=-1) prn = PACK( tr(:)%name, MASK = iGen == ig)
+    lg(:) = [(ANY(prn(:) == parent(iq)), iq=1, nq)]                  !--- Current generation tracers flag
     IF( ALL( .NOT. lg ) ) EXIT                                       !--- Empty current generation
-    ig = ig+1; WHERE(lg) tr(:)%iGeneration = ig
-  END DO
-  tr(:)%gen0Name = ancestor(tr)                                      !--- First generation ancestor name
-
-END SUBROUTINE setGeneration
+    ig = ig+1; WHERE(lg) iGen(:) = ig
+  END DO
+  tr%iGeneration = iGen; CALL addKey_mm('iGeneration', int2str(iGen(:)), tr(:)%keys)
+  CALL ancestor(tr, gen0)                                            !--- First generation ancestor name
+  tr%gen0Name    = gen0; CALL addKey_mm('gen0Name',    gen0,             tr(:)%keys)
+
+END FUNCTION setGeneration
 !==============================================================================================================================
 
@@ -581,5 +683,6 @@
   TYPE(trac_type), ALLOCATABLE :: ttr(:)
   INTEGER,   ALLOCATABLE ::  i0(:)
-  CHARACTER(LEN=maxlen)  :: nam, pha, trn
+  CHARACTER(LEN=maxlen)  :: nam, pha, tname
+  CHARACTER(LEN=maxlen), allocatable :: ph(:)
   CHARACTER(LEN=1) :: p
   INTEGER :: ip, np, iq, jq, nq, it, nt, nc, i, n
@@ -590,7 +693,7 @@
   DO iq = 1, nq                                                      !--- GET THE NUMBER OF TRACERS
     IF(tr(iq)%iGeneration /= 0) CYCLE                                !--- Only deal with generation 0 tracers
-    nc = COUNT(tr(:)%gen0Name==tr(iq)%name .AND. tr%iGeneration/=0)  !--- Number of childs of tr(iq)
-    tr(iq)%phase = fgetKey(iq, 'phases', tr(:)%keys)                 !--- Phases list      of tr(iq)
-    np = LEN_TRIM(tr(iq)%phase)                                      !--- Number of phases of tr(iq)
+    nc = COUNT(tr(:)%gen0Name==tr(iq)%name .AND. tr%iGeneration/=0)  !--- Number of children of tr(iq)
+    tr(iq)%phase = fgetKey(iq, 'phases', tr(:)%keys)                 !--- Phases list        of tr(iq)
+    np = LEN_TRIM(tr(iq)%phase)                                      !--- Number of phases   of tr(iq)
     nt = nt + (1+nc) * np                                            !--- Number of tracers after expansion
   END DO
@@ -609,15 +712,19 @@
       DO ip = 1, LEN_TRIM(pha)                                       !=== LOOP ON PHASES LISTS
         p = pha(ip:ip)
-        trn = TRIM(tr(iq)%name); nam = trn                           !--- Tracer name (regular case)
+        tname = TRIM(tr(iq)%name); nam = tname                       !--- Tracer name (regular case)
         IF(lTag) nam = TRIM(tr(iq)%parent)                           !--- Parent name (tagging case)
         IF(lExt) nam = addPhase(nam, p )                             !--- Phase extension needed
-        IF(lTag) nam = TRIM(nam)//'_'//TRIM(trn)                     !--- <parent>_<name> for tags
+        IF(lTag) nam = TRIM(nam)//'_'//TRIM(tname)                   !--- <parent>_<name> for tags
         ttr(it) = tr(iq)                                             !--- Same <key>=<val> pairs
         ttr(it)%name      = TRIM(nam)                                !--- Name with possibly phase suffix
         ttr(it)%keys%name = TRIM(nam)                                !--- Name inside the keys decriptor
         ttr(it)%phase     = p                                        !--- Single phase entry
+        CALL addKey('name', nam, ttr(it)%keys)
+        CALL addKey('phase', p,  ttr(it)%keys)
         IF(lExt .AND. tr(iq)%iGeneration>0) THEN
-          ttr(it)%parent   = addPhase(ttr(it)%parent,   p)
-          ttr(it)%gen0Name = addPhase(ttr(it)%gen0Name, p)
+          ttr(it)%parent   = addPhase(tr(iq)%parent,   p)
+          ttr(it)%gen0Name = addPhase(tr(iq)%gen0Name, p)
+          CALL addKey('parent',   ttr(it)%parent,   ttr(it)%keys)
+          CALL addKey('gen0Name', ttr(it)%gen0Name, ttr(it)%keys)
         END IF
         it = it+1
@@ -638,5 +745,5 @@
 !  * Put water at the beginning of the vector, in the "known_phases" order.
 !  * lGrowGen == T: in ascending generations numbers.
-!  * lGrowGen == F: tracer + its childs sorted by growing generation, one after the other.
+!  * lGrowGen == F: tracer + its children sorted by growing generation, one after the other.
 !   TO BE ADDED IF NECESSARY: HIGHER MOMENTS AT THE END
 !------------------------------------------------------------------------------------------------------------------------------
@@ -671,5 +778,5 @@
       ix(iq) = jq                                                    !--- Generation 0 ancestor index first
       iq = iq + 1                                                    !--- Next "iq" for next generations tracers
-      iy = strFind(tr(:)%gen0Name, TRIM(tr(jq)%name))                !--- Indexes of "tr(jq)" childs in "tr(:)"
+      iy = strFind(tr(:)%gen0Name, TRIM(tr(jq)%name))                !--- Indexes of "tr(jq)" children in "tr(:)"
       ng = MAXVAL(tr(iy)%iGeneration, MASK=.TRUE., DIM=1)            !--- Number of generations of the "tr(jq)" family
       DO ig = 1, ng                                                  !--- Loop   on generations of the "tr(jq)" family
@@ -683,4 +790,5 @@
 END SUBROUTINE sortTracers
 !==============================================================================================================================
+
 
 !==============================================================================================================================
@@ -793,5 +901,5 @@
   TYPE(trac_type), INTENT(INOUT) :: tr(:)
 
-  !--- Update %iqParent, %iqDescen, %nqDescen, %nqChilds
+  !--- Update %iqParent, %iqDescen, %nqDescen, %nqChildren
   CALL indexUpdate(tr)
 
@@ -808,5 +916,5 @@
   INTEGER :: idb, iq, nq
   INTEGER, ALLOCATABLE :: hadv(:), vadv(:)
-  CHARACTER(LEN=maxlen), ALLOCATABLE :: phas(:)
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: phas(:), prnt(:)
   TYPE(trac_type), POINTER :: tm(:)
   lerr = .FALSE.
@@ -816,14 +924,18 @@
   !--- BEWARE ! Can't use the "getKeyByName" functions yet.
   !             Names must first include the phases for tracers defined on multiple lines.
-  hadv = str2int([(fgetKey(iq, 'hadv',  tm(:)%keys, '10'), iq=1, nq)])
-  vadv = str2int([(fgetKey(iq, 'vadv',  tm(:)%keys, '10'), iq=1, nq)])
-  phas =         [(fgetKey(iq, 'phases',tm(:)%keys, 'g' ), iq=1, nq)]
+  hadv = str2int(fgetKey('hadv',  tm(:)%keys, '10'))
+  vadv = str2int(fgetKey('vadv',  tm(:)%keys, '10'))
+  prnt =         fgetKey('parent',tm(:)%keys,  '' )
+  IF(getKey('phases', phas, ky=tm(:)%keys)) phas = fGetKey('phase', tm(:)%keys, 'g')
   CALL msg(TRIM(message)//':', modname)
-  IF(ALL(tm(:)%parent == '')) THEN
-    IF(test(dispTable('iiiss',   ['iq    ','hadv  ','vadv  ','name  ','phase '], cat(tm%name, phas), &
+  IF(ALL(prnt == 'air')) THEN
+    IF(test(dispTable('iiiss',   ['iq    ','hadv  ','vadv  ','name  ','phase '],                   cat(tm%name,       phas),  &
+                 cat([(iq, iq=1, nq)], hadv, vadv),                 nColMax=maxTableWidth, nHead=2, sub=modname), lerr)) RETURN
+  ELSE IF(ALL(tm%iGeneration == -1)) THEN
+    IF(test(dispTable('iiisss', ['iq    ','hadv  ','vadv  ','name  ','parent','phase '],           cat(tm%name, prnt, phas),  &
                  cat([(iq, iq=1, nq)], hadv, vadv),                 nColMax=maxTableWidth, nHead=2, sub=modname), lerr)) RETURN
   ELSE
-    IF(test(dispTable('iiissis', ['iq    ','hadv  ','vadv  ','name  ','parent','igen  ','phase '], cat(tm%name, tm%parent, &
-      tm%phase), cat([(iq, iq=1, nq)], hadv, vadv, tm%iGeneration), nColMax=maxTableWidth, nHead=2, sub=modname), lerr)) RETURN
+    IF(test(dispTable('iiissis', ['iq    ','hadv  ','vadv  ','name  ','parent','igen  ','phase '], cat(tm%name, prnt, phas),  &
+                 cat([(iq, iq=1, nq)], hadv, vadv, tm%iGeneration), nColMax=maxTableWidth, nHead=2, sub=modname), lerr)) RETURN
   END IF
 END FUNCTION dispTraSection
@@ -884,7 +996,7 @@
 SUBROUTINE indexUpdate(tr)
   TYPE(trac_type), INTENT(INOUT) :: tr(:)
-  INTEGER :: iq, ig, ng, igen, ngen
-  INTEGER, ALLOCATABLE :: ix(:)
+  INTEGER :: iq, ig, ng, igen, ngen, ix(SIZE(tr))
   tr(:)%iqParent = strIdx( tr(:)%name, tr(:)%parent )                !--- Parent index
+  CALL addKey('iqParent', int2str(tr%iqParent), tr(:)%keys)
   ngen = MAXVAL(tr(:)%iGeneration, MASK=.TRUE.)
   DO iq = 1, SIZE(tr)
@@ -892,11 +1004,16 @@
     IF(ALLOCATED(tr(iq)%iqDescen)) DEALLOCATE(tr(iq)%iqDescen)
     ALLOCATE(tr(iq)%iqDescen(0))
-    ix = idxAncestor(tr, igen=ig)                                    !--- Ancestor of generation "ng" for each tr
+    CALL idxAncestor(tr, ix, ig)                                     !--- Ancestor of generation "ng" for each tr
     DO igen = ig+1, ngen
       tr(iq)%iqDescen = [tr(iq)%iqDescen, find(ix==iq .AND. tr%iGeneration==igen)]
       tr(iq)%nqDescen = SIZE(tr(iq)%iqDescen)
-      IF(igen == ig+1) tr(iq)%nqChilds=tr(iq)%nqDescen
+      IF(igen == ig+1) THEN
+        tr(iq)%nqChildren = tr(iq)%nqDescen
+        CALL addKey('nqChildren', int2str(tr(iq)%nqChildren), tr(iq)%keys)
+      END IF
     END DO
-  END DO
+    CALL addKey('iqDescen', strStack(int2str(tr(iq)%iqDescen)), tr(iq)%keys)
+  END DO
+  CALL addKey('nqDescen', int2str(tr(:)%nqDescen), tr(:)%keys)
 END SUBROUTINE indexUpdate
 !==============================================================================================================================
@@ -908,5 +1025,5 @@
 !===  * For each isotopes class, the <key>=<val> vector of each tracer is moved into the isotopes descriptor "isot"        ====
 !=== NOTES:                                                                                                                ====
-!===  * Most of the "isot" components have been defined in the calling routine (initIsotopes):                             ====
+!===  * Most of the "isot" components have been defined in the calling routine (readIsotopes):                             ====
 !===      parent,  nzone, zone(:),  niso, keys(:)%name,  ntiso, trac(:),  nphas, phas,  iqIsoPha(:,:),  itZonPhi(:,:)      ====
 !===  * Same syntax for isotopes file and "tracer.def": a tracers section contains one line for each of its isotopes       ====
@@ -916,5 +1033,5 @@
 !===  * The routine gives an error if a required isotope is not available in the database stored in "fnam"                 ====
 !==============================================================================================================================
-LOGICAL FUNCTION readIsotopesFile(fnam, isot) RESULT(lerr)
+LOGICAL FUNCTION readIsotopesFile_prv(fnam, isot) RESULT(lerr)
   CHARACTER(LEN=*),        INTENT(IN)    :: fnam                     !--- Input file name
   TYPE(isot_type), TARGET, INTENT(INOUT) :: isot(:)                  !--- Isotopes descriptors (field %parent must be defined!)
@@ -927,5 +1044,4 @@
   TYPE(trac_type),           POINTER ::   tt(:), t
   TYPE(dataBase_type),   ALLOCATABLE ::  tdb(:)
-  LOGICAL,               ALLOCATABLE :: liso(:)
   modname = 'readIsotopesFile'
 
@@ -953,14 +1069,12 @@
       is = strIdx(isot(iis)%keys(:)%name, t%name)                    !--- Index in "isot(iis)%keys(:)%name" of isotope "t%name"
       IF(is == 0) CYCLE
-      liso = reduceExpr(t%keys%val, vals)                            !--- Reduce expressions (for substituted variables)
-      IF(test(ANY(liso), lerr)) RETURN                               !--- Some non-numerical elements were found
-      isot(iis)%keys(is)%key = PACK(t%keys%key, MASK=.NOT.liso)
-      isot(iis)%keys(is)%val = PACK(  vals,     MASK=.NOT.liso)
+      IF(test(ANY(reduceExpr(t%keys%val, vals)), lerr)) RETURN       !--- Reduce expressions ; detect non-numerical elements
+      isot(iis)%keys(is)%key = t%keys%key
+      isot(iis)%keys(is)%val = vals
     END DO
 
     !--- CHECK FOR MISSING ISOTOPES (NO KEYS ALLOCATED)
-    liso = [( ALLOCATED(isot(iis)%keys(is)%key), is=1, SIZE(isot(iis)%keys) )]
-    IF(test(checkList(isot(iis)%keys(:)%name, .NOT.liso, &
-      'Check file "'//TRIM(fnam)//'" in section "'//TRIM(dBase(idb)%name)//'"', 'isotopes', 'missing'),lerr)) RETURN
+    IF(test(checkList(isot(iis)%keys(:)%name, .NOT.[( ALLOCATED(isot(iis)%keys(is)%key), is=1, SIZE(isot(iis)%keys) )], &
+      'Check file "'//TRIM(fnam)//'" in section "'//TRIM(dBase(idb)%name)//'"', 'isotopes', 'missing'), lerr)) RETURN
   END DO
 
@@ -975,7 +1089,34 @@
   CALL get_in('ok_iso_verif', isot(strIdx(isot%parent, 'H2O'))%check, .FALSE.)
 
-  lerr = dispIsotopes(isot, 'Isotopes parameters read from file "'//TRIM(fnam)//'"', modname)
-
-END FUNCTION readIsotopesFile
+  lerr = dispIsotopes()
+
+CONTAINS
+
+!------------------------------------------------------------------------------------------------------------------------------
+LOGICAL FUNCTION dispIsotopes() RESULT(lerr)
+  INTEGER :: ik, nk, ip, it, nt
+  CHARACTER(LEN=maxlen) :: prf
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: ttl(:), val(:,:)
+  CALL msg('Isotopes parameters read from file "'//TRIM(fnam)//'":', modname)
+  DO ip = 1, SIZE(isot)                                              !--- Loop on parents tracers
+    nk = SIZE(isot(ip)%keys(1)%key)                                  !--- Same keys for each isotope
+    nt = SIZE(isot(ip)%keys)                                         !--- Number of isotopes
+    prf = 'i'//REPEAT('s',nk+1)                                      !--- Profile for table printing
+    ALLOCATE(ttl(nk+2), val(nt,nk+1))
+    ttl(1:2) = ['it  ','name']; ttl(3:nk+2) = isot(ip)%keys(1)%key(:)!--- Titles line with keys names
+    val(:,1) = isot(ip)%keys(:)%name                                 !--- Values table 1st column: isotopes names  
+    DO ik = 1, nk
+      DO it = 1, nt
+        val(it,ik+1) = isot(ip)%keys(it)%val(ik)                     !--- Other columns: keys values
+      END DO
+    END DO
+    IF(test(fmsg('Problem with the table content', modname, dispTable(prf, ttl, val, &
+            cat([(it,it=1,nt)]), rFmt='(EN8.4)', nColMax=maxTableWidth, nHead=2, sub=modname)), lerr)) RETURN
+    DEALLOCATE(ttl, val)
+  END DO        
+END FUNCTION dispIsotopes
+!------------------------------------------------------------------------------------------------------------------------------
+
+END FUNCTION readIsotopesFile_prv
 !==============================================================================================================================
 
@@ -985,35 +1126,44 @@
 !===    * COMPUTE MOST OF THE RELATED QUANTITIES ("isot" COMPONENTS).                                                       ===
 !===    * COMPUTE FEW ISOTOPES-DEDICATED "trac" COMPONENTS                                                                  ===
-!===    * CALL readIsotopesFile TO GET PHYSICAL QUANTITIES (<key>=<val> PAIRS)                                              ===
+!===    * CALL readIsotopesFile_prv TO GET PHYSICAL QUANTITIES (<key>=<val> PAIRS)                                          ===
 !===      NOTE: THIS IS DONE HERE (IN A ROUTINE CALLED BY THE DYNAMIC), BECAUSE THE DYNAMIC NEEDS FEW PHYSICAL PARAMETERS.  ===
 !==============================================================================================================================
-LOGICAL FUNCTION initIsotopes(trac, isot) RESULT(lerr)
-  TYPE(trac_type), ALLOCATABLE, TARGET, INTENT(INOUT) :: trac(:)
-  TYPE(isot_type), ALLOCATABLE, TARGET, INTENT(INOUT) :: isot(:)
+LOGICAL FUNCTION readIsotopesFile(iNames) RESULT(lerr)
+  CHARACTER(LEN=maxlen), TARGET, OPTIONAL, INTENT(IN)  :: iNames(:)
   CHARACTER(LEN=maxlen), ALLOCATABLE :: p(:), str(:)                 !--- Temporary storage
-  CHARACTER(LEN=maxlen) :: iName
+  CHARACTER(LEN=maxlen) :: iName, modname
   CHARACTER(LEN=1)   :: ph                                           !--- Phase
-  INTEGER :: nbIso, ic, ip, iq, it, iz
+  INTEGER :: ic, ip, iq, it, iz
   LOGICAL, ALLOCATABLE :: ll(:)                                      !--- Mask
   TYPE(trac_type), POINTER   ::  t(:), t1
   TYPE(isot_type), POINTER   ::  i
   lerr = .FALSE.
-
-  t => trac
+  modname = 'readIsotopesFile'
+
+  t => tracers
 
   !--- GET FROM "tracers" THE FULL LIST OF AVAILABLE ISOTOPES CLASSES
   p = PACK(delPhase(t%parent), MASK = t%type=='tracer' .AND. t%iGeneration==1)
   CALL strReduce(p, nbIso)
-  ALLOCATE(isot(nbIso))
+
+  !--- CHECK WHETHER NEEDED ISOTOPES CLASSES "iNames" ARE AVAILABLE OR NOT
+  IF(PRESENT(iNames)) THEN
+    DO it = 1, SIZE(iNames)
+      IF(test(fmsg('No isotopes class "'//TRIM(iNames(it))//'" found among tracers', modname, ALL(p /= iNames(it))), lerr)) RETURN
+    END DO
+    p = iNames; nbIso = SIZE(p)
+  END IF
+  IF(ALLOCATED(isotopes)) DEALLOCATE(isotopes)
+  ALLOCATE(isotopes(nbIso))
 
   IF(nbIso==0) RETURN                                                !=== NO ISOTOPES: FINISHED
 
   !--- ISOTOPES RELATED VARIABLES ; NULL OR EMPTY IF NO ISOTOPES
-  isot(:)%parent = p
+  isotopes(:)%parent = p
   DO ic = 1, SIZE(p)                                                 !--- Loop on isotopes classes
-    i => isot(ic)
+    i => isotopes(ic)
     iname = i%parent                                                 !--- Current isotopes class name (parent tracer name)
 
-    !=== Isotopes childs of tracer "iname": mask, names, number (same for each phase of "iname")
+    !=== Isotopes children of tracer "iname": mask, names, number (same for each phase of "iname")
     ll = t(:)%type=='tracer' .AND. delPhase(t(:)%parent) == iname .AND. t(:)%phase == 'g'
     str = PACK(delPhase(t(:)%name), MASK = ll)                       !--- Effectively found isotopes of "iname"
@@ -1028,5 +1178,5 @@
     i%nzone = SIZE(i%zone)                                           !--- Tagging zones number for isotopes category "iname"
 
-    !=== Geographic tracers of the isotopes childs of tracer "iname" (same for each phase of "iname")
+    !=== Geographic tracers of the isotopes children of tracer "iname" (same for each phase of "iname")
     !    NOTE: One might like to create a similar variable for 2nd generation tagging tracers (tagging the gen1 tracers)
     str = PACK(delPhase(t(:)%name), MASK=ll)
@@ -1044,5 +1194,5 @@
     !=== Tables giving the index in a table of effectively found items for each dynamical tracer (1<=iq<=nqtot)
     DO iq = 1, SIZE(t)
-      t1 => trac(iq)
+      t1 => tracers(iq)
       IF(delPhase(t1%gen0Name)/=iname .OR. t1%iGeneration==0) CYCLE  !--- Only deal with tracers descending on "iname"
       t1%iso_iGroup = ic                                             !--- Isotopes family       idx in list "isotopes(:)%parent"
@@ -1055,5 +1205,5 @@
     !=== Table used to get iq (index in dyn array, size nqtot) from the isotope and phase indexes ; the full isotopes list
     !    (including tagging tracers) is sorted this way:  iso1, iso2, ..., iso1_zone1, iso2_zone1, ..., iso1_zoneN, iso2_zoneN
-    i%iqIsoPha = RESHAPE( [( (strIdx(t%name,  addPhase(i%trac(it),i%phase(ip:ip))),      it=1, i%ntiso), ip=1, i%nphas)], &
+    i%iqIsoPha = RESHAPE( [( (strIdx(t%name,  addPhase(i%trac(it),i%phase(ip:ip))),       it=1, i%ntiso), ip=1, i%nphas)], &
                          [i%ntiso, i%nphas] )
     !=== Table used to get ix (index in tagging tracers isotopes list, size ntiso) from the zone and isotope indexes
@@ -1062,37 +1212,82 @@
   END DO
 
-  !=== READ PHYSICAL PARAMETERS FROM "isotopes_params.def" FILE
-  lerr = readIsotopesFile('isotopes_params.def',isot)
-
-END FUNCTION initIsotopes
-!==============================================================================================================================
-
-
-!==============================================================================================================================
-LOGICAL FUNCTION dispIsotopes(ides, message, modname) RESULT(lerr)
-  TYPE(isot_type),  INTENT(IN) :: ides(:)                            !--- Isotopes descriptor vector
-  CHARACTER(LEN=*), INTENT(IN) :: message                            !--- Message to display
-  CHARACTER(LEN=*), INTENT(IN) :: modname                            !--- Calling subroutine name
-  INTEGER :: ik, nk, ip, it, nt
-  CHARACTER(LEN=maxlen) :: prf
-  CHARACTER(LEN=maxlen), ALLOCATABLE :: ttl(:), val(:,:)
-  CALL msg(TRIM(message)//':', modname)
-  DO ip = 1, SIZE(ides)                                              !--- Loop on parents tracers
-    nk = SIZE(ides(ip)%keys(1)%key)                                  !--- Same keys for each isotope
-    nt = SIZE(ides(ip)%keys)                                         !--- Number of isotopes
-    prf = 'i'//REPEAT('s',nk+1)                                      !--- Profile for table printing
-    ALLOCATE(ttl(nk+2), val(nt,nk+1))
-    ttl(1:2) = ['it  ','name']; ttl(3:nk+2) = ides(ip)%keys(1)%key(:)!--- Titles line with keys names
-    val(:,1) = ides(ip)%keys(:)%name                                 !--- Values table 1st column: isotopes names  
-    DO ik = 1, nk
-      DO it = 1, nt
-        val(it,ik+1) = ides(ip)%keys(it)%val(ik)                     !--- Other columns: keys values
-      END DO
+  !=== READ PHYSICAL PARAMETERS FROM isoFile FILE
+  IF(test(readIsotopesFile_prv(isoFile, isotopes), lerr)) RETURN
+
+  !=== CHECK CONSISTENCY
+  IF(test(testIsotopes(), lerr)) RETURN
+
+  !=== SELECT FIRST ISOTOPES CLASS OR, IF POSSIBLE, WATER CLASS
+  IF(.NOT.test(isoSelect(1, .TRUE.), lerr)) THEN; IF(isotope%parent == 'H2O') iH2O = ixIso; END IF
+
+CONTAINS
+
+!------------------------------------------------------------------------------------------------------------------------------
+LOGICAL FUNCTION testIsotopes() RESULT(lerr)     !--- MAKE SURE MEMBERS OF AN ISOTOPES FAMILY ARE PRESENT IN THE SAME PHASES
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: ix, it, ip, np, iz, nz
+  TYPE(isot_type), POINTER :: i
+  DO ix = 1, nbIso
+    i => isotopes(ix)
+    !--- Check whether each isotope and tagging isotopic tracer is present in the same number of phases
+    DO it = 1, i%ntiso
+      np = SUM([(COUNT(tracers(:)%name == addPhase(i%trac(it), i%phase(ip:ip))), ip=1, i%nphas)])
+      IF(test(fmsg(TRIM(int2str(np))//' phases instead of '//TRIM(int2str(i%nphas))//' for '//TRIM(i%trac(it)), &
+        modname, np /= i%nphas), lerr)) RETURN
     END DO
-    IF(test(fmsg('Problem with the table content', modname, dispTable(prf, ttl, val, &
-            cat([(it,it=1,nt)]), rFmt='(EN8.4)', nColMax=maxTableWidth, nHead=2, sub=modname)), lerr)) RETURN
-    DEALLOCATE(ttl, val)
-  END DO        
-END FUNCTION dispIsotopes
+    DO it = 1, i%niso
+      nz = SUM([(COUNT(i%trac == TRIM(i%trac(it))//'_'//i%zone(iz)), iz=1, i%nzone)])
+      IF(test(fmsg(TRIM(int2str(nz))//' tagging zones instead of '//TRIM(int2str(i%nzone))//' for '//TRIM(i%trac(it)), &
+        modname, nz /= i%nzone), lerr)) RETURN
+    END DO
+  END DO
+END FUNCTION testIsotopes
+!------------------------------------------------------------------------------------------------------------------------------
+
+END FUNCTION readIsotopesFile
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== THE ROUTINE isoSelect IS USED TO SWITCH FROM AN ISOTOPE FAMILY TO ANOTHER: ISOTOPES DEPENDENT PARAMETERS ARE UPDATED
+!     Single generic "isoSelect" routine, using the predefined index of the parent (fast version) or its name (first call).
+!==============================================================================================================================
+LOGICAL FUNCTION isoSelectByName(iName, lVerbose) RESULT(lerr)
+   IMPLICIT NONE
+   CHARACTER(LEN=*),  INTENT(IN) :: iName
+   LOGICAL, OPTIONAL, INTENT(IN) :: lVerbose
+   INTEGER :: iIso
+   LOGICAL :: lV
+   lV = .FALSE.; IF(PRESENT(lVerbose)) lV = lVerbose
+   iIso = strIdx(isotopes(:)%parent, iName)
+   IF(test(iIso == 0, lerr)) THEN
+      niso = 0; ntiso = 0; nzone = 0; nphas = 0; isoCheck=.FALSE.
+      CALL msg('no isotope family named "'//TRIM(iName)//'"', ll=lV)
+      RETURN
+   END IF
+   lerr = isoSelectByIndex(iIso, lV)
+END FUNCTION isoSelectByName
+!==============================================================================================================================
+LOGICAL FUNCTION isoSelectByIndex(iIso, lVerbose) RESULT(lerr)
+   IMPLICIT NONE
+   INTEGER,           INTENT(IN) :: iIso
+   LOGICAL, OPTIONAL, INTENT(IN) :: lVerbose
+   LOGICAL :: lV
+   lv = .FALSE.; IF(PRESENT(lVerbose)) lv = lVerbose
+   lerr = .FALSE.
+   IF(iIso == ixIso) RETURN                                          !--- Nothing to do if the index is already OK
+   lerr = iIso<=0 .OR. iIso>SIZE(isotopes)
+   CALL msg('Inconsistent isotopes family index '//TRIM(int2str(iIso))//': should be > 0 and <= '&
+          //TRIM(int2str(SIZE(isotopes)))//'"', ll = lerr .AND. lV)
+   IF(lerr) RETURN
+   ixIso = iIso                                                      !--- Update currently selected family index
+   isotope  => isotopes(ixIso)                                       !--- Select corresponding component
+   isoKeys  => isotope%keys;     niso     = isotope%niso
+   isoName  => isotope%trac;     ntiso    = isotope%ntiso
+   isoZone  => isotope%zone;     nzone    = isotope%nzone
+   isoPhas  => isotope%phase;    nphas    = isotope%nphas
+   itZonIso => isotope%itZonIso; isoCheck = isotope%check
+   iqIsoPha => isotope%iqIsoPha
+END FUNCTION isoSelectByIndex
 !==============================================================================================================================
 
@@ -1109,10 +1304,10 @@
   INTEGER :: iky, nky
   LOGICAL :: lo
-  lo=.FALSE.; IF(PRESENT(lOverWrite)) lo=lOverWrite
+  lo=.TRUE.; IF(PRESENT(lOverWrite)) lo=lOverWrite
   iky = strIdx(ky%key,key)
   IF(iky == 0) THEN
     nky = SIZE(ky%key)
     IF(nky == 0) THEN; ky%key = [key]; ky%val = [val]; ELSE; ky%key = [ky%key, key]; ky%val = [ky%val, val]; END IF
-  ELSE IF(lo) THEN                                                   !--- Overwriting
+  ELSE IF(lo) THEN
     ky%key(iky) = key; ky%val(iky) = val
   END IF
@@ -1125,8 +1320,15 @@
 !------------------------------------------------------------------------------------------------------------------------------
   INTEGER :: itr
-  LOGICAL :: lo
-  lo=.FALSE.; IF(PRESENT(lOverWrite)) lo=lOverWrite
-  DO itr = 1, SIZE(ky); CALL addKey_1(key, val, ky(itr), lo); END DO
+  DO itr = 1, SIZE(ky); CALL addKey_1(key, val, ky(itr), lOverWrite); END DO
 END SUBROUTINE addKey_m
+!==============================================================================================================================
+SUBROUTINE addKey_mm(key, val, ky, lOverWrite)
+  CHARACTER(LEN=*),  INTENT(IN)    :: key, val(:)
+  TYPE(keys_type),   INTENT(INOUT) :: ky(:)
+  LOGICAL, OPTIONAL, INTENT(IN)    :: lOverWrite
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: itr
+  DO itr = 1, SIZE(ky); CALL addKey_1(key, val(itr), ky(itr), lOverWrite); END DO
+END SUBROUTINE addKey_mm
 !==============================================================================================================================
 
@@ -1179,30 +1381,25 @@
 
 !==============================================================================================================================
-!=== getKey ROUTINE INITIALIZATION (TO BE EMBEDDED SOMEWHERE)  ================================================================
-!==============================================================================================================================
-SUBROUTINE getKey_init(tracers_, isotopes_)
-  TYPE(trac_type), OPTIONAL, INTENT(IN) ::  tracers_(:)
-  TYPE(isot_type), OPTIONAL, INTENT(IN) :: isotopes_(:)
-  IF(PRESENT( tracers_))  tracers =  tracers_
-  IF(PRESENT(isotopes_)) isotopes = isotopes_
-END SUBROUTINE getKey_init
-
-
-!==============================================================================================================================
 !================ GET THE VALUE OF A KEY FROM A "keys_type" DERIVED TYPE ; THE RESULT IS THE RETURNED VALUE ===================
 !==============================================================================================================================
-CHARACTER(LEN=maxlen) FUNCTION fgetKeyByIndex_s1(itr, keyn, ky, def_val) RESULT(val)
-  INTEGER,                    INTENT(IN) :: itr
-  CHARACTER(LEN=*),           INTENT(IN) :: keyn
-  TYPE(keys_type),            INTENT(IN) :: ky(:)
-  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: def_val
+CHARACTER(LEN=maxlen) FUNCTION fgetKeyIdx_s1(itr, keyn, ky, def_val, lerr) RESULT(val)
+  INTEGER,                    INTENT(IN)  :: itr
+  CHARACTER(LEN=*),           INTENT(IN)  :: keyn
+  TYPE(keys_type),            INTENT(IN)  :: ky(:)
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: def_val
+  LOGICAL,          OPTIONAL, INTENT(OUT) :: lerr
 !------------------------------------------------------------------------------------------------------------------------------
   INTEGER :: iky
-  iky = 0;  IF(itr >  0 .AND. itr <= SIZE(ky)) iky = strIdx(ky(itr)%key(:), keyn)
-  val = ''; IF(iky /= 0) val = ky(itr)%val(iky)                      !--- Key was found
-  IF(PRESENT(def_val) .AND. iky == 0) val = def_val                  !--- Default value from arguments
-END FUNCTION fgetKeyByIndex_s1
-!==============================================================================================================================
-CHARACTER(LEN=maxlen) FUNCTION fgetKeyByName_s1(tname, keyn, ky, def_val, lerr) RESULT(val)
+  LOGICAL :: ler
+  iky = 0; val = ''
+  IF(.NOT.test(itr <= 0 .OR. itr > SIZE(ky), ler)) iky = strIdx(ky(itr)%key(:), keyn)    !--- Correct index
+  IF(.NOT.test(iky == 0, ler))                     val = ky(itr)%val(iky)                !--- Found key
+  IF(iky == 0) THEN
+    IF(.NOT.test(.NOT.PRESENT(def_val), ler))      val = def_val                         !--- Default value
+  END IF
+  IF(PRESENT(lerr)) lerr = ler
+END FUNCTION fgetKeyIdx_s1
+!==============================================================================================================================
+CHARACTER(LEN=maxlen) FUNCTION fgetKeyNam_s1(tname, keyn, ky, def_val, lerr) RESULT(val)
   CHARACTER(LEN=*),           INTENT(IN)  :: tname, keyn
   TYPE(keys_type),            INTENT(IN)  :: ky(:)
@@ -1210,12 +1407,19 @@
   LOGICAL,          OPTIONAL, INTENT(OUT) :: lerr
 !------------------------------------------------------------------------------------------------------------------------------
-  INTEGER :: iky, itr
-  val = ''; iky = 0
-  itr = strIdx(ky(:)%name, tname)                                    !--- Get the index of the wanted tracer
-  IF(PRESENT(lerr)) lerr = itr==0; IF(itr == 0) RETURN
-  IF(itr >  0 .AND. itr <= SIZE(ky)) iky = strIdx(ky(itr)%key(:), keyn)
-  IF(iky /= 0) val = ky(itr)%val(iky)                                !--- Key was found
-  IF(PRESENT(def_val) .AND. iky == 0) val = def_val                  !--- Default value from arguments
-END FUNCTION fgetKeyByName_s1
+  val = fgetKeyIdx_s1(strIdx(ky(:)%name, tname), keyn, ky, def_val, lerr)
+END FUNCTION fgetKeyNam_s1
+!==============================================================================================================================
+FUNCTION fgetKey_sm(keyn, ky, def_val, lerr) RESULT(val)
+CHARACTER(LEN=maxlen),        ALLOCATABLE :: val(:)
+  CHARACTER(LEN=*),           INTENT(IN)  :: keyn
+  TYPE(keys_type),            INTENT(IN)  :: ky(:)
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: def_val
+  LOGICAL,          OPTIONAL, INTENT(OUT) :: lerr
+!------------------------------------------------------------------------------------------------------------------------------
+  LOGICAL :: ler(SIZE(ky))
+  INTEGER :: it
+  val = [(fgetKeyIdx_s1(it, keyn, ky, def_val, ler(it)), it = 1, SIZE(ky))]
+  IF(PRESENT(lerr)) lerr = ANY(ler)
+END FUNCTION fgetKey_sm
 !==============================================================================================================================
 
@@ -1236,43 +1440,68 @@
 !------------------------------------------------------------------------------------------------------------------------------
   CHARACTER(LEN=maxlen) :: tnam
-  INTEGER, ALLOCATABLE  :: is(:)
-  INTEGER :: i, itr
-  tnam = delPhase(strHead(tname,'_',.FALSE.))                        !--- Remove tag and phase
-  IF(PRESENT(ky)) THEN
-    val = fgetKeyByName_s1(tname, keyn, ky, lerr=lerr)               !--- "ky" and "tname"
-    IF(val /= '' .OR. lerr)      RETURN
-    val = fgetKeyByName_s1(tnam,  keyn, ky, lerr=lerr)               !--- "ky" and "tnam"
+  tnam = strHead(delPhase(tname),'_',.FALSE.)                                            !--- Remove tag and phase
+  IF(PRESENT(ky)) THEN                                                                   !=== KEY FROM "ky"
+               val = fgetKeyNam_s1(tname, keyn, ky,           lerr=lerr)                 !--- "ky" and "tname"
+    IF( lerr ) val = fgetKeyNam_s1(tnam,  keyn, ky,           lerr=lerr)                 !--- "ky" and "tnam"
   ELSE
-    IF(.NOT.ALLOCATED(tracers))  RETURN
-    val = fgetKeyByName_s1(tname, keyn, tracers(:)%keys, lerr=lerr)  !--- "tracers" and "tname"
-    IF(val /= ''.AND..NOT.lerr)  RETURN
-    IF(.NOT.ALLOCATED(isotopes)) RETURN
-    IF(SIZE(isotopes) == 0)      RETURN
-    !--- Search the "is" isotopes class index of the isotope named "tnam"
-    is = find([(ANY(isotopes(i)%keys(:)%name == tnam), i=1, SIZE(isotopes))])
-    IF(test(SIZE(is) == 0,lerr)) RETURN
-    val = fgetKeyByName_s1(tname, keyn, isotopes(is(1))%keys(:),lerr=lerr)!--- "isotopes" and "tnam"
+    IF(         .NOT.test(.NOT.ALLOCATED(tracers ), lerr)) lerr = SIZE(tracers ) == 0    !=== KEY FROM "tracers"
+    IF(.NOT.lerr) THEN
+               val = fgetKeyNam_s1(tname, keyn, tracers%keys, lerr=lerr)                 !--- "ky" and "tname"
+      IF(lerr) val = fgetKeyNam_s1(tnam,  keyn, tracers%keys, lerr=lerr)                 !--- "ky" and "tnam"
+    END IF
+    IF(lerr.AND..NOT.test(.NOT.ALLOCATED(isotopes), lerr)) lerr = SIZE(isotopes) == 0    !=== KEY FROM "isotope"
+    IF(.NOT.lerr) THEN
+               val = fgetKeyNam_s1(tname, keyn, isotope%keys, lerr=lerr)                 !--- "ky" and "tname"
+      IF(lerr) val = fgetKeyNam_s1(tnam,  keyn, isotope%keys, lerr=lerr)                 !--- "ky" and "tnam"
+    END IF
   END IF
 END FUNCTION getKeyByName_s1
 !==============================================================================================================================
-LOGICAL FUNCTION getKeyByName_sm(keyn, val, tname, ky) RESULT(lerr)
+LOGICAL FUNCTION getKeyByName_sm(keyn, val, tname, ky, nam) RESULT(lerr)
+  CHARACTER(LEN=*),                             INTENT(IN)  :: keyn
+  CHARACTER(LEN=maxlen), ALLOCATABLE,           INTENT(OUT) :: val(:)
+  CHARACTER(LEN=*),           TARGET, OPTIONAL, INTENT(IN)  :: tname(:)
+  TYPE(keys_type),            TARGET, OPTIONAL, INTENT(IN)  :: ky(:)
+  CHARACTER(LEN=maxlen), ALLOCATABLE, OPTIONAL, INTENT(OUT) :: nam(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: names(:)
+  TYPE(keys_type),       POINTER     ::  keys(:)
+  LOGICAL :: lk, lt, li, ll
+  INTEGER :: iq, nq
+
+  !--- DETERMINE THE DATABASE TO BE USED (ky, tracers or isotope)
+  lk = PRESENT(ky)
+  lt = .NOT.lk .AND. ALLOCATED(tracers);  IF(lt) lt = SIZE(tracers)  /= 0; IF(lt) lt = ANY(tracers(1)%keys%key(:) == keyn)
+  li = .NOT.lt .AND. ALLOCATED(isotopes); IF(li) li = SIZE(isotopes) /= 0; IF(li) li = ANY(isotope%keys(1)%key(:) == keyn)
+
+  IF(test(.NOT.ANY([lk,lt,li]), lerr)) RETURN
+  IF(lk) keys => ky(:)
+  IF(lt) keys => tracers(:)%keys
+  IF(li) keys => isotope%keys(:)
+
+  !--- DETERMINE THE NAMES
+  IF(PRESENT(tname)) THEN
+    ALLOCATE(names(SIZE(tname))); names(:) = tname(:)
+  ELSE
+    ALLOCATE(names(SIZE(keys)));  names(:) = keys(:)%name
+  END IF
+  nq = SIZE(names); ALLOCATE(val(nq)); IF(PRESENT(nam)) THEN; ALLOCATE(nam(nq)); nam(:) = names(:); END IF
+
+  !--- GET THE DATA
+  lerr = ANY([(getKeyByName_s1(keyn, val(iq), names(iq), keys(:)), iq=1, nq)])
+
+END FUNCTION getKeyByName_sm
+!==============================================================================================================================
+LOGICAL FUNCTION getKeyByName_s1m(keyn, val, tname, ky) RESULT(lerr)
   CHARACTER(LEN=*),                   INTENT(IN)  :: keyn
-  CHARACTER(LEN=maxlen), ALLOCATABLE, INTENT(OUT) ::   val(:)
-  CHARACTER(LEN=*),         OPTIONAL, INTENT(IN)  :: tname(:)
-  TYPE(keys_type),          OPTIONAL, INTENT(IN)  ::    ky(:)
-!------------------------------------------------------------------------------------------------------------------------------
-  TYPE(keys_type),           POINTER :: k(:)
-  CHARACTER(LEN=maxlen), ALLOCATABLE :: n(:)
-  INTEGER :: iq, nq
-  IF(test(.NOT.(PRESENT(tname).OR.PRESENT(ky)), lerr)) RETURN
-  IF(PRESENT(ky   )) nq = SIZE(ky%name)
-  IF(PRESENT(tname)) nq = SIZE(  tname)
-  ALLOCATE(val(nq))
-  IF(PRESENT(tname)) THEN
-    IF(     PRESENT(ky)) lerr = ANY([(getKeyByName_s1(keyn, val(iq), tname(iq),   ky), iq=1, nq)])
-    IF(.NOT.PRESENT(ky)) lerr = ANY([(getKeyByName_s1(keyn, val(iq), tname(iq)      ), iq=1, nq)])
-  ELSE;                  lerr = ANY([(getKeyByName_s1(keyn, val(iq), ky(iq)%name, ky), iq=1, nq)])
-  END IF
-END FUNCTION getKeyByName_sm
+  CHARACTER(LEN=maxlen), ALLOCATABLE, INTENT(OUT) :: val(:)
+  CHARACTER(LEN=*),                   INTENT(IN)  :: tname
+  TYPE(keys_type),          OPTIONAL, INTENT(IN)  :: ky(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen) :: sval
+  lerr = getKeyByName_s1(keyn, sval, tname, ky)
+  IF(test(fmsg('missing key "'//TRIM(keyn)//'" for tracer "'//TRIM(tname)//'"', modname, lerr), lerr)) RETURN
+  lerr = strParse(sval, ',', val)
+END FUNCTION getKeyByName_s1m
 !==============================================================================================================================
 LOGICAL FUNCTION getKeyByName_i1(keyn, val, tname, ky) RESULT(lerr)
@@ -1284,30 +1513,39 @@
   CHARACTER(LEN=maxlen) :: sval
   INTEGER :: ierr
-  IF(     PRESENT(ky)) lerr = getKeyByName_s1(keyn, sval, tname, ky)
-  IF(.NOT.PRESENT(ky)) lerr = getKeyByName_s1(keyn, sval, tname)
-  IF(test(fmsg('key "'//TRIM(keyn)//'" or tracer "'//TRIM(tname)//'" is missing',        modname, lerr), lerr)) RETURN
+  lerr = getKeyByName_s1(keyn, sval, tname, ky)
+  IF(test(fmsg('key "'//TRIM(keyn)//'" or tracer "'//TRIM(tname)//'" is missing', modname, lerr), lerr)) RETURN
   READ(sval, *, IOSTAT=ierr) val
+  IF(test(fmsg('key "'//TRIM(keyn)//'" of tracer "'//TRIM(tname)//'" is not an integer', modname, ierr/=0), lerr)) RETURN
+END FUNCTION getKeyByName_i1
+!==============================================================================================================================
+LOGICAL FUNCTION getKeyByName_im(keyn, val, tname, ky) RESULT(lerr)
+  CHARACTER(LEN=*),                   INTENT(IN)  :: keyn
+  INTEGER,               ALLOCATABLE, INTENT(OUT) ::   val(:)
+  CHARACTER(LEN=*), OPTIONAL, TARGET, INTENT(IN)  :: tname(:)
+  TYPE(keys_type),  OPTIONAL, TARGET, INTENT(IN)  ::    ky(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: sval(:), nam(:)
+  INTEGER :: ierr, iq
+  IF(test(getKeyByName_sm(keyn, sval, tname, ky, nam), lerr)) RETURN
+  ALLOCATE(val(SIZE(sval)))
+  DO iq = 1, SIZE(sval)                                              !--- CONVERT THE KEYS TO INTEGERS
+    READ(sval(iq), *, IOSTAT=ierr) val(iq)
+    IF(test(fmsg('key "'//TRIM(keyn)//'" of "'//TRIM(nam(iq))//'" is not an integer', modname, ierr/=0), lerr)) RETURN
+  END DO
+END FUNCTION getKeyByName_im
+!==============================================================================================================================
+LOGICAL FUNCTION getKeyByName_i1m(keyn, val, tname, ky) RESULT(lerr)
+  CHARACTER(LEN=*),           INTENT(IN)  :: keyn
+  INTEGER,       ALLOCATABLE, INTENT(OUT) :: val(:)
+  CHARACTER(LEN=*),           INTENT(IN)  :: tname
+  TYPE(keys_type),  OPTIONAL, INTENT(IN)  ::  ky(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: v(:)
+  INTEGER :: ierr, iq
+  IF(test(getKeyByName_s1m(keyn, v, tname, ky), lerr)) RETURN
+  ALLOCATE(val(SIZE(v)))
+  lerr = .FALSE.; DO iq=1, SIZE(v); READ(v(iq), *, IOSTAT=ierr) val(iq); lerr = lerr .OR. ierr /= 0; END DO
   IF(test(fmsg('key "'//TRIM(keyn)//'" of tracer "'//TRIM(tname)//'" is not an integer', modname, lerr), lerr)) RETURN
-END FUNCTION getKeyByName_i1
-!==============================================================================================================================
-LOGICAL FUNCTION getKeyByName_im(keyn, val, tname, ky) RESULT(lerr)
-  CHARACTER(LEN=*),           INTENT(IN)  :: keyn
-  INTEGER,       ALLOCATABLE, INTENT(OUT) ::   val(:)
-  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: tname(:)
-  TYPE(keys_type),  OPTIONAL, INTENT(IN)  ::    ky(:)
-!------------------------------------------------------------------------------------------------------------------------------
-  TYPE(keys_type),           POINTER :: k(:)
-  CHARACTER(LEN=maxlen), ALLOCATABLE :: n(:)
-  INTEGER :: iq, nq
-  IF(test(.NOT.(PRESENT(tname).OR.PRESENT(ky)), lerr)) RETURN
-  IF(PRESENT(ky   )) nq = SIZE(ky%name)
-  IF(PRESENT(tname)) nq = SIZE(  tname)
-  ALLOCATE(val(nq))
-  IF(PRESENT(tname)) THEN
-    IF(     PRESENT(ky)) lerr = ANY([(getKeyByName_i1(keyn, val(iq), tname(iq),   ky), iq=1, nq)])
-    IF(.NOT.PRESENT(ky)) lerr = ANY([(getKeyByName_i1(keyn, val(iq), tname(iq)      ), iq=1, nq)])
-  ELSE;                  lerr = ANY([(getKeyByName_i1(keyn, val(iq), ky(iq)%name, ky), iq=1, nq)])
-  END IF
-END FUNCTION getKeyByName_im
+END FUNCTION getKeyByName_i1m
 !==============================================================================================================================
 LOGICAL FUNCTION getKeyByName_r1(keyn, val, tname, ky) RESULT(lerr)
@@ -1319,9 +1557,8 @@
   CHARACTER(LEN=maxlen) :: sval
   INTEGER :: ierr
-  IF(     PRESENT(ky)) lerr = getKeyByName_s1(keyn, sval, tname, ky)
-  IF(.NOT.PRESENT(ky)) lerr = getKeyByName_s1(keyn, sval, tname)
+  lerr = getKeyByName_s1(keyn, sval, tname, ky)
   IF(test(fmsg('key "'//TRIM(keyn)//'" or tracer "'//TRIM(tname)//'" is missing',    modname, lerr), lerr)) RETURN
   READ(sval, *, IOSTAT=ierr) val
-  IF(test(fmsg('key "'//TRIM(keyn)//'" of tracer "'//TRIM(tname)//'" is not a real', modname, lerr), lerr)) RETURN
+  IF(test(fmsg('key "'//TRIM(keyn)//'" of tracer "'//TRIM(tname)//'" is not a real', modname, ierr/=0), lerr)) RETURN
 END FUNCTION getKeyByName_r1
 !==============================================================================================================================
@@ -1332,17 +1569,64 @@
   TYPE(keys_type),  OPTIONAL, INTENT(IN)  ::    ky(:)
 !------------------------------------------------------------------------------------------------------------------------------
-  TYPE(keys_type),           POINTER :: k(:)
-  CHARACTER(LEN=maxlen), ALLOCATABLE :: n(:)
-  INTEGER :: iq, nq
-  IF(test(.NOT.(PRESENT(tname).OR.PRESENT(ky)), lerr)) RETURN
-  IF(PRESENT(ky   )) nq = SIZE(ky%name)
-  IF(PRESENT(tname)) nq = SIZE(  tname)
-  ALLOCATE(val(nq))
-  IF(PRESENT(tname)) THEN
-    IF(     PRESENT(ky)) lerr = ANY([(getKeyByName_r1(keyn, val(iq), tname(iq),   ky), iq=1, nq)])
-    IF(.NOT.PRESENT(ky)) lerr = ANY([(getKeyByName_r1(keyn, val(iq), tname(iq)      ), iq=1, nq)])
-  ELSE;                  lerr = ANY([(getKeyByName_r1(keyn, val(iq), ky(iq)%name, ky), iq=1, nq)])
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: sval(:), nam(:)
+  INTEGER :: ierr, iq
+  IF(test(getKeyByName_sm(keyn, sval, tname, ky, nam), lerr)) RETURN
+  ALLOCATE(val(SIZE(sval)))
+  DO iq = 1, SIZE(sval)                                              !--- CONVERT THE KEYS TO INTEGERS
+    READ(sval(iq), *, IOSTAT=ierr) val(iq)
+    IF(test(fmsg('key "'//TRIM(keyn)//'" of "'//TRIM(nam(iq))//'" is not a real', modname, ierr/=0), lerr)) RETURN
+  END DO
+END FUNCTION getKeyByName_rm
+!==============================================================================================================================
+LOGICAL FUNCTION getKeyByName_r1m(keyn, val, tname, ky) RESULT(lerr)
+  CHARACTER(LEN=*),           INTENT(IN)  :: keyn
+  REAL,          ALLOCATABLE, INTENT(OUT) :: val(:)
+  CHARACTER(LEN=*),           INTENT(IN)  :: tname
+  TYPE(keys_type),  OPTIONAL, INTENT(IN)  ::  ky(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: v(:)
+  INTEGER :: ierr, iq
+  IF(     PRESENT(ky)) lerr = getKeyByName_s1m(keyn, v, tname, ky)
+  IF(.NOT.PRESENT(ky)) lerr = getKeyByName_s1m(keyn, v, tname)
+  ALLOCATE(val(SIZE(v)))
+  lerr = .FALSE.; DO iq=1, SIZE(v); READ(v(iq), *, IOSTAT=ierr) val(iq); lerr = lerr .OR. ierr /= 0; END DO
+  IF(test(fmsg('key "'//TRIM(keyn)//'" of tracer "'//TRIM(tname)//'" is not a real', modname, lerr), lerr)) RETURN
+END FUNCTION getKeyByName_r1m
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== ROUTINES TO GET OR PUT TRACERS AND ISOTOPES DATABASES, SINCE tracers AND isotopes ARE PRIVATE VARIABLES ==================
+!==============================================================================================================================
+SUBROUTINE setKeysDBase(tracers_, isotopes_, isotope_)
+  TYPE(trac_type), OPTIONAL, INTENT(IN) ::  tracers_(:)
+  TYPE(isot_type), OPTIONAL, INTENT(IN) :: isotopes_(:)
+  TYPE(isot_type), OPTIONAL, INTENT(IN) :: isotope_
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE(isot_type), ALLOCATABLE :: iso(:)
+  INTEGER :: ix, nbIso
+  IF(PRESENT( tracers_)) THEN;  tracers =  tracers_; ELSE; ALLOCATE( tracers(0)); END IF
+  IF(PRESENT(isotopes_)) THEN; isotopes = isotopes_; ELSE; ALLOCATE(isotopes(0)); END IF
+  IF(PRESENT(isotope_ )) THEN
+    ix = strIdx(isotopes(:)%parent, isotope%parent)
+    IF(ix /= 0) THEN
+      isotopes(ix) = isotope_
+    ELSE
+      nbIso = SIZE(isotopes); ALLOCATE(iso(nbIso+1)); iso(1:nbIso) = isotopes; iso(nbIso+1) = isotope_
+      CALL MOVE_ALLOC(FROM=iso, TO=isotopes)
+    END IF
   END IF
-END FUNCTION getKeyByName_rm
+END SUBROUTINE setKeysDBase
+!==============================================================================================================================
+SUBROUTINE getKeysDBase(tracers_, isotopes_, isotope_)
+  TYPE(trac_type), OPTIONAL, ALLOCATABLE, INTENT(OUT) ::  tracers_(:)
+  TYPE(isot_type), OPTIONAL, ALLOCATABLE, INTENT(OUT) :: isotopes_(:)
+  TYPE(isot_type), OPTIONAL,              INTENT(OUT) :: isotope_
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: ix
+  IF(PRESENT( tracers_)) THEN;  tracers_ =  tracers; ELSE; ALLOCATE( tracers_(0)); END IF
+  IF(PRESENT(isotopes_)) THEN; isotopes_ = isotopes; ELSE; ALLOCATE(isotopes_(0)); END IF
+  IF(PRESENT(isotope_ )) THEN; ix = strIdx(isotopes(:)%parent, isotope%parent); IF(ix /= 0) isotope_=isotopes(ix); END IF
+END SUBROUTINE getKeysDBase
 !==============================================================================================================================
 
@@ -1352,5 +1636,5 @@
 !==============================================================================================================================
 ELEMENTAL CHARACTER(LEN=maxlen) FUNCTION delPhase(s) RESULT(out)
-  CHARACTER(LEN=*), INTENT(IN)  :: s
+  CHARACTER(LEN=*), INTENT(IN) :: s
 !------------------------------------------------------------------------------------------------------------------------------
   INTEGER :: ix, ip, ns
@@ -1518,59 +1802,70 @@
 !=== GET THE NAME(S) OF THE ANCESTOR(S) OF TRACER(S) "tname" AT GENERATION "igen"  IN THE TRACERS DESCRIPTORS LIST "tr" =======
 !==============================================================================================================================
-CHARACTER(LEN=maxlen) FUNCTION ancestor_1(t, tname, igen) RESULT(out)
-  TYPE(trac_type),   INTENT(IN) :: t(:)
-  CHARACTER(LEN=*),  INTENT(IN) :: tname
-  INTEGER, OPTIONAL, INTENT(IN) :: igen
-!------------------------------------------------------------------------------------------------------------------------------
-  INTEGER :: ig, ix
-  ig = 0; IF(PRESENT(igen)) ig = igen
-  ix = idxAncestor_1(t, tname, ig)
+SUBROUTINE ancestor_1(t, out, tname, igen)
+  TYPE(trac_type),       INTENT(IN)  :: t(:)
+  CHARACTER(LEN=maxlen), INTENT(OUT) :: out
+  CHARACTER(LEN=*),      INTENT(IN)  :: tname
+  INTEGER,     OPTIONAL, INTENT(IN)  :: igen
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: ix
+  CALL idxAncestor_1(t, ix, tname, igen)
   out = ''; IF(ix /= 0) out = t(ix)%name
-END FUNCTION ancestor_1
-!==============================================================================================================================
-FUNCTION ancestor_m(t, tname, igen) RESULT(out)
-  CHARACTER(LEN=maxlen), ALLOCATABLE     ::   out(:)
-  TYPE(trac_type),            INTENT(IN) ::     t(:)
-  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: tname(:)
-  INTEGER,          OPTIONAL, INTENT(IN) :: igen
-!------------------------------------------------------------------------------------------------------------------------------
-  INTEGER, ALLOCATABLE :: ix(:)
+END SUBROUTINE ancestor_1
+!==============================================================================================================================
+SUBROUTINE ancestor_mt(t, out, tname, igen)
+  TYPE(trac_type),       INTENT(IN)  :: t(:)
+  CHARACTER(LEN=*),      INTENT(IN)  :: tname(:)
+  CHARACTER(LEN=maxlen), INTENT(OUT) :: out(SIZE(tname))
+  INTEGER,     OPTIONAL, INTENT(IN)  :: igen
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: ix(SIZE(tname))
+  CALL idxAncestor_mt(t, ix, tname, igen)
+  out(:) = ''; WHERE(ix /= 0) out = t(ix)%name
+END SUBROUTINE ancestor_mt
+!==============================================================================================================================
+SUBROUTINE ancestor_m(t, out, igen)
+  TYPE(trac_type),       INTENT(IN)  :: t(:)
+  CHARACTER(LEN=maxlen), INTENT(OUT) :: out(SIZE(t))
+  INTEGER,     OPTIONAL, INTENT(IN)  :: igen
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: ix(SIZE(t))
+  CALL idxAncestor_m(t, ix, igen)
+  out(:) = ''; WHERE(ix /= 0) out = t(ix)%name
+END SUBROUTINE ancestor_m
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== GET THE INDEX(ES) OF THE GENERATION "igen" ANCESTOR(S) OF "tname" (t%name IF UNSPECIFIED) IN THE "t" LIST ================
+!==============================================================================================================================
+SUBROUTINE idxAncestor_1(t, idx, tname, igen)
+  TYPE(trac_type),   INTENT(IN)  :: t(:)
+  INTEGER,           INTENT(OUT) :: idx
+  CHARACTER(LEN=*),  INTENT(IN)  :: tname
+  INTEGER, OPTIONAL, INTENT(IN)  :: igen
   INTEGER :: ig
   ig = 0; IF(PRESENT(igen)) ig = igen
-  IF(     PRESENT(tname)) ix = idxAncestor_m(t, tname,     ig)
-  IF(.NOT.PRESENT(tname)) ix = idxAncestor_m(t, t(:)%name, ig)
-  ALLOCATE(out(SIZE(ix))); out(:) = ''
-  WHERE(ix /= 0) out = t(ix)%name
-END FUNCTION ancestor_m
-!==============================================================================================================================
-
-
-!==============================================================================================================================
-!=== GET THE INDEX(ES) OF THE GENERATION "igen" ANCESTOR(S) OF "tname" (t%name IF UNSPECIFIED) IN THE "t" LIST ================
-!==============================================================================================================================
-INTEGER FUNCTION idxAncestor_1(t, tname, igen) RESULT(out)
-  TYPE(trac_type),   INTENT(IN) :: t(:)
-  CHARACTER(LEN=*),  INTENT(IN) :: tname
-  INTEGER, OPTIONAL, INTENT(IN) :: igen
-!------------------------------------------------------------------------------------------------------------------------------
-  INTEGER :: ig
-  ig = 0; IF(PRESENT(igen)) ig = igen
-  out = strIdx(t(:)%name, tname)
-  IF(out == 0)                 RETURN            !--- Tracer not found
-  IF(t(out)%iGeneration <= ig) RETURN            !--- Tracer has a lower generation number than asked generation 'igen"
-  DO WHILE(t(out)%iGeneration > ig); out = strIdx(t(:)%name, t(out)%parent); END DO
-END FUNCTION idxAncestor_1
-!==============================================================================================================================
-FUNCTION idxAncestor_m(t, tname, igen) RESULT(out)
-  INTEGER,          ALLOCATABLE          ::   out(:)
-  TYPE(trac_type),            INTENT(IN) ::     t(:)
-  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: tname(:)
-  INTEGER,          OPTIONAL, INTENT(IN) :: igen
-!------------------------------------------------------------------------------------------------------------------------------
-  INTEGER :: ig, ix
-  ig = 0; IF(PRESENT(igen)) ig = igen
-  IF(     PRESENT(tname)) out = [(idxAncestor_1(t, tname(ix),  ig), ix=1, SIZE(tname))]
-  IF(.NOT.PRESENT(tname)) out = [(idxAncestor_1(t, t(ix)%name, ig), ix=1, SIZE(t))]
-END FUNCTION idxAncestor_m
+  idx = strIdx(t(:)%name, tname)
+  IF(idx == 0)                 RETURN            !--- Tracer not found
+  IF(t(idx)%iGeneration <= ig) RETURN            !--- Tracer has a lower generation number than asked generation 'igen"
+  DO WHILE(t(idx)%iGeneration > ig); idx = strIdx(t(:)%name, t(idx)%parent); END DO
+END SUBROUTINE idxAncestor_1
+!------------------------------------------------------------------------------------------------------------------------------
+SUBROUTINE idxAncestor_mt(t, idx, tname, igen)
+  TYPE(trac_type),   INTENT(IN)  :: t(:)
+  CHARACTER(LEN=*),  INTENT(IN)  :: tname(:)
+  INTEGER,           INTENT(OUT) :: idx(SIZE(tname))
+  INTEGER, OPTIONAL, INTENT(IN)  :: igen
+  INTEGER :: ix
+  DO ix = 1, SIZE(tname); CALL idxAncestor_1(t, idx(ix), tname(ix), igen); END DO
+END SUBROUTINE idxAncestor_mt
+!------------------------------------------------------------------------------------------------------------------------------
+SUBROUTINE idxAncestor_m(t, idx, igen)
+  TYPE(trac_type),   INTENT(IN)  :: t(:)
+  INTEGER,           INTENT(OUT) :: idx(SIZE(t))
+  INTEGER, OPTIONAL, INTENT(IN)  :: igen
+  INTEGER :: ix
+  DO ix = 1, SIZE(t); CALL idxAncestor_1(t, idx(ix), t(ix)%name, igen); END DO
+END SUBROUTINE idxAncestor_m
 !==============================================================================================================================
 
Index: LMDZ6/trunk/libf/misc/strings_mod.F90
===================================================================
--- LMDZ6/trunk/libf/misc/strings_mod.F90	(revision 4319)
+++ LMDZ6/trunk/libf/misc/strings_mod.F90	(revision 4325)
@@ -1248,5 +1248,4 @@
   INTEGER      :: i, ni, io
 
-!  modname = 'reduceExpr_basic'
   lerr = .FALSE.
   IF(is_numeric(str)) THEN; val=TRIM(str); RETURN; END IF
Index: LMDZ6/trunk/libf/misc/trac_types_mod.F90
===================================================================
--- LMDZ6/trunk/libf/misc/trac_types_mod.F90	(revision 4319)
+++ 	(revision )
@@ -1,55 +1,0 @@
-MODULE trac_types_mod
-
-  USE strings_mod, ONLY: maxlen
-  PRIVATE
-
-!=== TRACERS DESCRIPTOR DERIVED TYPE AND ASSOCIATED ROUTINES INTERFACES =======================================================
-  PUBLIC :: trac_type, isot_type, keys_type
-!------------------------------------------------------------------------------------------------------------------------------
-  TYPE :: keys_type                                        !=== TYPE FOR A SET OF KEYS ASSOCIATED TO AN ELEMENT
-    CHARACTER(LEN=maxlen)              :: name             !--- Tracer name
-    CHARACTER(LEN=maxlen), ALLOCATABLE :: key(:)           !--- Keys string list
-    CHARACTER(LEN=maxlen), ALLOCATABLE :: val(:)           !--- Corresponding values string list
-  END TYPE keys_type
-!------------------------------------------------------------------------------------------------------------------------------
-  TYPE :: trac_type                                        !=== TYPE FOR A SINGLE TRACER NAMED "name"
-    CHARACTER(LEN=maxlen) :: name        = ''              !--- Name of the tracer
-    CHARACTER(LEN=maxlen) :: gen0Name    = ''              !--- First generation ancestor name
-    CHARACTER(LEN=maxlen) :: parent      = ''              !--- Parent name
-    CHARACTER(LEN=maxlen) :: longName    = ''              !--- Long name (with advection scheme suffix)
-    CHARACTER(LEN=maxlen) :: type        = 'tracer'        !--- Type  (so far: 'tracer' / 'tag')
-    CHARACTER(LEN=maxlen) :: phase       = 'g'             !--- Phase ('g'as / 'l'iquid / 's'olid)
-    CHARACTER(LEN=maxlen) :: component   = ''              !--- Coma-separated list of components (Ex: lmdz,inca)
-    INTEGER               :: iadv        = 10              !--- Advection scheme used
-    INTEGER               :: iGeneration = -1              !--- Generation number (>=0)
-    LOGICAL               :: isAdvected  = .FALSE.         !--- "true" tracers: iadv > 0.   COUNT(isAdvected )=nqtrue
-    LOGICAL               :: isInPhysics = .TRUE.          !--- "true" tracers: in tr_seri. COUNT(isInPhysics)=nqtottr
-    INTEGER               :: iqParent    = 0               !--- Parent index
-    INTEGER,  ALLOCATABLE :: iqDescen(:)                   !--- Descendants index (in growing generation order)
-    INTEGER               :: nqDescen    = 0               !--- Number of descendants (all generations)
-    INTEGER               :: nqChilds    = 0               !--- Number of childs    (first generation)
-    INTEGER               :: iso_iGroup  = 0               !--- Isotopes group index in isotopes(:)
-    INTEGER               :: iso_iName   = 0               !--- Isotope  name  index in isotopes(iso_iGroup)%trac(:)
-    INTEGER               :: iso_iZone   = 0               !--- Isotope  zone  index in isotopes(iso_iGroup)%zone(:)
-    INTEGER               :: iso_iPhase  = 0               !--- Isotope  phase index in isotopes(iso_iGroup)%phas
-    TYPE(keys_type)       :: keys                          !--- <key>=<val> pairs vector
-  END TYPE trac_type
-!------------------------------------------------------------------------------------------------------------------------------
-  TYPE :: isot_type                                        !=== TYPE FOR AN ISOTOPES FAMILY DESCENDING ON TRACER "parent"
-    CHARACTER(LEN=maxlen)              :: parent           !--- Isotopes family name (parent tracer name ; ex: H2O)
-    LOGICAL                            :: check=.FALSE.    !--- Triggering of the checking routines
-    TYPE(keys_type),       ALLOCATABLE :: keys(:)          !--- Isotopes keys/values pairs list     (length: niso)
-    CHARACTER(LEN=maxlen), ALLOCATABLE :: trac(:)          !--- Isotopes + tagging tracers list     (length: ntiso)
-    CHARACTER(LEN=maxlen), ALLOCATABLE :: zone(:)          !--- Geographic tagging zones names list (length: nzone)
-    CHARACTER(LEN=maxlen)              :: phase = 'g'      !--- Phases list: [g][l][s]              (length: nphas)
-    INTEGER                            :: niso  = 0        !--- Number of isotopes, excluding tagging tracers
-    INTEGER                            :: nzone = 0        !--- Number of geographic tagging zones
-    INTEGER                            :: ntiso = 0        !--- Number of isotopes, including tagging tracers
-    INTEGER                            :: nphas = 0        !--- Number phases
-    INTEGER,               ALLOCATABLE :: iqIsoPha(:,:)    !--- Idx in "tracers(1:nqtot)" = f(name(1:ntiso)),phas)
-                                                           !---        "iqIsoPha" former name: "iqiso"
-    INTEGER,               ALLOCATABLE :: itZonIso(:,:)    !--- Idx in "trac(1:ntiso)" = f(zone, name(1:niso))
-                                                           !---        "itZonIso" former name: "index_trac"
-  END TYPE isot_type
-
-END MODULE trac_types_mod
Index: LMDZ6/trunk/libf/phylmd/dyn1d/old_lmdz1d.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/dyn1d/old_lmdz1d.F90	(revision 4319)
+++ LMDZ6/trunk/libf/phylmd/dyn1d/old_lmdz1d.F90	(revision 4325)
@@ -613,5 +613,5 @@
       call init_dimphy1D(1,llm)
       call suphel
-      call infotrac_init 
+      call init_infotrac
 
       if (nqtot>nqmx) STOP 'Augmenter nqmx dans lmdz1d.F'
Index: LMDZ6/trunk/libf/phylmd/dyn1d/scm.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/dyn1d/scm.F90	(revision 4319)
+++ LMDZ6/trunk/libf/phylmd/dyn1d/scm.F90	(revision 4325)
@@ -393,5 +393,5 @@
       call init_dimphy1D(1,llm)
       call suphel
-      call infotrac_init 
+      call init_infotrac
 
       if (nqtot>nqmx) STOP 'Augmenter nqmx dans lmdz1d.F'
Index: LMDZ6/trunk/libf/phylmd/infotrac_phy.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/infotrac_phy.F90	(revision 4319)
+++ LMDZ6/trunk/libf/phylmd/infotrac_phy.F90	(revision 4325)
@@ -1,10 +1,10 @@
-
-! $Id: $
-
+!$Id: infotrac.F90 4301 2022-10-20 11:57:21Z dcugnet $
+!
 MODULE infotrac_phy
 
-   USE       strings_mod, ONLY: msg, maxlen, strStack, strHead, strParse, strIdx, int2str
-   USE readTracFiles_mod, ONLY: trac_type, isot_type, keys_type, delPhase, getKey, tnom_iso => newH2OIso
-
+   USE       strings_mod, ONLY: msg, fmsg, maxlen, cat, dispTable, int2str, bool2str, strStack, strParse, strIdx
+   USE readTracFiles_mod, ONLY: trac_type, nphas, readTracersFiles, tracers, setGeneration, itZonIso, nbIso, tran0, delPhase, &
+                        getKey, isot_type, nzone, readIsotopesFile, isotope, maxTableWidth, iqIsoPha, ntiso, ixIso, addPhase, &
+                   indexUpdate, isoSelect, niso,  testTracersFiles, isoPhas, isoZone, isoName, isoKeys, iH2O, isoCheck
    IMPLICIT NONE
 
@@ -22,5 +22,5 @@
 
    !=== FOR ISOTOPES: General
-   PUBLIC :: isotopes,  nbIso                              !--- Derived type, full isotopes families database + nb of families
+   PUBLIC :: isot_type, nbIso                              !--- Derived type, full isotopes families database + nb of families
    PUBLIC :: isoSelect, ixIso                              !--- Isotopes family selection tool + selected family index
    !=== FOR ISOTOPES: Specific to water
@@ -30,11 +30,9 @@
    PUBLIC :: isoName, isoZone, isoPhas                     !--- Isotopes and tagging zones names, phases
    PUBLIC :: niso,    nzone,   nphas,   ntiso              !---  " " numbers + isotopes & tagging tracers number
-   PUBLIC :: itZonIso                                      !--- iq = function(tagging zone idx, isotope idx)
-   PUBLIC :: iqIsoPha                                      !--- idx of tagging tracer in iName = function(isotope idx, phase idx)
+   PUBLIC :: itZonIso                                      !--- idx "it" (in "isoName(1:niso)") = function(tagging idx, isotope idx)
+   PUBLIC :: iqIsoPha                                      !--- idx "iq" (in "qx") = function(isotope idx, phase idx) + aliases
    PUBLIC :: isoCheck                                      !--- Run isotopes checking routines
    !=== FOR BOTH TRACERS AND ISOTOPES
    PUBLIC :: getKey                                        !--- Get a key from "tracers" or "isotope"
-
-   INTERFACE isoSelect; MODULE PROCEDURE isoSelectByIndex, isoSelectByName; END INTERFACE isoSelect
 
 !=== CONVENTIONS FOR TRACERS NUMBERS:
@@ -72,12 +70,10 @@
 !  | phase       | Phases list ("g"as / "l"iquid / "s"olid)             | /           | [g][l][s]              |
 !  | component   | Name(s) of the merged/cumulated section(s)           | /           | coma-separated names   |
-!  | iadv        | Advection scheme number                              | iadv        | 1-20,30 exc. 3-9,15,19 |
 !  | iGeneration | Generation (>=1)                                     | /           |                        |
-!  | isAdvected  | advected tracers flag (.TRUE. if iadv >= 0)          | /           | nqtrue  .TRUE. values  |
 !  | isInPhysics | tracers not extracted from the main table in physics | /           | nqtottr .TRUE. values  |
 !  | iqParent    | Index of the parent tracer                           | iqpere      | 1:nqtot                |
 !  | iqDescen    | Indexes of the childs       (all generations)        | iqfils      | 1:nqtot                |
 !  | nqDescen    | Number of the descendants   (all generations)        | nqdesc      | 1:nqtot                |
-!  | nqChilds    | Number of childs            (1st generation only)    | nqfils      | 1:nqtot                |
+!  | nqChildren  | Number of childs            (1st generation only)    | nqfils      | 1:nqtot                |
 !  | iso_iGroup  | Isotopes group index in isotopes(:)                  | /           | 1:nbIso                |
 !  | iso_iName   | Isotope  name  index in isotopes(iso_iGroup)%trac(:) | iso_indnum  | 1:niso                 |
@@ -102,36 +98,16 @@
 
    !=== DIMENSIONS OF THE TRACERS TABLES AND OTHER SCALAR VARIABLES
-   INTEGER,                 SAVE :: nqtot,  &                   !--- Tracers nb in dynamics (incl. higher moments + H2O)
-                                    nbtr,   &                   !--- Tracers nb in physics  (excl. higher moments + H2O)
-                                    nqo,    &                   !--- Number of water phases
-                                    nbIso,  &                   !--- Number of available isotopes family
-                                    nqtottr, &                  !--- Number of tracers passed to phytrac (TO BE DELETED ?)
-                                    nqCO2                       !--- Number of tracers of CO2  (ThL)
-   CHARACTER(LEN=maxlen),   SAVE :: type_trac                   !--- Keyword for tracers type
-   CHARACTER(LEN=maxlen),   SAVE, ALLOCATABLE :: types_trac(:)  !--- Keyword for tracers type
-!$OMP THREADPRIVATE(nqtot, nbtr, nqo, nbIso, nqtottr, nqCO2, type_trac, types_trac)
-
-   !=== DERIVED TYPES EMBEDDING MOST INFORMATIONS ABOUT TRACERS AND ISOTOPES FAMILIES
-   TYPE(trac_type), TARGET, SAVE, ALLOCATABLE ::  tracers(:)    !=== TRACERS DESCRIPTORS VECTOR
-   TYPE(isot_type), TARGET, SAVE, ALLOCATABLE :: isotopes(:)    !=== ISOTOPES PARAMETERS VECTOR
-!$OMP THREADPRIVATE(tracers, isotopes)
-
-   !=== ALIASES FOR CURRENTLY SELECTED ISOTOPES FAMILY OF VARIABLES EMBEDDED IN THE VECTOR "isotopes"
-   TYPE(isot_type),         SAVE, POINTER :: isotope            !--- CURRENTLY SELECTED ISOTOPES FAMILY DESCRIPTOR
-   INTEGER,                 SAVE          :: ixIso, iH2O        !--- Index of the selected isotopes family and H2O family
-   LOGICAL,                 SAVE          :: isoCheck           !--- Flag to trigger the checking routines
-   TYPE(keys_type),         SAVE, POINTER :: isoKeys(:)         !--- ONE SET OF KEYS FOR EACH ISOTOPE (LISTED IN isoName)
-   CHARACTER(LEN=maxlen),   SAVE, POINTER :: isoName(:),   &    !--- ISOTOPES NAMES FOR THE CURRENTLY SELECTED FAMILY
-                                             isoZone(:),   &    !--- TAGGING ZONES  FOR THE CURRENTLY SELECTED FAMILY
-                                             isoPhas            !--- USED PHASES    FOR THE CURRENTLY SELECTED FAMILY
-   INTEGER,                 SAVE          ::  niso, nzone, &    !--- NUMBER OF ISOTOPES, TAGGING ZONES AND PHASES
-                                             nphas, ntiso       !--- NUMBER OF PHASES AND ISOTOPES + ISOTOPIC TAGGING TRACERS
-   INTEGER,                 SAVE, POINTER ::itZonIso(:,:), &    !--- INDEX IN "isoTrac" AS f(tagging zone idx,  isotope idx)
-                                            iqIsoPha(:,:)       !--- INDEX IN "qx"      AS f(isotopic tracer idx, phase idx)
-!$OMP THREADPRIVATE(isotope, ixIso,iH2O, isoCheck, isoKeys, isoName,isoZone,isoPhas, niso,nzone,nphas,ntiso, itZonIso,iqIsoPha)
-
-   !=== VARIABLES FOR ISOTOPES INITIALIZATION AND FOR INCA
-   INTEGER,          SAVE,    ALLOCATABLE ::conv_flg(:),  &     !--- Convection     activation ; needed for INCA        (nbtr)
-                                             pbl_flg(:)         !--- Boundary layer activation ; needed for INCA        (nbtr)
+   INTEGER,               SAVE :: nqtot,  &                     !--- Tracers nb in dynamics (incl. higher moments + H2O)
+                                  nbtr,   &                     !--- Tracers nb in physics  (excl. higher moments + H2O)
+                                  nqo,    &                     !--- Number of water phases
+                                  nqtottr, &                    !--- Number of tracers passed to phytrac (TO BE DELETED ?)
+                                  nqCO2                         !--- Number of tracers of CO2  (ThL)
+   CHARACTER(LEN=maxlen), SAVE :: type_trac                     !--- Keyword for tracers type(s)
+   CHARACTER(LEN=maxlen), SAVE, ALLOCATABLE :: types_trac(:)    !--- Keyword for tracers type(s), parsed version
+!$OMP THREADPRIVATE(nqtot, nbtr, nqo, nqtottr, nqCO2, type_trac, types_trac)
+
+   !=== VARIABLES FOR INCA
+   INTEGER,               SAVE, ALLOCATABLE :: conv_flg(:), &   !--- Convection     activation ; needed for INCA        (nbtr)
+                                                pbl_flg(:)      !--- Boundary layer activation ; needed for INCA        (nbtr)
 !$OMP THREADPRIVATE(conv_flg, pbl_flg)
 
@@ -146,49 +122,313 @@
 CONTAINS
 
-SUBROUTINE init_infotrac_phy(type_trac_, tracers_, isotopes_, nqtottr_, nqCO2_, pbl_flg_, conv_flg_)
-
-   USE print_control_mod, ONLY: prt_level, lunout
-
+SUBROUTINE init_infotrac_phy
+   USE control_mod, ONLY: planet_type, config_inca
+   USE ioipsl_getin_p_mod, ONLY: getin_p
+#ifdef REPROBUS
+   USE CHEM_REP,    ONLY: Init_chem_rep_trac
+#endif
    IMPLICIT NONE
-   CHARACTER(LEN=*),INTENT(IN) :: type_trac_
-   TYPE(trac_type), INTENT(IN) ::  tracers_(:)
-   TYPE(isot_type), INTENT(IN) :: isotopes_(:)
-   INTEGER,         INTENT(IN) :: nqtottr_
-   INTEGER,         INTENT(IN) :: nqCO2_
-   INTEGER,         INTENT(IN) :: conv_flg_(:)
-   INTEGER,         INTENT(IN) ::  pbl_flg_(:)
-
-   INTEGER :: iq, ixt
+!==============================================================================================================================
+!
+!   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
+!   -------
+!
+!   Modifications:
+!   --------------
+!   05/94: F.Forget      Modif special traceur
+!   02/02: M-A Filiberti Lecture de traceur.def
+!   01/22: D. Cugnet     Nouveaux tracer.def et tracer_*.def + encapsulation (types trac_type et isot_type)
+!
+!   Objet:
+!   ------
+!   GCM LMD nouvelle grille
+!
+!==============================================================================================================================
+!   ... modification de l'integration de q ( 26/04/94 ) ....
+!------------------------------------------------------------------------------------------------------------------------------
+! Declarations:
+   INCLUDE "dimensions.h"
+   INCLUDE "iniprint.h"
+
+!------------------------------------------------------------------------------------------------------------------------------
+! Local variables
+   INTEGER, ALLOCATABLE :: hadv(:), vadv(:)                          !--- Horizontal/vertical transport scheme number
+#ifdef INCA
+   INTEGER, ALLOCATABLE :: had (:), hadv_inca(:), conv_flg_inca(:), &!--- Variables specific to INCA
+                           vad (:), vadv_inca(:),  pbl_flg_inca(:)
+   CHARACTER(LEN=8), ALLOCATABLE :: solsym_inca(:)                   !--- Tracers names for INCA
+   INTEGER :: nqINCA
+#endif
 #ifdef CPP_StratAer
    CHARACTER(LEN=maxlen), ALLOCATABLE :: tnames(:)
 #endif
-   CHARACTER(LEN=maxlen) :: modname="init_infotrac_phy"
-
-   type_trac = type_trac_
-   IF(strParse(type_trac, '|', types_trac)) CALL abort_physic(modname,'can''t parse "type_trac = '//TRIM(type_trac)//'"',1)
-   tracers   = tracers_
-   isotopes  = isotopes_
-   nqtottr   = nqtottr_
-   nqCO2     = nqCO2_
-   pbl_flg   =  pbl_flg_
-   conv_flg  = conv_flg_
-   nqtot     = SIZE(tracers_)
-   nqo       = COUNT(delPhase(tracers%name)=='H2O' .AND. tracers%iGeneration==0 .AND. tracers%component=='lmdz')
-   nbtr      = SIZE(conv_flg)
-   nbIso     = SIZE(isotopes_)
-
-   !=== Determine selected isotopes class related quantities:
-   !    ixIso, isotope, niso,isoKeys, ntiso,isoName, nzone,isoZone, nphas,isoPhas, itZonIso, iqIsoPha, isoCheck
-   IF(.NOT.isoSelect('H2O')) iH2O = ixIso
-   IF(prt_level > 1) THEN
-      CALL msg('nqtot   = '//TRIM(int2str(nqtot)),   modname)
-      CALL msg('nbtr    = '//TRIM(int2str(nbtr )),   modname)
-      CALL msg('nqo     = '//TRIM(int2str(nqo  )),   modname)
-      CALL msg('niso    = '//TRIM(int2str(niso )),   modname)
-      CALL msg('ntiso   = '//TRIM(int2str(ntiso)),   modname)
-      CALL msg('nqtottr = '//TRIM(int2str(nqtottr)), modname)
-      CALL msg('nqCO2   = '//TRIM(int2str(nqCO2)),   modname)
-   END IF
-
+   CHARACTER(LEN=2)      ::   suff(9)                                !--- Suffixes for schemes of order 3 or 4 (Prather)
+   CHARACTER(LEN=3)      :: descrq(30)                               !--- Advection scheme description tags
+   CHARACTER(LEN=maxlen) :: msg1                                     !--- String for messages
+   INTEGER :: fType                                                  !--- Tracers description file type ; 0: none
+                                                                     !--- 1/2/3: "traceur.def"/"tracer.def"/"tracer_*.def"
+   INTEGER :: nqtrue                                                 !--- Tracers nb from tracer.def (no higher order moments)
+   INTEGER :: iad                                                    !--- Advection scheme number
+   INTEGER :: ic, iq, jq, it, nt, im, nm, iz, k                      !--- Indexes and temporary variables
+   LOGICAL :: lerr, ll, lRepr, lInit
+   CHARACTER(LEN=1) :: p
+   TYPE(trac_type), ALLOCATABLE, TARGET :: ttr(:)
+   TYPE(trac_type), POINTER             :: t1, t(:)
+   INTEGER :: ierr
+
+   CHARACTER(LEN=*), PARAMETER :: modname="init_infotrac_phy"
+!------------------------------------------------------------------------------------------------------------------------------
+! Initialization :
+!------------------------------------------------------------------------------------------------------------------------------
+   suff          = ['x ','y ','z ','xx','xy','xz','yy','yz','zz']
+   descrq( 1: 2) = ['LMV','BAK']
+   descrq(10:20) = ['VL1','VLP','FH1','FH2','VLH','   ','PPM','PPS','PPP','   ','SLP']
+   descrq(30)    =  'PRA'
+
+   CALL getin_p('type_trac',type_trac)
+   CALL msg('type_trac = "'//TRIM(type_trac)//'"', modname)
+   IF(strParse(type_trac, '|', types_trac, n=nt)) CALL abort_gcm(modname,'can''t parse "type_trac = '//TRIM(type_trac)//'"',1)
+   lInit = .NOT.ALLOCATED(tracers)
+
+!##############################################################################################################################
+   IF(lInit) THEN                                                    !=== SKIPED IF ALREADY DONE IN dyn3d_common/infotrac  ####
+!##############################################################################################################################
+   !---------------------------------------------------------------------------------------------------------------------------
+   DO it = 1, nt                                                     !--- nt>1=> "type_trac": coma-separated keywords list
+   !---------------------------------------------------------------------------------------------------------------------------
+      !--- MESSAGE ABOUT THE CHOSEN CONFIGURATION
+      msg1 = 'For type_trac = "'//TRIM(types_trac(it))//'":'
+      SELECT CASE(types_trac(it))
+         CASE('inca'); CALL msg(TRIM(msg1)//' coupling with INCA chemistry model, config_inca='//config_inca, modname)
+         CASE('inco'); CALL msg(TRIM(msg1)//' coupling jointly with INCA and CO2 cycle',  modname)
+         CASE('repr'); CALL msg(TRIM(msg1)//' coupling with REPROBUS chemistry model',    modname)
+         CASE('co2i'); CALL msg(TRIM(msg1)//' you have chosen to run with CO2 cycle',     modname)
+         CASE('coag'); CALL msg(TRIM(msg1)//' tracers are treated for COAGULATION tests', modname)
+         CASE('lmdz'); CALL msg(TRIM(msg1)//' tracers are treated in LMDZ only',          modname)
+         CASE DEFAULT; CALL abort_gcm(modname,'type_trac='//TRIM(types_trac(it))//' not possible yet.',1)
+      END SELECT
+
+      !--- COHERENCE TEST BETWEEN "type_trac" AND "config_inca"
+      IF(ANY(['inca', 'inco'] == types_trac(it)) .AND. ALL(['aero', 'aeNP', 'chem'] /= config_inca)) &
+         CALL abort_gcm(modname, 'Incoherence between type_trac and config_inca. Please modify "run.def"',1)
+
+      !--- COHERENCE TEST BETWEEN "type_trac" AND PREPROCESSING KEYS
+      SELECT CASE(types_trac(it))
+         CASE('inca', 'inco')
+#ifndef INCA
+            CALL abort_gcm(modname, 'You must add cpp key INCA and compile with INCA code', 1)
+#endif
+         CASE('repr')
+#ifndef REPROBUS
+            CALL abort_gcm(modname, 'You must add cpp key REPROBUS and compile with REPROBUS code', 1)
+#endif
+         CASE('coag')
+#ifndef CPP_StratAer
+            CALL abort_gcm(modname, 'You must add cpp key StratAer and compile with StratAer code', 1)
+#endif
+      END SELECT
+
+   !---------------------------------------------------------------------------------------------------------------------------
+   END DO
+   !---------------------------------------------------------------------------------------------------------------------------
+
+!##############################################################################################################################
+   END IF
+!##############################################################################################################################
+
+   !--- DISABLE "config_inca" OPTION FOR A RUN WITHOUT "INCA" IF IT DIFFERS FROM "none"
+   IF(fmsg('Setting config_inca="none" as you do not couple with INCA model', &
+         modname, ALL(types_trac /= 'inca') .AND. ALL(types_trac /= 'inco') .AND. config_inca /= 'none')) config_inca = 'none'
+
+   nqCO2 = COUNT( [ANY(types_trac == 'inco') .OR. (ANY(types_trac == 'co2i') .AND. ANY(types_trac == 'inca'))] )
+
+!==============================================================================================================================
+! 1) Get the numbers of: true (first order only) tracers "nqtrue", water tracers "nqo" (vapor/liquid/solid)
+!==============================================================================================================================
+   lRepr = ANY(types_trac(:) == 'repr')
+!##############################################################################################################################
+   IF(lInit) THEN
+     IF(readTracersFiles(type_trac,  fType,  lRepr)) CALL abort_gcm(modname, 'problem with tracers file(s)',1)
+   ELSE
+     IF(testTracersFiles(modname, type_trac, fType)) CALL abort_gcm(modname, 'problem with tracers file(s)',1)
+   END IF
+!##############################################################################################################################
+
+   !---------------------------------------------------------------------------------------------------------------------------
+   IF(fType == 0) CALL abort_gcm(modname, 'Missing tracers file: "traceur.def", "tracer.def" or "tracer_<keyword>.def file.',1)
+   !---------------------------------------------------------------------------------------------------------------------------
+   IF(fType == 1 .AND. ANY(['inca','inco']==type_trac) .AND. lInit) THEN  !=== OLD STYLE INCA "traceur.def" (single type_trac)
+   !---------------------------------------------------------------------------------------------------------------------------
+#ifdef INCA
+      nqo = SIZE(tracers) - nqCO2
+      CALL Init_chem_inca_trac(nqINCA)                               !--- Get nqINCA from INCA
+      nbtr = nqINCA + nqCO2                                          !--- Number of tracers passed to phytrac
+      nqtrue = nbtr + nqo                                            !--- Total number of "true" tracers
+      IF(ALL([2,3] /= nqo)) CALL abort_gcm(modname, 'Only 2 or 3 water phases allowed ; found nqo='//TRIM(int2str(nqo)), 1)
+      ALLOCATE(hadv(nqtrue), hadv_inca(nqINCA), conv_flg_inca(nqINCA), solsym_inca(nqINCA))
+      ALLOCATE(vadv(nqtrue), vadv_inca(nqINCA),  pbl_flg_inca(nqINCA))
+      CALL init_transport(solsym_inca, conv_flg_inca, pbl_flg_inca, hadv_inca, vadv_inca)
+      ALLOCATE(ttr(nqtrue))
+      ttr(1:nqo+nqCO2)                    = tracers
+      ttr(1    :      nqo   )%component   = 'lmdz'
+      ttr(1+nqo:nqCO2+nqo   )%component   = 'co2i'
+      ttr(1+nqo+nqCO2:nqtrue)%component   = 'inca'
+      ttr(1+nqo      :nqtrue)%name        = [('CO2     ', k=1, nqCO2), solsym_inca]
+      ttr(1+nqo+nqCO2:nqtrue)%parent      = tran0
+      ttr(1+nqo+nqCO2:nqtrue)%phase       = 'g'
+      lerr = getKey('hadv', had, ky=tracers(:)%keys)
+      lerr = getKey('vadv', vad, ky=tracers(:)%keys)
+      hadv(1:nqo) = had(:); hadv(nqo+1:nqtrue) = hadv_inca
+      vadv(1:nqo) = vad(:); vadv(nqo+1:nqtrue) = vadv_inca
+      CALL MOVE_ALLOC(FROM=ttr, TO=tracers)
+      CALL setGeneration(tracers)                                    !--- SET FIELDS %iGeneration, %gen0Name
+      DEALLOCATE(had, hadv_inca, vad, vadv_inca, conv_flg_inca, pbl_flg_inca, solsym_inca)
+#endif
+   !---------------------------------------------------------------------------------------------------------------------------
+   ELSE                                                              !=== OTHER CASES (OLD OR NEW FORMAT, NO INCA MODULE)
+   !---------------------------------------------------------------------------------------------------------------------------
+      nqo    =        COUNT(delPhase(tracers(:)%name)     == 'H2O' &
+                               .AND. tracers(:)%component == 'lmdz') !--- Number of water phases
+      nqtrue = SIZE(tracers)                                         !--- Total number of "true" tracers
+      nbtr   = nqtrue-COUNT(delPhase(tracers(:)%gen0Name) == 'H2O' &
+                               .AND. tracers(:)%component == 'lmdz') !--- Number of tracers passed to phytrac
+#ifdef INCA
+      nqINCA = COUNT(tracers(:)%component == 'inca')
+#endif
+      lerr = getKey('hadv', hadv, ky=tracers(:)%keys)
+      lerr = getKey('vadv', vadv, ky=tracers(:)%keys)
+   !---------------------------------------------------------------------------------------------------------------------------
+   END IF
+   !---------------------------------------------------------------------------------------------------------------------------
+
+   !--- Transfert the number of tracers to Reprobus
+#ifdef REPROBUS
+   CALL Init_chem_rep_trac(nbtr, nqo, tracers(:)%name)
+#endif
+
+!##############################################################################################################################
+   IF(lInit) THEN                                                    !=== SKIPED IF ALREADY DONE IN dyn3d_common/infotrac  ####
+!##############################################################################################################################
+
+!==============================================================================================================================
+! 2) Calculate nqtot, number of tracers needed (greater if advection schemes 20 or 30 have been chosen).
+!==============================================================================================================================
+   DO iq = 1, nqtrue
+      IF( hadv(iq)<20 .OR. (ANY(hadv(iq)==[20,30]) .AND. hadv(iq)==vadv(iq)) ) CYCLE
+      WRITE(msg1,'("The choice hadv=",i0,", vadv=",i0,a)')hadv(iq),vadv(iq),' for "'//TRIM(tracers(iq)%name)//'" is not available'
+      CALL abort_gcm(modname, TRIM(msg1), 1)
+   END DO
+   nqtot =    COUNT( hadv< 20 .AND. vadv< 20 ) &                     !--- No additional tracer
+         +  4*COUNT( hadv==20 .AND. vadv==20 ) &                     !--- 3  additional tracers
+         + 10*COUNT( hadv==30 .AND. vadv==30 )                       !--- 9  additional tracers
+
+   !--- More tracers due to the choice of advection scheme => assign total number of tracers
+   IF( nqtot /= nqtrue ) THEN
+      CALL msg('The choice of advection scheme for one or more tracers makes it necessary to add tracers')
+      CALL msg('The number of true tracers is '//TRIM(int2str(nqtrue)))
+      CALL msg('The total number of tracers needed is '//TRIM(int2str(nqtot)))
+   END IF
+
+!==============================================================================================================================
+! 3) Determine the advection scheme ; needed to compute the full tracers list, the long names and nqtot.
+!==============================================================================================================================
+   ALLOCATE(ttr(nqtot))
+   jq = nqtrue+1; tracers(:)%iadv = -1
+   DO iq = 1, nqtrue
+      t1 => tracers(iq)
+
+      !--- VERIFY THE CHOICE OF ADVECTION SCHEME
+      iad = -1
+      IF(hadv(iq)     ==    vadv(iq)    ) iad = hadv(iq)
+      IF(hadv(iq)==10 .AND. vadv(iq)==16) iad = 11
+      WRITE(msg1,'("Bad choice of advection scheme for ",a,": hadv = ",i0,", vadv = ",i0)')TRIM(t1%name), hadv(iq), vadv(iq)
+      IF(iad == -1) CALL abort_gcm(modname, msg1, 1)
+
+      !--- SET FIELDS %longName, %isInPhysics
+      t1%longName   = t1%name; IF(iad > 0) t1%longName=TRIM(t1%name)//descrq(iad)
+      t1%isInPhysics= delPhase(t1%gen0Name) /= 'H2O' &
+                          .OR. t1%component /= 'lmdz' !=== OTHER EXCEPTIONS TO BE ADDED: CO2i, SURSATURATED WATER CLOUD...
+      ttr(iq)       = t1
+
+      !--- DEFINE THE HIGHER ORDER TRACERS, IF ANY
+      nm = 0
+      IF(iad == 20) nm = 3                                           !--- 2nd order scheme
+      IF(iad == 30) nm = 9                                           !--- 3rd order scheme
+      IF(nm == 0) CYCLE                                              !--- No higher moments
+      ttr(jq+1:jq+nm)             = t1
+      ttr(jq+1:jq+nm)%name        = [(TRIM(t1%name)    //'-'//TRIM(suff(im)), im=1, nm) ]
+      ttr(jq+1:jq+nm)%parent      = [(TRIM(t1%parent)  //'-'//TRIM(suff(im)), im=1, nm) ]
+      ttr(jq+1:jq+nm)%longName    = [(TRIM(t1%longName)//'-'//TRIM(suff(im)), im=1, nm) ]
+      jq = jq + nm
+   END DO
+   DEALLOCATE(hadv, vadv)
+   CALL MOVE_ALLOC(FROM=ttr, TO=tracers)
+
+   !--- SET FIELDS %iqParent, %nqChildren, %iGeneration, %iqDescen, %nqDescen
+   CALL indexUpdate(tracers)
+
+!##############################################################################################################################
+   END IF
+!##############################################################################################################################
+
+!##############################################################################################################################
+   IF(.NOT.lInit) THEN
+!##############################################################################################################################
+     nqtot = SIZE(tracers)
+!##############################################################################################################################
+   ELSE
+!##############################################################################################################################
+
+   !=== READ PHYSICAL PARAMETERS FOR ISOTOPES
+   niso = 0; nzone = 0; nphas = nqo; ntiso = 0; isoCheck = .FALSE.
+   IF(readIsotopesFile()) CALL abort_gcm(modname, 'Problem when reading isotopes parameters', 1)
+
+!##############################################################################################################################
+   END IF
+!##############################################################################################################################
+   !--- Convection / boundary layer activation for all tracers
+   ALLOCATE(conv_flg(nbtr)); conv_flg(1:nbtr) = 1
+   ALLOCATE( pbl_flg(nbtr));  pbl_flg(1:nbtr) = 1
+
+   !--- Note: nqtottr can differ from nbtr when nmom/=0
+   nqtottr = nqtot - COUNT(delPhase(tracers%gen0Name) == 'H2O' .AND. tracers%component == 'lmdz')
+   IF(COUNT(tracers%iso_iName == 0) - COUNT(delPhase(tracers%name) == 'H2O' .AND. tracers%component == 'lmdz') /= nqtottr) &
+      CALL abort_gcm(modname, 'pb dans le calcul de nqtottr', 1)
+
+   !=== DISPLAY THE RESULTS
+!   IF(prt_level > 1) THEN
+      CALL msg('nqo    = '//TRIM(int2str(nqo)),    modname)
+      CALL msg('nbtr   = '//TRIM(int2str(nbtr)),   modname)
+      CALL msg('nqtrue = '//TRIM(int2str(nqtrue)), modname)
+      CALL msg('nqtot  = '//TRIM(int2str(nqtot)),  modname)
+      CALL msg('niso   = '//TRIM(int2str(niso)),   modname)
+      CALL msg('ntiso  = '//TRIM(int2str(ntiso)),  modname)
+#ifdef INCA
+      CALL msg('nqCO2  = '//TRIM(int2str(nqCO2)),  modname)
+      CALL msg('nqINCA = '//TRIM(int2str(nqINCA)), modname)
+#endif
+!   END IF
+   t => tracers
+   CALL msg('Information stored in infotrac_phy :', modname)
+   IF(dispTable('issssssssiiiiiiii', &
+      ['iq    ', 'name  ', 'lName ', 'gen0N ', 'parent', 'type  ', 'phase ', 'compon', 'isPhy ',           &
+                 'iGen  ', 'iqPar ', 'nqDes ', 'nqChld', 'iGroup', 'iName ', 'iZone ', 'iPhase'],          &
+      cat(t%name, t%longName, t%gen0Name, t%parent, t%type, t%phase, t%component, bool2str(t%isInPhysics)),&
+      cat([(iq, iq=1, nqtot)], t%iGeneration, t%iqParent, t%nqDescen, t%nqChildren, t%iso_iGroup,          &
+                  t%iso_iName, t%iso_iZone, t%iso_iPhase), nColMax=maxTableWidth, nHead=2, sub=modname))   &
+      CALL abort_gcm(modname, "problem with the tracers table content", 1)
+   IF(niso > 0) THEN
+      CALL msg('Where, for isotopes family "'//TRIM(isotope%parent)//'":', modname)
+      CALL msg('  isoKeys%name = '//strStack(isoKeys%name), modname)
+      CALL msg('  isoName = '//strStack(isoName),      modname)
+      CALL msg('  isoZone = '//strStack(isoZone),      modname)
+      CALL msg('  isoPhas = '//TRIM(isoPhas),          modname)
+   ELSE
+      CALL msg('No isotopes identified.', modname)
+   END IF
+
+#ifdef ISOVERIF
+   CALL msg('iso_iName = '//strStack(int2str(PACK(tracers(:)%iso_iName, MASK=tracers(:)%iso_iGroup==iH2O))), modname)
+#endif
 #ifdef CPP_StratAer
    IF (ANY(types_trac == 'coag')) THEN
@@ -210,55 +450,7 @@
    END IF
 #endif
-#ifdef ISOVERIF
-   CALL msg('iso_iName = '//strStack(int2str(PACK(tracers(:)%iso_iName, MASK=tracers(:)%iso_iGroup==iH2O))), modname)
-#endif
+   CALL msg('end', modname)
 
 END SUBROUTINE init_infotrac_phy
 
-
-!==============================================================================================================================
-!=== THE ROUTINE isoSelect IS USED TO SWITCH FROM AN ISOTOPE FAMILY TO ANOTHER: ISOTOPES DEPENDENT PARAMETERS ARE UPDATED
-!     Single generic "isoSelect" routine, using the predefined index of the parent (fast version) or its name (first call).
-!==============================================================================================================================
-LOGICAL FUNCTION isoSelectByName(iName, lVerbose) RESULT(lerr)
-   IMPLICIT NONE
-   CHARACTER(LEN=*),  INTENT(IN)  :: iName
-   LOGICAL, OPTIONAL, INTENT(IN) :: lVerbose
-   INTEGER :: iIso
-   LOGICAL :: lV
-   lV = .FALSE.; IF(PRESENT(lVerbose)) lV = lVerbose
-   iIso = strIdx(isotopes(:)%parent, iName)
-   lerr = iIso == 0
-   IF(lerr) THEN
-      niso = 0; ntiso = 0; nzone=0; nphas=nqo; isoCheck=.FALSE.
-      CALL msg('no isotope family named "'//TRIM(iName)//'"', ll=lV)
-      RETURN
-   END IF
-   lerr = isoSelectByIndex(iIso, lV)
-END FUNCTION isoSelectByName
-!==============================================================================================================================
-LOGICAL FUNCTION isoSelectByIndex(iIso, lVerbose) RESULT(lerr)
-   IMPLICIT NONE
-   INTEGER,           INTENT(IN) :: iIso
-   LOGICAL, OPTIONAL, INTENT(IN) :: lVerbose
-   LOGICAL :: lv
-   lv = .FALSE.; IF(PRESENT(lVerbose)) lv = lVerbose
-   lerr = .FALSE.
-   IF(iIso == ixIso) RETURN                                      !--- Nothing to do if the index is already OK
-   lerr = iIso<=0 .OR. iIso>nbIso
-   CALL msg('Inconsistent isotopes family index '//TRIM(int2str(iIso))//': should be > 0 and <= '//TRIM(int2str(nbIso))//'"',&
-            ll=lerr .AND. lV)
-   IF(lerr) RETURN
-   ixIso = iIso                                                  !--- Update currently selected family index
-   isotope  => isotopes(ixIso)                                   !--- Select corresponding component
-   isoKeys  => isotope%keys;     niso     = isotope%niso
-   isoName  => isotope%trac;     ntiso    = isotope%ntiso
-   isoZone  => isotope%zone;     nzone    = isotope%nzone
-   isoPhas  => isotope%phase;    nphas    = isotope%nphas
-   itZonIso => isotope%itZonIso; isoCheck = isotope%check
-   iqIsoPha => isotope%iqIsoPha
-END FUNCTION isoSelectByIndex
-!==============================================================================================================================
-
-
 END MODULE infotrac_phy
Index: LMDZ6/trunk/libf/phylmd/traclmdz_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/traclmdz_mod.F90	(revision 4319)
+++ LMDZ6/trunk/libf/phylmd/traclmdz_mod.F90	(revision 4325)
@@ -175,5 +175,5 @@
     it = 0
     DO iq = 1, nqtot
-       IF(.NOT.(tracers(iq)%isAdvected .AND. tracers(iq)%isInPhysics)) CYCLE
+       IF(.NOT.(tracers(iq)%isInPhysics)) CYCLE
        it = it+1
        SELECT CASE(strLower(tracers(iq)%name))
Index: LMDZ6/trunk/libf/phylmdiso/isotopes_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/isotopes_mod.F90	(revision 4319)
+++ LMDZ6/trunk/libf/phylmdiso/isotopes_mod.F90	(revision 4325)
@@ -133,5 +133,4 @@
 
 SUBROUTINE iso_init()
-   USE ioipsl_getin_p_mod, ONLY: getin_p
    USE infotrac_phy,       ONLY: ntiso, niso, getKey
     USE strings_mod,       ONLY: maxlen
@@ -181,10 +180,5 @@
    iso_HTO = strIdx(isoName, 'H[3]HO');  CALL msg('iso_HTO='//int2str(iso_HTO), modname)
 
-   ! initialisation
-   ! lecture des parametres isotopiques:
-   ! pour que ca marche en openMP, il faut utiliser getin_p. Car le getin ne peut
-   ! etre appele que par un thread a la fois, et ca pose tout un tas de problemes,
-   ! d'ou tout un tas de magouilles comme dans conf_phys_m. A terme, tout le monde
-   ! lira par getin_p.
+   !--- Initialiaation: reading the isotopic parameters.
    CALL get_in('lambda',     lambda_sursat, 0.004)
    CALL get_in('thumxt1',    thumxt1,       0.75*1.2)
@@ -339,16 +333,16 @@
    USE mod_phys_lmdz_omp_data, ONLY :  is_omp_root
    USE mod_phys_lmdz_transfert_para, ONLY : bcast
-   CHARACTER(LEN=*),  INTENT(IN)    :: nam
-   CHARACTER(LEN=*),  INTENT(INOUT) :: val
-   CHARACTER(LEN=*),  INTENT(IN)    :: def
-   LOGICAL, OPTIONAL, INTENT(IN)    :: lDisp
+   CHARACTER(LEN=*),           INTENT(IN)    :: nam
+   CHARACTER(LEN=*),           INTENT(INOUT) :: val
+   CHARACTER(LEN=*), OPTIONAL, INTENT(IN)    :: def
+   LOGICAL,          OPTIONAL, INTENT(IN)    :: lDisp
    LOGICAL :: lD
 !$OMP BARRIER
    IF(is_mpi_root.AND.is_omp_root) THEN
-      val=def; CALL getin(nam,val)
+      IF(PRESENT(def)) val=def; CALL getin(nam,val)
       lD=.TRUE.; IF(PRESENT(lDisp)) lD=lDisp
       IF(lD) CALL msg(TRIM(nam)//' = '//TRIM(val))
-   END IF
-   CALL bcast(val)
+  END IF
+  CALL bcast(val)
 END SUBROUTINE getinp_s
 
@@ -360,14 +354,14 @@
    CHARACTER(LEN=*),  INTENT(IN)    :: nam
    INTEGER,           INTENT(INOUT) :: val
-   INTEGER,           INTENT(IN)    :: def
+   INTEGER, OPTIONAL, INTENT(IN)    :: def
    LOGICAL, OPTIONAL, INTENT(IN)    :: lDisp
    LOGICAL :: lD
 !$OMP BARRIER
    IF(is_mpi_root.AND.is_omp_root) THEN
-      val=def; CALL getin(nam,val)
+      IF(PRESENT(def)) val=def; CALL getin(nam,val)
       lD=.TRUE.; IF(PRESENT(lDisp)) lD=lDisp
       IF(lD) CALL msg(TRIM(nam)//' = '//TRIM(int2str(val)))
-   END IF
-   CALL bcast(val)
+  END IF
+  CALL bcast(val)
 END SUBROUTINE getinp_i
 
@@ -379,14 +373,14 @@
    CHARACTER(LEN=*),  INTENT(IN)    :: nam
    REAL,              INTENT(INOUT) :: val
-   REAL,              INTENT(IN)    :: def
+   REAL,    OPTIONAL, INTENT(IN)    :: def
    LOGICAL, OPTIONAL, INTENT(IN)    :: lDisp
    LOGICAL :: lD
 !$OMP BARRIER
    IF(is_mpi_root.AND.is_omp_root) THEN
-      Val=def; CALL getin(nam,val)
+      IF(PRESENT(def)) val=def; CALL getin(nam,val)
       lD=.TRUE.; IF(PRESENT(lDisp)) lD=lDisp
       IF(lD) CALL msg(TRIM(nam)//' = '//TRIM(real2str(val)))
-   ENDIF
-   CALL bcast(val)
+  END IF
+  CALL bcast(val)
 END SUBROUTINE getinp_r
 
@@ -398,14 +392,14 @@
    CHARACTER(LEN=*),  INTENT(IN)    :: nam
    LOGICAL,           INTENT(INOUT) :: val
-   LOGICAL,           INTENT(IN)    :: def
+   LOGICAL, OPTIONAL, INTENT(IN)    :: def
    LOGICAL, OPTIONAL, INTENT(IN)    :: lDisp
    LOGICAL :: lD
 !$OMP BARRIER
    IF(is_mpi_root.AND.is_omp_root) THEN
-      val=def; CALL getin(nam,val)
+      IF(PRESENT(def)) val=def; CALL getin(nam,val)
       lD=.TRUE.; IF(PRESENT(lDisp)) lD=lDisp
       IF(lD) CALL msg(TRIM(nam)//' = '//TRIM(bool2str(val)))
-   END IF
-   CALL bcast(val)
+  END IF
+  CALL bcast(val)
 END SUBROUTINE getinp_l
 
Index: LMDZ6/trunk/libf/phylmdiso/isotopes_routines_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/isotopes_routines_mod.F90	(revision 4319)
+++ LMDZ6/trunk/libf/phylmdiso/isotopes_routines_mod.F90	(revision 4325)
@@ -16532,4 +16532,7 @@
  USE isotrac_mod, ONLY: strtrac,initialisation_isotrac,index_iso, &
 &       index_zone,izone_init
+ USE readTracFiles_mod, ONLY: newH2Oiso, oldH2Oiso
+ USE strings_mod, ONLY: strIdx, strHead, strTail
+
 #endif
         implicit none
@@ -16563,5 +16566,5 @@
         CHARACTER*5 str5
         real xmin,xmax    
-        CHARACTER*50 outiso  
+        CHARACTER*50 outiso, oldIso
         integer lnblnk
         LOGICAL :: found,phyetat0_get,phyetat0_srf
@@ -16585,9 +16588,11 @@
    do ixt=1,ntraciso
 
-      outiso=TRIM(isoName(ixt))
-      i = INDEX(outiso, '_', .TRUE.)
-      outiso = outiso(1:i-1)//outiso(i+1:LEN_TRIM(outiso))
-      write(*,*) 'phyiso_etat0_fichier 16621: ixt,outiso=',ixt,TRIM(outiso)
-
+      outiso = isoName(ixt)
+      k = strIdx(newH2Oiso, strHead(outiso, '_'))
+      oldIso = outiso; IF(k /= 0) oldIso = oldH2Oiso(k)
+      IF(INDEX(outiso, '_') /= 0) THEN
+        outiso = TRIM(outiso)//TRIM(strTail(outiso, '_'))
+        oldIso = TRIM(oldIso)//TRIM(strTail(outiso, '_'))
+      END IF
            
       ! on lit seulement si ixt<=niso ou si on initialise les traceurs d'après
@@ -16598,30 +16603,38 @@
 
       found=phyetat0_srf(1,iso_tmp_lonsrf,"XTSNOW"//TRIM(outiso),"Surface snow",0.)
-      if (.NOT.found) CALL abort_physic('isotopes_routines_mod', &
-                            'phyiso_etat0_fichier 16581: variable isotopique not found',1)
+      if (.NOT.found.AND.k/=0) found=phyetat0_srf(1,iso_tmp_lonsrf,"XTSNOW"//TRIM(oldIso),"Surface snow",0.)
+      if (.NOT.found) CALL abort_physic('isotopes_routines_mod', 'phyiso_etat0_fichier 16581: variable isotopique not found',1)
       xtsnow(ixt,:,:)=iso_tmp_lonsrf(:,:)
       
       found=phyetat0_srf(1,iso_tmp_lonsrf,"XTEVAP"//TRIM(outiso),"evaporation",0.)
+      if (.NOT.found.AND.k/=0) found=phyetat0_srf(1,iso_tmp_lonsrf,"XTEVAP"//TRIM(oldIso),"evaporation",0.)
       fxtevap(ixt,:,:)=iso_tmp_lonsrf(:,:)
 
       found=phyetat0_get(1,iso_tmp,"xtrain_f"//TRIM(outiso),"xrain fall",0.)
+      if (.NOT.found.AND.k/=0) found=phyetat0_get(1,iso_tmp,"xtrain_f"//TRIM(oldIso),"xrain fall",0.)
       xtrain_fall(ixt,:)=iso_tmp(:)
 
       found=phyetat0_get(1,iso_tmp,"xtsnow_f"//TRIM(outiso),"snow fall",0.)
+      if (.NOT.found.AND.k/=0) found=phyetat0_get(1,iso_tmp,"xtsnow_f"//TRIM(oldIso),"snow fall",0.)
       xtsnow_fall(ixt,:)=iso_tmp(:)
 
       found=phyetat0_get(klev,iso_tmp_lonlev,"XTANCIEN"//TRIM(outiso),"QANCIEN",0.)
+      if (.NOT.found.AND.k/=0) found=phyetat0_get(1,iso_tmp_lonlev,"XTANCIEN"//TRIM(oldIso),"QANCIEN",0.)
       xt_ancien(ixt,:,:)=iso_tmp_lonlev(:,:)
 
       found=phyetat0_get(klev,iso_tmp_lonlev,"XTLANCIEN"//TRIM(outiso),"QLANCIEN",0.)
+      if (.NOT.found.AND.k/=0) found=phyetat0_get(1,iso_tmp_lonlev,"XTLANCIEN"//TRIM(oldIso),"QLANCIEN",0.)
       xtl_ancien(ixt,:,:)=iso_tmp_lonlev(:,:)
 
       found=phyetat0_get(klev,iso_tmp_lonlev,"XTSANCIEN"//TRIM(outiso),"QSANCIEN",0.)
+      if (.NOT.found.AND.k/=0) found=phyetat0_get(1,iso_tmp_lonlev,"XTSANCIEN"//TRIM(oldIso),"QSANCIEN",0.)
       xts_ancien(ixt,:,:)=iso_tmp_lonlev(:,:)
 
       found=phyetat0_get(1,iso_tmp,"XTRUNOFFLIC0"//TRIM(outiso),"RUNOFFLIC0",0.)  
+      if (.NOT.found.AND.k/=0) found=phyetat0_get(1,iso_tmp,"XTRUNOFFLIC0"//TRIM(oldIso),"RUNOFFLIC0",0.)
       xtrun_off_lic_0(ixt,:)=iso_tmp(:)
 
       found=phyetat0_get(klev,iso_tmp_lonlev,"WAKE_DELTAXT"//TRIM(outiso),"Delta hum. wake/env",0.)  
+      if (.NOT.found.AND.k/=0) found=phyetat0_get(1,iso_tmp_lonlev,"WAKE_DELTAXT"//TRIM(oldIso),"Delta hum. wake/env",0.)
       wake_deltaxt(ixt,:,:)=iso_tmp_lonlev(:,:)
 
@@ -16666,7 +16679,9 @@
        if (ixt.le.niso) then
         found=phyetat0_get(1,iso_tmp,"XTSOL"//TRIM(outiso),"Surface hmidity / bucket",0.)  
+        if (.NOT.found.AND.k/=0) found=phyetat0_get(1,iso_tmp,"XTSOL"//TRIM(oldIso),"Surface hmidity / bucket",0.)
         xtsol(ixt,:)=iso_tmp(:)
 
         found=phyetat0_get(1,iso_tmp,"Rland_ice"//TRIM(outiso),"R land ice",0.)
+        if (.NOT.found.AND.k/=0) found=phyetat0_get(1,iso_tmp,"Rland_ice"//TRIM(oldIso),"R land ice",0.)
         Rland_ice(ixt,:)=iso_tmp(:)
 
Index: LMDZ6/trunk/libf/phylmdiso/isotrac_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/isotrac_mod.F90	(revision 4319)
+++ LMDZ6/trunk/libf/phylmdiso/isotrac_mod.F90	(revision 4325)
@@ -1,104 +1,27 @@
 #ifdef ISO
 #ifdef ISOTRAC
-! $Id: $
 
 MODULE isotrac_mod
-use infotrac_phy, ONLY: niso,ntiso,ntraceurs_zone=>nzone
-use isotopes_mod, only: ridicule
-
-IMPLICIT NONE
-SAVE
-
-! contient toutes les variables traceurs isotopiques + les routines specifiquement
-! traceurs isotopiques
-
-      real ridicule_trac
-      parameter (ridicule_trac=ridicule*1e4)
-
-integer, save ::  option_traceurs 
-integer, save ::  ntraceurs_zone_opt ! ntraceurs_zone propre à l'option
-! on vérifie que ça correspond bien à ntraceurs_zone d'infotrac
-integer, save ::  ntraceurs_zoneOR
-!$OMP THREADPRIVATE(option_traceurs,ntraceurs_zone_opt,ntraceurs_zoneOR)
-integer, save ::  initialisation_isotrac
-                ! 1 pour idéalisé
-                ! 0 pour lecture dans fichier
-!$OMP THREADPRIVATE(initialisation_isotrac)
-
-! variables spécifiques aux différentes options, mais necessaires au
-! calcul du nombre de zones de traceurs
-! si option=3
-integer, save :: use_bassin_atlantic
-!$OMP THREADPRIVATE(use_bassin_atlantic)
-integer, save :: use_bassin_medit
-!$OMP THREADPRIVATE(use_bassin_medit)
-integer, save :: use_bassin_indian
-!$OMP THREADPRIVATE(use_bassin_indian)
-integer, save :: use_bassin_austral
-!$OMP THREADPRIVATE(use_bassin_austral)
-integer, save :: use_bassin_pacific
-!$OMP THREADPRIVATE(use_bassin_pacific)
-integer, save :: use_bassin_merarabie
-!$OMP THREADPRIVATE(use_bassin_merarabie)
-integer, save :: use_bassin_golfebengale
-!$OMP THREADPRIVATE(use_bassin_golfebengale)
-integer, save :: use_bassin_indiansud
-!$OMP THREADPRIVATE(use_bassin_indiansud)
-integer, save :: use_bassin_tropics
-!$OMP THREADPRIVATE(use_bassin_tropics)
-integer, save :: use_bassin_midlats
-!$OMP THREADPRIVATE(use_bassin_midlats)
-integer, save :: use_bassin_hauteslats
-!$OMP THREADPRIVATE(use_bassin_hauteslats)
-integer, save :: bassin_atlantic
-!$OMP THREADPRIVATE(bassin_atlantic)
-integer, save :: bassin_medit
-!$OMP THREADPRIVATE(bassin_medit)
-integer, save :: bassin_indian
-!$OMP THREADPRIVATE(bassin_indian)
-integer, save :: bassin_austral
-!$OMP THREADPRIVATE(bassin_austral)
-integer, save :: bassin_pacific
-!$OMP THREADPRIVATE(bassin_pacific)
-integer, save :: bassin_merarabie
-!$OMP THREADPRIVATE(bassin_merarabie)
-integer, save :: bassin_golfebengale
-!$OMP THREADPRIVATE(bassin_golfebengale)
-integer, save :: bassin_indiansud
-!$OMP THREADPRIVATE(bassin_indiansud)
-integer, save :: bassin_tropics
-!$OMP THREADPRIVATE(bassin_tropics)
-integer, save :: bassin_midlats
-!$OMP THREADPRIVATE(bassin_midlats)
-integer, save :: bassin_hauteslats
-!$OMP THREADPRIVATE(bassin_hauteslats)
-! si option=4
-integer nzone_temp
-parameter (nzone_temp=1)
-real, save :: zone_temp1,zone_tempf,zone_tempa  
-!$OMP THREADPRIVATE(zone_temp1,zone_tempf,zone_tempa)
-! si option 14
-integer nzone_lat
-parameter (nzone_lat=4)
-integer nzone_pres
-parameter (nzone_pres=3)
-real, save :: zone_pres1,zone_presf,zone_presa 
-!$OMP THREADPRIVATE(zone_pres1,zone_presf,zone_presa) 
-real, save :: dlattag,lattag_min
-!$OMP THREADPRIVATE(dlattag,lattag_min)
-
-
-! option 1: on trace evap ocean et continent séparement  
+  USE infotrac_phy,      ONLY: niso, ntiso, nzone
+  USE readTracFiles_mod, ONLY: delPhase
+  USE isotopes_mod,      ONLY: ridicule, get_in
+
+  IMPLICIT NONE
+  SAVE
+
+!=== CONTENT: ALL THE ISOTOPIC TRACERS RELATED VARIABLES ===
+!
+! option 1: on trace evap ocean et continent separement  
 ! option 2: on trace evap ocean, continent et evap precip
-! option 3: on trace evap différents bassins océaniques 
-!       + continents + résidu
-!       attention, choisir dans ce cas les bassins océaniques
+! option 3: on trace evap differents bassins oceaniques 
+!       + continents + residu
+!       attention, choisir dans ce cas les bassins oceaniques
 !       dans iso_traceurs_opt3F90.h 
-! option 4: tracage par température minimale
-!       dans ce cas, on définit des bins dans iso_traceurs_opt4.h
-! option 5: pour AMMA: on taggue résidu/AEJ/flux mousson/Harmattan
+! option 4: tracage par temperature minimale
+!       dans ce cas, on definit des bins dans iso_traceurs_opt4.h
+! option 5: pour AMMA: on taggue residu/AEJ/flux mousson/Harmattan
 ! option 6: taggage des ddfts
-! option 7: pour Sandrine: taggage de la vapeur à 700hPa pour omega500<-20 TODO
-! option 8: pour Sandrine: taggage de la vapeur entre 950 et 800hPa, omega de 0 à 25 hPa et de l'évaoration en omega<-20. TODO
+! option 7: pour Sandrine: taggage de la vapeur a 700hPa pour omega500<-20 TODO
+! option 8: pour Sandrine: taggage de la vapeur entre 950 et 800hPa, omega de 0 a 25 hPa et de l'evaoration en omega<-20. TODO
 ! option 9: taggage du condensat et de la revap precip
 ! option 10: taggage evap oce, transpiration et evaporation
@@ -107,910 +30,649 @@
 ! option 12: taggage evap oce, sol nu, canop et reste evap cont. 
 ! A utiliser quand on couple avec ORCHIDEE
-! option 13: taggage température minimale + revap precip
-! option 14: taggage lat et altitude de dernière saturation (niveaux de pression) + evap surf
+! option 13: taggage temperature minimale + revap precip
+! option 14: taggage lat et altitude de derniere saturation (niveaux de pression) + evap surf
 ! otion 15: taggage irrigation
 ! option 16: taggage precip selon saisons et fonte neige: seulement pour ORCHIDEE
-! option 17: taggage température minimum de condensation directement dans la convection et la cond LS, + evap sfc, condensat et precipitation
+! option 17: taggage temperature minimum de condensation directement dans la convection et la cond LS, + evap sfc, condensat et precipitation
 ! option 18: idem 17 mais on tague qsmin au lieu de Tmin
 ! option 19: on tag vap residuelle, vap residuelle dans ddfts, sfc, cond, rev
 ! option 20: on taggue vapeur tropicale vs vapeur extratropicale
 ! option 21: taggage de 2 boites 3D: extratropiques (>35°) et UT tropicale (15-15°, > 500hPa)
-! option 22: tagage de la vapeur proccessée dans les zones très convectives
+! option 22: tagage de la vapeur proccessee dans les zones tres convectives
                 
-        ! ces variables sont initialisées dans traceurs_init
+   !--- nzone_opt (value of nzone for the selected option) must be equal to nzone as defined in onfotrac
+   REAL, PARAMETER :: ridicule_trac = ridicule * 1e4
+   INTEGER, SAVE :: option_traceurs, nzone_opt, nzoneOR
+!$OMP THREADPRIVATE(option_traceurs,nzone_opt,nzoneOR)
+   INTEGER, SAVE :: initialisation_isotrac
+!$OMP THREADPRIVATE(initialisation_isotrac)      
+                ! 1 pour idealise
+                ! 0 pour lecture dans fichier
+
+   !=== VARIABLES SPECIFIC TO THE SELECTED OPTION, BUT NEEDED FOR THE COMPUTATION OF THE NUMBER OF ZONES ; TO BE INITIALIZED IN traceurs_init
+
+   !--- option 3
+   LOGICAL, SAVE :: use_bassin_Austral, use_bassin_Atlantic, use_bassin_MidLats, use_bassin_SouthIndian, use_bassin_MerArabie
+!$OMP THREADPRIVATE(use_bassin_Austral, use_bassin_Atlantic, use_bassin_MidLats, use_bassin_SouthIndian, use_bassin_MerArabie)
+   INTEGER, SAVE ::     bassin_Austral,     bassin_Atlantic,     bassin_MidLats,     bassin_SouthIndian,     bassin_MerArabie
+!$OMP THREADPRIVATE(    bassin_Austral,     bassin_Atlantic,     bassin_MidLats,     bassin_SouthIndian,     bassin_MerArabie)
+   LOGICAL, SAVE :: use_bassin_Pacific, use_bassin_Indian,   use_bassin_Tropics, use_bassin_BengalGolf,  use_bassin_HighLats, use_bassin_Medit
+!$OMP THREADPRIVATE(use_bassin_Pacific, use_bassin_Indian,   use_bassin_Tropics, use_bassin_BengalGolf,  use_bassin_HighLats, use_bassin_Medit)
+   INTEGER, SAVE ::     bassin_Pacific,     bassin_Indian,       bassin_Tropics,     bassin_BengalGolf,      bassin_HighLats,     bassin_Medit
+!$OMP THREADPRIVATE(    bassin_Pacific,     bassin_Indian,       bassin_Tropics,     bassin_BengalGolf,      bassin_HighLats,     bassin_Medit)
+
+   !--- option 4
+   INTEGER, PARAMETER :: nzone_temp = 1
+   REAL,   SAVE ::  zone_temp1, zone_tempf, zone_tempa  
+!$OMP THREADPRIVATE(zone_temp1, zone_tempf, zone_tempa)
+   REAL,   SAVE ::  zone_temp(nzone_temp-1)
+!$OMP THREADPRIVATE(zone_temp)
+
+   !--- option 5
+   INTEGER, SAVE :: izone_aej, izone_harmattan, izone_mousson
+!$OMP THREADPRIVATE(izone_aej, izone_harmattan, izone_mousson)
+
+   !--- option 6
+   INTEGER, SAVE :: izone_ddft
+!$OMP THREADPRIVATE(izone_ddft)
+
+   !--- option 10
+   INTEGER, SAVE :: izone_contfrac
+!$OMP THREADPRIVATE(izone_contfrac)
+
+   !--- option 12	
+   INTEGER, SAVE :: izone_contcanop
+!$OMP THREADPRIVATE(izone_contcanop)
+
+   !--- option 13
+   INTEGER, PARAMETER :: nzone_pres = 3
+   REAL, SAVE ::  zone_pres(nzone_pres-1)
+!$OMP THREADPRIVATE(zone_pres)
+
+   !--- option 14
+   INTEGER, PARAMETER :: nzone_lat = 4
+   REAL,    SAVE :: zone_pres1, zone_presf, zone_presa 
+!$OMP THREADPRIVATE(zone_pres1, zone_presf, zone_presa) 
+   REAL,    SAVE :: dlattag, lattag_min, zone_lat(nzone_lat-1)
+!$OMP THREADPRIVATE(dlattag, lattag_min, zone_lat)
+
+   !--- option 15
+   INTEGER, SAVE :: izone_irrig
+!$OMP THREADPRIVATE(izone_irrig)
+
+   !--- option 17
+   REAL,    SAVE :: seuil_tag_tmin, seuil_tag_tmin_ls
+!$OMP THREADPRIVATE(seuil_tag_tmin, seuil_tag_tmin_ls)
+  INTEGER,  SAVE :: option_seuil_tag_tmin
+!$OMP THREADPRIVATE(option_seuil_tag_tmin)
+
+   !--- option 20
+   INTEGER, SAVE :: izone_trop, izone_extra
+!$OMP THREADPRIVATE(izone_trop, izone_extra)
+   REAL,    SAVE :: lim_tag20 
+!$OMP THREADPRIVATE(lim_tag20)
+
+   !--- option 21: on garde izone_trop, izone_extra  
+
+   !--- option 22
+   INTEGER, SAVE :: izone_conv_BT, izone_conv_UT
+!$OMP THREADPRIVATE(izone_conv_BT, izone_conv_UT)
+   REAL,    SAVE :: lim_precip_tag22 
+!$OMP THREADPRIVATE(lim_precip_tag22)
+
         
-integer, ALLOCATABLE, DIMENSION(:), save :: index_iso
-!$OMP THREADPRIVATE(index_iso)
-integer, ALLOCATABLE, DIMENSION(:), save ::  index_zone
-!$OMP THREADPRIVATE(index_zone)
-integer, ALLOCATABLE, DIMENSION(:,:), save ::  itZonIso_loc ! il y a déjà un itZonIso dans infotrac: vérifier que c'est le même
-!$OMP THREADPRIVATE(itZonIso_loc)
-character*3, ALLOCATABLE, DIMENSION(:), save :: strtrac
-!$OMP THREADPRIVATE(strtrac)
-! -> tout ça passe maintenant par infotrac
-
-integer, ALLOCATABLE, DIMENSION(:), save :: bassin_map
-integer, ALLOCATABLE, DIMENSION(:,:), save :: boite_map
-!$OMP THREADPRIVATE(bassin_map,boite_map)
-
-
-        ! traitement recyclage et evap
-integer, save :: izone_cont ! pour le recyclage continental
-!$OMP THREADPRIVATE(izone_cont)
-integer, save :: izone_oce ! pour l'océan
-!$OMP THREADPRIVATE(izone_oce)
-integer, save :: izone_poubelle ! pour les petits résidus numériques
+  INTEGER, ALLOCATABLE, SAVE :: index_iso(:), index_zone(:), itZonIso_loc(:,:)
+!$OMP THREADPRIVATE(            index_iso,    index_zone,    itZonIso_loc)
+  CHARACTER(LEN=3), ALLOCATABLE :: strtrac(:)
+!$OMP THREADPRIVATE(               strtrac)
+  INTEGER, ALLOCATABLE, SAVE :: bassin_map(:), boite_map(:,:)
+!$OMP THREADPRIVATE(            bassin_map,    boite_map)
+
+   !=== RECYCLING AND EVAPORATION TREATMENT
+   INTEGER, SAVE :: izone_cont, izone_oce        !--- For land and ocean recycling
+!$OMP THREADPRIVATE(izone_cont, izone_oce)
+   INTEGER, SAVE :: izone_poubelle               !--- For small numerical residues
 !$OMP THREADPRIVATE(izone_poubelle)
-integer, save :: izone_init ! pour l'initialisation par défaut
+   INTEGER, SAVE :: izone_init                   !--- For default initialization
 !$OMP THREADPRIVATE(izone_init)
-integer, save :: izone_revap ! pour l'évap des gouttes
+   INTEGER, SAVE :: izone_revap                  !--- For droplets evaporation
 !$OMP THREADPRIVATE(izone_revap)
-integer, save :: option_revap
-!$OMP THREADPRIVATE(option_revap)
-integer, save :: option_tmin
-!$OMP THREADPRIVATE(option_tmin)
-integer, save :: option_cond
-!$OMP THREADPRIVATE(option_cond)
-integer, save :: izone_cond
-!$OMP THREADPRIVATE(izone_cond)
-real evap_franche
-parameter (evap_franche=1e-6) ! en kg/m2/s
-
-! specifique à option 4:
-real, save ::  zone_temp(nzone_temp-1)
-!$OMP THREADPRIVATE(zone_temp)
-! si option 5
-integer, save :: izone_aej
-!$OMP THREADPRIVATE(izone_aej)
-integer, save :: izone_harmattan
-!$OMP THREADPRIVATE(izone_harmattan)
-integer, save :: izone_mousson
-!$OMP THREADPRIVATE(izone_mousson)
-! si option 6
-integer, save :: izone_ddft
-!$OMP THREADPRIVATE(izone_ddft)
-! si option 10
-integer, save :: izone_contfrac
-!$OMP THREADPRIVATE(izone_contfrac)
-! si option 12	
-integer, save :: izone_contcanop
-!$OMP THREADPRIVATE(izone_contcanop)
-! specifique à option 13:
-real, save ::  zone_pres(nzone_pres-1)
-!$OMP THREADPRIVATE(zone_pres)
-! si option 14
-real, save ::  zone_lat(nzone_lat-1)
-!$OMP THREADPRIVATE(zone_lat)
-! si option 15
-integer, save :: izone_irrig
-!$OMP THREADPRIVATE(izone_irrig)
-! si option 17
-real, save ::  seuil_tag_tmin
-!$OMP THREADPRIVATE(seuil_tag_tmin)
-real, save ::  seuil_tag_tmin_ls
-!$OMP THREADPRIVATE(seuil_tag_tmin_ls)
-integer, save :: option_seuil_tag_tmin
-!$OMP THREADPRIVATE(option_seuil_tag_tmin)
-! si option 20
-integer, save :: izone_trop,izone_extra
-real, save ::  lim_tag20 
-!$OMP THREADPRIVATE(izone_trop,izone_extra,lim_tag20)
-! si option 21: on garde izone_trop,izone_extra  
-! si opt 22
-integer, save :: izone_conv_BT,izone_conv_UT
-real, save ::  lim_precip_tag22 
-!$OMP THREADPRIVATE(izone_conv_BT,izone_conv_UT,lim_precip_tag22)
-
+   INTEGER, SAVE :: option_revap, option_tmin, option_cond, izone_cond
+!$OMP THREADPRIVATE(option_revap, option_tmin, option_cond, izone_cond)
+   REAL, PARAMETER :: evap_franche = 1e-6        !--- In kg/m2/s
 
 CONTAINS
 
-      subroutine iso_traceurs_init()
-
-      use IOIPSL ! getin
-      USE infotrac_phy, ONLY: itZonIso
-      USE isotopes_mod, ONLY: iso_eau,ntracisoOR,initialisation_iso
-      USE dimphy, only: klon,klev
-
-        implicit none
-
-
-        ! définition de quelles zones et quelles isotopes représentent
-        ! les traceurs
-
-        ! inputs, outputs
-        ! ! c'est les variables dans traceurs.h qui sont modifiées
-
-        ! locals
-        integer itrac,izone,ixt,k
-        integer izone_pres,izone_lat
-        character*2 strz,strz_preslat
-        character*1 strz_pres,strz_lat
-        integer ntraceurs_zone_opt
-
-        ! vérifier que on a bien l'eau comme traceurs
-        if (iso_eau.eq.0) then
-            write(*,*) 'traceurs_init 18: isotrac ne marche que si ', &
-     &            'on met l''eau comme isotope'
-            stop
-        endif
-
-        ! initialiser
-        option_traceurs=0
-        initialisation_isotrac=0
-
-        ! allouer
-        allocate (index_iso(ntiso))
-        allocate (index_zone(ntiso))
-        allocate (itZonIso_loc(ntraceurs_zone,niso))
-        allocate (strtrac(ntraceurs_zone))
-        allocate (bassin_map(klon))
-        allocate (boite_map(klon,klev))
-
-        if (initialisation_iso.eq.0) then
-          call getin('initialisation_isotrac',initialisation_isotrac)
-          write(*,*) 'initialisation_isotrac=',initialisation_isotrac
-        endif !if (initialisation_iso.eq.0) then
-
-        ! lire l'option de traçage
-        call getin('option_traceurs',option_traceurs)
-        write(*,*) 'option_traceurs=',option_traceurs
-
-        ! cas général: pas de traceurs dans ORCHIDEE
-        ntracisoOR=niso
-
-        ! partie à éditer ! pour définir les différentes zones
-        if (option_traceurs.eq.1) then
-          ! on trace continents/ocean  
-
-          ntraceurs_zone_opt=2
-          izone_cont=1
-          izone_oce=2          
-          izone_poubelle=2 ! zone où on met les flux non physiques, de
-                ! réajustement
-          izone_init=2 ! zone d'initialisation par défaut          
-          option_revap=0
-          option_tmin=0
-          izone_revap=0
-          option_cond=0
-
-          strtrac(izone_cont)='con'
-          strtrac(izone_oce)='oce'
-
-        elseif (option_traceurs.eq.2) then
-          ! on trace continent/ ocean/reevap des gouttes
-
-          ntraceurs_zone_opt=3
-          izone_cont=1
-          izone_oce=2
-          izone_poubelle=2 ! zone où on met les flux non physiques, de
-                ! réajustement
-          izone_init=2 ! zone d'initialisation par défaut
-          option_revap=1
-          option_tmin=0
-          izone_revap=3
-          option_cond=0
-
-          strtrac(izone_cont)='con'
-          strtrac(izone_oce)='oce'
-          strtrac(izone_revap)='rev'
-          
-
-        else if (option_traceurs.eq.3) then
-            ! on trace des bassins océaniques + un résidu. On ne trace
-            ! pas l'évap des gouttes à part
-            ! le résidu est la dernère dimension
-            
-          ! lire les use_bassin
-          call getin('use_bassin_atlantic',use_bassin_atlantic)     
-          call getin('use_bassin_medit',use_bassin_medit)      
-          call getin('use_bassin_indian',use_bassin_indian)     
-          call getin('use_bassin_austral',use_bassin_austral)     
-          call getin('use_bassin_pacific',use_bassin_pacific)     
-          call getin('use_bassin_merarabie',use_bassin_merarabie)     
-          call getin('use_bassin_golfebengale',use_bassin_golfebengale)     
-          call getin('use_bassin_indiansud',use_bassin_indiansud)     
-          call getin('use_bassin_tropics',use_bassin_tropics)     
-          call getin('use_bassin_midlats',use_bassin_midlats)     
-          call getin('use_bassin_hauteslats',use_bassin_hauteslats)
-
-          write(*,*) 'use_bassin_atlantic=' ,use_bassin_atlantic  
-          write(*,*) 'use_bassin_medit=' ,use_bassin_medit
-          write(*,*) 'use_bassin_indian=' ,use_bassin_indian
-          write(*,*) 'use_bassin_austral=' ,use_bassin_austral
-          write(*,*) 'use_bassin_merarabie=' ,use_bassin_merarabie
-          write(*,*) 'use_bassin_golfebengale=' ,use_bassin_golfebengale
-          write(*,*) 'use_bassin_indiansud=' ,use_bassin_indiansud
-          write(*,*) 'use_bassin_tropics=' ,use_bassin_tropics
-          write(*,*) 'use_bassin_midlats=' ,use_bassin_midlats
-          write(*,*) 'use_bassin_hauteslats=' ,use_bassin_hauteslats
-
-        
-          ntraceurs_zone_opt=2 &
-     &                   +use_bassin_atlantic &
-     &                   +use_bassin_medit &
-     &                   +use_bassin_indian &
-     &                   +use_bassin_austral &
-     &                   +use_bassin_pacific &
-     &                   +use_bassin_merarabie &
-     &                   +use_bassin_golfebengale &
-     &                   +use_bassin_indiansud &
-     &                   +use_bassin_tropics &
-     &                   +use_bassin_midlats &
-     &                   +use_bassin_hauteslats
-
-          izone_cont=ntraceurs_zone
-          izone_oce=0 ! pas de sens car séparée en bassins          
-          izone_poubelle=ntraceurs_zone-1 ! zone où on met les flux non physiques, de
-                ! réajustement
-          izone_init=ntraceurs_zone-1 ! zone d'initialisation par défaut
-          option_revap=0 ! on ne trace pas les gouttes 
-          option_tmin=0
-          izone_revap=0 ! pas de sens car on taggue pas les gouttes séparemment  
-          option_cond=0
-
-          ! si on a use_bassin_indian, on n'a pas le découpage détaillé
-          ! de l'indian:
+   SUBROUTINE iso_traceurs_init()
+
+   USE infotrac_phy, ONLY: itZonIso, isoName, isoZone
+   USE isotopes_mod, ONLY: iso_eau, ntracisoOR, initialisation_iso
+   USE dimphy,       ONLY: klon, klev
+   USE  strings_mod, ONLY: int2str, strStack, strTail, strHead, fmsg
+
+   IMPLICIT NONE
+   ! Define which zones and isotopes correspond to isotopic tagging tracers
+   ! Modify traceurs.h variables
+   INTEGER :: izone, ixt, k
+   INTEGER :: izone_pres, izone_lat
+   INTEGER :: nzone_opt
+
+   IF(fmsg("traceurs_init 18: isotrac ne marche que si on met l'eau comme isotope", 'iso_traceurs_init', iso_eau==0)) STOP
+
+   !--- Initialize
+   option_traceurs = 0
+   initialisation_isotrac = 0
+
+   !--- Allocate
+   ALLOCATE(index_iso (ntiso))
+   ALLOCATE(index_zone(ntiso))
+   ALLOCATE(itZonIso_loc(nzone,niso))
+   ALLOCATE(strtrac(nzone))
+   ALLOCATE(bassin_map(klon))
+   ALLOCATE( boite_map(klon,klev))
+
+   IF(initialisation_iso == 0) CALL get_in('initialisation_isotrac', initialisation_isotrac)
+
+   !--- Read tracing option
+   CALL get_in('option_traceurs', option_traceurs)
+
+   !--- Genral case: no traceurs in ORCHIDEE
+   ntracisoOR=niso
+
+   ! partie a editer ! pour definir les differentes zones
+   SELECT CASE(option_traceurs)
+      !========================================================================================================================
+      CASE(1)      !=== TRACING LAND/OCEAN
+      !========================================================================================================================
+         nzone_opt=2
+         izone_cont=1
+         izone_oce=2
+         izone_poubelle=2    ! zone ou on met les flux non physiques, de reajustement
+         izone_init=2        ! zone d'initialisation par defaut          
+         option_revap=0
+         option_tmin=0
+         izone_revap=0
+         option_cond=0
+         strtrac(izone_cont) = 'con'
+         strtrac(izone_oce)  = 'oce'
+      !========================================================================================================================
+      CASE(2)      !=== TRACING LAND/OCEAN/DROPLETS REEVAPORATION
+      !========================================================================================================================
+         nzone_opt=3
+         izone_cont=1
+         izone_oce=2
+         izone_poubelle=2    ! zone ou on met les flux non physiques, de reajustement
+         izone_init=2        ! zone d'initialisation par defaut          
+         option_revap=1
+         option_tmin=0
+         izone_revap=3
+         option_cond=0
+         strtrac(izone_cont) = 'con'
+         strtrac(izone_oce)  = 'oce'
+         strtrac(izone_revap)= 'rev'
+      !========================================================================================================================
+      CASE(3)      !=== TRACING OCEANS BASINS + RESIDUE (LAST DIMENSION). NO DROPLETS EVAPORATION TRACING.
+      !========================================================================================================================
+         ! lire les use_bassin
+         CALL get_in('use_bassin_Atlantic',   use_bassin_Atlantic)
+         CALL get_in('use_bassin_Medit',      use_bassin_Medit)
+         CALL get_in('use_bassin_Indian',     use_bassin_Indian)
+         CALL get_in('use_bassin_Austral',    use_bassin_Austral)
+         CALL get_in('use_bassin_Pacific',    use_bassin_Pacific)
+         CALL get_in('use_bassin_MerArabie',  use_bassin_MerArabie)
+         CALL get_in('use_bassin_BengalGolf', use_bassin_BengalGolf)
+         CALL get_in('use_bassin_SouthIndian',use_bassin_SouthIndian)
+         CALL get_in('use_bassin_Tropics',    use_bassin_Tropics)
+         CALL get_in('use_bassin_Midlats',    use_bassin_Midlats)
+         CALL get_in('use_bassin_HighLats',   use_bassin_HighLats)
+         nzone_opt  =  2  +  COUNT([use_bassin_Atlantic, use_bassin_Medit,     use_bassin_Indian,     &
+            use_bassin_Austral,     use_bassin_Pacific,  use_bassin_MerArabie, use_bassin_BengalGolf, &
+            use_bassin_SouthIndian, use_bassin_Tropics,  use_bassin_Midlats,   use_bassin_HighLats])
+         izone_cont=nzone
+         izone_oce=0             ! pas de sens car separee en bassins          
+         izone_poubelle=nzone-1  ! zone ou on met les flux non physiques, de reajustement
+         izone_init=nzone-1      ! zone d'initialisation par defaut
+         option_revap=0          ! on ne trace pas les gouttes 
+         option_tmin=0
+         izone_revap=0           ! pas de sens car on taggue pas les gouttes separemment  
+         option_cond=0
 #ifdef ISOVERIF
-          if (use_bassin_indian.eq.1) then
-!              call iso_verif_egalite(float(use_bassin_merarabie), &
-!     &            0.0,'iso_traceurs_init 73: revoir def des bassins')
-               if ((use_bassin_merarabie.ne.0).or. &
-      &            (use_bassin_indiansud.ne.0).or. &
-      &            (use_bassin_golfebengale.ne.0)) then
-                write(*,*) 'traceurs_init 73'
-                stop
-               endif 
-!              call iso_verif_egalite(float(use_bassin_golfebengale), &
-!     &            0.0,'iso_traceurs_init 73: revoir def des bassins')
-!              call iso_verif_egalite(float(use_bassin_indiansud), &
-!     &            0.0,'iso_traceurs_init 73: revoir def des bassins')
-          endif
+         IF(use_bassin_Indian) THEN   !=== NON COMPATIBLE WITH A DETAILED INDIAN CUTTING
+            IF(use_bassin_MerArabie .OR. use_bassin_SouthIndian .OR. use_bassin_BengalGolf) THEN
+               WRITE(*,*) 'traceurs_init 73'; STOP
+            END IF
+!           CALL iso_verif_egalite(float(use_bassin_MerArabie),   0.0, 'iso_traceurs_init 73: revoir def des bassins')
+!           CALL iso_verif_egalite(float(use_bassin_BengalGolf),  0.0, 'iso_traceurs_init 73: revoir def des bassins')
+!           CALL iso_verif_egalite(float(use_bassin_SouthIndian), 0.0, 'iso_traceurs_init 73: revoir def des bassins')
+         END IF
 #endif    
-          
-          bassin_atlantic= max(use_bassin_atlantic,1)
-          bassin_medit=max(use_bassin_atlantic &
-     &           +use_bassin_medit,1)
-          bassin_indian=max(use_bassin_atlantic &
-     &           +use_bassin_medit &
-     &           +use_bassin_indian,1)
-          bassin_austral=max(use_bassin_atlantic &
-     &           +use_bassin_medit &
-     &           +use_bassin_indian &
-     &           +use_bassin_austral,1)
-          bassin_pacific=max(use_bassin_atlantic &
-     &           +use_bassin_medit &
-     &           +use_bassin_indian &
-     &           +use_bassin_austral &
-     &           +use_bassin_pacific,1)
-          bassin_merarabie=max(use_bassin_atlantic &
-     &           +use_bassin_medit &
-     &           +use_bassin_indian &
-     &           +use_bassin_austral &
-     &           +use_bassin_pacific &
-     &           +use_bassin_merarabie,1)
-          bassin_golfebengale=max(use_bassin_atlantic& 
-     &           +use_bassin_medit &
-     &           +use_bassin_indian & 
-     &           +use_bassin_austral &
-     &           +use_bassin_pacific &
-     &           +use_bassin_merarabie &
-     &           +use_bassin_golfebengale,1)
-          bassin_indiansud=max(use_bassin_atlantic &
-     &           +use_bassin_medit &
-     &           +use_bassin_indian &
-     &           +use_bassin_austral &
-     &           +use_bassin_pacific &
-     &           +use_bassin_merarabie &
-     &           +use_bassin_golfebengale &
-     &           +use_bassin_indiansud,1) 
-          bassin_tropics=max(use_bassin_atlantic &
-     &                       +use_bassin_medit &
-     &                       +use_bassin_indian &
-     &                       +use_bassin_austral &
-     &                       +use_bassin_pacific &
-     &                       +use_bassin_merarabie &
-     &                       +use_bassin_golfebengale &
-     &                       +use_bassin_indiansud, &
-     &                       +use_bassin_tropics,1) 
-          bassin_midlats=max(use_bassin_atlantic &
-     &                       +use_bassin_medit &
-     &                       +use_bassin_indian &
-     &                       +use_bassin_austral &
-     &                       +use_bassin_pacific &
-     &                       +use_bassin_merarabie &
-     &                       +use_bassin_golfebengale &
-     &                       +use_bassin_indiansud &
-     &                       +use_bassin_tropics &
-     &                       +use_bassin_midlats,1) 
-          bassin_hauteslats=max(use_bassin_atlantic & 
-     &                       +use_bassin_medit &
-     &                       +use_bassin_indian &
-     &                       +use_bassin_austral &
-     &                       +use_bassin_pacific &
-     &                       +use_bassin_merarabie &
-     &                       +use_bassin_golfebengale &
-     &                       +use_bassin_indiansud &
-     &                       +use_bassin_tropics &
-     &                       +use_bassin_midlats &
-     &                       +use_bassin_hauteslats,1) 
-
-          write(*,*) 'bassin_atlantic=' ,bassin_atlantic  
-          write(*,*) 'bassin_medit=' ,bassin_medit
-          write(*,*) 'bassin_indian=' ,bassin_indian
-          write(*,*) 'bassin_austral=' ,bassin_austral
-          write(*,*) 'bassin_merarabie=' ,bassin_merarabie
-          write(*,*) 'bassin_golfebengale=' ,bassin_golfebengale
-          write(*,*) 'bassin_indiansud=' ,bassin_indiansud
-          write(*,*) 'bassin_tropics=' ,bassin_tropics
-          write(*,*) 'bassin_midlats=' ,bassin_midlats
-          write(*,*) 'bassin_hauteslats=' ,bassin_hauteslats
-
-          if (use_bassin_atlantic.eq.1) then
-            strtrac(bassin_atlantic)='atl'
-          endif
-          if (use_bassin_medit.eq.1) then
-            strtrac(bassin_medit)='med'
-          endif
-          if (use_bassin_indian.eq.1) then
-            strtrac(bassin_indian)='ind'
-          endif
-          if (use_bassin_austral.eq.1) then
-            strtrac(bassin_austral)='aus'
-          endif
-          if (use_bassin_pacific.eq.1) then
-            strtrac(bassin_pacific)='pac'
-          endif
-          if (use_bassin_merarabie.eq.1) then
-            strtrac(bassin_merarabie)='ara'
-          endif
-          if (use_bassin_golfebengale.eq.1) then
-            strtrac(bassin_golfebengale)='ben'
-          endif
-          if (use_bassin_indiansud.eq.1) then
-            strtrac(bassin_indiansud)='ins'
-          endif
-          if (use_bassin_tropics.eq.1) then
-            strtrac(bassin_tropics)='tro'
-          endif
-          if (use_bassin_midlats.eq.1) then
-            strtrac(bassin_midlats)='mid'
-          endif
-          if (use_bassin_hauteslats.eq.1) then
-            strtrac(bassin_hauteslats)='hau'
-          endif
-          strtrac(ntraceurs_zone-1)='res'
-          strtrac(ntraceurs_zone)='con'
-
-        else if (option_traceurs.eq.4) then
-          ! on trace les température minimales vécues
-          ! comme dans article sur LdG sauf pas de revap
-            
-          zone_temp1=293.0 ! en K
-!          zone_tempf=223.0 ! en K
-          zone_tempf=243.0 ! en K
- ! courbure de la relation entre l'indice et la température: 0 pour linéaire, <0 pour plus de détal en bas
-
+         bassin_Atlantic   = 1
+         bassin_Medit      = bassin_Atlantic    + COUNT([use_bassin_Medit]);       WRITE(*,*) 'bassin_Atlantic    =' ,bassin_Atlantic
+         bassin_Indian     = bassin_Medit       + COUNT([use_bassin_Indian]);      WRITE(*,*) 'bassin_Medit       =' ,bassin_Medit
+         bassin_Austral    = bassin_Indian      + COUNT([use_bassin_Austral]);     WRITE(*,*) 'bassin_Indian      =' ,bassin_Indian
+         bassin_Pacific    = bassin_Austral     + COUNT([use_bassin_Pacific]);     WRITE(*,*) 'bassin_Austral     =' ,bassin_Austral
+         bassin_MerArabie  = bassin_Pacific     + COUNT([use_bassin_MerArabie]);   WRITE(*,*) 'bassin_MerArabie   =' ,bassin_MerArabie
+         bassin_BengalGolf = bassin_MerArabie   + COUNT([use_bassin_BengalGolf]);  WRITE(*,*) 'bassin_BengalGolf  =' ,bassin_BengalGolf
+         bassin_SouthIndian= bassin_BengalGolf  + COUNT([use_bassin_SouthIndian]); WRITE(*,*) 'bassin_SouthIndian =' ,bassin_SouthIndian
+         bassin_Tropics    = bassin_SouthIndian + COUNT([use_bassin_Tropics]);     WRITE(*,*) 'bassin_Tropics     =' ,bassin_Tropics
+         bassin_MidLats    = bassin_Tropics     + COUNT([use_bassin_MidLats]);     WRITE(*,*) 'bassin_MidLats     =' ,bassin_MidLats
+         bassin_HighLats   = bassin_MidLats     + COUNT([use_bassin_HighLats]);    WRITE(*,*) 'bassin_HighLats    =' ,bassin_HighLats
+         IF(use_bassin_atlantic   ) strtrac(bassin_atlantic)   = 'atl'
+         IF(use_bassin_medit      ) strtrac(bassin_medit)      = 'med'
+         IF(use_bassin_indian     ) strtrac(bassin_indian)     = 'ind'
+         IF(use_bassin_austral    ) strtrac(bassin_austral)    = 'aus'
+         IF(use_bassin_pacific    ) strtrac(bassin_pacific)    = 'pac'
+         IF(use_bassin_merarabie  ) strtrac(bassin_merarabie)  = 'ara'
+         IF(use_bassin_BengalGolf ) strtrac(bassin_BengalGolf) = 'ben'
+         IF(use_bassin_SouthIndian) strtrac(bassin_SouthIndian)= 'ins'
+         IF(use_bassin_tropics    ) strtrac(bassin_tropics)    = 'tro'
+         IF(use_bassin_midlats    ) strtrac(bassin_midlats)    = 'mid'
+         IF(use_bassin_HighLats   ) strtrac(bassin_HighLats)   = 'hau'
+         strtrac(nzone-1)='res'
+         strtrac(nzone)='con'
+      !========================================================================================================================
+      CASE(4)      !=== TRACING MINIMAL EXPERIENCED TEMPERATURE AS IN THE ARTICLE ON LfG, EXCEPT NO REVAPORATION
+      !========================================================================================================================
+         zone_temp1 = 293.0  ! en K
+!        zone_tempf = 223.0  ! en K
+         zone_tempf = 243.0  ! en K
+        ! courbure de la relation entre l'indice et la temperature: 0 pour lineaire, <0 pour plus de detal en bas
         ! zone 1: >= zone_temp1
-        ! zone 2 à 4: intermédiaire, 
+        ! zone 2 a 4: intermediaire, 
         ! zone 5: <zone_tempf
-        
-          ntraceurs_zone_opt=nzone_temp+1
-
-          zone_tempa=-4.0 ! en K
-          izone_cont=ntraceurs_zone
-          izone_oce=ntraceurs_zone  
-          izone_poubelle=ntraceurs_zone
-          izone_init=ntraceurs_zone ! zone d'initialisation par défaut
-          option_revap=0 
-          option_tmin=0  
-          izone_revap=0
-          option_cond=0
-          do izone=1,nzone_temp
-            write(strz,'(i2.2)') izone
-            strtrac(izone)='t'//strz
-            write(*,*) 'izone,strz,strtrac=',izone,strz,strtrac(izone)
-          enddo
-          strtrac(izone_poubelle)='pou'
-
-          ! initialisation des zones de tempéarture
-          do izone=1,nzone_temp-1
-            zone_temp(izone)=zone_temp1+float(izone-1) &
-     &                      *(zone_tempa*float(izone-nzone_temp+1) &
-     &                      +(zone_tempf-zone_temp1)/float(nzone_temp-2))
-          enddo
-          write(*,*) 'iso_trac_init 183: zone_temp=',zone_temp          
-
-        elseif (option_traceurs.eq.5) then
-          ! on trace AEJ/flux de mousson/Harmattan
-!          write(*,*) 'iso_traceurs_init 129'
-
-          ntraceurs_zone_opt=4
-          izone_cont=1
-          izone_oce=1
-          izone_poubelle=1 ! zone où on met les flux non physiques, de
-                ! réajustement
-          izone_init=1 ! zone d'initialisation par défaut
-          option_revap=0
-          option_tmin=0
-          izone_revap=0
-          izone_aej=2
-          izone_mousson=3
-          izone_harmattan=4 
-          option_cond=0
-
-          strtrac(izone_poubelle)='res'
-          strtrac(izone_aej)='aej'
-          strtrac(izone_mousson)='mou'
-          strtrac(izone_harmattan)='sah'
-
-        elseif (option_traceurs.eq.6) then
-          ! on trace les ddfts
-
-          ntraceurs_zone_opt=2
-          izone_cont=1
-          izone_oce=1
-          izone_poubelle=1 ! zone où on met les flux non physiques, de
-                ! réajustement
-          izone_init=1 ! zone d'initialisation par défaut
-          option_revap=0
-          option_tmin=0
-          izone_revap=0
-          izone_ddft=2
-          option_cond=0
-
-          strtrac(izone_poubelle)='res'
-          strtrac(izone_ddft)='dft'
-
-        elseif (option_traceurs.eq.9) then
-          ! on trace le condensat
-
-          ntraceurs_zone_opt=3
-          izone_cont=1
-          izone_oce=1
-          izone_poubelle=1 ! zone où on met les flux non physiques, de
-                ! réajustement
-          izone_init=1 ! zone d'initialisation par défaut
-          option_revap=1
-          option_tmin=0
-          izone_revap=2
-          izone_cond=3
-          option_cond=1
-
-          ! 1 par défaut pour colorier à la fois condensat LS et
-          ! condensat convectif. Mais on peut mettre 2 si on ne veut que
-          ! collorier que le condensat convectif.
-          call getin('option_cond',option_cond)
-          write(*,*) 'option_cond=',option_cond
-
-          strtrac(izone_poubelle)='res'
-          strtrac(izone_cond)='con'
-          strtrac(izone_revap)='rev'
-
-        elseif (option_traceurs.eq.10) then
-          ! on trace l'évap venant de ocean/continent no frac/continent frac
-          !  utilse seulement si couplé avec ORCHIDEE
-#ifdef CPP_VEGET
-#else
-          write(*,*) 'iso_traceurs_init 219: option_traceurs=10 ', &
-     &                      'inutile si on ne couple pas avec ORCHIDEE'
-          stop
+         nzone_opt=nzone_temp+1
+         zone_tempa=-4.0     ! en K
+         izone_cont=nzone
+         izone_oce=nzone  
+         izone_poubelle=nzone
+         izone_init=nzone    ! zone d'initialisation par defaut
+         option_revap=0 
+         option_tmin=0  
+         izone_revap=0
+         option_cond=0
+         DO izone=1,nzone_temp
+            strtrac(izone) = 't'//TRIM(int2str(izone))
+            WRITE(*,*) 'izone, strtrac=', izone, strtrac(izone)
+         END DO
+         strtrac(izone_poubelle)='pou'
+         ! Initialization of temperatures zones
+         DO izone=1,nzone_temp-1
+            zone_temp(izone) = zone_temp1+float(izone-1)            &
+                            * (zone_tempa*float(izone-nzone_temp+1) &
+                            + (zone_tempf-zone_temp1)/float(nzone_temp-2))
+         END DO
+         WRITE(*,*) 'iso_trac_init 183: zone_temp=', zone_temp
+      !========================================================================================================================
+      CASE(5)      !=== TRACING AEJ/MOONSOON FLUX/Harmattan
+      !========================================================================================================================
+!        WRITE*,*) 'iso_traceurs_init 129'
+         nzone_opt=4
+         izone_cont=1
+         izone_oce=1
+         izone_poubelle=1    ! zone ou on met les flux non physiques, de reajustement
+         izone_init=1        ! zone d'initialisation par defaut          
+         option_revap=0
+         option_tmin=0
+         izone_revap=0
+         izone_aej=2
+         izone_mousson=3
+         izone_harmattan=4 
+         option_cond=0
+         strtrac(izone_poubelle) = 'res'
+         strtrac(izone_aej)      = 'aej'
+         strtrac(izone_mousson)  = 'mou'
+         strtrac(izone_harmattan)= 'sah'
+      !========================================================================================================================
+      CASE(6)      !=== TRACING DDFTS
+      !========================================================================================================================
+         nzone_opt=2
+         izone_cont=1
+         izone_oce=1
+         izone_poubelle=1    ! zone ou on met les flux non physiques, de reajustement
+         izone_init=1        ! zone d'initialisation par defaut          
+         option_revap=0
+         option_tmin=0
+         izone_revap=0
+         izone_ddft=2
+         option_cond=0
+         strtrac(izone_poubelle)='res'
+         strtrac(izone_ddft)='dft'
+      !========================================================================================================================
+      CASE(9)      !=== TRACING CONDENSATION
+      !========================================================================================================================
+         nzone_opt=3
+         izone_cont=1
+         izone_oce=1
+         izone_poubelle=1    ! zone ou on met les flux non physiques, de reajustement
+         izone_init=1        ! zone d'initialisation par defaut          
+         option_revap=1
+         option_tmin=0
+         izone_revap=2
+         izone_cond=3
+         option_cond=1
+         ! 1 par defaut pour colorier a la fois condensat LS et condensat convectif.
+         ! Mais on peut mettre 2 si on ne veut que colorier que le condensat convectif.
+         CALL get_in('option_cond',option_cond)
+         strtrac(izone_poubelle)='res'
+         strtrac(izone_cond)='con'
+         strtrac(izone_revap)='rev'
+      !========================================================================================================================
+      CASE(10)     !=== TRACING EVAPORATION FROM OCEAN/LAND, NON FRAC/LAND FRAC ; ONLY WHEN COUPLED WITH ORCHIDEE
+      !========================================================================================================================
+#ifndef CPP_VEGET
+         WRITE(*,*) 'iso_traceurs_init 219: option_traceurs=10 inutile si on ne couple pas avec ORCHIDEE'; STOP
 #endif          
-
-          ntraceurs_zone_opt=3
-          izone_cont=1 ! sous-entendu non fractionnant
-          izone_oce=2
-          izone_poubelle=2 ! zone où on met les flux non physiques, de
-                ! réajustement
-          izone_init=2 ! zone d'initialisation par défaut
-          option_revap=0
-          option_tmin=0
-          izone_revap=0
-          izone_contfrac=3
-          izone_contcanop=3
-          izone_irrig=0
-          option_cond=0
-
-          strtrac(izone_oce)='oce'
-          strtrac(izone_cont)='con'  
-          strtrac(izone_contfrac)='enu'  ! evap sol nu
-
-        elseif (option_traceurs.eq.11) then
-          ! on trace reevap des gouttes et le reste
-
-          ntraceurs_zone_opt=2
-          izone_cont=1
-          izone_oce=1
-          izone_poubelle=1 ! zone où on met les flux non physiques, de
-                ! réajustement
-          izone_init=1 ! zone d'initialisation par défaut
-          option_revap=1
-          option_tmin=0
-          izone_revap=2
-          izone_irrig=0
-          option_cond=0
-
-          strtrac(izone_poubelle)='res'
-          strtrac(izone_revap)='rev'
-
-        elseif (option_traceurs.eq.12) then
-          ! on trace evap du sol nu, evap de la canopée, reste de l'evap cont et
-          ! evap oce
-#ifdef CPP_VEGET
-#else
-          write(*,*) 'iso_traceurs_init 257: option_traceurs=10 ', &
-     &                      'inutile si on ne couple pas avec ORCHIDEE'
-          stop
+         nzone_opt=3
+         izone_cont=1        ! sous-entendu non fractionnant
+         izone_oce=2
+         izone_poubelle=2    ! zone ou on met les flux non physiques, de reajustement
+         izone_init=2        ! zone d'initialisation par defaut
+         option_revap=0
+         option_tmin=0
+         izone_revap=0
+         izone_contfrac=3
+         izone_contcanop=3
+         izone_irrig=0
+         option_cond=0
+         strtrac(izone_oce)='oce'
+         strtrac(izone_cont)='con'  
+         strtrac(izone_contfrac)='enu'  ! evap sol nu
+      !========================================================================================================================
+      CASE(11)     !=== TRACING DROPLETS REEVAPORATION + REST
+      !========================================================================================================================
+         nzone_opt=2
+         izone_cont=1
+         izone_oce=1
+         izone_poubelle=1    ! zone ou on met les flux non physiques, de reajustement
+         izone_init=1        ! zone d'initialisation par defaut
+         option_revap=1
+         option_tmin=0
+         izone_revap=2
+         izone_irrig=0
+         option_cond=0
+         strtrac(izone_poubelle)='res'
+         strtrac(izone_revap)='rev'
+      !========================================================================================================================
+      CASE(12)     !=== TRACING NAKED GROUND EVAPORATION, CANOPY EVAPORATION, REST OF LAND EVAPORATION AND OCEAN EVAPORATION
+      !========================================================================================================================
+#ifndef CPP_VEGET
+         WRITE(*,*) 'iso_traceurs_init 257: option_traceurs=10 inutile si on ne couple pas avec ORCHIDEE'; STOP
 #endif           
-
-          ntraceurs_zone_opt=2
-          izone_cont=1
-          izone_oce=2
-          izone_poubelle=2 ! zone où on met les flux non physiques, de
-                ! réajustement
-          izone_init=2 ! zone d'initialisation par défaut
-          option_revap=0
-          option_tmin=0
-          izone_revap=0
-          izone_contfrac=3
-          izone_contcanop=4
-          izone_irrig=0   
-          option_cond=0
-
-          strtrac(izone_oce)='oce'
-          strtrac(izone_cont)='con'
-          strtrac(izone_contfrac)='enu'  ! evap sol nu
-          strtrac(izone_contcanop)='eca'  ! evap canop
-
-       else if (option_traceurs.eq.13) then
-          ! on trace les température minimales vécues + la revap
-          ! comme dans article sur LdG
-            
-        zone_temp1=293.0         ! en K        
-!        parameter (zone_tempf=223.0) ! en K
-        zone_tempf=243.0 ! en K
-        zone_tempa=-4.0 ! courbure de la relation entre l'indice et la température: 0 pour linéaire, <0 pour plus de détal en bas
-
-        ! zone 1: >= zone_temp1
-        ! zone 2 à 4: intermédiaire, 
-        ! zone 5: <zone_tempf
-        
-          ntraceurs_zone_opt=nzone_temp+1
-          
-          izone_cont=1
-          izone_oce=1  
-          izone_poubelle=1
-          izone_init=1 ! zone d'initialisation par défaut
-          option_revap=1   
-          option_tmin=0
-          izone_revap=ntraceurs_zone
-          izone_irrig=0
-          option_cond=0
-          do izone=1,nzone_temp
-            write(strz,'(i2.2)') izone
-            strtrac(izone)='t'//strz
-            write(*,*) 'izone,strz,strtrac=',izone,strz,strtrac(izone)
-          enddo
-          strtrac(izone_revap)='rev'
-
-          ! initialisation des zones de tempéarture
-          do izone=1,nzone_temp-1
-            zone_temp(izone)=zone_temp1+float(izone-1) &
-     &                      *(zone_tempa*float(izone-nzone_temp+1) &
-     &                      +(zone_tempf-zone_temp1)/float(nzone_temp-2))
-          enddo
-          write(*,*) 'zone_temp=',zone_temp
-
-       else if (option_traceurs.eq.14) then
-          ! on trace les pres et lat de dernière saturation définies
-          ! comme rh>90%
-            
-        zone_pres1=600.0*100.0 ! en Pa       
-        zone_presf=300.0*100.0 ! en Pa
-        zone_presa=0.0 ! courbure de la relation entre l'indice et la température: 0 pour linéaire, <0 pour plus de détal en bas
-
-        lattag_min=10.0 ! en degrès
-        dlattag=15.0
-
-        ! zone 1: >= zone_pres1
-        ! zone 2 à 4: intermédiaire, 
-        ! zone 5: <zone_presf
-        
-         ntraceurs_zone_opt=nzone_pres*nzone_lat+1          
-          izone_cont=ntraceurs_zone
-          izone_oce=ntraceurs_zone
-          izone_poubelle=ntraceurs_zone
-          izone_init=ntraceurs_zone ! zone d'initialisation par défaut
-          option_revap=0  
-          option_tmin=0 
-          izone_revap=0
-          izone_irrig=0
-          option_cond=0
-          do izone_pres=1,nzone_pres
-           do izone_lat=1,nzone_lat
-            write(strz_pres,'(i1.1)') izone_pres
-            write(strz_lat,'(i1.1)') izone_lat
-            strz_preslat=strz_pres//strz_lat
-            izone=izone_lat+(izone_pres-1)*nzone_lat
-            strtrac(izone)='t'//strz_preslat
-            write(*,*) 'izone_pres,izone_lat,strtrac=', &
-     &                        izone_pres,izone_lat,izone,strtrac(izone)
-           enddo !do izone_lat=1,nzone_lat
-          enddo !do izone_pres=1,nzone_pres
-          strtrac(ntraceurs_zone)='sfc'
-
-          ! initialisation des zones de tempéarture
-          do izone=1,nzone_pres-1
-            zone_pres(izone)=zone_pres1+float(izone-1) &
-     &                      *(zone_presa*float(izone-nzone_pres+1) &
-     &                      +(zone_presf-zone_pres1)/float(nzone_pres-2))
-          enddo !do izone=1,nzone_pres-1
-          write(*,*) 'traceurs_init 332: zone_pres=',zone_pres
-!          stop
-!
-       elseif (option_traceurs.eq.15) then
-          ! on trace l'irrigation dans ORCHIDEE
-#ifdef CPP_VEGET
-#else
-          write(*,*) 'iso_traceurs_init 257: option_traceurs=15 ', &
-     &                      'inutile si on ne couple pas avec ORCHIDEE'
-          stop
+         nzone_opt=2
+         izone_cont=1
+         izone_oce=2
+         izone_poubelle=2    ! zone ou on met les flux non physiques, de reajustement
+         izone_init=2        ! zone d'initialisation par defaut
+         option_revap=0
+         option_tmin=0
+         izone_revap=0
+         izone_contfrac=3
+         izone_contcanop=4
+         izone_irrig=0   
+         option_cond=0
+         strtrac(izone_oce)='oce'
+         strtrac(izone_cont)='con'
+         strtrac(izone_contfrac)='enu' ! evap sol nu
+         strtrac(izone_contcanop)='eca'! evap canop
+      !========================================================================================================================
+      CASE(13)     !=== TRACING MINIMUM EXPERIENCED TEMPERATIRES + REEVAPORATION AS IN THE ARTICLE ON LdG
+      !========================================================================================================================
+         zone_temp1=293.0    ! en K        
+!        zone_tempf=223.0    ! en K
+         zone_tempf=243.0    ! en K
+         zone_tempa=-4.0     ! courbure de la relation entre l'indice et la temperature: 0 pour lineaire, <0 pour plus de detal en bas
+         ! zone 1: >= zone_temp1
+         ! zone 2 a 4: intermediaire, 
+         ! zone 5: <zone_tempf
+         nzone_opt=nzone_temp+1
+         izone_cont=1
+         izone_oce=1  
+         izone_poubelle=1
+         izone_init=1        ! zone d'initialisation par defaut
+         option_revap=1   
+         option_tmin=0
+         izone_revap=nzone
+         izone_irrig=0
+         option_cond=0
+         DO izone=1,nzone_temp
+            strtrac(izone) = 't'//TRIM(int2str(izone))
+            WRITE(*,*) 'izone, strtrac = ', izone, strtrac(izone)
+         END DO
+         strtrac(izone_revap)='rev'
+         ! initialisation des zones de tempearture
+         DO izone=1,nzone_temp-1
+            zone_temp(izone) = zone_temp1+float(izone-1) &
+                             *(zone_tempa*float(izone-nzone_temp+1) &
+                             +(zone_tempf-zone_temp1)/float(nzone_temp-2))
+         END DO
+         WRITE(*,*) 'zone_temp=',zone_temp
+      !========================================================================================================================
+      CASE(14)     !=== TRACING PRES AND LAT OF LAST SATURATION DEFINED AS rh>90%
+      !========================================================================================================================
+         zone_pres1=600.0*100.0   ! en Pa       
+         zone_presf=300.0*100.0   ! en Pa
+         zone_presa=0.0           ! courbure de la relation entre l'indice et la temperature: 0 pour lineaire
+         lattag_min=10.0          ! en degres
+         dlattag=15.0
+         ! zone 1: >= zone_pres1
+         ! zone 2 a 4: intermediaire, 
+         ! zone 5: <zone_presf
+         nzone_opt=nzone_pres*nzone_lat+1          
+         izone_cont=nzone
+         izone_oce=nzone
+         izone_poubelle=nzone
+         izone_init=nzone         ! zone d'initialisation par defaut
+         option_revap=0  
+         option_tmin=0 
+         izone_revap=0
+         izone_irrig=0
+         option_cond=0
+         DO izone_pres=1,nzone_pres
+            DO izone_lat=1,nzone_lat
+               izone=izone_lat+(izone_pres-1)*nzone_lat
+               strtrac(izone) = 't'//TRIM(int2str(izone_pres))//TRIM(int2str(izone_lat))
+               write(*,*) 'izone_pres, izone_lat, izone, strtrac = ',izone_pres, izone_lat, izone, strtrac(izone)
+            END DO
+         END DO
+         strtrac(nzone)='sfc'
+         ! initialisation des zones de temperature
+         DO izone=1,nzone_pres-1
+            zone_pres(izone) = zone_pres1+float(izone-1) &
+                             *(zone_presa*float(izone-nzone_pres+1) &
+                             +(zone_presf-zone_pres1)/float(nzone_pres-2))
+         END DO
+         WRITE(*,*) 'traceurs_init 332: zone_pres=',zone_pres
+      !========================================================================================================================
+      CASE(15)     !=== TRACING IRRIGATION IN ORCHIDEE
+      !========================================================================================================================
+#ifndef CPP_VEGET
+         WRITE(*,*) 'iso_traceurs_init 257: option_traceurs=15 inutile si on ne couple pas avec ORCHIDEE'; STOP
 #endif
-
-          ntraceurs_zone_opt=1
-          izone_cont=1
-          izone_oce=1
-          izone_poubelle=1 ! zone où on met les flux non physiques, de
-                ! réajustement
-          izone_init=1 ! zone d'initialisation par défaut
-          option_revap=0
-          option_tmin=0
-          izone_revap=0
-          izone_contfrac=0
-          izone_contcanop=0 
-          izone_irrig=2 
-          option_cond=0
-          
-          strtrac(izone_poubelle)='res'
-          strtrac(izone_irrig)='irrig'
-
-          ! dans ce cas particulier, il y a des traceurs dans ORCHIDEE
-          ntracisoOR=ntiso
-
-        else if ((option_traceurs.eq.17).or. &
-     &           (option_traceurs.eq.18)) then
-          ! on trace les température minimales vécues
-          ! comme dans article sur LdG sauf pas de revap
-            
-        zone_temp1=12.0e-3 ! en kg/kg        
-        zone_tempf=0.2e-3 ! en kg/kg
-        zone_tempa=1.2e-3 ! courbure de la relation entre l'indice et la température: 0 pour linéaire, <0 pour plus de détail en bas
-
-!	parameter (zone_temp1=14.0e-3) ! en kg/kg        
-!   	parameter (zone_tempf=0.2e-3) ! en kg/kg
-!	parameter (zone_tempa=0.5e-3)       
-
-!        parameter (zone_temp1=10.0e-3) ! en kg/kg
-!   	parameter (zone_tempf=0.5e-3) ! en kg/kg
-!	parameter (zone_tempa=0.5e-3)
-
-        ! zone 1: >= zone_temp1
-        ! zone 2 à 4: intermédiaire, 
-        ! zone 5: <zone_tempf
-        
-        ntraceurs_zone_opt=nzone_temp+3
-        
-          izone_cont=nzone_temp+1
-          izone_oce=nzone_temp+1 
-          izone_poubelle=nzone_temp+1
-          izone_init=nzone_temp+1 ! zone d'initialisation par défaut
-          option_revap=1  
-          option_tmin=1 
-          option_cond=1
-
-          izone_revap=nzone_temp+3
-          izone_cond=nzone_temp+2
-          do izone=1,nzone_temp
-            write(strz,'(i2.2)') izone
-            strtrac(izone)='t'//strz
-            write(*,*) 'izone,strz,strtrac=',izone,strz,strtrac(izone)
-          enddo !do izone=1,nzone_temp
-          strtrac(izone_poubelle)='sfc'
-          strtrac(izone_cond)='con'
-          strtrac(izone_revap)='rev'
-
-          ! initialisation des zones de tempéarture
-          do izone=1,nzone_temp-1
-            zone_temp(izone)=zone_temp1+float(izone-1) &
-     &                      *(zone_tempa*float(izone-nzone_temp+1) &
-     &             +(zone_tempf-zone_temp1)/float(nzone_temp-2))
-          enddo
-         write(*,*) 'zone_temp1,zone_tempf,zone_tempa=', &
-     &              zone_temp1,zone_tempf,zone_tempa
-          write(*,*) 'zone_temp=',zone_temp
-!          stop          
-
-        else if (option_traceurs.eq.19) then
-
-        zone_temp1=12.0e-3 ! en kg/kg        
-        zone_tempf=0.2e-3 ! en kg/kg
-        zone_tempa=1.2e-3 ! courbure de la relation entre l'indice et la température: 0 pour linéaire, <0 pour plus de détail en bas
-
-!	parameter (zone_temp1=14.0e-3) ! en kg/kg        
-!   	parameter (zone_tempf=0.2e-3) ! en kg/kg
-!	parameter (zone_tempa=0.5e-3)       
-
-!        parameter (zone_temp1=10.0e-3) ! en kg/kg
-!   	parameter (zone_tempf=0.5e-3) ! en kg/kg
-!	parameter (zone_tempa=0.5e-3)
-
-        ! zone 1: >= zone_temp1
-        ! zone 2 à 4: intermédiaire, 
-        ! zone 5: <zone_tempf
-        
-        ntraceurs_zone_opt=nzone_temp+4
-        
-          izone_cont=nzone_temp+1
-          izone_oce=nzone_temp+1 
-          izone_poubelle=nzone_temp+1
-          if (option_seuil_tag_tmin.eq.1) then
-            izone_init=nzone_temp+1 ! zone d'initialisation par défaut
-          else
+         nzone_opt=1
+         izone_cont=1
+         izone_oce=1
+         izone_poubelle=1    ! zone ou on met les flux non physiques, de reajustement
+         izone_init=1        ! zone d'initialisation par defaut
+         option_revap=0
+         option_tmin=0
+         izone_revap=0
+         izone_contfrac=0
+         izone_contcanop=0 
+         izone_irrig=2 
+         option_cond=0
+         strtrac(izone_poubelle)='res'
+         strtrac(izone_irrig)='irrig'
+         ! dans ce cas particulier, il y a des traceurs dans ORCHIDEE
+         ntracisoOR=ntiso
+      !========================================================================================================================
+      CASE(17,18)  !=== TRACING MINIMAL EXPERIENCES TEMPERATURES AS IN THE ARTICLE ABOUT LdG, BUT NO EVAPORATION
+      !========================================================================================================================
+         zone_temp1=12.0e-3  ! en kg/kg        
+         zone_tempf=0.2e-3   ! en kg/kg
+         zone_tempa=1.2e-3   ! courbure de la relation entre l'indice et la temperature: 0 pour lineaire
+!	 zone_temp1=14.0e-3  ! en kg/kg        
+!   	 zone_tempf=0.2e-3   ! en kg/kg
+!	 zone_tempa=0.5e-3       
+!        zone_temp1=10.0e-3  ! en kg/kg
+!   	 zone_tempf=0.5e-3   ! en kg/kg
+!	 zone_tempa=0.5e-3
+         ! zone 1: >= zone_temp1
+         ! zone 2 a 4: intermediaire, 
+         ! zone 5: <zone_tempf
+         nzone_opt=nzone_temp+3
+         izone_cont=nzone_temp+1
+         izone_oce=nzone_temp+1 
+         izone_poubelle=nzone_temp+1
+         izone_init=nzone_temp+1 ! zone d'initialisation par defaut
+         option_revap=1  
+         option_tmin=1 
+         option_cond=1
+         izone_revap=nzone_temp+3
+         izone_cond=nzone_temp+2
+         DO izone=1,nzone_temp
+            strtrac(izone) = 't'//TRIM(int2str(izone))
+            WRITE(*,*) 'izone, strtrac = ', izone, strtrac(izone)
+         END DO !do izone=1,nzone_temp
+         strtrac(izone_poubelle)='sfc'
+         strtrac(izone_cond)='con'
+         strtrac(izone_revap)='rev'
+         ! initialisation des zones de tempearture
+         DO izone=1,nzone_temp-1
+            zone_temp(izone) = zone_temp1+float(izone-1) &
+                             *(zone_tempa*float(izone-nzone_temp+1) &
+                             +(zone_tempf-zone_temp1)/float(nzone_temp-2))
+         END DO
+         WRITE(*,*) 'zone_temp1,zone_tempf,zone_tempa=',zone_temp1,zone_tempf,zone_tempa
+         WRITE(*,*) 'zone_temp=',zone_temp
+!        STOP          
+      !========================================================================================================================
+      CASE(19)     !=== TRACING TROPICAL AND EXTRATROPICAL VAPOUR
+      !========================================================================================================================
+         zone_temp1=12.0e-3  ! en kg/kg        
+         zone_tempf=0.2e-3   ! en kg/kg
+         zone_tempa=1.2e-3   ! courbure de la relation entre l'indice et la temperature: 0 pour lineaire, <0 pour plus de detail en bas
+!        zone_temp1=14.0e-3  ! en kg/kg        
+!        zone_tempf=0.2e-3   ! en kg/kg
+!        zone_tempa=0.5e-3
+!        zone_temp1=10.0e-3  ! en kg/kg        
+!        zone_tempf=0.5e-3   ! en kg/kg
+!        zone_tempa=0.5e-3
+         ! zone 1: >= zone_temp1
+         ! zone 2 a 4: intermediaire, 
+         ! zone 5: <zone_tempf
+         nzone_opt=nzone_temp+4
+         izone_cont=nzone_temp+1
+         izone_oce=nzone_temp+1 
+         izone_poubelle=nzone_temp+1
+         IF(option_seuil_tag_tmin == 1) THEN
+            izone_init=nzone_temp+1 ! zone d'initialisation par defaut
+         ELSE
             izone_init=nzone_temp
-          endif
-          option_revap=1   
-          izone_revap=nzone_temp+3
-          izone_cond=nzone_temp+2
-          izone_ddft=nzone_temp+4
-          option_tmin=1          
-          option_cond=1
-          do izone=1,nzone_temp
-            write(strz,'(i2.2)') izone
-            strtrac(izone)='t'//strz
-            write(*,*) 'izone,strz,strtrac=',izone,strz,strtrac(izone)
-          enddo !do izone=1,nzone_temp
-          strtrac(izone_poubelle)='sfc'
-          strtrac(izone_cond)='con'
-          strtrac(izone_revap)='rev'
-          strtrac(izone_ddft)='dft'
-
-        elseif (option_traceurs.eq.20) then
-          ! on vapeur tropical/extractropicale/recyclage extractropical
-          ! pour comprendre controles humidité et isotopes subtropicaux.        
-          
-          lim_tag20=35.0
-          call getin('lim_tag20',lim_tag20)
-          write(*,*) 'lim_tag20=',lim_tag20
-
-          ntraceurs_zone_opt=3
-          izone_cont=1
-          izone_oce=1
-          izone_poubelle=2 ! zone où on met les flux non physiques, de
-                ! réajustement
-          izone_init=2 ! zone d'initialisation par défaut
-          option_revap=0
-          option_tmin=0
-          izone_revap=0
-          izone_trop=2
-          izone_extra=3
-
-          strtrac(izone_trop)='tro' ! vapeur tropicale
-          strtrac(izone_extra)='ext' ! vapeur extractropicale evaporée 
-                ! dans les tropiques
-          strtrac(izone_cont)='rec' ! recyclage
-
-        elseif (option_traceurs.eq.21) then
-          ! on trace 2 boites 3D: UT tropicale et extratropiques 
-          ! fonctionnement similaire à option 5 pour taggage des zones
-          ! AMMA
-!          write(*,*) 'iso_traceurs_init 129'
-
-          ntraceurs_zone_opt=3
-          izone_cont=1
-          izone_oce=1
-          izone_poubelle=1 ! zone où on met les flux non physiques, de
-                ! réajustement
-          izone_init=1 ! zone d'initialisation par défaut
-          option_revap=0
-          option_tmin=0
-          izone_revap=0
-          izone_trop=2
-          izone_extra=3
-          option_cond=0
-
-          strtrac(izone_poubelle)='res'
-          strtrac(izone_trop)='tro'
-          strtrac(izone_extra)='ext'
-
-        elseif (option_traceurs.eq.22) then
-          ! on trace la vapeur qui a été processée dans zones de
-          ! convections à 3 niveaux: BT, MT et UT
-
-          lim_precip_tag22=20.0
-          call getin('lim_precip_tag22',lim_precip_tag22)
-          write(*,*) 'lim_precip_tag22=',lim_precip_tag22
-
-          ntraceurs_zone_opt=3
-          izone_cont=1
-          izone_oce=1
-          izone_poubelle=1 ! zone où on met les flux non physiques, de
-                ! réajustement
-          izone_init=1 ! zone d'initialisation par défaut
-          option_revap=0
-          option_tmin=0
-          izone_revap=0
-          izone_conv_BT=2
-          izone_conv_UT=3
-          option_cond=0
-
-          strtrac(izone_poubelle)='res'
-          strtrac(izone_conv_BT)='cbt'
-          strtrac(izone_conv_UT)='cut'
-
-        else
-            write(*,*) 'traceurs_init 36: option pas encore prévue'
-            stop
-        endif
-
-        
-          if (ntraceurs_zone_opt.ne.ntraceurs_zone) then
-                write(*,*) 'ntraceurs_zone_opt,ntraceurs_zone=', &
-                        & ntraceurs_zone_opt,ntraceurs_zone
-                call abort_physic ('isotrac_mod','ntraceurs_zone incoherent',1)
-          endif
-
-       
-        ! seuil sur le taux de condensation
-        if (option_tmin.eq.1) then
-          seuil_tag_tmin=0.01
-          call getin('seuil_tag_tmin',seuil_tag_tmin)
-          write(*,*) 'seuil_tag_tmin=',seuil_tag_tmin
-
-          seuil_tag_tmin_ls=seuil_tag_tmin
-          call getin('seuil_tag_tmin_ls',seuil_tag_tmin_ls)
-          write(*,*) 'seuil_tag_tmin_ls=',seuil_tag_tmin_ls
-
-          option_seuil_tag_tmin=1
-          call getin('option_seuil_tag_tmin',option_seuil_tag_tmin)
-          write(*,*) 'option_seuil_tag_tmin=',option_seuil_tag_tmin
-        endif
-
-
-        do ixt=1,niso
-           index_zone(ixt)=0
-           index_iso(ixt)=ixt
-        enddo
-        itrac=niso        
-        do izone=1,ntraceurs_zone
-          do ixt=1,niso
-            itrac=itrac+1
-            index_zone(itrac)=izone
-            index_iso(itrac)=ixt
-            itZonIso_loc(izone,ixt)=itrac
-            if (itZonIso(izone,ixt).ne.itZonIso_loc(izone,ixt)) then
-                write(*,*) 'isotrac 989: izone,ixt,itrac=',izone,ixt,itrac
-                CALL abort_physic ('isotrac','isotrac 989',1)
-            endif
-          enddo
-        enddo
+         END IF
+         option_revap=1   
+         izone_revap=nzone_temp+3
+         izone_cond=nzone_temp+2
+         izone_ddft=nzone_temp+4
+         option_tmin=1          
+         option_cond=1
+         DO izone=1,nzone_temp
+            strtrac(izone) = 't'//TRIM(int2str(izone))
+            WRITE(*,*) 'izone, strtrac = ', izone, strtrac(izone)
+         END DO
+         strtrac(izone_poubelle)='sfc'
+         strtrac(izone_cond)='con'
+         strtrac(izone_revap)='rev'
+         strtrac(izone_ddft)='dft'
+      !========================================================================================================================
+      CASE(20)     !=== TRACING TROPICAL/EXTRATROPICAL/EXTRATROPICAL RECYCLING TO STUDY HUMIDITY AND SUBTROPICAL ISOTOPES CONTROL
+      !========================================================================================================================
+         CALL get_in('lim_tag20', lim_tag20, 35.0)
+         nzone_opt=3
+         izone_cont=1
+         izone_oce=1
+         izone_poubelle=2    ! zone ou on met les flux non physiques, de reajustement
+         izone_init=2        ! zone d'initialisation par defaut
+         option_revap=0
+         option_tmin=0
+         izone_revap=0
+         izone_trop=2
+         izone_extra=3
+         strtrac(izone_trop)='tro'     ! tropical vapour
+         strtrac(izone_extra)='ext'    ! extratropical vapour evaporated in the tropics
+         strtrac(izone_cont)='rec'     ! recycling
+      !========================================================================================================================
+      CASE(21)     !=== TRACING TWO 3D BOXES: TROPICAL UT AND EXTRATROPICS ; SIMILAR TO 5 FOR AMMA ZONES TAGGING
+      !========================================================================================================================
+!        WRITE(*,*) 'iso_traceurs_init 129'
+         nzone_opt=3
+         izone_cont=1
+         izone_oce=1
+         izone_poubelle=1    ! zone ou on met les flux non physiques, de reajustement
+         izone_init=1        ! zone d'initialisation par defaut
+         option_revap=0
+         option_tmin=0
+         izone_revap=0
+         izone_trop=2
+         izone_extra=3
+         option_cond=0
+         strtrac(izone_poubelle)='res'
+         strtrac(izone_trop)='tro'
+         strtrac(izone_extra)='ext'
+      !========================================================================================================================
+      CASE(22)     !=== TRACING WATER VAPOUR PROCESSED IN THE 3-LEVELS SCONVECTION ZONES BT, MT AND UT
+      !========================================================================================================================
+         CALL get_in('lim_precip_tag22', lim_precip_tag22, 20.0)
+         nzone_opt=3
+         izone_cont=1
+         izone_oce=1
+         izone_poubelle=1    ! zone ou on met les flux non physiques, de reajustement
+         izone_init=1        ! zone d'initialisation par defaut
+         option_revap=0
+         option_tmin=0
+         izone_revap=0
+         izone_conv_BT=2
+         izone_conv_UT=3
+         option_cond=0
+         strtrac(izone_poubelle)='res'
+         strtrac(izone_conv_BT)='cbt'
+         strtrac(izone_conv_UT)='cut'
+      CASE DEFAULT
+         WRITE(*,*) 'traceurs_init 36: option pas encore prevue' ; STOP
+   END SELECT
+
+   IF(nzone_opt /= nzone) THEN
+      WRITE(*,*) 'nzone_opt, nzone=', nzone_opt, nzone
+      CALL abort_physic ('isotrac_mod','nzone incoherent',1)
+   END IF
+
+   !--- Condensation rate threshold
+   IF(option_tmin == 1) THEN
+      seuil_tag_tmin = 0.01
+      CALL get_in('seuil_tag_tmin',        seuil_tag_tmin,        0.01)
+      CALL get_in('seuil_tag_tmin_ls',     seuil_tag_tmin_ls,     seuil_tag_tmin)
+      CALL get_in('option_seuil_tag_tmin', option_seuil_tag_tmin, 1)
+   END IF
+   DO ixt=1,niso
+      index_zone(ixt)=0
+      index_iso(ixt)=ixt
+   END DO
+
+   index_zone = [(INDEX(isoZone, strTail(         isoName(ixt) ,'_')), ixt=1, ntiso)]
+   index_iso  = [(INDEX(isoName, strHead(delPhase(isoName(ixt)),'_')), ixt=1, ntiso)]
+   itZonIso_loc = itZonIso(:,:)
 #ifdef ISOVERIF
-!        call iso_verif_egalite(float(itrac),float(ntiso), &
-!     &           'traceurs_init 50')
-        if (itrac.ne.ntiso) then
-          write(*,*) 'traceurs_init 50'
-          stop
-        endif
-     
-        write(*,*) 'traceurs_init 65: bilan de l''init:'
-        write(*,*) 'index_zone=',index_zone(1:ntiso)
-        write(*,*) 'index_iso=',index_iso(1:ntiso)
-        write(*,*) 'itZonIso=',itZonIso(1:ntraceurs_zone,1:niso)
-        do izone=1,ntraceurs_zone
-          write(*,*) 'strtrac(',izone,')=',strtrac(izone)
-        enddo !do izone=1,ntraceurs_zone
-        write(*,*) 'ntracisoOR=',ntracisoOR
+   WRITE(*,*) 'traceurs_init 65: bilan de l''init:'
+   WRITE(*,*) 'index_zone = '//TRIM(strStack(int2str(index_zone(1:ntiso))))
+   WRITE(*,*) 'index_iso  = '//TRIM(strStack(int2str(index_iso (1:ntiso))))
+   DO izone=1,nzone
+      WRITE(*,*)'itZonIso('//TRIM(int2str(izone))//',:) = '//strStack(int2str(itZonIso(izone,:)))
+   END DO
+   DO izone=1,nzone
+      WRITE(*,*)'strtrac('//TRIM(int2str(izone))//',:) = '//TRIM(strtrac(izone))
+   END DO
+   WRITE(*,*) 'ntracisoOR=',ntracisoOR
 #endif  
 
-        end subroutine iso_traceurs_init
-
+END SUBROUTINE iso_traceurs_init
 
 END MODULE isotrac_mod
Index: LMDZ6/trunk/libf/phylmdiso/isotrac_routines_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/isotrac_routines_mod.F90	(revision 4319)
+++ LMDZ6/trunk/libf/phylmdiso/isotrac_routines_mod.F90	(revision 4325)
@@ -1113,10 +1113,10 @@
       USE isotrac_mod, only: use_bassin_atlantic,use_bassin_medit, &
 &       use_bassin_indian,use_bassin_austral,use_bassin_pacific, &
-&       use_bassin_merarabie,use_bassin_golfebengale,use_bassin_indiansud, &
-&       use_bassin_tropics,use_bassin_midlats,use_bassin_hauteslats, &
+&       use_bassin_MerArabie,use_bassin_BengalGolf,use_bassin_SouthIndian, &
+&       use_bassin_tropics,use_bassin_midlats,use_bassin_HighLats, &
 &       bassin_atlantic,bassin_medit, &
 &       bassin_indian,bassin_austral,bassin_pacific, &
-&       bassin_merarabie,bassin_golfebengale,bassin_indiansud, &
-&       bassin_tropics,bassin_midlats,bassin_hauteslats
+&       bassin_MerArabie,bassin_BengalGolf,bassin_SouthIndian, &
+&       bassin_tropics,bassin_midlats,bassin_HighLats
       implicit none
       ! répond true si lat,lon se trouve dans le bassin numéroté bassin
@@ -1137,6 +1137,5 @@
       write(*,*) 'is_in_basin 84: entree,bassin=',bassin
 #endif
-      if ((use_bassin_atlantic.eq.1).and. &
-     &           (bassin.eq.bassin_atlantic)) then
+      if (use_bassin_atlantic .and. bassin==bassin_atlantic) then
 #ifdef ISOVERIF           
           write(*,*) 'bassin Atlantique?'
@@ -1169,6 +1168,5 @@
           endif
 
-      else if ((use_bassin_medit.eq.1).and. &
-     &           (bassin.eq.bassin_medit)) then
+      else if (use_bassin_medit .and. bassin==bassin_medit) then
 #ifdef ISOVERIF           
           write(*,*) 'bassin Medit?'
@@ -1183,6 +1181,5 @@
           endif
 
-      else if ((use_bassin_indian.eq.1).and. &
-     &           (bassin.eq.bassin_indian)) then
+      else if (use_bassin_indian .and. bassin==bassin_indian) then
 #ifdef ISOVERIF           
           write(*,*) 'bassin indian?'
@@ -1199,6 +1196,5 @@
           endif   
 
-      else if ((use_bassin_indiansud.eq.1).and. &
-     &           (bassin.eq.bassin_indiansud)) then
+      else if (use_bassin_SouthIndian .and. bassin==bassin_SouthIndian) then
 #ifdef ISOVERIF           
           write(*,*) 'bassin indian hemisphere Sud?'
@@ -1209,6 +1205,5 @@
           endif
           
-      else if ((use_bassin_merarabie.eq.1).and. &
-     &           (bassin.eq.bassin_merarabie)) then
+      else if (use_bassin_MerArabie .and. bassin==bassin_MerArabie) then
 #ifdef ISOVERIF           
           write(*,*) 'bassin Mer d''Arabie?'
@@ -1219,6 +1214,5 @@
           endif 
 
-      else if ((use_bassin_golfebengale.eq.1).and. &
-     &           (bassin.eq.bassin_golfebengale)) then
+      else if (use_bassin_BengalGolf .and. bassin==bassin_BengalGolf) then
 #ifdef ISOVERIF           
           write(*,*) 'bassin Golfe du Bengale?'
@@ -1229,6 +1223,5 @@
           endif         
 
-      else if ((use_bassin_pacific.eq.1).and. &
-     &           (bassin.eq.bassin_pacific)) then
+      else if (use_bassin_pacific .and. bassin==bassin_pacific) then
 #ifdef ISOVERIF           
           write(*,*) 'bassin Pacific?'
@@ -1267,6 +1260,5 @@
           endif
 
-      else if ((use_bassin_austral.eq.1).and. &
-     &           (bassin.eq.bassin_austral)) then  
+      else if (use_bassin_austral .and. bassin==bassin_austral) then  
 #ifdef ISOVERIF           
           write(*,*) 'bassin austral?'
@@ -1277,6 +1269,5 @@
           endif  
 
-      else if ((use_bassin_hauteslats.eq.1).and. &
-     &           (bassin.eq.bassin_hauteslats)) then  
+      else if (use_bassin_HighLats .and. bassin==bassin_HighLats) then  
 #ifdef ISOVERIF           
           write(*,*) 'bassin hautes lats?'
@@ -1287,6 +1278,5 @@
           endif
 
-      else if ((use_bassin_tropics.eq.1).and. &
-     &           (bassin.eq.bassin_tropics)) then  
+      else if (use_bassin_tropics .and. bassin==bassin_tropics) then  
 #ifdef ISOVERIF           
           write(*,*) 'bassin tropics?'
@@ -1297,6 +1287,5 @@
           endif
 
-       else if ((use_bassin_midlats.eq.1).and. &
-     &           (bassin.eq.bassin_midlats)) then  
+       else if (use_bassin_midlats .and. bassin==bassin_midlats) then  
 #ifdef ISOVERIF           
           write(*,*) 'bassin mid lats?'
@@ -1314,14 +1303,14 @@
           write(*,*) 'bassin_indian=' ,bassin_indian
           write(*,*) 'bassin_austral=' ,bassin_austral
-          write(*,*) 'bassin_merarabie=' ,bassin_merarabie
-          write(*,*) 'bassin_golfebengale=' ,bassin_golfebengale
-          write(*,*) 'bassin_indiansud=' ,bassin_indiansud
+          write(*,*) 'bassin_MerArabie=' ,bassin_MerArabie
+          write(*,*) 'bassin_BengalGolf=' ,bassin_BengalGolf
+          write(*,*) 'bassin_SouthIndian=' ,bassin_SouthIndian
           write(*,*) 'use_bassin_atlantic=' ,use_bassin_atlantic  
           write(*,*) 'use_bassin_medit=' ,use_bassin_medit
           write(*,*) 'use_bassin_indian=' ,use_bassin_indian
           write(*,*) 'use_bassin_austral=' ,use_bassin_austral
-          write(*,*) 'use_bassin_merarabie=' ,use_bassin_merarabie
-          write(*,*) 'use_bassin_golfebengale=' ,use_bassin_golfebengale
-          write(*,*) 'use_bassin_indiansud=' ,use_bassin_indiansud
+          write(*,*) 'use_bassin_MerArabie=' ,use_bassin_MerArabie
+          write(*,*) 'use_bassin_BengalGolf=' ,use_bassin_BengalGolf
+          write(*,*) 'use_bassin_SouthIndian=' ,use_bassin_SouthIndian
           stop
       endif
