Index: LMDZ6/branches/LMDZ_ECRad/libf/misc/readTracFiles_mod.f90
===================================================================
--- LMDZ6/branches/LMDZ_ECRad/libf/misc/readTracFiles_mod.f90	(revision 4203)
+++ LMDZ6/branches/LMDZ_ECRad/libf/misc/readTracFiles_mod.f90	(revision 4482)
@@ -1,7 +1,6 @@
 MODULE readTracFiles_mod
 
-  USE strings_mod,    ONLY: msg, testFile,  strFind, strStack, strReduce,  strHead, strCount, find, fmsg, reduceExpr, &
-             removeComment, cat, checkList, str2int, strParse, strReplace, strTail, strIdx, maxlen, test, dispTable, get_in
-  USE trac_types_mod, ONLY: trac_type, isot_type, keys_type
+  USE strings_mod,    ONLY: msg, find, get_in, str2int, dispTable, strHead,  strReduce,  strFind, strStack, strIdx, &
+       test, removeComment, cat, fmsg, maxlen, int2str, checkList, strParse, strReplace, strTail, strCount, reduceExpr
 
   IMPLICIT NONE
@@ -9,34 +8,103 @@
   PRIVATE
 
-  PUBLIC :: initIsotopes, maxlen, trac_type, isot_type, keys_type
-  PUBLIC :: readTracersFiles, indexUpdate, setGeneration             !--- TOOLS ASSOCIATED TO TRACERS  DESCRIPTORS
-  PUBLIC :: readIsotopesFile                                         !--- TOOLS ASSOCIATED TO ISOTOPES DESCRIPTORS
-  PUBLIC :: getKey_init, getKey, fGetKey, setDirectKeys              !--- GET/SET KEYS FROM/TO tracers & isotopes
-
-  PUBLIC :: addPhase, new2oldName,  getPhase, &                      !--- FUNCTIONS RELATED TO THE PHASES
-            delPhase, old2newName, getiPhase, &                      !--- + ASSOCIATED VARIABLES
-            known_phases, old_phases, phases_sep, phases_names, nphases
-
-  PUBLIC :: oldH2OIso, newH2OIso                                     !--- NEEDED FOR BACKWARD COMPATIBILITY (OLD traceur.def)
-
-  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, fGetKeys, addKey, 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, 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               :: iGeneration = -1              !--- Generation number (>=0)
+    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)
+    TYPE(keys_type)       :: keys                          !--- <key>=<val> pairs vector
+    INTEGER               :: iadv        = 10              !--- Advection scheme used
+    LOGICAL               :: isAdvected  = .FALSE.         !--- "true" tracers: iadv > 0.   COUNT(isAdvected )=nqtrue
+    LOGICAL               :: isInPhysics = .TRUE.          !--- "true" tracers: in tr_seri. COUNT(isInPhysics)=nqtottr
+    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
+  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
+    MODULE PROCEDURE getKeyByName_s1, getKeyByName_s1m, getKeyByName_sm, getKey_sm, &
+                     getKeyByName_i1, getKeyByName_i1m, getKeyByName_im, getKey_im, &
+                     getKeyByName_r1, getKeyByName_r1m, getKeyByName_rm, getKey_rm, &
+                     getKeyByName_l1, getKeyByName_l1m, getKeyByName_lm, getKey_lm
   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;        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; END INTERFACE addKey!,      addKey_m,     addKey_mm;     END INTERFACE addKey
   INTERFACE    addPhase;   MODULE PROCEDURE   addPhase_s1,   addPhase_sm,   addPhase_i1,   addPhase_im; END INTERFACE addPhase
-  INTERFACE old2newName;   MODULE PROCEDURE old2newName_1, old2newName_m; END INTERFACE old2newName
-  INTERFACE new2oldName;   MODULE PROCEDURE new2oldName_1, new2oldName_m; END INTERFACE new2oldName
 !------------------------------------------------------------------------------------------------------------------------------
 
@@ -45,24 +113,42 @@
 
   !--- SOME PARAMETERS THAT ARE NOT LIKELY TO CHANGE OFTEN
-  CHARACTER(LEN=maxlen), SAVE      :: tran0        = 'air'           !--- Default transporting fluid
-  CHARACTER(LEN=maxlen), PARAMETER :: old_phases   = 'vlir'          !--- Old phases for water (no separator)
-  CHARACTER(LEN=maxlen), PARAMETER :: known_phases = 'glsr'          !--- Known phases initials
-  INTEGER,               PARAMETER :: nphases=LEN_TRIM(known_phases) !--- Number of phases
-  CHARACTER(LEN=maxlen), SAVE      :: phases_names(nphases) &        !--- Known phases names
-                                = ['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
-
-  !--- KEPT JUST TO MANAGE OLD WATER ISOTOPES NAMES
-  !--- Apart from that context, on limitaion on isotopes names (as long as they have a corresponding line in isotopes_params.def)
-  CHARACTER(LEN=maxlen), SAVE :: oldH2OIso(5) = ['eau',     'HDO',     'O18',     'O17',     'HTO'    ]
-  CHARACTER(LEN=maxlen), SAVE :: newH2OIso(5) = ['H2[16]O', 'H[2]HO ', 'H2[18]O', 'H2[17]O', 'H[3]HO ']
-
-  !=== LOCAL COPIES OF THE TRACERS AND ISOTOPES DESCRIPTORS, USED BY getKey (INITIALIZED IN getKey_init)
+  CHARACTER(LEN=maxlen), SAVE      :: tran0        = 'air'      !--- Default transporting fluid
+  CHARACTER(LEN=maxlen), PARAMETER :: old_phases   = 'vlir'     !--- Old phases for water (no separator)
+  CHARACTER(LEN=maxlen), PARAMETER :: known_phases = 'glsr'     !--- Known phases initials
+  INTEGER, PARAMETER :: nphases = LEN_TRIM(known_phases)        !--- Number of phases
+  CHARACTER(LEN=maxlen), SAVE      :: phases_names(nphases) &   !--- Known phases names
+                                             = ['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
+  CHARACTER(LEN=maxlen), SAVE :: oldH2OIso(5) = ['eau',   'HDO',   'O18',   'O17',   'HTO'  ]
+  CHARACTER(LEN=maxlen), SAVE :: newH2OIso(5) = ['H216O', 'HDO  ', 'H218O', 'H217O', 'HTO  ']
+
+  !--- CORRESPONDANCE BETWEEN OLD AND NEW HNO3 RELATED SPECIES NAMES
+  CHARACTER(LEN=maxlen), SAVE ::   oldHNO3(2) = ['HNO3_g ', 'HNO3   ']
+  CHARACTER(LEN=maxlen), SAVE ::   newHNO3(2) = ['HNO3   ', 'HNO3tot']
+
+  !=== 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
 
@@ -72,5 +158,5 @@
 !==============================================================================================================================
 !=== READ ONE OR SEVERAL TRACER FILES AND FILL A "tr" TRACERS DESCRIPTOR DERIVED TYPE.
-!=== THE RETURN VALUE fType DEPENDS ON WHAT IS FOUND:
+!=== THE RETURNED VALUE fType DEPENDS ON WHAT IS FOUND:
 !===  0: NO ADEQUATE FILE FOUND ; DEFAULT VALUES MUST BE USED
 !===  1: AN "OLD STYLE" TRACERS FILE "traceur.def":
@@ -93,47 +179,27 @@
 !     * 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, fType, tracs) RESULT(lerr)
-!------------------------------------------------------------------------------------------------------------------------------
-  CHARACTER(LEN=*),             INTENT(IN)  :: type_trac              !--- List of components used
-  INTEGER,                      INTENT(OUT) :: fType                  !--- Type of input file found
-  TYPE(trac_type), ALLOCATABLE, INTENT(OUT) :: tracs(:)
+LOGICAL FUNCTION readTracersFiles(type_trac, lRepr) RESULT(lerr)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=*),  INTENT(IN)  :: type_trac                        !--- List of components used
+  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
-  INTEGER               :: is, nsec, ierr, it, ntrac, ns, ip, ix
-  LOGICAL, ALLOCATABLE  :: ll(:), lGen3(:)
+  CHARACTER(LEN=maxlen) :: str, fname, tname, pname, cname
+  INTEGER               :: nsec, ierr, it, ntrac, ns, ip, ix, fType
+  LOGICAL :: lRep
+  TYPE(keys_type), POINTER :: k
 !------------------------------------------------------------------------------------------------------------------------------
   lerr = .FALSE.
   modname = 'readTracersFiles'
   IF(.NOT.ALLOCATED(dBase)) ALLOCATE(dBase(0))
-
-  !--- 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
-  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)
+  lRep=.FALSE.; IF(PRESENT(lRepr)) lRep = lRepr
+
+  !--- Required sections + corresponding files names (new style single section case) for tests
+  IF(test(testTracersFiles(modname, type_trac, fType, .FALSE., trac_files, sections), lerr)) RETURN
+  nsec = SIZE(sections)
 
   !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SELECT CASE(fType)                         !--- Set %name, %genOName, %parent, %type, %phase, %component, %iGeneration, %keys
   !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-    CASE(1)                                                               !=== OLD FORMAT "traceur.def"
+    CASE(1)                                                          !=== OLD FORMAT "traceur.def"
     !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       !--- OPEN THE "traceur.def" FILE
@@ -145,27 +211,54 @@
 
       !--- 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
         IF(test(fmsg('Invalid format for "' //TRIM(fname)//'"', modname, ierr>0), lerr)) RETURN
         IF(test(fmsg('Not enough lines in "'//TRIM(fname)//'"', modname, ierr<0), lerr)) RETURN
-        ll = strParse(str, ' ', s, n=ns)
+        lerr = strParse(str, ' ', s, ns)
         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)
-        tracs(it)%name   = old2newName(s(3), ip)                     !--- Set %name:   name   of the tracer
-        tracs(it)%parent = tran0                                     !--- Default transporting fluid name
-        IF(ns == 4) tracs(it)%parent = old2newName(s(4))             !--- Set %parent: parent of the tracer
-        tracs(it)%phase = known_phases(ip:ip)                        !--- Set %phase:  tracer phase (default: "g"azeous)
-        tracs(it)%component = TRIM(type_trac)                        !--- Set %component: model component name
-        tracs(it)%keys%key = ['hadv', 'vadv']                        !--- Set %keys%key
-        tracs(it)%keys%val = s(1:2)                                  !--- Set %keys%val
+        k => tracers(it)%keys
+
+        !=== NAME OF THE TRACER
+        tname = old2newH2O(s(3), ip)
+        ix = strIdx(oldHNO3, s(3))
+        IF(ix /= 0 .AND. lRep) tname = newHNO3(ix)                   !--- Exception for HNO3 (REPROBUS ONLY)
+        tracers(it)%name = tname                                     !--- Set %name
+        CALL addKey_1('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'
+        tracers(it)%component = cname                                !--- Set %component
+        CALL addKey_1('component', cname, k)                         !--- Set the name of the model component
+
+        !=== NAME OF THE PARENT
+        pname = tran0                                                !--- Default name: default transporting fluid (air)
+        IF(ns == 4) THEN
+          pname = old2newH2O(s(4))
+          ix = strIdx(oldHNO3, s(4))
+          IF(ix /= 0 .AND. lRep) pname = newHNO3(ix)                 !--- Exception for HNO3 (REPROBUS ONLY)
+        END IF
+        tracers(it)%parent = pname                                   !--- Set %parent
+        CALL addKey_1('parent', pname, k)
+
+        !=== PHASE AND ADVECTION SCHEMES NUMBERS
+        tracers(it)%phase = known_phases(ip:ip)                      !--- Set %phase:  tracer phase (default: "g"azeous)
+        CALL addKey_1('phase', known_phases(ip:ip), k)               !--- Set the phase  of the tracer (default: "g"azeous)
+        CALL addKey_1('hadv', s(1),  k)                              !--- Set the horizontal advection schemes number
+        CALL addKey_1('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'
+      DO it=1,ntrac
+        CALL addKey_1('type', tracers(it)%type, tracers(it)%keys)    !--- Set the type of tracer
+      END DO
+      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
@@ -175,19 +268,57 @@
   END SELECT
   !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
   IF(ALL([2,3] /= fType)) RETURN
 
   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
-  WHERE(tracs%gen0Name(1:3) /= 'H2O') tracs%isInPhysics=.TRUE.       !--- Set %isInPhysics: passed to physics
-  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, lDisp, tracf, sects) RESULT(lerr)
+  CHARACTER(LEN=*),                             INTENT(IN)  :: modname, type_trac
+  INTEGER,                                      INTENT(OUT) :: fType
+  LOGICAL,                            OPTIONAL, INTENT(IN)  :: lDisp
+  CHARACTER(LEN=maxlen), ALLOCATABLE, OPTIONAL, INTENT(OUT) :: tracf(:), sects(:)
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: trac_files(:), sections(:)
+  LOGICAL, ALLOCATABLE :: ll(:)
+  LOGICAL :: lD, lFound
+  INTEGER :: is, nsec
+  lD = .FALSE.; IF(PRESENT(lDisp)) lD = lDisp
+  lerr = .FALSE.
+
+  !--- PARSE "type_trac" LIST AND DETERMINE THE TRACERS FILES NAMES (FOR CASE 3: MULTIPLE FILES, SINGLE SECTION PER FILE)
+  !--- If type_trac is a scalar (case 1), "sections" and "trac_files" are not usable, but are meaningless for case 1 anyway.
+  IF(test(strParse(type_trac, '|', sections,  n=nsec), lerr)) RETURN !--- Parse "type_trac" list
+  IF(PRESENT(sects)) sects = sections
+  ALLOCATE(trac_files(nsec), ll(nsec))
+  DO is=1, nsec
+     trac_files(is) = 'tracer_'//TRIM(sections(is))//'.def'
+     INQUIRE(FILE=TRIM(trac_files(is)), EXIST=ll(is))
+  END DO
+  IF(PRESENT(tracf)) tracf = trac_files
+  fType = 0
+  INQUIRE(FILE='traceur.def', EXIST=lFound); IF(lFound)  fType = 1   !--- OLD STYLE FILE
+  INQUIRE(FILE='tracer.def',  EXIST=lFound); IF(lFound)  fType = 2   !--- NEW STYLE ; SINGLE  FILE, SEVERAL SECTIONS
+                                             IF(ALL(ll)) fType = 3   !--- NEW STYLE ; SEVERAL FILES, SINGLE SECTION USED
+  IF(.NOT.lD) RETURN                                                 !--- NO CHECKING/DISPLAY NEEDED: JUST GET type_trac,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
+
+  !--- 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
 !==============================================================================================================================
 
@@ -206,5 +337,4 @@
   INTEGER,  ALLOCATABLE :: ndb(:)                                    !--- Number of sections for each file
   INTEGER,  ALLOCATABLE :: ixf(:)                                    !--- File index for each section of the expanded list
-  LOGICAL,  ALLOCATABLE :: lTg(:)                                    !--- Tagging tracers mask
   CHARACTER(LEN=maxlen) :: fnm, snm
   INTEGER               :: idb, i
@@ -224,5 +354,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
@@ -246,5 +376,5 @@
   CHARACTER(LEN=maxlen), ALLOCATABLE :: sec(:)
   INTEGER,               ALLOCATABLE ::  ix(:)
-  INTEGER :: n0, idb, ndb, i, j
+  INTEGER :: n0, idb, ndb
   LOGICAL :: ll
 !------------------------------------------------------------------------------------------------------------------------------
@@ -299,5 +429,5 @@
       CALL MOVE_ALLOC(FROM=tdb, TO=dBase)
     ELSE                                                             !=== TRACER LINE
-      ll = strParse(str,' ', keys = s, vals = v, n = n)              !--- Parse <key>=<val> pairs
+      ll = strParse(str,' ', s, n, v)                                !--- Parse <key>=<val> pairs
       tt = dBase(ndb)%trac(:)
       tmp%name = s(1); tmp%keys = keys_type(s(1), s(2:n), v(2:n))    !--- Set %name and %keys
@@ -330,6 +460,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"
@@ -381,5 +511,5 @@
   CHARACTER(LEN=maxlen), ALLOCATABLE :: ta(:), pa(:)
   CHARACTER(LEN=maxlen) :: msg1, modname
-  INTEGER :: it, nt, iq, nq, itr, ntr, ipr, npr, i
+  INTEGER :: it, nt, iq, nq, itr, ntr, ipr, npr
   LOGICAL :: ll
   modname = 'expandSection'
@@ -394,4 +524,6 @@
     tr(it)%type      = fgetKey(it, 'type'  , tr(:)%keys, 'tracer')
     tr(it)%component = sname
+!   CALL addKey_m('component', sname, tr(:)%keys)
+    DO iq=1,SIZE(tr); CALL addKey_1('component', sname, tr(iq)%keys); END DO
 
     !--- Determine the number of tracers and parents ; coherence checking
@@ -409,5 +541,4 @@
   END DO
   !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-  CALL delKey(['parent','type  '], tr)
 
   ALLOCATE(ttr(nq))
@@ -416,20 +547,18 @@
   DO it = 1, nt                                                      !=== EXPAND TRACERS AND PARENTS NAMES LISTS
   !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-    ll = strParse(tr(it)%name,   ',', ta, n=ntr)                     !--- Number of tracers
-    ll = strParse(tr(it)%parent, ',', pa, n=npr)                     !--- Number of parents
+    ll = strParse(tr(it)%name,   ',', ta, ntr)                       !--- Number of tracers
+    ll = strParse(tr(it)%parent, ',', pa, npr)                       !--- Number of parents
     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_1('name',      ta(itr),          ttr(iq)%keys)
+        ttr(iq)%parent    = TRIM(pa(ipr));    CALL addKey_1('parent',    pa(ipr),          ttr(iq)%keys)
+        ttr(iq)%type      = tr(it)%type;      CALL addKey_1('type',      tr(it)%type,      ttr(iq)%keys)
+        ttr(iq)%component = tr(it)%component; CALL addKey_1('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
@@ -441,33 +570,33 @@
 !==============================================================================================================================
 
-!==============================================================================================================================
-SUBROUTINE setGeneration(tr)
+
+!==============================================================================================================================
+LOGICAL FUNCTION setGeneration(tr) RESULT(lerr)
 !------------------------------------------------------------------------------------------------------------------------------
 ! Purpose: Determine, for each tracer of "tr(:)":
 !   * %iGeneration: the generation number 
 !   * %gen0Name:    the generation 0 ancestor name
+!          Check also for orphan tracers (tracers not descending on "tran0").
 !------------------------------------------------------------------------------------------------------------------------------
   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
-  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
-
-  !=== 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( 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
-!==============================================================================================================================
+  INTEGER                            :: iq, jq, ig
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: parent(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen) :: modname
+  modname = 'setGeneration'
+  IF(test(fmsg('missing "parent" attribute', modname, getKey('parent', parent, ky=tr(:)%keys)), lerr)) RETURN
+  DO iq = 1, SIZE(tr)
+    jq = iq; ig = 0
+    DO WHILE(parent(jq) /= tran0)
+      jq = strIdx(tr(:)%name, parent(jq))
+      IF(test(fmsg('Orphan tracer "'//TRIM(tr(iq)%name)//'"', modname, jq == 0), lerr)) RETURN
+      ig = ig + 1
+    END DO
+    tr(iq)%gen0Name = tr(jq)%name; CALL addKey_1('gen0Name',    tr(iq)%gen0Name,   tr(iq)%keys)
+    tr(iq)%iGeneration = ig;       CALL addKey_1('iGeneration', TRIM(int2str(ig)), tr(iq)%keys)
+  END DO
+END FUNCTION setGeneration
+!==============================================================================================================================
+
 
 !==============================================================================================================================
@@ -503,4 +632,5 @@
 END FUNCTION checkTracers
 !==============================================================================================================================
+
 
 !==============================================================================================================================
@@ -543,4 +673,5 @@
 !==============================================================================================================================
 
+
 !==============================================================================================================================
 SUBROUTINE expandPhases(tr)
@@ -552,8 +683,8 @@
   TYPE(trac_type), ALLOCATABLE :: ttr(:)
   INTEGER,   ALLOCATABLE ::  i0(:)
-  CHARACTER(LEN=maxlen)  :: nam, pha, trn
+  CHARACTER(LEN=maxlen)  :: nam, pha, tname
   CHARACTER(LEN=1) :: p
   INTEGER :: ip, np, iq, jq, nq, it, nt, nc, i, n
-  LOGICAL :: lTg, lEx
+  LOGICAL :: lTag, lExt
 !------------------------------------------------------------------------------------------------------------------------------
   nq = SIZE(tr, DIM=1)
@@ -561,7 +692,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
@@ -569,9 +700,9 @@
   it = 1                                                             !--- Current "ttr(:)" index
   DO iq = 1, nq                                                      !--- Loop on "tr(:)" indexes
-    lTg = tr(iq)%type=='tag'                                         !--- Current tracer is a tag
+    lTag = tr(iq)%type=='tag'                                        !--- Current tracer is a tag
     i0 = strFind(tr(:)%name, TRIM(tr(iq)%gen0Name), n)               !--- Indexes of first generation ancestor copies
     np = SUM([( LEN_TRIM(tr(i0(i))%phase),i=1,n )], 1)               !--- Number of phases for current tracer tr(iq)
-    lEx = np>1                                                       !--- Phase suffix only required if phases number is > 1
-    IF(lTg) lEx = lEx .AND. tr(iq)%iGeneration>0                     !--- No phase suffix for generation 0 tags
+    lExt = np>1                                                      !--- Phase suffix only required if phases number is > 1
+    IF(lTag) lExt = lExt .AND. tr(iq)%iGeneration>0                  !--- No phase suffix for generation 0 tags
     DO i=1,n                                                         !=== LOOP ON GENERATION 0 ANCESTORS
       jq = i0(i)                                                     !--- tr(jq): ith tracer with same gen 0 ancestor as tr(iq)
@@ -580,15 +711,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)
-        IF(lTg) nam = TRIM(tr(iq)%parent)                            !--- Parent name (tagging case)
-        IF(lEx) nam = addPhase(nam, p )                              !--- Phase extension needed
-        IF(lTg) nam = TRIM(nam)//'_'//TRIM(trn)                      !--- <parent>_<name> for tags
+        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(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
-        IF(lEx .AND. tr(iq)%iGeneration>0) THEN
-          ttr(it)%parent   = addPhase(ttr(it)%parent,   p)
-          ttr(it)%gen0Name = addPhase(ttr(it)%gen0Name, p)
+        CALL addKey_1('name', nam, ttr(it)%keys)
+        CALL addKey_1('phase', p,  ttr(it)%keys)
+        IF(lExt .AND. tr(iq)%iGeneration>0) THEN
+          ttr(it)%parent   = addPhase(tr(iq)%parent,   p)
+          ttr(it)%gen0Name = addPhase(tr(iq)%gen0Name, p)
+          CALL addKey_1('parent',   ttr(it)%parent,   ttr(it)%keys)
+          CALL addKey_1('gen0Name', ttr(it)%gen0Name, ttr(it)%keys)
         END IF
         it = it+1
@@ -603,4 +738,5 @@
 !==============================================================================================================================
 
+
 !==============================================================================================================================
 SUBROUTINE sortTracers(tr)
@@ -609,10 +745,13 @@
 !  * 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
 !------------------------------------------------------------------------------------------------------------------------------
   TYPE(trac_type), ALLOCATABLE, INTENT(INOUT) :: tr(:)         !--- Tracer derived type vector
-  INTEGER,         ALLOCATABLE :: iy(:), iz(:)
-  INTEGER :: ig, ng, iq, jq, ip, nq, n, ix(SIZE(tr)), k
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE(trac_type), ALLOCATABLE        :: tr2(:)
+  INTEGER,         ALLOCATABLE        :: iy(:), iz(:)
+  INTEGER                             :: ig, ng, iq, jq, ip, nq, n, ix(SIZE(tr)), k
+!  tr2 is introduced in order to cope with a bug in gfortran 4.8.5 compiler
 !------------------------------------------------------------------------------------------------------------------------------
   nq = SIZE(tr)
@@ -620,5 +759,6 @@
     iq = strIdx(tr(:)%name, addPhase('H2O', ip))
     IF(iq == 0) CYCLE
-    tr = [tr(iq), tr(1:iq-1), tr(iq+1:nq)]
+    tr2 = tr(:)
+    tr = [tr2(iq), tr2(1:iq-1), tr2(iq+1:nq)]
   END DO
   IF(lSortByGen) THEN
@@ -637,5 +777,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
@@ -649,4 +789,5 @@
 END SUBROUTINE sortTracers
 !==============================================================================================================================
+
 
 !==============================================================================================================================
@@ -724,5 +865,5 @@
   TYPE(dataBase_type),  TARGET, INTENT(IN)  :: sections(:)
   TYPE(trac_type), ALLOCATABLE, INTENT(OUT) ::       tr(:)
-  TYPE(trac_type), POINTER     :: t1(:), t2(:)
+  TYPE(trac_type), POINTER     :: t(:)
   INTEGER,   ALLOCATABLE :: nt(:)
   CHARACTER(LEN=maxlen)  :: tnam, tnam_new
@@ -735,10 +876,10 @@
   DO is=1, nsec                                                      !=== LOOP ON SECTIONS
   !----------------------------------------------------------------------------------------------------------------------------
-    t1 => sections(is)%trac(:)
+    t => sections(is)%trac(:)
     !--------------------------------------------------------------------------------------------------------------------------
     DO iq=1, nt(is)                                                  !=== LOOP ON TRACERS
     !--------------------------------------------------------------------------------------------------------------------------
-      tnam = TRIM(t1(iq)%name)                                       !--- Original name
-      IF(COUNT(t1%name == tnam) == 1) CYCLE                          !--- Current tracer is not duplicated: finished
+      tnam = TRIM(t(iq)%name)                                        !--- Original name
+      IF(COUNT(t%name == tnam) == 1) CYCLE                           !--- Current tracer is not duplicated: finished
       tnam_new = TRIM(tnam)//'_'//TRIM(sections(is)%name)            !--- Same with section extension
       nq = SUM(nt(1:is-1))                                           !--- Number of tracers in previous sections
@@ -759,5 +900,5 @@
   TYPE(trac_type), INTENT(INOUT) :: tr(:)
 
-  !--- Update %iqParent, %iqDescen, %nqDescen, %nqChilds
+  !--- Update %iqParent, %iqDescen, %nqDescen, %nqChildren
   CALL indexUpdate(tr)
 
@@ -774,5 +915,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.
@@ -782,17 +923,20 @@
   !--- 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(fgetKeys('hadv',  tm(:)%keys, '10'))
+  vadv = str2int(fgetKeys('vadv',  tm(:)%keys, '10'))
+  prnt =         fgetKeys('parent',tm(:)%keys,  '' )
+  IF(getKey('phases', phas, ky=tm(:)%keys)) phas = fGetKeys('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
-!==============================================================================================================================
 !==============================================================================================================================
 
@@ -809,5 +953,5 @@
   out => NULL(); IF(it /= 0) out => t(it)
 END FUNCTION aliasTracer
-!------------------------------------------------------------------------------------------------------------------------------
+!==============================================================================================================================
 
 
@@ -830,5 +974,5 @@
   CALL indexUpdate(out)
 END FUNCTION trSubset_Name
-!------------------------------------------------------------------------------------------------------------------------------
+!==============================================================================================================================
 
 
@@ -843,5 +987,5 @@
   CALL indexUpdate(out)
 END FUNCTION trSubset_gen0Name
-!------------------------------------------------------------------------------------------------------------------------------
+!==============================================================================================================================
 
 
@@ -851,7 +995,7 @@
 SUBROUTINE indexUpdate(tr)
   TYPE(trac_type), INTENT(INOUT) :: tr(:)
-  INTEGER :: iq, ig, ng, igen, ngen
-  INTEGER, ALLOCATABLE :: ix(:)
+  INTEGER :: iq, ig, igen, ngen, ix(SIZE(tr))
   tr(:)%iqParent = strIdx( tr(:)%name, tr(:)%parent )                !--- Parent index
+  DO iq = 1, SIZE(tr); CALL addKey_1('iqParent', int2str(tr(iq)%iqParent), tr(iq)%keys); END DO
   ngen = MAXVAL(tr(:)%iGeneration, MASK=.TRUE.)
   DO iq = 1, SIZE(tr)
@@ -859,13 +1003,18 @@
     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_1('nqChildren', int2str(tr(iq)%nqChildren), tr(iq)%keys)
+      END IF
     END DO
+    CALL addKey_1('iqDescen', strStack(int2str(tr(iq)%iqDescen)), tr(iq)%keys)
+    CALL addKey_1('nqDescen',          int2str(tr(iq)%nqDescen) , tr(iq)%keys)
   END DO
 END SUBROUTINE indexUpdate
-!------------------------------------------------------------------------------------------------------------------------------
+!==============================================================================================================================
   
   
@@ -875,5 +1024,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       ====
@@ -883,20 +1032,18 @@
 !===  * 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 "prnt" must be defined !)
-  INTEGER :: ik, is, it, idb, nk0, i, iis
-  INTEGER :: nk, ns, nt, ndb, nb0, i0
-  CHARACTER(LEN=maxlen), POINTER     :: k(:), v(:), k0(:), v0(:)
+  TYPE(isot_type), TARGET, INTENT(INOUT) :: isot(:)                  !--- Isotopes descriptors (field %parent must be defined!)
+  LOGICAL :: lFound
+  INTEGER :: is, iis, it, idb, ndb, nb0
   CHARACTER(LEN=maxlen), ALLOCATABLE :: vals(:)
-  CHARACTER(LEN=maxlen)              :: val, modname
-  TYPE(keys_type),           POINTER ::   ky(:)
+  CHARACTER(LEN=maxlen)              :: modname
   TYPE(trac_type),           POINTER ::   tt(:), t
   TYPE(dataBase_type),   ALLOCATABLE ::  tdb(:)
-  LOGICAL,               ALLOCATABLE :: liso(:)
   modname = 'readIsotopesFile'
 
   !--- THE INPUT FILE MUST BE PRESENT
-  IF(test(fmsg('Missing isotopes parameters file "'//TRIM(fnam)//'"', modname, testFile(fnam)),lerr)) RETURN
+  INQUIRE(FILE=TRIM(fnam), EXIST=lFound); lerr = .NOT.lFound
+  IF(fmsg('Missing isotopes parameters file "'//TRIM(fnam)//'"', modname, lerr)) RETURN
 
   !--- READ THE FILE SECTIONS, ONE EACH PARENT TRACER
@@ -905,5 +1052,5 @@
   ndb = SIZE(dBase, DIM=1)                                           !--- Current database size
   DO idb = nb0, ndb
-   iis = idb-nb0+1
+    iis = idb-nb0+1
 
     !--- GET FEW GLOBAL KEYS FROM "def" FILES AND ADD THEM TO THE 'params' SECTION
@@ -920,14 +1067,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
 
@@ -942,8 +1087,36 @@
   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
+!==============================================================================================================================
+
 
 !==============================================================================================================================
@@ -951,34 +1124,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
-
-  p = PACK(delPhase(t%parent), MASK = t%type=='tracer' .AND. t%iGeneration==1) !--- Parents of generation 1 isotopes
+  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"
@@ -989,9 +1172,9 @@
     !=== Geographic tagging tracers descending on tracer "iname": mask, names, number
     ll = t(:)%type=='tag'    .AND. delPhase(t(:)%gen0Name) == iname .AND. t(:)%iGeneration == 2
-    i%zone = PACK(strTail(t(:)%name,'_'), MASK = ll)                 !--- Tagging zones names  for isotopes category "iname"
+    i%zone = PACK(strTail(t(:)%name,'_',.TRUE.), MASK = ll)          !--- Tagging zones names  for isotopes category "iname"
     CALL strReduce(i%zone)
     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)
@@ -1009,9 +1192,9 @@
     !=== 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"
-      t1%iso_iName  = strIdx(i%trac, delPhase(strHead(t1%name,'_'))) !--- Current isotope       idx in effective isotopes list
-      t1%iso_iZone  = strIdx(i%zone,          strTail(t1%name,'_') ) !--- Current isotope zone  idx in effective zones    list
+      t1%iso_iName  = strIdx(i%trac, strHead(delPhase(t1%name),'_',.TRUE.)) !--- Current isotope       idx in effective isotopes list
+      t1%iso_iZone  = strIdx(i%zone,          strTail(t1%name, '_',.TRUE.)) !--- Current isotope zone  idx in effective zones    list
       t1%iso_iPhase =  INDEX(i%phase,TRIM(t1%phase))                 !--- Current isotope phase idx in effective phases   list
       IF(t1%iGeneration /= 2) t1%iso_iZone = 0                       !--- Skip possible generation 1 tagging tracers
@@ -1020,5 +1203,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
@@ -1027,59 +1210,108 @@
   END DO
 
-  !=== READ PHYSICAL PARAMETERS FROM "isotopes_params.def" FILE
-  !    DONE HERE, AND NOT ONLY IN "infotrac_phy", BECAUSE SOME PHYSICAL PARAMS ARE NEEDED FOR RESTARTS (tnat AND alpha_ideal)
-  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) = ['iq  ','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
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== ADD THE <key>=<val> PAIR TO THE "ky[(:)]" KEY[S] DESCRIPTOR[S] OR THE <key>=<val(:)> PAIRS TO THE "ky(:)" KEYS DESCRIPTORS
 !==============================================================================================================================
 SUBROUTINE addKey_1(key, val, ky, lOverWrite)
-!------------------------------------------------------------------------------------------------------------------------------
-! Purpose: Add the <key>=<val> pair in the "ky" keys descriptor.
-!------------------------------------------------------------------------------------------------------------------------------
   CHARACTER(LEN=*),  INTENT(IN)    :: key, val
   TYPE(keys_type),   INTENT(INOUT) :: ky
   LOGICAL, OPTIONAL, INTENT(IN)    :: lOverWrite
+!------------------------------------------------------------------------------------------------------------------------------
   CHARACTER(LEN=maxlen), ALLOCATABLE :: k(:), v(:)
   INTEGER :: iky, nky
   LOGICAL :: lo
-!------------------------------------------------------------------------------------------------------------------------------
-  lo=.FALSE.; IF(PRESENT(lOverWrite)) lo=lOverWrite
+  lo=.TRUE.; IF(PRESENT(lOverWrite)) lo=lOverWrite
+  IF(.NOT.ALLOCATED(ky%key)) THEN
+    ALLOCATE(ky%key(1)); ky%key(1)=key
+    ALLOCATE(ky%val(1)); ky%val(1)=val
+    RETURN
+  END IF
   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
+    ALLOCATE(k(nky+1)); k(1:nky) = ky%key; k(nky+1) = key; ky%key = k
+    ALLOCATE(v(nky+1)); v(1:nky) = ky%val; v(nky+1) = val; ky%val = v
+  ELSE IF(lo) THEN
     ky%key(iky) = key; ky%val(iky) = val
   END IF
@@ -1087,23 +1319,32 @@
 !==============================================================================================================================
 SUBROUTINE addKey_m(key, val, ky, lOverWrite)
-!------------------------------------------------------------------------------------------------------------------------------
-! Purpose: Add the <key>=<val> pair in all the components of the "ky" keys descriptor.
-!------------------------------------------------------------------------------------------------------------------------------
   CHARACTER(LEN=*),  INTENT(IN)    :: key, val
   TYPE(keys_type),   INTENT(INOUT) :: ky(:)
   LOGICAL, OPTIONAL, INTENT(IN)    :: lOverWrite
+!------------------------------------------------------------------------------------------------------------------------------
   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
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== OVERWRITE THE KEYS OF THE TRACER NAMED "tr0" WITH THE VALUES FOUND IN THE *.def FILES, IF ANY. ===========================
+!==============================================================================================================================
 SUBROUTINE addKeysFromDef(t, tr0)
-!------------------------------------------------------------------------------------------------------------------------------
-! Purpose: The values of the keys of the tracer named "tr0" are overwritten by the values found in the *.def files, if any.
-!------------------------------------------------------------------------------------------------------------------------------
   TYPE(trac_type), ALLOCATABLE, INTENT(INOUT) :: t(:)
   CHARACTER(LEN=*),             INTENT(IN)    :: tr0
+!------------------------------------------------------------------------------------------------------------------------------
   CHARACTER(LEN=maxlen) :: val
   INTEGER               :: ik, jd
@@ -1116,16 +1357,17 @@
 END SUBROUTINE addKeysFromDef
 !==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== REMOVE THE KEYS NAMED "keyn(:)" FROM EITHER THE "itr"th OR ALL THE KEYS DESCRIPTORS OF "ky(:)" ===========================
+!==============================================================================================================================
 SUBROUTINE delKey_1(itr, keyn, ky)
-!------------------------------------------------------------------------------------------------------------------------------
-! Purpose: Internal routine.
-!   Remove <key>=<val> pairs in the "itr"th component of the "ky" keys descriptor.
-!------------------------------------------------------------------------------------------------------------------------------
   INTEGER,          INTENT(IN)    :: itr
   CHARACTER(LEN=*), INTENT(IN)    :: keyn(:)
   TYPE(trac_type),  INTENT(INOUT) :: ky(:)
+!------------------------------------------------------------------------------------------------------------------------------
   CHARACTER(LEN=maxlen), ALLOCATABLE :: k(:), v(:)
   LOGICAL,               ALLOCATABLE :: ll(:)
   INTEGER :: iky
-!------------------------------------------------------------------------------------------------------------------------------
   IF(itr<=0 .OR. itr>SIZE(ky, DIM=1)) RETURN                          !--- Index is out of range
   ll = [( ALL(keyn/=ky(itr)%keys%key(iky)), iky=1, SIZE(ky(itr)%keys%key) )]
@@ -1135,12 +1377,8 @@
 !==============================================================================================================================
 SUBROUTINE delKey(keyn, ky)
-!------------------------------------------------------------------------------------------------------------------------------
-! Purpose: Internal routine.
-!   Remove <key>=<val> pairs in all the components of the "t" tracers descriptor.
-!------------------------------------------------------------------------------------------------------------------------------
   CHARACTER(LEN=*), INTENT(IN)    :: keyn(:)
   TYPE(trac_type),  INTENT(INOUT) :: ky(:)
+!------------------------------------------------------------------------------------------------------------------------------
   INTEGER :: iky
-!------------------------------------------------------------------------------------------------------------------------------
   DO iky = 1, SIZE(ky); CALL delKey_1(iky, keyn, ky); END DO
 END SUBROUTINE delKey
@@ -1149,34 +1387,25 @@
 
 !==============================================================================================================================
-!=== PUBLIC ROUTINES: GET A KEY FROM A <key>=<val> LIST ; VECTORS, TRACER AND DATABASE VERSIONS ===============================
-!=== BEWARE !!! IF THE "ky" ARGUMENT IS NOT PRESENT, THEN THE VARIABLES "tracers" AND "isotopes" ARE USED. ====================
-!===     THEY ARE LOCAL TO THIS MODULE, SO MUST MUST BE INITIALIZED FIRST USING the "getKey_init" ROUTINE  ====================
-!==============================================================================================================================
-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
-!==============================================================================================================================
-CHARACTER(LEN=maxlen) FUNCTION fgetKeyByIndex_s1(itr, keyn, ky, def_val) RESULT(val)
-!------------------------------------------------------------------------------------------------------------------------------
-! Purpose: Internal function ; get a string formatted key value (returned argument) from its name and the tracer index.
-!------------------------------------------------------------------------------------------------------------------------------
-  INTEGER,                    INTENT(IN) :: itr
-  CHARACTER(LEN=*),           INTENT(IN) :: keyn
-  TYPE(keys_type),            INTENT(IN) :: ky(:)
-  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: def_val
+!================ GET THE VALUE OF A KEY FROM A "keys_type" DERIVED TYPE ; THE RESULT IS THE RETURNED VALUE ===================
+!==============================================================================================================================
+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)
-!------------------------------------------------------------------------------------------------------------------------------
-! Purpose: Internal function ; get a string formatted key value (returned argument) from its name and the tracer name.
-!------------------------------------------------------------------------------------------------------------------------------
+  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(:)
@@ -1184,61 +1413,116 @@
   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 fgetKeys(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 fgetKeys
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!========== GET THE VALUE OF A KEY FROM A "keys_type" DERIVED TYPE ; THE RETURNED VALUE IS THE ERROR CODE        ==============
+!==========  The key "keyn" is searched in: 1)           "ky(:)%name" (if given)                                 ==============
+!==========                                 2)      "tracers(:)%name"                                            ==============
+!==========                                 3) "isotope%keys(:)%name"                                            ==============
+!==========  for the tracer[s] "tname[(:)]" (if given) or all the available tracers from the used set otherwise. ==============
+!==========  The type of the returned value(s) can be string, integer or real, scalar or vector                  ==============
 !==============================================================================================================================
 LOGICAL FUNCTION getKeyByName_s1(keyn, val, tname, ky) RESULT(lerr)
-  !--- Purpose: Get the value of the key named "keyn" for the tracer named "tnam".
-  !     * "ky" unspecified: try in "tracers" for "tnam" with phase and tagging suffixes, then in "isotopes" without.
-  !     * "ky"   specified: try in "ky"      for "tnam" with phase and tagging suffixes, then without.
-  !    The returned error code is always .FALSE.: an empty string is returned when the key hasn't been found.
   CHARACTER(LEN=*),          INTENT(IN)  :: keyn
   CHARACTER(LEN=maxlen),     INTENT(OUT) :: val
   CHARACTER(LEN=*),          INTENT(IN)  :: tname
   TYPE(keys_type), OPTIONAL, INTENT(IN)  :: ky(:)
+!------------------------------------------------------------------------------------------------------------------------------
   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),'_',.TRUE.)                                             !--- Remove phase and tag
+  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_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(:)
+  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_sm(keyn, val, tname, ky, nam) RESULT(lerr)
+  CHARACTER(LEN=*),                             INTENT(IN)  :: keyn
+  CHARACTER(LEN=maxlen),           ALLOCATABLE, INTENT(OUT) :: val(:)
+  CHARACTER(LEN=*),                             INTENT(IN)  :: tname(:)
+  TYPE(keys_type),       OPTIONAL, TARGET,      INTENT(IN)  :: ky(:)
+  CHARACTER(LEN=maxlen), OPTIONAL, ALLOCATABLE, INTENT(OUT) :: nam(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE(keys_type), POINTER ::  keys(:)
+  LOGICAL :: lk, lt, li
   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)
+
+  !--- 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)
+
+  !--- LINK "keys" TO THE RIGHT DATABASE
+  IF(test(.NOT.ANY([lk,lt,li]), lerr)) RETURN
+  IF(lk) keys => ky(:)
+  IF(lt) keys => tracers(:)%keys
+  IF(li) keys => isotope%keys(:)
+
+  !--- GET THE DATA
+  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)])
+  lerr = ANY([(getKeyByName_s1(keyn, val(iq), tname(iq), keys(:)), iq=1, nq)])
+  IF(PRESENT(nam)) nam = tname(:)
+
+END FUNCTION getKeyByName_sm
+!==============================================================================================================================
+LOGICAL FUNCTION getKey_sm(keyn, val, ky, nam) RESULT(lerr)
+  CHARACTER(LEN=*),                             INTENT(IN)  :: keyn
+  CHARACTER(LEN=maxlen),           ALLOCATABLE, INTENT(OUT) :: val(:)
+  TYPE(keys_type),       OPTIONAL,              INTENT(IN)  :: ky(:)
+  CHARACTER(LEN=maxlen), OPTIONAL, ALLOCATABLE, INTENT(OUT) :: nam(:)
+!------------------------------------------------------------------------------------------------------------------------------
+! Note: if names are repeated, getKeyByName_sm can't be used ; this routine, using indexes, must be used instead.
+  IF(PRESENT(ky)) THEN                                                                   !=== KEY FROM "ky"
+    val = fgetKeys(keyn, ky, lerr=lerr)
+    IF(PRESENT(nam)) nam = ky(:)%name
+  ELSE
+    IF(         .NOT.test(.NOT.ALLOCATED(tracers ), lerr)) lerr = SIZE(tracers ) == 0    !=== KEY FROM "tracers"
+    IF(.NOT.lerr) val = fgetKeys(keyn, tracers%keys, lerr=lerr)
+    IF(.NOT.lerr.AND.PRESENT(nam)) nam = tracers(:)%keys%name
+    IF(lerr.AND..NOT.test(.NOT.ALLOCATED(isotopes), lerr)) lerr = SIZE(isotopes) == 0    !=== KEY FROM "isotope"
+    IF(.NOT.lerr) val = fgetKeys(keyn, isotope%keys, lerr=lerr)
+    IF(.NOT.lerr.AND.PRESENT(nam)) nam = isotope%keys(:)%name
   END IF
-END FUNCTION getKeyByName_sm
+END FUNCTION getKey_sm
 !==============================================================================================================================
 LOGICAL FUNCTION getKeyByName_i1(keyn, val, tname, ky) RESULT(lerr)
@@ -1247,31 +1531,61 @@
   CHARACTER(LEN=*),          INTENT(IN)  :: tname
   TYPE(keys_type), OPTIONAL, INTENT(IN)  :: ky(:)
+!------------------------------------------------------------------------------------------------------------------------------
   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_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 :: sval(:)
+  INTEGER :: ierr, iq, nq
+  IF(test(getKeyByName_s1m(keyn, sval, tname, ky), lerr)) RETURN
+  nq = SIZE(sval); ALLOCATE(val(nq))
+  lerr = .FALSE.; DO iq=1, nq; READ(sval(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_i1m
+!==============================================================================================================================
+LOGICAL FUNCTION getKeyByName_im(keyn, val, tname, ky, nam) 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), OPTIONAL, ALLOCATABLE, INTENT(OUT) ::   nam(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: sval(:), names(:)
+  INTEGER :: ierr, iq, nq
+  IF(test(getKeyByName_sm(keyn, sval, tname, ky, names), lerr)) RETURN
+  nq = SIZE(sval); ALLOCATE(val(nq))
+  DO iq = 1, nq                                                      !--- CONVERT THE KEYS TO INTEGERS
+    READ(sval(iq), *, IOSTAT=ierr) val(iq)
+    IF(test(fmsg('key "'//TRIM(keyn)//'" of tracer "'//TRIM(names(iq))//'" is not an integer', modname, ierr/=0), lerr)) RETURN
+  END DO
+  IF(PRESENT(nam)) nam = names(:)
 END FUNCTION getKeyByName_im
+!==============================================================================================================================
+LOGICAL FUNCTION getKey_im(keyn, val, ky, nam) RESULT(lerr)
+  CHARACTER(LEN=*),                             INTENT(IN)  :: keyn
+  INTEGER,                         ALLOCATABLE, INTENT(OUT) ::  val(:)
+  TYPE(keys_type),       OPTIONAL,              INTENT(IN)  ::   ky(:)
+  CHARACTER(LEN=maxlen), OPTIONAL, ALLOCATABLE, INTENT(OUT) ::  nam(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: sval(:), names(:)
+  INTEGER :: ierr, iq, nq
+  IF(test(getKey_sm(keyn, sval, ky, names), lerr)) RETURN
+  nq = SIZE(sval); ALLOCATE(val(nq))
+  DO iq = 1, nq
+    READ(sval(iq), *, IOSTAT=ierr) val(iq)
+    IF(test(fmsg('key "'//TRIM(keyn)//'" of tracer "'//TRIM(names(iq))//'" is not an integer', modname, ierr/=0), lerr)) RETURN
+  END DO
+  IF(PRESENT(nam)) nam = names
+END FUNCTION getKey_im
 !==============================================================================================================================
 LOGICAL FUNCTION getKeyByName_r1(keyn, val, tname, ky) RESULT(lerr)
@@ -1280,31 +1594,153 @@
   CHARACTER(LEN=*),          INTENT(IN)  :: tname
   TYPE(keys_type), OPTIONAL, INTENT(IN)  :: ky(:)
+!------------------------------------------------------------------------------------------------------------------------------
   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
 !==============================================================================================================================
-LOGICAL FUNCTION getKeyByName_rm(keyn, val, tname, ky) RESULT(lerr)
+LOGICAL FUNCTION getKeyByName_r1m(keyn, val, tname, ky) RESULT(lerr)
   CHARACTER(LEN=*),           INTENT(IN)  :: keyn
-  REAL,          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(:)
+  REAL,          ALLOCATABLE, INTENT(OUT) :: val(:)
+  CHARACTER(LEN=*),           INTENT(IN)  :: tname
+  TYPE(keys_type),  OPTIONAL, INTENT(IN)  ::  ky(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: sval(:)
+  INTEGER :: ierr, iq, nq
+  IF(test(getKeyByName_s1m(keyn, sval, tname, ky), lerr)) RETURN
+  nq = SIZE(sval); ALLOCATE(val(nq))
+  lerr = .FALSE.; DO iq=1, nq; READ(sval(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 vector of reals', modname, lerr), lerr)) RETURN
+END FUNCTION getKeyByName_r1m
+!==============================================================================================================================
+LOGICAL FUNCTION getKeyByName_rm(keyn, val, tname, ky, nam) 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), OPTIONAL, ALLOCATABLE, INTENT(OUT) ::   nam(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: sval(:), names(:)
+  INTEGER :: ierr, iq, nq
+  IF(test(getKeyByName_sm(keyn, sval, tname, ky, names), lerr)) RETURN
+  nq = SIZE(sval); ALLOCATE(val(nq))
+  DO iq = 1, nq                                                      !--- CONVERT THE KEYS TO INTEGERS
+    READ(sval(iq), *, IOSTAT=ierr) val(iq)
+    IF(test(fmsg('key "'//TRIM(keyn)//'" of tracer "'//TRIM(names(iq))//'" is not a real', modname, ierr/=0), lerr)) RETURN
+  END DO
+  IF(PRESENT(nam)) nam = names
+END FUNCTION getKeyByName_rm
+!==============================================================================================================================
+LOGICAL FUNCTION getKey_rm(keyn, val, ky, nam) RESULT(lerr)
+  CHARACTER(LEN=*),                             INTENT(IN)  :: keyn
+  REAL,                            ALLOCATABLE, INTENT(OUT) ::  val(:)
+  TYPE(keys_type),       OPTIONAL,              INTENT(IN)  ::   ky(:)
+  CHARACTER(LEN=maxlen), OPTIONAL, ALLOCATABLE, INTENT(OUT) ::  nam(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: sval(:), names(:)
+  INTEGER :: ierr, iq, nq
+  IF(test(getKey_sm(keyn, sval, ky, names), lerr)) RETURN
+  nq = SIZE(sval); ALLOCATE(val(nq))
+  DO iq = 1, nq                                                      !--- CONVERT THE KEYS TO INTEGERS
+    READ(sval(iq), *, IOSTAT=ierr) val(iq)
+    IF(test(fmsg('key "'//TRIM(keyn)//'" of tracer "'//TRIM(names(iq))//'" is not a real', modname, ierr/=0), lerr)) RETURN
+  END DO
+  IF(PRESENT(nam)) nam = names
+END FUNCTION getKey_rm
+!==============================================================================================================================
+LOGICAL FUNCTION getKeyByName_l1(keyn, val, tname, ky) RESULT(lerr)
+  USE strings_mod, ONLY: str2bool
+  CHARACTER(LEN=*),          INTENT(IN)  :: keyn
+  LOGICAL,                   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('key "'//TRIM(keyn)//'" or tracer "'//TRIM(tname)//'" is missing', modname, lerr), lerr)) RETURN
+  val = str2bool(sval)
+END FUNCTION getKeyByName_l1
+!==============================================================================================================================
+LOGICAL FUNCTION getKeyByName_l1m(keyn, val, tname, ky) RESULT(lerr)
+  USE strings_mod, ONLY: str2bool
+  CHARACTER(LEN=*),           INTENT(IN)  :: keyn
+  LOGICAL,       ALLOCATABLE, INTENT(OUT) :: val(:)
+  CHARACTER(LEN=*),           INTENT(IN)  :: tname
+  TYPE(keys_type),  OPTIONAL, INTENT(IN)  ::  ky(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: sval(:)
   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)])
+  IF(test(getKeyByName_s1m(keyn, sval, tname, ky), lerr)) RETURN
+  nq = SIZE(sval); ALLOCATE(val(nq))
+  lerr = .FALSE.; DO iq=1, nq; val(iq)=str2bool(sval(iq)); END DO
+END FUNCTION getKeyByName_l1m
+!==============================================================================================================================
+LOGICAL FUNCTION getKeyByName_lm(keyn, val, tname, ky, nam) RESULT(lerr)
+  USE strings_mod, ONLY: str2bool
+  CHARACTER(LEN=*),                             INTENT(IN)  ::  keyn
+  LOGICAL,                         ALLOCATABLE, INTENT(OUT) ::   val(:)
+  CHARACTER(LEN=*),                             INTENT(IN)  :: tname(:)
+  TYPE(keys_type),       OPTIONAL,              INTENT(IN)  ::    ky(:)
+  CHARACTER(LEN=maxlen), OPTIONAL, ALLOCATABLE, INTENT(OUT) ::   nam(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: sval(:)
+  INTEGER :: iq, nq
+  IF(test(getKeyByName_sm(keyn, sval, tname, ky, nam), lerr)) RETURN
+  nq = SIZE(sval); ALLOCATE(val(nq))
+  lerr = .FALSE.; DO iq=1, nq; val(iq)=str2bool(sval(iq)); END DO
+END FUNCTION getKeyByName_lm
+!==============================================================================================================================
+LOGICAL FUNCTION getKey_lm(keyn, val, ky, nam) RESULT(lerr)
+  USE strings_mod, ONLY: str2bool
+  CHARACTER(LEN=*),                             INTENT(IN)  :: keyn
+  LOGICAL,                         ALLOCATABLE, INTENT(OUT) ::  val(:)
+  TYPE(keys_type),       OPTIONAL,              INTENT(IN)  ::   ky(:)
+  CHARACTER(LEN=maxlen), OPTIONAL, ALLOCATABLE, INTENT(OUT) ::  nam(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: sval(:)
+  INTEGER :: iq, nq
+  IF(test(getKey_sm(keyn, sval, ky, nam), lerr)) RETURN
+  nq = SIZE(sval); ALLOCATE(val(nq))
+  lerr = .FALSE.; DO iq=1, nq; val(iq)=str2bool(sval(iq)); END DO
+END FUNCTION getKey_lm
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== 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
 !==============================================================================================================================
 
@@ -1315,28 +1751,21 @@
 ELEMENTAL CHARACTER(LEN=maxlen) FUNCTION delPhase(s) RESULT(out)
   CHARACTER(LEN=*), INTENT(IN) :: s
-  INTEGER :: l, i, ix
-  CHARACTER(LEN=maxlen) :: sh, st
-  out = s
-  IF(s == '') RETURN                                                           !--- Empty string: nothing to do
-
-  !--- Special case: old phases for water, no phases separator
-  i = INDEX(s,'_'); sh = s; IF(i/=0) sh=s(1:i-1); st='H2O'; IF(i/=0) st='H2O_'//s(i+1:LEN_TRIM(s))
-  IF(ANY([('H2O'//old_phases(ix:ix), ix=1, nphases)] == sh)) THEN; out=st; RETURN; END IF
-
-  !--- Index of found phase in "known_phases"
-  ix = MAXLOC( [(i, i=1,nphases)], MASK=[( INDEX(s, phases_sep//known_phases(i:i))/=0, i=1, nphases)], DIM=1 )
-  IF(ix == 0) RETURN                                                           !--- No phase pattern found
-  i = INDEX(s, phases_sep//known_phases(ix:ix))                                !--- Index of <sep><pha> pattern in "str"
-  l = LEN_TRIM(s)
-  IF(i == l-1) THEN                                                            !--- <var><sep><pha>       => return <var>
-    out = s(1:l-2)
-  ELSE IF(s(i+2:i+2) == '_') THEN                                              !--- <var><sep><pha>_<tag> => return <var>_<tag>
-    out = s(1:i-1)//s(i+2:l)
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: ix, ip, ns
+  out = s; ns = LEN_TRIM(s)
+  IF(s == '' .OR. ns<=2) RETURN                                                !--- Empty string or LEN(name)<=2: nothing to do
+  IF(s(1:3)=='H2O' .AND. INDEX(old_phases,s(4:4))/=0 .AND. (ns==4 .OR. s(5:5)=='_')) THEN
+    out='H2O'//s(5:ns)                                                         !--- H2O<phase>[_<iso>][_<tag>]
+  ELSE IF(s(ns-1:ns-1)==phases_sep .AND. INDEX(known_phases,s(ns:ns))/=0) THEN
+    out = s(1:ns-2); RETURN                                                    !--- <var><phase_sep><phase>
+  ELSE; DO ip = 1, nphases; ix = INDEX(s, phases_sep//known_phases(ip:ip)//'_'); IF(ix /= 0) EXIT; END DO
+    IF(ip /= nphases+1) out = s(1:ix-1)//s(ix+2:ns)                            !--- <var><phase_sep><phase>_<tag>
   END IF
 END FUNCTION delPhase
-!------------------------------------------------------------------------------------------------------------------------------
+!==============================================================================================================================
 CHARACTER(LEN=maxlen) FUNCTION addPhase_s1(s,pha) RESULT(out)
   CHARACTER(LEN=*),           INTENT(IN) :: s
   CHARACTER(LEN=1),           INTENT(IN) :: pha
+!------------------------------------------------------------------------------------------------------------------------------
   INTEGER :: l, i
   out = s
@@ -1347,24 +1776,26 @@
   IF(i /= 0) out = s(1:i-1)//phases_sep//pha//'_'//s(i+1:l)                    !--- <var>_<tag> => return <var><sep><pha>_<tag>
 END FUNCTION addPhase_s1
-!------------------------------------------------------------------------------------------------------------------------------
+!==============================================================================================================================
 FUNCTION addPhase_sm(s,pha) RESULT(out)
   CHARACTER(LEN=*),           INTENT(IN) :: s(:)
   CHARACTER(LEN=1),           INTENT(IN) :: pha
   CHARACTER(LEN=maxlen),     ALLOCATABLE :: out(:)
+!------------------------------------------------------------------------------------------------------------------------------
   INTEGER :: k
   out = [( addPhase_s1(s(k), pha), k=1, SIZE(s) )]
 END FUNCTION addPhase_sm
-!------------------------------------------------------------------------------------------------------------------------------
+!==============================================================================================================================
 CHARACTER(LEN=maxlen) FUNCTION addPhase_i1(s,ipha,phases) RESULT(out)
   CHARACTER(LEN=*),           INTENT(IN) :: s
   INTEGER,                    INTENT(IN) :: ipha
   CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: phases
+!------------------------------------------------------------------------------------------------------------------------------
   out = s
   IF(s == '') RETURN                                                           !--- Empty string: nothing to do
-  IF(ipha==0) RETURN                                                           !--- Null index: no phase to add
+  IF(ipha == 0 .OR. ipha > nphases) RETURN                                     !--- Absurd index: no phase to add
   IF(     PRESENT(phases)) out = addPhase_s1(s,       phases(ipha:ipha))
   IF(.NOT.PRESENT(phases)) out = addPhase_s1(s, known_phases(ipha:ipha))
 END FUNCTION addPhase_i1
-!------------------------------------------------------------------------------------------------------------------------------
+!==============================================================================================================================
 FUNCTION addPhase_im(s,ipha,phases) RESULT(out)
   CHARACTER(LEN=*),           INTENT(IN) :: s(:)
@@ -1372,9 +1803,10 @@
   CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: phases
   CHARACTER(LEN=maxlen),     ALLOCATABLE :: out(:)
+!------------------------------------------------------------------------------------------------------------------------------
   INTEGER :: k
   IF(     PRESENT(phases)) out = [( addPhase_i1(s(k), ipha,       phases), k=1, SIZE(s) )]
   IF(.NOT.PRESENT(phases)) out = [( addPhase_i1(s(k), ipha, known_phases), k=1, SIZE(s) )]
 END FUNCTION addPhase_im
-!------------------------------------------------------------------------------------------------------------------------------
+!==============================================================================================================================
 
 
@@ -1385,15 +1817,17 @@
   CHARACTER(LEN=*),           INTENT(IN)  :: tname
   CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: phases
+!------------------------------------------------------------------------------------------------------------------------------
   CHARACTER(LEN=maxlen) :: phase
   IF(     PRESENT(phases)) phase = getPhase(tname,       phases, iPhase)
   IF(.NOT.PRESENT(phases)) phase = getPhase(tname, known_phases, iPhase)
 END FUNCTION getiPhase
-!------------------------------------------------------------------------------------------------------------------------------
+!==============================================================================================================================
 CHARACTER(LEN=1) FUNCTION getPhase(tname, phases, iPhase) RESULT(phase)
   CHARACTER(LEN=*),           INTENT(IN)  :: tname
   CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: phases
   INTEGER,          OPTIONAL, INTENT(OUT) :: iPhase
+!------------------------------------------------------------------------------------------------------------------------------
   INTEGER :: ip
-  phase = TRIM(strHead(strTail(tname, phases_sep, .TRUE.), '_', .TRUE.))
+  phase = TRIM(strHead(strTail(tname, phases_sep), '_', .TRUE.))     !--- <nam><sep><pha>[_<tag>] -> <pha>[_<tag>] -> <pha>
   IF(     PRESENT(phases)) ip = INDEX(      phases, phase)
   IF(.NOT.PRESENT(phases)) ip = INDEX(known_phases, phase)
@@ -1401,75 +1835,76 @@
   IF(PRESENT(iPhase)) iPhase = ip
 END FUNCTION getPhase
-!------------------------------------------------------------------------------------------------------------------------------
-
-
-!------------------------------------------------------------------------------------------------------------------------------
-CHARACTER(LEN=maxlen) FUNCTION old2newName_1(oldName, iPhase) RESULT(newName)
-  !--- Convert an old style name into a new one.
-  !    Only usable with old style "traceur.def" files, in which only water isotopes are allowed.
-  !    In these files, H2O descendants names are: H2O<phase>[_<isotope>][_<tag>], with:
-  !    phase = v, l or i ; isotope = eau, HDO, O18, O17 or HTO.
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!============ CONVERT WATER-DERIVED NAMES FROM FORMER TO CURRENT CONVENTION ; OTHER NAMES ARE LEFT UNTOUCHED ==================
+!======= NAMES STRUCTURE: H2O[<phase>][_<isotope>][_<tag>] (<phase> from "old_phases", <isotope> from "oldH2OIso") ============
+!==============================================================================================================================
+CHARACTER(LEN=maxlen) FUNCTION old2newH2O_1(oldName, iPhase) RESULT(newName)
   CHARACTER(LEN=*),  INTENT(IN)  :: oldName
   INTEGER, OPTIONAL, INTENT(OUT) :: iPhase
+!------------------------------------------------------------------------------------------------------------------------------
   CHARACTER(LEN=maxlen), ALLOCATABLE :: tmp(:)
-  INTEGER :: ix, ip, it, nt
-  LOGICAL :: lerr, lH2O
+  INTEGER :: ix, ip, nt
+  LOGICAL :: lerr
   newName = oldName
   IF(PRESENT(iPhase)) iPhase = 1                                               !--- Default: gaseous phase
-  lH2O=.FALSE.
-  IF(LEN_TRIM(oldName) > 3) THEN
-    lH2O = oldName(1:3)=='H2O' .AND. INDEX(old_phases,oldName(4:4))/=0         !--- H2O<phase>*,  with phase=="v", "l", "i" or "r"
-    IF(LEN_TRIM(oldName) > 4) lH2O = lH2O .AND. oldName(5:5) == '_'            !--- H2O<phase>_*, with phase=="v", "l", "i" or "r"
+  lerr = strParse(oldName, '_', tmp, nt)                                       !--- Parsing: 1 up to 3 elements.
+  ip = strIdx( [('H2O'//old_phases(ix:ix), ix=1, nphases)], tmp(1) )           !--- Phase index
+  IF(ip /= 0 .AND. PRESENT(iPhase)) iPhase = ip                                !--- Returning phase index
+  IF(ip == 0 .AND. tmp(1) /= 'H2O')   RETURN                                   !--- Not an old-style water-related species
+  IF(nt == 1) THEN
+    newName = addPhase('H2O',ip)                                               !=== WATER WITH OR WITHOUT PHASE
+  ELSE
+    ix = strIdx(oldH2OIso, tmp(2))                                             !--- Index in the known isotopes list
+    IF(ix /= 0) newName = newH2OIso(ix)                                        !--- Move to new isotope name
+    IF(ip /= 0) newName = addPhase(newName, ip)                                !--- Add phase to isotope name
+    IF(nt == 3) newName = TRIM(newName)//'_'//TRIM(tmp(3))                     !=== WATER ISOTOPE OR TAGGING TRACER
   END IF
-  IF(.NOT.lH2O) RETURN
-  IF(LEN_TRIM(oldName)>3) THEN; IF(INDEX(old_Phases,oldName(4:4))==0) RETURN; END IF 
-  lerr = strParse(oldName, '_', tmp, n=nt)
-  ip = strIdx([('H2O'//old_phases(ip:ip), ip=1, nphases)], tmp(1))             !--- Phase index (/=0 if any)
-  IF(PRESENT(iPhase)) iPhase = ip
-  newName = addPhase('H2O', ip)                                                !--- Water
-  IF(nt == 1) RETURN                                                           !--- Water: finished
-  ix = strIdx(oldH2OIso, tmp(2))                                               !--- Index in the known isotopes list
-  IF(ix == 0) newName = addPhase(tmp(2),        ip)                            !--- Not an isotope
-  IF(ix /= 0) newName = addPhase(newH2OIso(ix), ip)                            !--- Isotope
-  IF(nt == 3) newName = TRIM(newName)//'_'//TRIM(tmp(3))                       !--- Tagging tracer
-END FUNCTION old2newName_1
-!------------------------------------------------------------------------------------------------------------------------------
-FUNCTION old2newName_m(oldName, iPhase) RESULT(newName)
-  CHARACTER(LEN=*),  INTENT(IN)  :: oldName(:)
+END FUNCTION old2newH2O_1
+!==============================================================================================================================
+FUNCTION old2newH2O_m(oldName) RESULT(newName)
+  CHARACTER(LEN=*), INTENT(IN) :: oldName(:)
+  CHARACTER(LEN=maxlen)        :: newName(SIZE(oldName))
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: i
+  newName = [(old2newH2O_1(oldName(i)), i=1, SIZE(oldName))]
+END FUNCTION old2newH2O_m
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!============ CONVERT WATER-DERIVED NAMES FROM CURRENT TO FORMER CONVENTION ; OTHER NAMES ARE LEFT UNTOUCHED ==================
+!==== NAMES STRUCTURE: <var>[<phase_sep><phase>][_<tag>] (<phase> from "known_phases", <var> = 'H2O' or from "newH2OIso") =====
+!==============================================================================================================================
+CHARACTER(LEN=maxlen) FUNCTION new2oldH2O_1(newName, iPhase) RESULT(oldName)
+  CHARACTER(LEN=*),  INTENT(IN)  :: newName
   INTEGER, OPTIONAL, INTENT(OUT) :: iPhase
-  CHARACTER(LEN=maxlen)          :: newName(SIZE(oldName))
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: ix, ip
+  CHARACTER(LEN=maxlen) :: var
+  oldName = newName
+  ip = getiPhase(newName)                                                      !--- Phase index
+  IF(PRESENT(iPhase)) iPhase = MAX(ip, 1)                                      !--- Return phase index ; default: 1 (gazeous)
+  var = TRIM(strHead(newName, phases_sep, .TRUE.))                             !--- Variable without phase and tag
+  ix = strIdx(newH2OIso, var)                                                  !--- Index in the known H2O isotopes list
+  IF(ix == 0 .AND. var /= 'H2O') RETURN                                        !--- Neither H2O nor an H2O isotope => finished
+  oldName = 'H2O'
+  IF(ip /= 0) oldName = TRIM(oldName)//old_phases(ip:ip)                       !--- Add phase if needed
+  IF(ix /= 0) oldName = TRIM(oldName)//'_'//oldH2OIso(ix)                      !--- H2O isotope name
+  IF(newName /= addPhase(var, ip)) &
+    oldName = TRIM(oldName)//strTail(newName, '_', .TRUE.)                     !--- Add the tag suffix
+  IF(ip == 0 .AND. ix /= 0) oldName = strTail(oldName, '_')                    !--- Isotope with no phase: remove 'H2O_' prefix
+END FUNCTION new2oldH2O_1
+!==============================================================================================================================
+FUNCTION new2oldH2O_m(newName) RESULT(oldName)
+  CHARACTER(LEN=*), INTENT(IN) :: newName(:)
+  CHARACTER(LEN=maxlen)        :: oldName(SIZE(newName))
+!------------------------------------------------------------------------------------------------------------------------------
   INTEGER :: i
-  newName = [(old2newName_1(oldName(i), iPhase), i=1, SIZE(oldName))]
-END FUNCTION old2newName_m
-!------------------------------------------------------------------------------------------------------------------------------
-
-!------------------------------------------------------------------------------------------------------------------------------
-CHARACTER(LEN=maxlen) FUNCTION new2oldName_1(newName, iPhase) RESULT(oldName)
-  !--- Convert a new style name into an old one.
-  !    Only convertable names are water descendants names H2O_<phase>, <isotope>_<phase>, <isotope>_<phase>_<tag>, with:
-  !    phase = g, l or s ; isotope = H2[16]O, H[2]O, H2<[18]O, H2[17]O or H[3]O.
-  CHARACTER(LEN=*),  INTENT(IN)    :: newName
-  INTEGER, OPTIONAL, INTENT(OUT)   :: iPhase
-  INTEGER :: ix, ip, it, nt
-  LOGICAL :: lH2O
-  CHARACTER(LEN=maxlen) :: tag
-  ix = strIdx([(addPhase('H2O',ip), ip=1, nphases)], newName)                  !--- Phase index for H2O_<phase>
-  IF(ix /= 0) THEN; oldName = 'H2O'//old_phases(ix:ix); RETURN; END IF         !--- H2O_<phase> case
-  ix = strIdx(newH2OIso, strHead(newName, phases_sep, .TRUE.))                 !--- Isotope index
-  IF(ix == 0) THEN; oldName = newName;                  RETURN; END IF         !--- Not a water descendant
-  ip = getiPhase(newName)                                                      !--- Phase index
-  oldName = TRIM(oldH2OIso(ix))//old_phases(ip:ip)                             !--- <isotope>_<phase>
-  tag = strTail(delPhase(newName), TRIM(newH2OIso(ix)))                        !--- Get "_<tag>" if any
-  IF(tag /= delPhase(newName) .AND. tag /= '') oldName = TRIM(oldName)//tag    !--- Tagging tracer
-END FUNCTION new2oldName_1
-!------------------------------------------------------------------------------------------------------------------------------
-FUNCTION new2oldName_m(newName, iPhase) RESULT(oldName)
-  CHARACTER(LEN=*),  INTENT(IN)  :: newName(:)
-  INTEGER, OPTIONAL, INTENT(OUT) :: iPhase
-  CHARACTER(LEN=maxlen)          :: oldName(SIZE(newName))
-  INTEGER :: i
-  oldName = [(new2oldName_1(newName(i), iPhase), i=1, SIZE(newName))]
-END FUNCTION new2oldName_m
-!------------------------------------------------------------------------------------------------------------------------------
+  oldName = [(new2oldH2O_1(newName(i)), i=1, SIZE(newName))]
+END FUNCTION new2oldH2O_m
+!==============================================================================================================================
 
 
@@ -1477,56 +1912,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 ANCESTOR(S) OF TRACER(S) "tname" AT GENERATION "igen"  IN THE TRACERS DESCRIPTORS LIST "tr" =====
-!==============================================================================================================================
-INTEGER FUNCTION idxAncestor_1(t, tname, igen) RESULT(out)
-! Return the name of the generation "igen" (>=0) ancestor of "tname"
-  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/branches/LMDZ_ECRad/libf/misc/strings_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_ECRad/libf/misc/strings_mod.F90	(revision 4203)
+++ LMDZ6/branches/LMDZ_ECRad/libf/misc/strings_mod.F90	(revision 4482)
@@ -10,5 +10,5 @@
   PUBLIC :: is_numeric, bool2str, int2str, real2str, dble2str
   PUBLIC :: reduceExpr, str2bool, str2int, str2real, str2dble
-  PUBLIC :: addQuotes, testFile, checkList, removeComment, test
+  PUBLIC :: addQuotes, checkList, removeComment, test
 
   INTERFACE get_in;     MODULE PROCEDURE getin_s,  getin_i,  getin_r,  getin_l;  END INTERFACE get_in
@@ -21,5 +21,4 @@
   INTERFACE strIdx;     MODULE PROCEDURE     strIdx_1,                 strIdx_m; END INTERFACE strIdx
   INTERFACE strCount;   MODULE PROCEDURE  strCount_m1, strCount_11, strCount_1m; END INTERFACE strCount
-  INTERFACE strParse;   MODULE PROCEDURE   strParse_1,               strParse_m; END INTERFACE strParse
   INTERFACE strReplace; MODULE PROCEDURE strReplace_1,             strReplace_m; END INTERFACE strReplace
   INTERFACE cat;        MODULE PROCEDURE   horzcat_s1,  horzcat_i1,  horzcat_r1, &
@@ -30,5 +29,4 @@
   INTERFACE reduceExpr;   MODULE PROCEDURE   reduceExpr_1,   reduceExpr_m; END INTERFACE reduceExpr
   INTERFACE addQuotes;    MODULE PROCEDURE    addQuotes_1,    addQuotes_m; END INTERFACE addQuotes
-  INTERFACE testFile;     MODULE PROCEDURE     testFile_1,     testFile_m; END INTERFACE testFile
 
   INTEGER, PARAMETER :: maxlen    = 256                    !--- Standard maximum length for strings
@@ -38,4 +36,5 @@
 CONTAINS
 
+!==============================================================================================================================
 LOGICAL FUNCTION test(lcond, lout) RESULT(lerr)
   LOGICAL, INTENT(IN)  :: lcond
@@ -43,9 +42,12 @@
   lerr = lcond; lout = lcond
 END FUNCTION test
+!==============================================================================================================================
+
 
 !==============================================================================================================================
 SUBROUTINE init_printout(lunout_, prt_level_)
   INTEGER, INTENT(IN) :: lunout_, prt_level_
-  lunout = lunout_
+  lunout    = lunout_
+  prt_level = prt_level_
 END SUBROUTINE init_printout
 !==============================================================================================================================
@@ -105,4 +107,5 @@
   LOGICAL,          OPTIONAL, INTENT(IN) :: ll
   INTEGER,          OPTIONAL, INTENT(IN) :: unit
+!------------------------------------------------------------------------------------------------------------------------------
   CHARACTER(LEN=maxlen) :: subn
   INTEGER :: unt
@@ -121,4 +124,5 @@
   INTEGER,          OPTIONAL, INTENT(IN) :: unit
   INTEGER,          OPTIONAL, INTENT(IN) :: nmax
+!------------------------------------------------------------------------------------------------------------------------------
   CHARACTER(LEN=maxlen), ALLOCATABLE :: s(:)
   CHARACTER(LEN=maxlen) :: subn
@@ -138,4 +142,5 @@
   LOGICAL,          OPTIONAL, INTENT(IN) :: ll
   INTEGER,          OPTIONAL, INTENT(IN) :: unit
+!------------------------------------------------------------------------------------------------------------------------------
   CHARACTER(LEN=maxlen) :: subn
   INTEGER :: unt
@@ -152,4 +157,5 @@
   INTEGER,          OPTIONAL, INTENT(IN) :: unit
   INTEGER,          OPTIONAL, INTENT(IN)  :: nmax
+!------------------------------------------------------------------------------------------------------------------------------
   CHARACTER(LEN=maxlen) :: subn
   INTEGER :: unt, nmx
@@ -187,69 +193,73 @@
 
 !==============================================================================================================================
-!=== Extract the substring in front of the last (first if lFirst==TRUE) occurrence of separator "sep" in "str" ================
+!=== Extract the substring in front of the first (last if lBackward==TRUE) occurrence of "sep" in "str"        ================
 !=== Examples for str='a_b_c' and sep='_' and the bash command with the same effect:                           ================
 !===    * strHead(..,.FALSE.) = 'a'           ${str%%$sep*}                                                    ================
 !===    * strHead(..,.TRUE.)  = 'a_b'         ${str%$sep*}                                                     ================
 !==============================================================================================================================
-CHARACTER(LEN=maxlen) FUNCTION strHead_1(str,sep,lFirst) RESULT(out)
+CHARACTER(LEN=maxlen) FUNCTION strHead_1(str, sep, lBackward) RESULT(out)
   CHARACTER(LEN=*),           INTENT(IN) :: str
   CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: sep
-  LOGICAL,          OPTIONAL, INTENT(IN) :: lFirst
-  LOGICAL :: lf
-  lf = .FALSE.; IF(PRESENT(lFirst)) lf=lFirst
+  LOGICAL,          OPTIONAL, INTENT(IN) :: lBackward
+!------------------------------------------------------------------------------------------------------------------------------
   IF(PRESENT(sep)) THEN
-    out = str(1:INDEX(str,sep,.NOT.lf)-1)
+    IF(     PRESENT(lBackWard)) out = str(1:INDEX(str,sep,lBackWard)-1)
+    IF(.NOT.PRESENT(lBackWard)) out = str(1:INDEX(str,sep)-1)
   ELSE
-    out = str(1:INDEX(str,'/',.NOT.lf)-1)
+    IF(     PRESENT(lBackWard)) out = str(1:INDEX(str,'/',lBackWard)-1)
+    IF(.NOT.PRESENT(lBackWard)) out = str(1:INDEX(str,'/')-1)
   END IF
   IF(out == '') out = str
 END FUNCTION strHead_1
 !==============================================================================================================================
-FUNCTION strHead_m(str,sep,lFirst) RESULT(out)
+FUNCTION strHead_m(str, sep, lBackward) RESULT(out)
   CHARACTER(LEN=maxlen),     ALLOCATABLE :: out(:)
   CHARACTER(LEN=*),           INTENT(IN) :: str(:)
   CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: sep
-  LOGICAL,          OPTIONAL, INTENT(IN) :: lFirst
-  LOGICAL :: lf
+  LOGICAL,          OPTIONAL, INTENT(IN) :: lBackward
+!------------------------------------------------------------------------------------------------------------------------------
   INTEGER :: k
-  lf = .FALSE.; IF(PRESENT(lFirst)) lf=lFirst
   IF(PRESENT(sep)) THEN
-    out = [(strHead_1(str(k), sep,   lf), k=1, SIZE(str))]
+    IF(     PRESENT(lBackWard)) out = [(strHead_1(str(k), sep, lBackWard), k=1, SIZE(str))]
+    IF(.NOT.PRESENT(lBackWard)) out = [(strHead_1(str(k), sep),            k=1, SIZE(str))]
   ELSE
-    out = [(strHead_1(str(k), lFirst=lf), k=1, SIZE(str))]
+    IF(     PRESENT(lBackWard)) out = [(strHead_1(str(k), '/', lBackWard), k=1, SIZE(str))]
+    IF(.NOT.PRESENT(lBackWard)) out = [(strHead_1(str(k), '/'),            k=1, SIZE(str))]
   END IF
 END FUNCTION strHead_m
 !==============================================================================================================================
-!=== Extract the substring following the last (first if lFirst==TRUE) occurrence of separator "sep" in "str"   ================
+!=== Extract the substring following the first (last if lBackward==TRUE) occurrence of "sep" in "str"          ================
 !=== Examples for str='a_b_c' and sep='_' and the bash command with the same effect:                           ================
-!===    * strTail(..,.FALSE.) = 'c'           ${str#*$sep}                                                     ================
-!===    * strTail(..,.TRUE.)  = 'b_c'         ${str##*$sep}                                                    ================
-!==============================================================================================================================
-CHARACTER(LEN=maxlen) FUNCTION strTail_1(str,sep,lFirst) RESULT(out)
+!===    * strTail(str, '_', .FALSE.) = 'b_c'         ${str#*$sep}                                              ================
+!===    * strTail(str, '_', .TRUE.)  = 'c'           ${str##*$sep}                                             ================
+!==============================================================================================================================
+CHARACTER(LEN=maxlen) FUNCTION strTail_1(str, sep, lBackWard) RESULT(out)
   CHARACTER(LEN=*),           INTENT(IN) :: str
   CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: sep
-  LOGICAL,          OPTIONAL, INTENT(IN) :: lFirst
-  LOGICAL :: lf
-  lf = .FALSE.; IF(PRESENT(lFirst)) lf=lFirst
+  LOGICAL,          OPTIONAL, INTENT(IN) :: lBackWard
+!------------------------------------------------------------------------------------------------------------------------------
   IF(PRESENT(sep)) THEN
-    out = str(INDEX(str,sep,.NOT.lf)+LEN(sep):LEN_TRIM(str))
+    IF(     PRESENT(lBackWard)) out = str(INDEX(str,sep,lBackWard)+LEN(sep):LEN_TRIM(str))
+    IF(.NOT.PRESENT(lBackWard)) out = str(INDEX(str,sep)          +LEN(sep):LEN_TRIM(str))
   ELSE
-    out = str(INDEX(str,'/',.NOT.lf)+1:LEN_TRIM(str))
+    IF(     PRESENT(lBackWard)) out = str(INDEX(str,'/',lBackWard)+1:LEN_TRIM(str))
+    IF(.NOT.PRESENT(lBackWard)) out = str(INDEX(str,'/')          +1:LEN_TRIM(str))
   END IF
   IF(out == '') out = str
 END FUNCTION strTail_1
 !==============================================================================================================================
-FUNCTION strTail_m(str,sep,lFirst) RESULT(out)
+FUNCTION strTail_m(str, sep, lBackWard) RESULT(out)
   CHARACTER(LEN=maxlen),     ALLOCATABLE :: out(:)
   CHARACTER(LEN=*),           INTENT(IN) :: str(:)
   CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: sep
-  LOGICAL,          OPTIONAL, INTENT(IN) :: lFirst
-  LOGICAL :: lf
+  LOGICAL,          OPTIONAL, INTENT(IN) :: lBackWard
+!------------------------------------------------------------------------------------------------------------------------------
   INTEGER :: k
-  lf = .FALSE.; IF(PRESENT(lFirst)) lf=lFirst
   IF(PRESENT(sep)) THEN
-    out = [(strTail_1(str(k), sep,   lf), k=1, SIZE(str))]
+    IF(     PRESENT(lBackWard)) out = [(strTail_1(str(k), sep, lBackWard), k=1, SIZE(str))]
+    IF(.NOT.PRESENT(lBackWard)) out = [(strTail_1(str(k), sep),            k=1, SIZE(str))]
   ELSE
-    out = [(strTail_1(str(k), lFirst=lf), k=1, SIZE(str))]
+    IF(     PRESENT(lBackWard)) out = [(strTail_1(str(k), '/', lBackWard), k=1, SIZE(str))]
+    IF(.NOT.PRESENT(lBackWard)) out = [(strTail_1(str(k), '/'),            k=1, SIZE(str))]
   END IF
 END FUNCTION strTail_m
@@ -265,4 +275,5 @@
   CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: sep
   LOGICAL,          OPTIONAL, INTENT(IN) :: mask(:)
+!------------------------------------------------------------------------------------------------------------------------------
   CHARACTER(LEN=:), ALLOCATABLE :: s
   INTEGER :: is, i0
@@ -285,4 +296,5 @@
   CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: sep
   INTEGER,          OPTIONAL, INTENT(IN) :: nmax
+!------------------------------------------------------------------------------------------------------------------------------
   CHARACTER(LEN=maxlen), ALLOCATABLE :: t(:)
   CHARACTER(LEN=maxlen) :: sp
@@ -339,4 +351,5 @@
   CHARACTER(LEN=*), ALLOCATABLE, INTENT(INOUT) :: str(:)
   INTEGER,          OPTIONAL,    INTENT(OUT)   :: nb
+!------------------------------------------------------------------------------------------------------------------------------
   CHARACTER(LEN=maxlen), ALLOCATABLE :: s1(:)
   INTEGER :: k, n, n1
@@ -355,4 +368,5 @@
   CHARACTER(LEN=*),   ALLOCATABLE, INTENT(INOUT) :: str1(:)
   CHARACTER(LEN=*),                INTENT(IN)    :: str2(:)
+!------------------------------------------------------------------------------------------------------------------------------
   CHARACTER(LEN=maxlen), ALLOCATABLE :: s1(:), s2(:)
   INTEGER :: k
@@ -374,5 +388,6 @@
 
 !==============================================================================================================================
-!=== GET THE INDEX OF THE FIRST APPEARANCE IN THE STRING VECTOR "str(:)" OF THE STRING(s) "s" =================================
+!=== GET THE INDEX OF THE FIRST APPEARANCE IN THE STRING VECTOR "str(:)" OF THE STRING(s) "s[(:)]" ============================
+!=== OPTIONALY: GET THE NUMBER OF FOUND ELEMENTS "n". NB: UNFOUND => INDEX=0                       ============================
 !==============================================================================================================================
 INTEGER FUNCTION strIdx_1(str, s) RESULT(out)
@@ -386,4 +401,5 @@
   INTEGER, OPTIONAL, INTENT(OUT) :: n
   INTEGER,           ALLOCATABLE :: out(:)
+!------------------------------------------------------------------------------------------------------------------------------
   INTEGER :: k
   out = [(strIdx_1(str(:), s(k)), k=1, SIZE(s))]
@@ -400,4 +416,5 @@
   INTEGER, OPTIONAL, INTENT(OUT) :: n
   INTEGER,           ALLOCATABLE :: out(:)
+!------------------------------------------------------------------------------------------------------------------------------
   INTEGER :: k
   out = PACK( [(k, k=1, SIZE(str(:), DIM=1))], MASK = str(:) == s )
@@ -409,4 +426,5 @@
   INTEGER, OPTIONAL, INTENT(OUT) :: n
   INTEGER,           ALLOCATABLE :: out(:)
+!------------------------------------------------------------------------------------------------------------------------------
   INTEGER :: k
   out = PACK( [(k, k=1, SIZE(i(:), DIM=1))], MASK = i(:) == j )
@@ -418,15 +436,9 @@
   INTEGER, OPTIONAL, INTENT(OUT) :: n
   INTEGER,           ALLOCATABLE :: out(:)
+!------------------------------------------------------------------------------------------------------------------------------
   INTEGER :: k
   out = PACK( [(k, k=1, SIZE(l(:), DIM=1))], MASK = l(:) )
   IF(PRESENT(n)) n = SIZE(out(:), DIM=1)
 END FUNCTION find_boo
-!==============================================================================================================================
-
-
-
-!==============================================================================================================================
-!=== GET 1ST APPEARANCE INDEX OF EACH ELEMENT OF "t(:)" IN "s(:)" (UNFOUND: INDEX=0) ==========================================
-!=== OPTIONALY: GET THE NUMBER OF FOUND ELEMENTS "n"                                 ==========================================
 !==============================================================================================================================
 
@@ -444,54 +456,41 @@
   INTEGER,           INTENT(OUT) :: idel                             !--- Index of the identified delimiter (0 if idx==0)
   LOGICAL, OPTIONAL, INTENT(IN)  :: lSc                              !--- Care about nbs with front sign or in scient. notation
-
-  INTEGER               :: idx0                                      !--- Used to display an identified non-numeric string
-  INTEGER, ALLOCATABLE  :: ii(:)
-  LOGICAL               :: ll, ls
-  CHARACTER(LEN=maxlen) :: d
-!  modname = 'strIdx'
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: idx0                                                    !--- Used to display an identified non-numeric string
   lerr = .FALSE.
-  idx = strIdx1(rawList, del, ibeg, idel)                            !--- del(idel) appears in "rawList" at position idx
+  idx = strIdx1(rawList, del, ibeg, idel)                            !--- idx/=0: del(idel) is at position "idx" in "rawList" 
   IF(.NOT.PRESENT(lSc))               RETURN                         !--- No need to check exceptions for numbers => finished
   IF(.NOT.        lSc )               RETURN                         !--- No need to check exceptions for numbers => finished
-  IF(idx == 0) THEN                                                  !--- No element of "del" in "rawList":
-    lerr = .NOT.is_numeric(rawList(ibeg:))                           !--- String must be a number
-    IF(lerr) idx = LEN_TRIM(rawList); RETURN                         !--- Update idx => rawList(ibeg:idx-1) is the whole string
-  END IF
-  idx0 = idx
-  IF(test(idx==1.AND.INDEX('+-',del(idel))/=0, lerr)) RETURN         !--- Front separator different from +/-: error
-  IF(idx/=1.AND.is_numeric(rawList(ibeg:idx-1)))      RETURN         !--- The input string tail is a valid number
-  idx = strIdx1(rawList, del, idx+1, idel)                           !---   => TO THE NEXT DELIMITER
+
+  !=== No delimiter found: the whole string must be a valid number
+  IF(idx == 0) THEN                                                  !--- No element of "del" in "rawList"
+    lerr = .NOT.is_numeric(rawList(ibeg:LEN_TRIM(rawList)))          !--- String must be a number
+    IF(lerr) idx = LEN_TRIM(rawList); RETURN                         !--- Set idx so that rawList(ibeg:idx-1) = whole string
+  END IF
+
+  IF(test(idx == 1 .AND. INDEX('+-',del(idel)) /= 0, lerr)) RETURN   !--- The front delimiter is different from +/-: error
+  IF(     idx /= 1 .AND. is_numeric(rawList(ibeg:idx-1)))   RETURN   !--- The input string head is a valid number
+
+  !=== The string part in front of the 1st delimiter is not a valid number: search for next delimiter index "idx"
+  idx0 = idx ; idx = strIdx1(rawList, del, idx+1, idel)              !--- Keep start index because idx is recycled
   IF(idx == 0) THEN
-    lerr = .NOT.is_numeric(rawList(ibeg:))                           !--- No delimiter detected: string must be a number
+    lerr = .NOT.is_numeric(rawList(ibeg:LEN_TRIM(rawList)))          !--- No other delimiter: whole string must be a valid numb
     IF(lerr) idx = idx0; RETURN
   END IF
-  idx0 = idx
-  IF(is_numeric(rawList(ibeg:idx-1)))                 RETURN         !--- The input string tail is a valid number
-  IF(test(          INDEX('eE',rawList(idx-1:idx-1)) /= 0  &         !--- Sole possible exception: scientific notation: E+/-
-               .OR. INDEX('+-',del(idel)) /=0, lerr)) RETURN
-  idx = strIdx1(rawList, del, idx+1, idel)                           !---   => TO THE NEXT DELIMITER
-  IF(idx == 0) THEN
-    lerr = .NOT.is_numeric(rawList(ibeg:))                           !--- No separator detected: string must be a number
-    IF(lerr) idx = idx0; RETURN
-  END IF
   lerr = .NOT.is_numeric(rawList(ibeg:idx-1))
+
 CONTAINS
 
-INTEGER FUNCTION strIdx1(str, del, ib, id) RESULT(idx)
-  CHARACTER(LEN=*),  INTENT(IN)  :: str
-  CHARACTER(LEN=*),  INTENT(IN)  :: del(:)
+!------------------------------------------------------------------------------------------------------------------------------
+INTEGER FUNCTION strIdx1(str, del, ib, id) RESULT(i)
+!--- Get the index of the first appereance of one of the delimiters "del(:)" in "str" starting from position "ib".
+!--- "id" is the index in "del(:)" of the first delimiter found.
+  IMPLICIT NONE
+  CHARACTER(LEN=*),  INTENT(IN)  :: str, del(:)
   INTEGER,           INTENT(IN)  :: ib
   INTEGER,           INTENT(OUT) :: id
-
-  INTEGER              :: nd, ns, i
-  INTEGER, ALLOCATABLE :: ii(:)
-
-  nd  = SIZE(del)                                                    !--- Number of separators
-  ns  = LEN_TRIM(str)                                                !--- Length of the raw chain
-  ii  = [(INDEX( str(ib:ns), del(i) ), i = 1, nd)]                   !--- Determine the next separator start index
-  id  =  MINLOC( ii, MASK = ii /= 0, DIM = 1 )                       !--- Current delimiter index in the "delimiter(:)" list
-  idx = 0
-  IF(ANY(ii /= 0)) idx = MINVAL( ii, MASK = ii /= 0 ) + ib - 1       !--- Index in "str(1:ns)" of the delimiter first character
-  IF(idx == 0) id = 0
+!------------------------------------------------------------------------------------------------------------------------------
+  DO i = ib, LEN_TRIM(str); id = strIdx(del, str(i:i)); IF(id /= 0) EXIT; END DO
+  IF(i > LEN_TRIM(str)) THEN; i = 0; id = 0; END IF
 END FUNCTION strIdx1
 
@@ -501,42 +500,30 @@
 
 !==============================================================================================================================
-!=== Return the index of first appearance of "del" in "str" starting from index "ib"
-!==============================================================================================================================
-INTEGER FUNCTION strIndex(str, del, ib) RESULT(idx)
-  CHARACTER(LEN=*),  INTENT(IN)  :: str
-  CHARACTER(LEN=*),  INTENT(IN)  :: del
-  INTEGER,           INTENT(IN)  :: ib
-  idx  = INDEX( str(ib:LEN_TRIM(str)), del ) + ib -1
-END FUNCTION strIndex
-!==============================================================================================================================
-
-
-!==============================================================================================================================
 !=== Count the number of elements separated by "delimiter" in list "rawList". =================================================
 !==============================================================================================================================
-LOGICAL FUNCTION strCount_11(rawList, delimiter, nb, lSc) RESULT(out)
+LOGICAL FUNCTION strCount_11(rawList, delimiter, nb, lSc) RESULT(lerr)
   CHARACTER(LEN=*),  INTENT(IN)  :: rawList
   CHARACTER(LEN=*),  INTENT(IN)  :: delimiter
   INTEGER,           INTENT(OUT) :: nb
   LOGICAL, OPTIONAL, INTENT(IN)  :: lSc
+!------------------------------------------------------------------------------------------------------------------------------
   LOGICAL :: ll
   ll = .FALSE.; IF(PRESENT(lSc)) ll = lSc
-  out = strCount_1m(rawList, [delimiter], nb, ll)
+  lerr = strCount_1m(rawList, [delimiter], nb, ll)
 END FUNCTION strCount_11
 !==============================================================================================================================
-LOGICAL FUNCTION strCount_m1(rawList, delimiter, nb, lSc) RESULT(out)
+LOGICAL FUNCTION strCount_m1(rawList, delimiter, nb, lSc) RESULT(lerr)
   CHARACTER(LEN=*),     INTENT(IN)  :: rawList(:)
   CHARACTER(LEN=*),     INTENT(IN)  :: delimiter
   INTEGER, ALLOCATABLE, INTENT(OUT) :: nb(:)
   LOGICAL,    OPTIONAL, INTENT(IN)  :: lSc
-
+!------------------------------------------------------------------------------------------------------------------------------
   LOGICAL :: ll
   INTEGER :: id
-
   ll  = .FALSE.; IF(PRESENT(lSc)) ll = lSc .AND. INDEX('+-', delimiter) /= 0
-  out = .TRUE.
+  lerr = .TRUE.
   ALLOCATE(nb(SIZE(rawList)))
   DO id = 1, SIZE(rawList)
-    out = out .AND. strCount_1m(rawList(id), [delimiter], nb(id), ll)
+    lerr = lerr .AND. strCount_1m(rawList(id), [delimiter], nb(id), ll)
   END DO
 END FUNCTION strCount_m1
@@ -547,8 +534,8 @@
   INTEGER,           INTENT(OUT) :: nb
   LOGICAL, OPTIONAL, INTENT(IN)  :: lSc
+!------------------------------------------------------------------------------------------------------------------------------
   INTEGER              :: ib, ie, jd, nr
   LOGICAL              :: ll
   CHARACTER(LEN=1024)  :: r
-!  modname = 'strCount'
   lerr = .FALSE.
   ll   = .FALSE.; IF(PRESENT(lSc)) ll = lSc
@@ -559,5 +546,5 @@
     lerr = strIdx_prv(r, delimiter, ib, ie, jd, ll)
     IF(fmsg('"'//TRIM(r(ib:ie-1))//'" is not numeric', ll=lerr)) RETURN
-    IF(jd == 0) EXIT
+    IF(ie == 0 .OR. jd == 0) EXIT
     ib = ie + LEN(delimiter(jd))
     DO WHILE(r(ib:ib) == ' ' .AND. ib < nr); ib = ib + 1; END DO     !--- Skip spaces before next chain
@@ -572,31 +559,74 @@
 !===          Corresponding "vals" remains empty if the element does not contain "=" sign. ====================================
 !==============================================================================================================================
-LOGICAL FUNCTION strParse_1(rawList, delimiter, keys, lSc, vals, n) RESULT(lerr)
+LOGICAL FUNCTION strParse(rawList, delimiter, keys, n, vals) RESULT(lerr)
   CHARACTER(LEN=*),                             INTENT(IN)  :: rawList, delimiter
   CHARACTER(LEN=maxlen), ALLOCATABLE,           INTENT(OUT) :: keys(:)
-  LOGICAL,                            OPTIONAL, INTENT(IN)  :: lSc      !--- Take care about numbers in scientific notation
+  INTEGER,                            OPTIONAL, INTENT(OUT) :: n
   CHARACTER(LEN=maxlen), ALLOCATABLE, OPTIONAL, INTENT(OUT) :: vals(:)
-  INTEGER,                            OPTIONAL, INTENT(OUT) :: n
-  LOGICAL :: ll
-!  modname = 'strParse'
-  ll = .FALSE.; IF(PRESENT(lSc)) ll = lSc
-  IF(.NOT.PRESENT(vals)) lerr = strParse_m(rawList, [delimiter], keys, ll)
-  IF(     PRESENT(vals)) lerr = strParse_m(rawList, [delimiter], keys, ll, vals)
-  IF(PRESENT(n)) n = SIZE(keys)
-END FUNCTION strParse_1
-!==============================================================================================================================
-LOGICAL FUNCTION strParse_m(rawList, delimiter, keys, lSc, vals, n, id) RESULT(lerr)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=1024) :: r
+  INTEGER :: nr, nk
+  lerr = .FALSE.
+  r  = TRIM(ADJUSTL(rawList))
+  nr = LEN_TRIM(r); IF(nr == 0) THEN; keys = ['']; RETURN; END IF
+  CALL strParse_prv(nk)                                              !--- COUNT THE ELEMENTS
+  ALLOCATE(keys(nk))
+  IF(PRESENT(vals)) THEN
+    ALLOCATE(vals(nk)); CALL strParse_prv(nk, keys, vals)            !--- PARSE THE KEYS
+  ELSE
+    CALL strParse_prv(nk, keys)                                      !--- PARSE THE KEYS
+  END IF
+  IF(PRESENT(n)) n = nk
+
+CONTAINS
+
+!------------------------------------------------------------------------------------------------------------------------------
+SUBROUTINE strParse_prv(nkeys, keys, vals)
+!--- * Get the number of elements after parsing ("nkeys" only is present)
+!--- * Parse the <key>=<val> pairs and store result in "keys" and "vals" (already allocated)
+  IMPLICIT NONE
+  INTEGER,                         INTENT(OUT) :: nkeys
+  CHARACTER(LEN=maxlen), OPTIONAL, INTENT(OUT) :: keys(:)
+  CHARACTER(LEN=maxlen), OPTIONAL, INTENT(OUT) :: vals(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: ib, ie
+  nkeys = 1; ib = 1
+  DO
+    ie = INDEX(rawList(ib:nr), delimiter)+ib-1                       !--- Determine the next separator start index
+    IF(ie == ib-1) EXIT
+    IF(PRESENT(keys)) keys(nkeys) = r(ib:ie-1)                       !--- Get the ikth key
+    IF(PRESENT(vals)) CALL parseKeys(keys(nkeys), vals(nkeys))       !--- Parse the ikth <key>=<val> pair
+    ib = ie + LEN(delimiter)
+    DO WHILE(r(ib:ib) == ' ' .AND. ib < nr); ib = ib + 1; END DO     !--- Skip spaces before next chain
+    nkeys = nkeys+1
+  END DO
+  IF(PRESENT(keys)) keys(nkeys) = r(ib:nr)                           !--- Get the last key
+  IF(PRESENT(vals)) CALL parseKeys(keys(nkeys), vals(nkeys))         !--- Parse the last <key>=<val> pair
+END SUBROUTINE strParse_prv
+
+!------------------------------------------------------------------------------------------------------------------------------
+SUBROUTINE parseKeys(key, val)
+  CHARACTER(LEN=*), INTENT(INOUT) :: key
+  CHARACTER(LEN=*), INTENT(OUT)   :: val
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: ix
+  ix = INDEX(key, '='); IF(ix == 0) RETURN                           !--- First "=" index in "key"
+  val = ADJUSTL(key(ix+1:LEN_TRIM(key)))
+  key = ADJUSTL(key(1:ix-1))
+END SUBROUTINE parseKeys
+
+END FUNCTION strParse
+!==============================================================================================================================
+LOGICAL FUNCTION strParse_m(rawList, delimiter, keys, n, vals, lSc, id) RESULT(lerr)
   CHARACTER(LEN=*),                             INTENT(IN)  :: rawList, delimiter(:)
   CHARACTER(LEN=maxlen),           ALLOCATABLE, INTENT(OUT) :: keys(:)  !--- Parsed keys vector
+  INTEGER,               OPTIONAL,              INTENT(OUT) :: n        !--- Length of the parsed vector
+  CHARACTER(LEN=maxlen), OPTIONAL, ALLOCATABLE, INTENT(OUT) :: vals(:)  !--- Values for <name>=<value> keys
   LOGICAL,               OPTIONAL,              INTENT(IN)  :: lSc      !--- Take care about numbers in scientific notation
-  CHARACTER(LEN=maxlen), OPTIONAL, ALLOCATABLE, INTENT(OUT) :: vals(:)  !--- Values for <name>=<value> keys
-  INTEGER,               OPTIONAL,              INTENT(OUT) :: n        !--- Length of the parsed vector
   INTEGER,               OPTIONAL, ALLOCATABLE, INTENT(OUT) :: id(:)    !--- Indexes of the separators in "delimiter(:)" vector
-
+!------------------------------------------------------------------------------------------------------------------------------
   CHARACTER(LEN=1024) :: r
   INTEGER :: nr, ik, nk, ib, ie, jd
   LOGICAL :: ll
-
-!  modname = 'strParse'
   ll = .FALSE.; IF(PRESENT(lSc)) ll = lSc
   IF(test(fmsg("Couldn't parse list: non-numerical strings were found", ll=strCount_1m(rawList, delimiter, nk, ll)),lerr)) RETURN
@@ -625,7 +655,9 @@
 CONTAINS
 
+!------------------------------------------------------------------------------------------------------------------------------
 SUBROUTINE parseKeys(key, val)
   CHARACTER(LEN=*), INTENT(INOUT) :: key
   CHARACTER(LEN=*), INTENT(OUT)   :: val
+!------------------------------------------------------------------------------------------------------------------------------
   INTEGER :: ix
   ix = INDEX(key, '='); IF(ix == 0) RETURN                           !--- First "=" index in "key"
@@ -645,9 +677,7 @@
   CHARACTER(LEN=*),  INTENT(IN)    :: key, val   !--- "key" will be replaced by "val"
   LOGICAL, OPTIONAL, INTENT(IN)    :: lsurr      !--- TRUE => key must be surrounded by special characters to be substituted
-
-  CHARACTER(LEN=1024) :: s, t
+!------------------------------------------------------------------------------------------------------------------------------
   INTEGER :: i0, ix, nk, ns
   LOGICAL :: lsur, lb, le
-
   lsur = .FALSE.; IF(PRESENT(lsurr)) lsur = lsurr
   nk = LEN_TRIM(key)
@@ -688,4 +718,5 @@
   CHARACTER(LEN=*), OPTIONAL, TARGET, INTENT(IN) :: s1, s2, s3, s4, s5, s6, s7, s8, s9
   CHARACTER(LEN=maxlen), ALLOCATABLE :: out(:)
+!------------------------------------------------------------------------------------------------------------------------------
   CHARACTER(LEN=maxlen), POINTER     :: s
   LOGICAL :: lv(10)
@@ -707,4 +738,5 @@
   CHARACTER(LEN=*), OPTIONAL, TARGET, DIMENSION(:), INTENT(IN) :: s1, s2, s3, s4, s5, s6, s7, s8, s9
   CHARACTER(LEN=maxlen), ALLOCATABLE :: out(:,:)
+!------------------------------------------------------------------------------------------------------------------------------
   CHARACTER(LEN=maxlen), POINTER     :: s(:)
   LOGICAL :: lv(10)
@@ -729,4 +761,5 @@
   INTEGER, OPTIONAL, TARGET, INTENT(IN) :: i1, i2, i3, i4, i5, i6, i7, i8, i9
   INTEGER, ALLOCATABLE :: out(:)
+!------------------------------------------------------------------------------------------------------------------------------
   INTEGER, POINTER     :: i
   LOGICAL :: lv(10)
@@ -748,4 +781,5 @@
   INTEGER, OPTIONAL, TARGET, DIMENSION(:), INTENT(IN) :: i1, i2, i3, i4, i5, i6, i7, i8, i9
   INTEGER, ALLOCATABLE :: out(:,:)
+!------------------------------------------------------------------------------------------------------------------------------
   INTEGER, POINTER     :: i(:)
   LOGICAL :: lv(10)
@@ -770,4 +804,5 @@
   REAL, OPTIONAL, TARGET, INTENT(IN) :: r1, r2, r3, r4, r5, r6, r7, r8, r9
   REAL, ALLOCATABLE :: out(:)
+!------------------------------------------------------------------------------------------------------------------------------
   REAL, POINTER     :: r
   LOGICAL :: lv(10)
@@ -789,4 +824,5 @@
   REAL, OPTIONAL, TARGET, DIMENSION(:), INTENT(IN) :: r1, r2, r3, r4, r5, r6, r7, r8, r9
   REAL, ALLOCATABLE :: out(:,:)
+!------------------------------------------------------------------------------------------------------------------------------
   REAL, POINTER     :: r(:)
   LOGICAL :: lv(10)
@@ -811,4 +847,5 @@
   DOUBLE PRECISION, OPTIONAL, TARGET, INTENT(IN) :: d1, d2, d3, d4, d5, d6, d7, d8, d9
   DOUBLE PRECISION, ALLOCATABLE :: out(:)
+!------------------------------------------------------------------------------------------------------------------------------
   DOUBLE PRECISION, POINTER     :: d
   LOGICAL :: lv(10)
@@ -830,4 +867,5 @@
   DOUBLE PRECISION, OPTIONAL, TARGET, DIMENSION(:), INTENT(IN) :: d1, d2, d3, d4, d5, d6, d7, d8, d9
   DOUBLE PRECISION, ALLOCATABLE :: out(:,:)
+!------------------------------------------------------------------------------------------------------------------------------
   DOUBLE PRECISION, POINTER     :: d(:)
   LOGICAL :: lv(10)
@@ -869,5 +907,5 @@
   INTEGER,          OPTIONAL, INTENT(IN)  :: unit          !--- Output unit (default: screen)
   CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: sub           !--- Subroutine name
-
+!------------------------------------------------------------------------------------------------------------------------------
   CHARACTER(LEN=2048) :: row
   CHARACTER(LEN=maxlen)  :: rFm, el, subn
@@ -879,5 +917,4 @@
   INTEGER, PARAMETER   :: nm=1                             !--- Space between values & columns
   LOGICAL :: ls, li, lr
-
   subn = '';    IF(PRESENT(sub)) subn = sub
   rFm = '*';    IF(PRESENT(rFmt)) rFm = rFmt               !--- Specified format for reals
@@ -954,5 +991,5 @@
       nr = LEN_TRIM(row)-1                                           !--- Final separator removed
       CALL msg(row(1:nr), subn, unit=unt)
-      IF(ir /= 1) CYCLE                                              !--- Titles are underlined
+      IF(ir /= 1) CYCLE                                              !--- Titles only are underlined
       row=''; DO ic=1,nHd; row=TRIM(row)//REPEAT('-',n(ic))//'+'; END DO
       DO ic = ib,ncmx(it); row=TRIM(row)//REPEAT('-',n(ic))//'+'; END DO
@@ -975,5 +1012,5 @@
   CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: rFmt          !--- Format for reals
   LOGICAL,          OPTIONAL, INTENT(IN)  :: llast         !--- Last variable: no final ','
-
+!------------------------------------------------------------------------------------------------------------------------------
   CHARACTER(LEN=maxlen)  :: rFm, el
   CHARACTER(LEN=maxlen), ALLOCATABLE :: d(:,:)
@@ -983,6 +1020,4 @@
   INTEGER, ALLOCATABLE :: n(:)
   LOGICAL :: ls, li, lr, la
-
-!  modname = 'dispNamelist'
   rFm = '*';    IF(PRESENT(rFmt)) rFm = rFmt               !--- Specified format for reals
   ls = PRESENT(s); li = PRESENT(i); lr = PRESENT(r)
@@ -1056,5 +1091,4 @@
   REAL,                       INTENT(IN)  ::  a(:)                   !--- Linearized array of values
   INTEGER,                    INTENT(IN)  ::  n(:)                   !--- Profile before linearization
-
   CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: err_msg, nam(:), subn   !--- Error message, variables and calling subroutine names
   INTEGER,          OPTIONAL, INTENT(IN)  :: nRowMax                 !--- Maximum number of lines to display    (default: all)
@@ -1062,4 +1096,5 @@
   INTEGER,          OPTIONAL, INTENT(IN)  :: nHead                   !--- Number of front columns to duplicate  (default: 1)
   INTEGER,          OPTIONAL, INTENT(IN)  :: unit                    !--- Output unit                           (def: lunout)
+!------------------------------------------------------------------------------------------------------------------------------
   CHARACTER(LEN=maxlen),      ALLOCATABLE :: ttl(:)
   LOGICAL,                    ALLOCATABLE :: m(:)
@@ -1139,10 +1174,9 @@
   INTEGER,          OPTIONAL, INTENT(IN)  :: nHead                   !--- Number of front columns to duplicate  (default: 1)
   INTEGER,          OPTIONAL, INTENT(IN)  :: unit                    !--- Output unit                           (def: lunout)
-
+!------------------------------------------------------------------------------------------------------------------------------
   CHARACTER(LEN=maxlen)                   :: mes, sub, fm='(f12.9)', prf
   CHARACTER(LEN=maxlen),      ALLOCATABLE :: ttl(:), vnm(:)
-  LOGICAL,                    ALLOCATABLE :: m(:)
   INTEGER,                    ALLOCATABLE :: ki(:), kj(:), kl(:)
-  INTEGER                                 :: i, j, k, rk, ib, ie, itr, nm, nv, unt, nRmx, nCmx, nHd, rk1
+  INTEGER                                 :: i, j, k, rk, nv, unt, nRmx, nCmx, nHd
   REAL,                       ALLOCATABLE :: val(:,:)
 
@@ -1160,5 +1194,5 @@
   lerr= SIZE(a,1) /= PRODUCT(n); IF(fmsg('profile "n" does not match "a" and "ll"', sub, lerr, unt)) RETURN
 
-  SELECT CASE(rk1)                                                   !--- Indices list
+  SELECT CASE(rk)                                                   !--- Indices list
     CASE(0); IF(ll(1)) THEN; WRITE(unt,'(a,", ",a," = ",2f12.9)')TRIM(vnm(1)),TRIM(vnm(2)),a(1,1),a(1,2); RETURN; END IF
     CASE(1); ki = [  (i,i=1,n(1)) ]
@@ -1189,5 +1223,5 @@
   CHARACTER(LEN=*),      INTENT(IN)  :: str
   CHARACTER(LEN=maxlen), INTENT(OUT) :: val
-
+!------------------------------------------------------------------------------------------------------------------------------
   CHARACTER(LEN=maxlen)              :: v
   CHARACTER(LEN=1024)                :: s, vv
@@ -1196,5 +1230,4 @@
   INTEGER :: nl, nn, i, j, im, ix
   LOGICAL :: ll
-!  modname = 'reduceExpr_1'
   s = str
 
@@ -1243,15 +1276,13 @@
   CHARACTER(LEN=maxlen), ALLOCATABLE :: ky(:)
   CHARACTER(LEN=1),      ALLOCATABLE :: op(:)
-
+!------------------------------------------------------------------------------------------------------------------------------
   CHARACTER(LEN=1024) :: s
   DOUBLE PRECISION :: v, vm, vp
   INTEGER      :: i, ni, io
-
-!  modname = 'reduceExpr_basic'
   lerr = .FALSE.
   IF(is_numeric(str)) THEN; val=TRIM(str); RETURN; END IF
   op = ['^','/','*','+','-']                                                   !--- List of recognized operations
   s = str
-  IF(test(strParse_m(s, op, ky, .TRUE., id = id), lerr)) RETURN                !--- Parse the values
+  IF(test(strParse_m(s, op, ky, lSc=.TRUE., id = id), lerr)) RETURN            !--- Parse the values
   vl = str2dble(ky)                                                            !--- Conversion to doubles
   lerr = ANY(vl >= HUGE(1.d0))
@@ -1262,10 +1293,10 @@
       IF(id(i) /= io) CYCLE                                                    !--- Current found operator is not op(io)
       vm = vl(i); vp = vl(i+1)                                                 !--- Couple of values used for current operation
-      SELECT CASE(io)                                                          !--- Perform operation on the two values
-        CASE(1); v = vm**vp  !--- ^
-        CASE(2); v = vm/vp   !--- /
-        CASE(3); v = vm*vp   !--- *
-        CASE(4); v = vm+vp   !--- +
-        CASE(5); v = vm-vp   !--- +
+      SELECT CASE(op(io))                                                          !--- Perform operation on the two values
+        CASE('^'); v = vm**vp
+        CASE('/'); v = vm/vp
+        CASE('*'); v = vm*vp
+        CASE('+'); v = vm+vp
+        CASE('-'); v = vm-vp
       END SELECT
       IF(i == ni) THEN; vl = [vl(1:ni-1), v]; ELSE; vl = [vl(1:i-1), v, vl(i+2:ni+1)]; END IF
@@ -1283,4 +1314,5 @@
   CHARACTER(LEN=*),                   INTENT(IN)  :: str(:)
   CHARACTER(LEN=maxlen), ALLOCATABLE, INTENT(OUT) :: val(:)
+!------------------------------------------------------------------------------------------------------------------------------
   INTEGER :: i
   ALLOCATE(lerr(SIZE(str)),val(SIZE(str)))
@@ -1302,5 +1334,4 @@
   READ(str,fmt,IOSTAT=e) x
   out = e==0 .AND. INDEX('Ee',str(LEN_TRIM(str):LEN_TRIM(str)))==0
-  IF(str == '') out = .FALSE.
 END FUNCTION is_numeric
 !==============================================================================================================================
@@ -1348,4 +1379,5 @@
   INTEGER,           INTENT(IN) :: i
   INTEGER, OPTIONAL, INTENT(IN) :: nDigits
+!------------------------------------------------------------------------------------------------------------------------------
   WRITE(out,*)i
   out = ADJUSTL(out)
@@ -1357,4 +1389,5 @@
   REAL,                       INTENT(IN) :: r
   CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: fmt
+!------------------------------------------------------------------------------------------------------------------------------
   IF(     PRESENT(fmt)) WRITE(out,fmt)r
   IF(.NOT.PRESENT(fmt)) WRITE(out, * )r
@@ -1365,4 +1398,5 @@
   DOUBLE PRECISION,           INTENT(IN) :: d
   CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: fmt
+!------------------------------------------------------------------------------------------------------------------------------
   IF(     PRESENT(fmt)) WRITE(out,fmt)d
   IF(.NOT.PRESENT(fmt)) WRITE(out, * )d
@@ -1385,4 +1419,5 @@
 !==============================================================================================================================
 
+
 !==============================================================================================================================
 FUNCTION addQuotes_1(s) RESULT(out)
@@ -1395,4 +1430,5 @@
   CHARACTER(LEN=*), INTENT(IN)  :: s(:)
   CHARACTER(LEN=:), ALLOCATABLE :: out(:)
+!------------------------------------------------------------------------------------------------------------------------------
   INTEGER :: k, n
   n = MAXVAL(LEN_TRIM(s), MASK=.TRUE.)
@@ -1406,27 +1442,9 @@
   CHARACTER(LEN=*), INTENT(IN) :: s
   CHARACTER(LEN=1) :: b, e
+!------------------------------------------------------------------------------------------------------------------------------
   out = .TRUE.; IF(TRIM(s) == '') RETURN
   b = s(1:1); e = s(MAX(1,LEN_TRIM(s)):MAX(1,LEN_TRIM(s)))
   out = .NOT.is_numeric(s) .AND. (b /= "'" .OR. e /= "'") .AND. ( b /= '"' .OR. e /= '"')
 END FUNCTION needQuotes
-!==============================================================================================================================
-
-
-!==============================================================================================================================
-!=== TEST WHETHER A FILE IS PRESENT OR NOT ====================================================================================
-!==============================================================================================================================
-LOGICAL FUNCTION testFile_1(fname) RESULT(out)
-  CHARACTER(LEN=*), INTENT(IN) :: fname
-  INTEGER :: ierr
-  OPEN(90, FILE=fname, FORM='formatted', STATUS='old', IOSTAT=ierr); CLOSE(99)
-  out = ierr/=0
-END FUNCTION testFile_1
-!==============================================================================================================================
-FUNCTION testFile_m(fname) RESULT(out)
-  LOGICAL,         ALLOCATABLE ::   out(:)
-  CHARACTER(LEN=*), INTENT(IN) :: fname(:)
-  INTEGER :: k
-  out = [(testFile_1(fname(k)), k=1, SIZE(fname))]
-END FUNCTION testFile_m
 !==============================================================================================================================
 
@@ -1442,4 +1460,5 @@
   CHARACTER(LEN=*),   INTENT(IN)  :: message, items, reason
   INTEGER,  OPTIONAL, INTENT(IN)  :: nmax
+!------------------------------------------------------------------------------------------------------------------------------
   CHARACTER(LEN=maxlen), ALLOCATABLE :: s(:)
   INTEGER :: i, nmx
@@ -1464,4 +1483,3 @@
 
 
-
 END MODULE strings_mod
Index: LMDZ6/branches/LMDZ_ECRad/libf/misc/trac_types_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_ECRad/libf/misc/trac_types_mod.F90	(revision 4203)
+++ 	(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/branches/LMDZ_ECRad/libf/misc/wxios.F90
===================================================================
--- LMDZ6/branches/LMDZ_ECRad/libf/misc/wxios.F90	(revision 4203)
+++ LMDZ6/branches/LMDZ_ECRad/libf/misc/wxios.F90	(revision 4482)
@@ -308,4 +308,8 @@
         IF (xios_is_valid_axis("axis_lat")) THEN
            CALL xios_set_axis_attr( "axis_lat", n_glo=nbp_lat, n=jj_nb, begin=jj_begin-1, value=io_lat(jj_begin:jj_end))
+        ENDIF
+        IF (xios_is_valid_axis("axis_lat_greordered")) THEN
+           CALL xios_set_axis_attr( "axis_lat_greordered", n_glo=nbp_lat, n=jj_nb, begin=jj_begin-1, &
+                                    value=io_lat(jj_begin:jj_end)*(-1.))
         ENDIF
 
