Index: /readTracFiles_mod.f90
===================================================================
--- /readTracFiles_mod.f90	(revision 24)
+++ /readTracFiles_mod.f90	(revision 25)
@@ -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,61 +1802,73 @@
 !=== 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
 !==============================================================================================================================
 
 
 END MODULE readTracFiles_mod
+
Index: /trac_types_mod.F90
===================================================================
--- /trac_types_mod.F90	(revision 24)
+++ /trac_types_mod.F90	(revision 25)
@@ -28,5 +28,5 @@
     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               :: 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(:)
