Index: LMDZ6/branches/LMDZ-tracers/libf/misc/readTracFiles_mod.f90
===================================================================
--- LMDZ6/branches/LMDZ-tracers/libf/misc/readTracFiles_mod.f90	(revision 3852)
+++ LMDZ6/branches/LMDZ-tracers/libf/misc/readTracFiles_mod.f90	(revision 3852)
@@ -0,0 +1,1259 @@
+MODULE readTracFiles_mod
+
+  USE strings_mod,    ONLY: msg, testFile,  strFind, strStack, strCount,   strHead, removeComment, dispTable, fmsg, &
+                            cat, checkList, strIdx,  strParse, strReplace, strTail,  reduceExpr, modname, find, test
+  USE trac_types_mod, ONLY : tra, iso, db, kys
+
+  IMPLICIT NONE
+
+  PRIVATE
+
+  PUBLIC :: readTracersFiles, aliasTracer, tracersSubset, indexUpdate     !--- TOOLS ASSOCIATED TO TRACERS  DESCRIPTORS
+  PUBLIC :: readIsotopesFile                                              !--- TOOLS ASSOCIATED TO ISOTOPES DESCRIPTORS
+  PUBLIC :: getKey_init, getKey, setDirectKeys                            !--- FUNCTIONS TO GET KEYS FROM tracers & isotopes
+
+  PUBLIC :: known_phases, nphases, phases_names, delPhase, addPhase       !--- PHASES RELATED VARIABLES AND ROUTINES
+
+  PUBLIC :: tran0, idxAncestor, ancestor                                  !--- GEN 0 TRACER + TOOLS FOR GENERATIONS
+
+!------------------------------------------------------------------------------------------------------------------------------
+  INTERFACE getKey
+    MODULE PROCEDURE getKeyByName_s1, getKeyByName_i1, getKeyByName_r1, getKeyByName_sm, getKeyByName_im, getKeyByName_rm
+  END INTERFACE getKey
+!------------------------------------------------------------------------------------------------------------------------------
+  INTERFACE tracersSubset; MODULE PROCEDURE trSubset_Indx, trSubset_Name, trSubset_Nam1; 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    addPhase;   MODULE PROCEDURE    addPhase_1,    addPhase_m; END INTERFACE    addPhase
+!------------------------------------------------------------------------------------------------------------------------------
+
+  !=== MAIN DATABASE: files sections descriptors
+  TYPE(db), SAVE, ALLOCATABLE, TARGET :: dBase(:)
+
+  !--- SOME PARAMETERS THAT ARE NOT LIKELY TO CHANGE OFTEN
+  CHARACTER(LEN=256), SAVE      :: tran0        = 'air'              !--- Default transporting fluid
+  CHARACTER(LEN=256), PARAMETER :: known_phases = 'gls'              !--- Known phases initials
+  INTEGER,            PARAMETER :: nphases = LEN_TRIM(known_phases)  !--- Number of phases
+  CHARACTER(LEN=256), SAVE      :: phases_names(nphases) &           !--- Known phases names
+                                 = ['gaseous', 'liquid ', 'solid  ']
+  LOGICAL,            SAVE      :: tracs_merge = .TRUE.              !--- Merge/stack tracers lists
+  LOGICAL,            SAVE      :: lSortByGen  = .TRUE.              !--- Sort by growing generation
+
+  !=== LOCAL COPIES OF THE TRACERS AND ISOTOPES DESCRIPTORS, USED BY getKey (INITIALIZED IN getKey_init)
+  TYPE(tra), ALLOCATABLE, TARGET, SAVE ::  tracers(:)
+  TYPE(iso), ALLOCATABLE, TARGET, SAVE :: isotopes(:)
+
+
+CONTAINS
+
+!==============================================================================================================================
+!==============================================================================================================================
+!=== READ ONE OR SEVERAL TRACER FILES AND FILL A "tr" TRACERS DESCRIPTOR DERIVED TYPE.
+!=== THE RETURN 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":
+!===    First line: <nb tracers>     Other lines: <hadv> <vadv> <tracer name> [<parent name>]
+!===  2: A  "NEW STYLE" TRACERS FILE  "tracer.def" WITH SEVERAL SECTIONS.
+!===  3: SEVERAL  "  "  TRACERS FILES "tracer_<component>.def" WITH A SINGLE SECTION IN EACH.
+!=== REMARKS:
+!===  * EACH SECTION BEGINS WITH A "&<section name> LINE
+!===  * DEFAULT VALUES FOR ALL THE SECTIONS OF THE FILE ARE DEFINED IN THE SPECIAL SECTION "&default"
+!===  * EACH SECTION LINE HAS THE STRUCTURE:  <name(s)>  <key1>=<value1> <key2>=<value2> ...
+!===  * SO FAR, THE DEFINED KEYS ARE: parent, phases, hadv, vadv, type
+!===  * <name> AND <parent> CAN BE LISTS OF COMA-SEPARATED TRACERS ; THE ROUTINE EXPAND THESE FACTORIZATIONS.
+!=== FUNCTION RETURN VALUE "lerr" IS FALSE IN CASE SOMETHING WENT WRONG.
+!=== ABOUT THE KEYS:
+!     * The "keys" component (of type kys) is in principle enough to store everything we could need.
+!     But some variables are stored as direct-access keys to make the code more readable and because they are used often.
+!     * Most of the direct-access keys are set in this module, but some are not (lnam, iadv and itr for example).
+!     * Some of the direct-access keys must be updated (using the routine "setDirectKeys") is a subset of "tracers(:)"
+!     is extracted: the indexes are no longer valid for a subset (examples: tracers(:)iprnt or tracers(:)%ichld).
+!     * 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(tra), ALLOCATABLE, INTENT(OUT) :: tracs(:)
+  CHARACTER(LEN=256),     ALLOCATABLE ::  s(:), sections(:), trac_files(:)
+  CHARACTER(LEN=256)   :: str, fname, mesg
+  INTEGER              :: is, nsec, ierr, it, ntrac, ns
+  LOGICAL, ALLOCATABLE :: ll(:), lGen3(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  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
+
+  !--- TELLS WHAT WAS IS ABOUT TO BE USED
+  IF(test(fmsg(fType==0, 'No adequate tracers description file(s) found ; default values will be used'), lerr)) RETURN
+  CALL msg(fType==1, 'Trying to read old-style tracers description file "traceur.def"')
+  CALL msg(fType==2, 'Trying to read the new style multi-sections tracers description file "tracer.def"')
+  CALL msg(fType==3, 'Trying to read the new style single section tracers description files "tracer_*.def"')
+
+  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+  IF(fType==1) THEN                                                  !=== OLD FORMAT "traceur.def"
+  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+    !--- OPEN THE "traceur.def" FILE
+    OPEN(90, FILE="traceur.def", FORM='formatted', STATUS='old', IOSTAT=ierr)
+
+    !--- GET THE TRACERS NUMBER
+    READ(90,'(i3)',IOSTAT=ierr)ntrac                                 !--- Number of lines/tracers
+    IF(test(fmsg(ierr /= 0, 'Invalid format for "'//TRIM(fname)//'"'), lerr)) RETURN
+
+    !--- READ THE REMAINING LINES: <hadv> <vadv> <tracer> [<transporting fluid>]
+    ALLOCATE(tracs(ntrac))
+    DO it=1,ntrac                                                    !=== READ RAW DATA: loop on the line/tracer number
+      READ(90,'(a)',IOSTAT=ierr) str
+      IF(test(fmsg(ierr>0, 'Invalid format for "' //TRIM(fname)//'"'), lerr)) RETURN
+      IF(test(fmsg(ierr<0, 'Not enough lines in "'//TRIM(fname)//'"'), lerr)) RETURN
+      ll = strParse(str, ' ', s, n=ns)
+      tracs(it)%keys%key = ['hadv', 'vadv']
+      tracs(it)%keys%val = s(1:2)
+      CALL msg(ns == 3 .AND. it == 1, 'This file is for air tracers only')
+      CALL msg(ns == 4 .AND. it == 1, 'This files specifies the transporting fluid')
+      SELECT CASE(s(3))                                              !--- name and phase of the tracer
+        CASE('H2Ov'); tracs(it)%name = 'H2O-g'; tracs(it)%phas = 'g'
+        CASE('H2Ol'); tracs(it)%name = 'H2O-l'; tracs(it)%phas = 'l'
+        CASE('H2Oi'); tracs(it)%name = 'H2O-s'; tracs(it)%phas = 's'
+        CASE DEFAULT; tracs(it)%name = s(3)   ; tracs(it)%phas = 'g'
+      END SELECT
+      tracs(it)%prnt = tran0                                         !--- Default transporting fluid: Air
+      IF(ns == 4) tracs(it)%prnt = s(4)                              !--- Transporting fluid name
+    END DO
+    CLOSE(90)
+
+    lGen3 = tracs%igen==3
+    CALL setGeneration(tracs)                                        !--- Determine tracs(:)%igen values
+    IF(test(checkTracers(tracs,      fname,fname),lerr)) RETURN      !--- Detect orphans and check phases
+    IF(test(checkUnique (tracs,lGen3,fname,fname),lerr)) RETURN      !--- Detect repeated tracers
+    CALL sortTracers  (tracs)                                        !--- Sort the tracers
+    CALL setDirectKeys(tracs)                                        !--- Set the directly accessible keys
+    tracs(:)%keys%name = tracs(:)%name                               !--- Copy tracers names in keys components
+    RETURN
+  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+  END IF
+  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+  !=== USING NEW FORMAT TRACERS DESCRIPTION FILES WITH POSSIBLY SEVERAL SECTIONS
+  CALL msg(nsec > 1 .AND.      tracs_merge, 'The multiple required sections will be MERGED.')
+  CALL msg(nsec > 1 .AND. .NOT.tracs_merge, 'The multiple required sections will be CUMULATED.')
+
+  !=== FEED THE DATABASE WITH THE RAW CONTENT OF THE FILE
+  SELECT CASE(fType)
+    CASE(2); lerr = feedDBase(["tracer.def"],[type_trac])            !--- Single file, coma-separated sections list
+    CASE(3); lerr = feedDBase(trac_files,     sections  )            !--- Multiple files, one section name each file
+  END SELECT
+  IF(lerr) RETURN
+  IF(     tracs_merge) lerr = mergeTracers(dBase, tracs)             !--- MERGE    THE COMPONENTS OF THE DATABASE
+  IF(.NOT.tracs_merge) lerr = cumulTracers(dBase, tracs)             !--- CUMULATE THE COMPONENTS OF THE DATABASE
+  IF(lerr) RETURN
+  CALL setDirectKeys(tracs)                                          !--- Set the directly accessible keys
+
+END FUNCTION readTracersFiles
+!==============================================================================================================================
+
+!==============================================================================================================================
+LOGICAL FUNCTION feedDBase(fnames, snames) RESULT(lerr)
+! Purpose: Read the sections "snames(is)" (coma-separated list) from each "fnames(is)"
+!   file and create the corresponding tracers set descriptors in the database "dBase":
+! * dBase(id)%name                : section name
+! * dBase(id)%trac(:)%name        : tracers names
+! * dBase(id)%trac(it)%keys%key(:): names  of keys associated to tracer dBase(id)%trac(it)%name
+! * dBase(id)%trac(it)%keys%val(:): values of keys associated to tracer dBase(id)%trac(it)%name
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=*), INTENT(IN)  :: fnames(:)                         !--- Files names
+  CHARACTER(LEN=*), INTENT(IN)  :: snames(:)                         !--- Coma-deparated list of sections (one list each file)
+  INTEGER,  ALLOCATABLE :: ndb(:)                                    !--- Nuber 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=256)    :: fnm, snm
+  INTEGER               :: idb, i
+  LOGICAL :: ll
+!------------------------------------------------------------------------------------------------------------------------------
+
+  !=== READ THE REQUIRED SECTIONS
+  ll = strCount(snames, ',', ndb)                                    !--- Number of sections for each file
+  ALLOCATE(ixf(SUM(ndb)))
+  DO i=1, SIZE(fnames)
+    IF(test(readSections(fnames(i), snames(i), 'default'), lerr)) RETURN
+    ixf(1+SUM(ndb(1:i-1)):SUM(ndb(1:i))) = i                         !--- File index for each section of the expanded list
+  END DO
+  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+  DO idb=1,SIZE(dBase)                                               !--- LOOP ON THE LOADED SECTIONS
+  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+    fnm = fnames(ixf(idb)); snm = dBase(idb)%name
+    IF(test(expandSection(dBase(idb)%trac,   snm, fnm),lerr)) RETURN !--- EXPAND NAMES AND PARENTS LISTS
+    CALL setGeneration   (dBase(idb)%trac)                           !--- DETERMINE GENERATION NUMBER
+    IF(test(checkTracers (dBase(idb)%trac,   snm, fnm),lerr)) RETURN !--- CHECK ORPHANS AND PHASES
+    lTg = dBase(idb)%trac(:)%type == 'tag'                           !--- Flag for tagging tracers
+    IF(test(checkUnique  (dBase(idb)%trac,lTg,snm,fnm),lerr)) RETURN !--- CHECK TRACERS UNIQUENESS
+    CALL expandPhases    (dBase(idb)%trac)                           !--- EXPAND THE PHASES
+    CALL sortTracers     (dBase(idb)%trac)                           !--- SORT TRACERS
+  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+  END DO
+  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+  !=== DISPLAY BASIC INFORMATION
+  lerr = ANY([(dispTraSection('Expanded list for section "'//TRIM(dBase(idb)%name)//'"',dBase(idb)%name), idb=1, SIZE(dBase))])
+END FUNCTION feedDBase
+!------------------------------------------------------------------------------------------------------------------------------
+
+!------------------------------------------------------------------------------------------------------------------------------
+LOGICAL FUNCTION readSections(fnam,snam,defName) RESULT(lerr)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=*),           INTENT(IN) :: fnam                     !--- File name
+  CHARACTER(LEN=*),           INTENT(IN) :: snam                     !--- Coma-separated sections list
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: defName                  !--- Special section (default values) name
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE(db),           ALLOCATABLE :: tdb(:)
+  CHARACTER(LEN=256), ALLOCATABLE :: sec(:)
+  INTEGER,            ALLOCATABLE ::  ix(:)
+  INTEGER :: n0, idb, ndb
+  LOGICAL :: ll
+!------------------------------------------------------------------------------------------------------------------------------
+  n0 = SIZE(dBase) + 1                                               !--- Index for next entry in the database
+  CALL readSections_all()                                            !--- Read all the sections of file "fnam"
+  ndb= SIZE(dBase)                                                   !--- Current number of sections in the database
+  IF(PRESENT(defName)) THEN                                          !--- Add default values to all the tracers
+    DO idb=n0,ndb; CALL addDefault(dBase(idb)%trac, defName); END DO !--- and remove the virtual tracer "defName"
+  END IF
+  ll = strParse(snam, ',', keys = sec)                               !--- Requested sections names
+  ix = strIdx(dBase(:)%name, sec(:))                                 !--- Indexes of requested sections in database
+  IF(test(checkList(sec, ix == 0, 'In file "'//TRIM(fnam)//'"','section(s)', 'missing'), lerr)) RETURN
+  tdb = dBase(:); dBase = [tdb(1:n0-1),tdb(PACK(ix, MASK=ix/=0))]    !--- Keep requested sections only
+
+CONTAINS
+
+!------------------------------------------------------------------------------------------------------------------------------
+SUBROUTINE readSections_all()
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=256), ALLOCATABLE ::  s(:), v(:)
+  TYPE(tra),          ALLOCATABLE :: tt(:)
+  TYPE(tra)           :: tmp
+  CHARACTER(LEN=1024) :: str
+  CHARACTER(LEN=256)  :: secn
+  INTEGER             :: ierr, n
+!------------------------------------------------------------------------------------------------------------------------------
+  IF(.NOT.ALLOCATED(dBase)) ALLOCATE(dBase(0))
+  OPEN(90, FILE=fnam, FORM='formatted', STATUS='old')
+  DO; READ(90,'(a)', IOSTAT=ierr)str
+    IF(ierr    /= 0 ) EXIT                                           !--- Finished: error or end of file
+    IF(str(1:1)=='#') CYCLE                                          !--- Skip comments lines
+    CALL removeComment(str)                                          !--- Skip comments at the end of a line
+    IF(str     == '') CYCLE                                          !--- Skip empty line (probably at the end of the file)
+    IF(str(1:1)=='&') THEN                                           !=== SECTION HEADER LINE
+      ndb  = SIZE(dBase)                                             !--- Number of sections so far
+      secn = str(2:LEN_TRIM(str))//' '                               !--- Current section name
+      IF(ANY(dBase(:)%name == secn)) CYCLE                           !--- Already known section
+      IF(secn(1:7) == 'version')     CYCLE                           !--- Skip the "version" special section
+      ndb = ndb + 1                                                  !--- Extend database
+      ALLOCATE(tdb(ndb))
+      tdb(1:ndb-1)  = dBase
+      tdb(ndb)%name = secn
+      ALLOCATE(tdb(ndb)%trac(0))
+      CALL MOVE_ALLOC(FROM=tdb, TO=dBase)
+    ELSE                                                             !=== TRACER LINE
+      ll = strParse(str,' ', keys = s, vals = v, n = n)              !--- Parse <key>=<val> pairs
+      tt = dBase(ndb)%trac(:)
+      tmp%name = s(1); tmp%keys = kys(s(1), s(2:n), v(2:n))
+      dBase(ndb)%trac = [tt(:), tmp]
+      DEALLOCATE(tt)
+!      dBase(ndb)%trac = [dBase(ndb)%trac(:), tra(name=s(1), comp=secn, keys=kys(s(1), s(2:n), v(2:n)))]
+    END IF
+  END DO
+  CLOSE(90)
+
+END SUBROUTINE readSections_all
+!------------------------------------------------------------------------------------------------------------------------------
+
+END FUNCTION readSections
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+SUBROUTINE addDefault(t, defName)
+!------------------------------------------------------------------------------------------------------------------------------
+! Purpose: Add the keys from virtual tracer named "defName" (if any) and remove this virtual tracer.
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE(tra), ALLOCATABLE, TARGET, INTENT(INOUT) :: t(:)
+  CHARACTER(LEN=*),               INTENT(IN)    :: defName
+  INTEGER :: i0, it, k
+  TYPE(kys), POINTER     :: k0
+  TYPE(tra), ALLOCATABLE :: tt(:)
+  i0 = strIdx(t(:)%name, defName)
+  IF(i0 == 0) RETURN
+  k0 => t(i0)%keys
+  DO k = 1, SIZE(k0%key)                                             !--- Loop on the keys of the tracer named "defName"
+    CALL addKey_tra(TRIM(k0%key(k)), TRIM(k0%val(k)), t)             !--- Add key to all the tracers (no overwriting)
+  END DO
+  tt = [t(1:i0-1),t(i0+1:SIZE(t))]; CALL MOVE_ALLOC(FROM=tt, TO=t)   !--- Remove the virtual tracer named "defName"
+END SUBROUTINE addDefault
+!==============================================================================================================================
+
+!==============================================================================================================================
+SUBROUTINE subDefault(t, defName, lSubLocal)
+!------------------------------------------------------------------------------------------------------------------------------
+! Purpose: Substitute the keys from virtual tracer named "defName" (if any) and remove this virtual tracer.
+!          Substitute the keys locally (for the current tracer) if the flag "lSubLocal" is .TRUE.
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE(tra), ALLOCATABLE, TARGET, INTENT(INOUT) :: t(:)
+  CHARACTER(LEN=*),               INTENT(IN)    :: defName
+  LOGICAL,                        INTENT(IN)    :: lSubLocal
+  INTEGER :: i0, it, ik
+  TYPE(kys), POINTER     :: k0, ky
+  TYPE(tra), ALLOCATABLE :: tt(:)
+  i0 = strIdx(t(:)%name, defName)
+  IF(i0 == 0) RETURN
+  k0 => t(i0)%keys 
+  DO it = 1, SIZE(t); IF(it == i0) CYCLE                             !--- Loop on the tracers
+    ky => t(it)%keys
+
+    !--- Substitute in the values of <key>=<val> pairs the keys defined in the virtual tracer "defName"
+    DO ik = 1, SIZE(k0%key); CALL strReplace(ky%val, k0%key(ik), k0%val(ik), .TRUE.); END DO
+
+    IF(.NOT.lSubLocal) CYCLE
+    !--- Substitute in the values of <key>=<val> pairs the keys defined locally (in the current tracer)
+    DO ik = 1, SIZE(ky%key); CALL strReplace(ky%val, ky%key(ik), ky%val(ik), .TRUE.); END DO
+  END DO
+  tt = [t(1:i0-1),t(i0+1:SIZE(t))]; CALL MOVE_ALLOC(FROM=tt, TO=t)   !--- Remove the virtual tracer named "defName"
+
+END SUBROUTINE subDefault
+!==============================================================================================================================
+
+
+
+
+!==============================================================================================================================
+LOGICAL FUNCTION expandSection(tr, sname, fname) RESULT(lerr)
+!------------------------------------------------------------------------------------------------------------------------------
+! Purpose: Expand tracers and parents lists in the tracers descriptor "tra".
+! Note:  * The following keys are expanded, so are accessible explicitely using "%" operator: "parent" "type".
+!        * Default values are provided for these keys because they are necessary.
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE(tra),     ALLOCATABLE, INTENT(INOUT) :: tr(:)                 !--- Tracer derived type vector
+  CHARACTER(LEN=*),           INTENT(IN)    :: sname
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)    :: fname
+  TYPE(tra),          ALLOCATABLE :: ttr(:)
+  CHARACTER(LEN=256), ALLOCATABLE :: ta(:), pa(:)
+  CHARACTER(LEN=256) :: msg1
+  INTEGER :: it, nt, iq, nq, itr, ntr, ipr, npr, i
+  LOGICAL :: ll
+  lerr = .FALSE.
+  nt = SIZE(tr)
+  nq = 0
+  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+  DO it = 1, nt    !=== GET TRACERS NB AFTER EXPANSION + NEEDED KEYS (name, parent, type)
+  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+    !--- Extract useful keys: parent name, type, component name
+    tr(it)%prnt = fgetKey(it, 'parent', tr(:)%keys,  tran0  )
+    tr(it)%type = fgetKey(it, 'type'  , tr(:)%keys, 'tracer')
+    tr(it)%comp = sname
+
+    !--- Determine the number of tracers and parents ; coherence checking
+    ll = strCount(tr(it)%name, ',', ntr)
+    ll = strCount(tr(it)%prnt, ',', npr)
+
+    !--- Tagging tracers only can have multiple parents
+    IF(test(npr/=1 .AND. TRIM(tr(it)%type)/='tag', lerr)) THEN
+      msg1 = 'Check section "'//TRIM(sname)//'"'
+      IF(PRESENT(fname)) msg1=TRIM(msg1)//' in file "'//TRIM(fname)//'"'
+      CALL msg(TRIM(msg1)//': "'//TRIM(tr(it)%name)//'" has several parents but is not a tag'); RETURN
+    END IF
+    nq = nq + ntr*npr                 
+  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+  END DO
+  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+  CALL delKey(['parent','type  '], tr)
+
+  ALLOCATE(ttr(nq))
+  iq = 1
+  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+  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)%prnt, ',', pa, n=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 = ta(itr); ttr(i)%prnt = pa(ipr)
+        ttr(i)%keys = kys(ta(itr), tr(it)%keys%key, tr(it)%keys%val)
+      END DO
+    END DO
+    ttr(iq:iq+ntr*npr-1)%type = tr(it)%type                          !--- Duplicating type
+    iq = iq + ntr*npr
+  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+  END DO
+  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+  DEALLOCATE(ta,pa)
+  CALL MOVE_ALLOC(FROM=ttr, TO=tr)
+
+END FUNCTION expandSection
+!==============================================================================================================================
+
+!==============================================================================================================================
+SUBROUTINE setGeneration(tr)
+!------------------------------------------------------------------------------------------------------------------------------
+! Purpose: Determine, for each tracer of "tr(:)":
+!   * the generation number 
+!   * the first generation ancestor name
+!------------------------------------------------------------------------------------------------------------------------------
+! Arguments:
+  TYPE(tra), INTENT(INOUT) :: tr(:)                                  !--- Tracer derived type vector
+!------------------------------------------------------------------------------------------------------------------------------
+! Local variables:
+  INTEGER :: iq, nq, ig
+  LOGICAL,            ALLOCATABLE :: lg(:)
+  CHARACTER(LEN=256), ALLOCATABLE :: prn(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  tr(:)%igen = 0                                                     !--- error if 0
+  nq = SIZE(tr, DIM=1)                                               !--- Number of tracers lines
+  lg = tr(:)%prnt == tran0                                           !--- First generation tracers flag
+  WHERE(lg) tr(:)%igen = 1                                           !--- First generation tracers
+
+  !=== Determine generation for each tracer
+  ig=0; prn = [tran0]
+  DO                                                                 !--- Update current generation flag
+    IF(ig/=0) prn = PACK( tr(:)%name, MASK=tr(:)%igen == ig)
+    lg(:) = [(ANY(prn(:) == tr(iq)%prnt), iq=1, nq)]                 !--- Current generation tracers flag
+    IF( ALL( .NOT. lg ) ) EXIT                                       !--- Empty current generation
+    ig = ig+1; WHERE(lg) tr(:)%igen = ig
+  END DO
+  tr(:)%nam1 = ancestor(tr)                                          !--- First generation ancestor name
+
+END SUBROUTINE setGeneration
+!==============================================================================================================================
+
+!==============================================================================================================================
+LOGICAL FUNCTION checkTracers(tr, sname, fname) RESULT(lerr)
+!------------------------------------------------------------------------------------------------------------------------------
+! Purpose:
+!   * check for orphan tracers (without known parent)
+!   * check wether the phases are known or not ("g"aseous, "l"iquid or "s"olid so far)
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE(tra),                  INTENT(IN) :: tr(:)                    !--- Tracer derived type vector
+  CHARACTER(LEN=*),           INTENT(IN) :: sname                    !--- Section name
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: fname                    !--- File name
+  CHARACTER(LEN=256) :: mesg
+  CHARACTER(LEN=256) :: bp(SIZE(tr, DIM=1)), pha
+  CHARACTER(LEN=1) :: p
+  INTEGER :: ip, np, iq, nq
+!------------------------------------------------------------------------------------------------------------------------------
+  nq = SIZE(tr,DIM=1)                                                !--- Number of tracers lines
+  mesg = 'Check section "'//TRIM(sname)//'"'
+  IF(PRESENT(fname)) mesg=TRIM(mesg)//' in file "'//TRIM(fname)//'"'
+
+  !=== CHECK FOR ORPHAN TRACERS
+  IF(test(checkList(tr%name, tr%igen==0, mesg, 'tracers', 'orphan'), lerr)) RETURN
+
+  !=== CHECK PHASES
+  DO iq=1,nq; IF(tr(iq)%igen/=1) CYCLE                               !--- Generation 1 only is checked
+    pha = fgetKey(iq, 'phases', tr(:)%keys, 'g')                     !--- Phases
+    np = LEN_TRIM(pha); bp(iq)=' '
+    DO ip=1,np; p = pha(ip:ip); IF(INDEX(known_phases,p)==0) bp(iq) = TRIM(bp(iq))//p; END DO
+    IF(TRIM(bp(iq)) /= '') bp(iq) = TRIM(tr(iq)%name)//': '//TRIM(bp(iq))
+  END DO
+  lerr = checkList(bp, tr%igen==1 .AND. bp/='', mesg, 'tracers phases', 'unknown')
+END FUNCTION checkTracers
+!==============================================================================================================================
+
+!==============================================================================================================================
+LOGICAL FUNCTION checkUnique(tr, lTag, sname, fname) RESULT(lerr)
+!------------------------------------------------------------------------------------------------------------------------------
+! Purpose: Make sure that tracers are not repeated.
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE(tra),                  INTENT(IN) :: tr(:)                    !--- Tracer derived type vector
+  LOGICAL,                    INTENT(IN) :: lTag(:)                  !--- Tagging tracer flag
+  CHARACTER(LEN=*),           INTENT(IN) :: sname                    !--- Section name
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: fname                    !--- File name
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: ip, np, iq, nq, k
+  LOGICAL, ALLOCATABLE :: ll(:)
+  CHARACTER(LEN=256) :: mesg, tnam, tdup(SIZE(tr,DIM=1))
+  CHARACTER(LEN=1)   :: p
+!------------------------------------------------------------------------------------------------------------------------------
+  mesg = 'Check section "'//TRIM(sname)//'"'
+  IF(PRESENT(fname)) mesg=TRIM(mesg)//' in file "'//TRIM(fname)//'"'
+  nq=SIZE(tr,DIM=1); lerr=.FALSE.                                    !--- Number of lines ; error flag
+  tdup(:) = ''
+  DO iq=1,nq; IF(lTag(iq)) CYCLE                                     !--- Tags can be repeated
+    tnam = tr(iq)%name
+    ll = tr(:)%name==tnam                                            !--- Mask for current tracer name
+    IF(COUNT(ll)==1 ) CYCLE                                          !--- Tracer is not repeated
+    IF(tr(iq)%igen>1) THEN
+      tdup(iq) = tnam                                                !--- gen>1: MUST be unique
+    ELSE
+      DO ip=1,nphases; p=known_phases(ip:ip)                         !--- Loop on known phases
+        !--- Number of appearances of the current tracer with known phase "p"
+        np = COUNT( PACK( [(INDEX(fgetKey(k, 'phases', tr(:)%keys, 'g'),p), k=1, nq)] /=0 , MASK=ll ) )
+        IF(np <=1) CYCLE
+        tdup(iq) = TRIM(tdup(iq))//TRIM(phases_names(ip))
+        IF(ANY(tdup(1:iq-1) == tdup(iq))) tdup(iq)=''                !--- Avoid repeating same messages
+      END DO
+    END IF
+    IF(tdup(iq) /= '') tdup(iq)=TRIM(tnam)//' in '//TRIM(tdup(iq))//' phase(s)'
+  END DO
+  lerr = checkList(tdup, tdup/='', mesg, 'tracers', 'duplicated')
+END FUNCTION checkUnique
+!==============================================================================================================================
+
+!==============================================================================================================================
+SUBROUTINE expandPhases(tr)
+!------------------------------------------------------------------------------------------------------------------------------
+! Purpose: Expand the phases in the tracers descriptor "tr".
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE(tra), ALLOCATABLE, INTENT(INOUT) :: tr(:)                     !--- Tracer derived type vector
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE(tra), ALLOCATABLE :: ttr(:)
+  INTEGER,  ALLOCATABLE ::  i0(:)
+  CHARACTER(LEN=256)    :: nam, pha, trn
+  INTEGER :: ip, np, iq, jq, nq, it, nt, nc, i, n
+  LOGICAL :: lTg, lEx
+!------------------------------------------------------------------------------------------------------------------------------
+  nq = SIZE(tr, DIM=1)
+  nt = 0
+  DO iq = 1, nq                                                      !--- GET THE NUMBER OF TRACERS
+    IF(tr(iq)%igen /= 1) CYCLE
+    nc = COUNT(tr(:)%nam1==tr(iq)%name .AND. tr%igen/=1)             !--- Number of childs of tr(iq)
+    tr(iq)%phas = fgetKey(iq, 'phases', tr(:)%keys)                  !--- Phases list      of tr(iq)
+    np = LEN_TRIM(tr(iq)%phas)                                       !--- Number of phases of tr(iq)
+    nt = nt + (1+nc) * np                                            !--- Number of tracers after expansion
+  END DO
+  ALLOCATE(ttr(nt))
+  it = 1                                                             !--- Current "ttr(:)" index
+  DO iq = 1, nq                                                      !--- Loop on "tr(:)" indexes
+    lTg = tr(iq)%type=='tag'                                         !--- Current tracer is a tag
+    i0 = strFind(tr(:)%name, tr(iq)%nam1, n)                         !--- Indexes of first generation ancestor copies
+    np = SUM( [( LEN_TRIM(tr(i0(i))%phas),i=1,n )],1)                !--- Number of phases for current tracer tr(iq)
+    lEx = np>1                                                       !--- Need of a phase suffix
+    IF(lTg) lEx=lEx.AND.tr(iq)%igen>1                                !--- No phase suffix for first generation tags
+    DO i=1,n                                                         !=== LOOP ON FIRST GENERATION ANCESTORS
+      jq=i0(i)                                                       !--- tr(jq): ith copy of 1st gen. ancestor of tr(iq)
+      IF(tr(iq)%igen==1) jq=iq                                       !--- Generation 1: current tracer phases only
+      pha = tr(jq)%phas                                              !--- Phases list for tr(jq)
+      DO ip=1,LEN_TRIM(pha)                                          !=== LOOP ON PHASES LISTS
+        trn=TRIM(tr(iq)%name); nam=trn                               !--- Tracer name (regular case)
+        IF(lTg) nam = TRIM(tr(iq)%prnt)                              !--- Parent name (tagging case)
+        IF(lEx) nam = TRIM(nam)//'-'//pha(ip:ip)                     !--- Phase extension needed
+        IF(lTg) nam = TRIM(nam)//'_'//TRIM(trn)                      !--- <parent>_<name> for tags
+        ttr(it) = tr(iq)                                             !--- Same <key>=<val> pairs
+        ttr(it)%name = nam                                           !--- Name with possibly phase suffix
+        ttr(it)%keys%name = nam                                      !--- Name inside the keys decriptor
+        ttr(it)%phas = pha(ip:ip)                                    !--- Single phase entry
+        IF(lEx.AND.tr(iq)%igen>1) THEN
+          ttr(it)%prnt = TRIM(ttr(it)%prnt)//'-'//pha(ip:ip)
+          ttr(it)%nam1 = TRIM(ttr(it)%nam1)//'-'//pha(ip:ip)
+        END IF
+        it=it+1
+      END DO
+      IF(tr(iq)%igen==1) EXIT                                        !--- Break phase loop for gen 1
+    END DO
+  END DO
+  CALL MOVE_ALLOC(FROM=ttr, TO=tr)
+  CALL delKey(['phases'],tr)                                         !--- Remove few keys entries
+
+END SUBROUTINE expandPhases
+!==============================================================================================================================
+
+!==============================================================================================================================
+SUBROUTINE sortTracers(tr)
+!------------------------------------------------------------------------------------------------------------------------------
+! Purpose: Sort tracers:
+!  * lGrowGen == T: in ascending generations numbers.
+!  * lGrowGen == F: tracer + its childs sorted by growing generation, one after the other.
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE(tra), ALLOCATABLE, INTENT(INOUT) :: tr(:)                     !--- Tracer derived type vector
+  INTEGER :: ig, ng, iq, jq, n, ix(SIZE(tr)), k
+  INTEGER, ALLOCATABLE :: iy(:), iz(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  iq = 1
+  IF(lSortByGen) THEN
+    ng = MAXVAL(tr(:)%igen, MASK=.TRUE., DIM=1)                      !--- Number of generations
+    DO ig = 0, ng                                                    !--- Loop on generations
+      iy = PACK([(k, k=1, SIZE(tr))], MASK=tr(:)%igen==ig)           !--- Generation ig tracers indexes
+      n = SIZE(iy)
+      ix(iq:iq+n-1) = iy                                             !--- Stack growing generations idxs
+      iq = iq + n
+    END DO
+  ELSE
+    DO jq = 1, SIZE(tr,DIM=1)                                        !--- Loop on first generation tracers
+      IF(tr(jq)%igen /= 1) CYCLE                                     !--- Skip generations >= 1
+      ix(iq) = jq                                                    !--- First generation ancestor index first
+      iq = iq + 1
+      iy = strFind(tr(:)%nam1, tr(jq)%name)                          !--- Indexes of "tr(jq)" childs in "tr(:)"
+      ng = MAXVAL(tr(iy)%igen, MASK=.TRUE., DIM=1)                   !--- Generations number of the "tr(jq)" family
+      DO ig = 2, ng                                                  !--- Loop   on generations for the tr(jq) family
+        iz = find(tr(iy)%igen, ig, n)                                !--- Indexes of the tracers "tr(iy(:))" of generation "ig"
+        ix(iq:iq+n-1) = iy(iz)                                       !--- Same indexes in "tr(:)"
+        iq = iq + n
+      END DO
+    END DO
+  END IF
+  tr = tr(ix)                                                        !--- Reorder the tracers
+END SUBROUTINE sortTracers
+!==============================================================================================================================
+
+!==============================================================================================================================
+LOGICAL FUNCTION mergeTracers(sections, tr) RESULT(lerr)
+  TYPE(db),  TARGET,      INTENT(IN)  :: sections(:)
+  TYPE(tra), ALLOCATABLE, INTENT(OUT) ::       tr(:)
+  TYPE(tra), POINTER     ::   t1(:),   t2(:)
+  INTEGER,  ALLOCATABLE :: ixct(:), ixck(:)
+  INTEGER :: is, k1, k2, nk2, i1, i2, nt2
+  CHARACTER(LEN=256) :: s1, v1, v2, tnam, knam
+  lerr = .FALSE.
+  t1 => sections(1)%trac(:)                                          !--- Alias: first tracers section
+  tr = t1
+  !----------------------------------------------------------------------------------------------------------------------------
+  DO is=2,SIZE(sections)                                             !=== SEVERAL SECTIONS: MERGE THEM
+  !----------------------------------------------------------------------------------------------------------------------------
+    t2  => sections(is)%trac(:)                                      !--- Alias: current tracers section
+    nt2  = SIZE(t2(:), DIM=1)                                        !--- Number of tracers in section
+    ixct = strIdx(t1(:)%name, t2(:)%name)                            !--- Indexes of common tracers
+    tr = [tr, PACK(t2, MASK= ixct==0)]                               !--- Append with new tracers
+    IF( ALL(ixct == 0) ) CYCLE                                       !--- No common tracers => done
+    CALL msg('Tracers defined in previous sections and duplicated in "'//TRIM(sections(is)%name)//'":')
+    CALL msg( t1(PACK(ixct, MASK = ixct/=0))%name, nmax=128 )        !--- Display duplicates (the 128 first at most)
+    !--------------------------------------------------------------------------------------------------------------------------
+    DO i2=1,nt2; tnam = t2(i2)%name                                  !=== LOOP ON COMMON TRACERS
+    !--------------------------------------------------------------------------------------------------------------------------
+      i1 = ixct(i2); IF(i1 == 0) CYCLE                               !--- Idx in t1(:) ; skip new tracers
+
+      !=== CHECK WETHER ESSENTIAL KEYS ARE IDENTICAL OR NOT
+      s1=' of "'//TRIM(tnam)//'" in "'//TRIM(sections(is)%name)//'" not matching previous value'
+      
+      IF(test(fmsg(t1(i1)%prnt /= t2(i2)%prnt, 'Parent name'//TRIM(s1)), lerr)) RETURN
+      IF(test(fmsg(t1(i1)%type /= t2(i2)%type, 'Type'       //TRIM(s1)), lerr)) RETURN
+      IF(test(fmsg(t1(i1)%igen /= t2(i2)%igen, 'Generation' //TRIM(s1)), lerr)) RETURN
+
+      !=== APPEND <key>=<val> PAIRS NOT PREVIOULSLY DEFINED
+      nk2  = SIZE(t2(i2)%keys%key(:))                                !--- Keys number in current section
+      ixck = strIdx(t1(i1)%keys%key(:), t2(i2)%keys%key(:))          !--- Common keys indexes
+
+      !=== APPEND NEW KEYS
+      tr(i1)%keys%key = [ tr(i1)%keys%key, PACK(tr(i2)%keys%key, MASK = ixck==0)]
+      tr(i1)%keys%val = [ tr(i1)%keys%val, PACK(tr(i2)%keys%val, MASK = ixck==0)]
+
+      !--- KEEP TRACK OF THE COMPONENTS NAMES
+      tr(i1)%comp = TRIM(tr(i1)%comp)//','//TRIM(tr(i2)%comp)
+
+      !--- SELECT COMMON TRACERS WITH DIFFERING KEYS VALUES (PREVIOUS VALUE IS KEPT)
+      DO k2=1,nk2
+        k1 = ixck(k2); IF(k1 == 0) CYCLE
+        IF(t1(i1)%keys%val(k1) == t2(i2)%keys%val(k2)) ixck(k2)=0
+      END DO
+      IF(ALL(ixck==0)) CYCLE                                         !--- No identical keys with /=values
+
+      !--- DISPLAY INFORMATION: OLD VALUES ARE KEPT FOR THE KEYS FOUND IN PREVIOUS AND CURRENT SECTIONS
+      CALL msg('Key(s)'//TRIM(s1))
+      DO k2 = 1, nk2                                                 !--- Loop on keys found in both t1(:) and t2(:)
+        knam = t2(i2)%keys%key(k2)                                   !--- Name of the current key
+        k1 = ixck(k2)                                                !--- Corresponding index in t1(:)
+        IF(k1 == 0) CYCLE                                            !--- New keys are skipped
+        v1 = t1(i1)%keys%val(k1); v2 = t2(i2)%keys%val(k2)           !--- Key values in t1(:) and t2(:)
+        CALL msg(' * '//TRIM(knam)//'='//TRIM(v2)//' ; previous value kept:'//TRIM(v1))
+      END DO
+      !------------------------------------------------------------------------------------------------------------------------
+    END DO
+    !--------------------------------------------------------------------------------------------------------------------------
+  END DO
+  CALL sortTracers(tr)
+
+END FUNCTION mergeTracers
+!==============================================================================================================================
+
+!==============================================================================================================================
+LOGICAL FUNCTION cumulTracers(sections, tr) RESULT(lerr)
+  TYPE(db),  TARGET,      INTENT(IN)  :: sections(:)
+  TYPE(tra), ALLOCATABLE, INTENT(OUT) ::       tr(:)
+  TYPE(tra), POINTER     :: t1(:), t2(:)
+  INTEGER,   ALLOCATABLE :: nt(:)
+  CHARACTER(LEN=256)     :: tnam, tnam_new
+  INTEGER :: iq, nq, is, ns, nsec
+  lerr = .FALSE.                                                     !--- Can't fail ; kept to match "mergeTracer" interface.
+  nsec =  SIZE(sections)
+  tr = [(      sections(is)%trac(:) , is=1, nsec )]                  !--- Concatenated tracers vector
+  nt = [( SIZE(sections(is)%trac(:)), is=1, nsec )]                  !--- Number of tracers in each section
+  !----------------------------------------------------------------------------------------------------------------------------
+  DO is=1, nsec                                                      !=== LOOP ON SECTIONS
+  !----------------------------------------------------------------------------------------------------------------------------
+    t1 => sections(is)%trac(:)
+    !--------------------------------------------------------------------------------------------------------------------------
+    DO iq=1, nt(is)                                                  !=== LOOP ON TRACERS
+    !--------------------------------------------------------------------------------------------------------------------------
+      tnam = t1(iq)%name                                             !--- Original name
+      IF(COUNT(t1%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
+      ns = nt(is)                                                    !--- Number of tracers in the current section
+      tr(iq + nq)%name = tnam_new                                    !--- Modify tracer name
+      WHERE(tr(1+nq:ns+nq)%prnt==tnam) tr(1+nq:ns+nq)%prnt=tnam_new  !--- Modify parent name
+    !--------------------------------------------------------------------------------------------------------------------------
+    END DO
+  !----------------------------------------------------------------------------------------------------------------------------
+  END DO
+  !----------------------------------------------------------------------------------------------------------------------------
+  CALL sortTracers(tr)
+END FUNCTION cumulTracers
+!==============================================================================================================================
+
+!==============================================================================================================================
+SUBROUTINE setDirectKeys(tr)
+  TYPE(tra), INTENT(INOUT) :: tr(:)
+  CALL indexUpdate(tr)                                               !--- Update iparnt and idesc indexes vectors
+!  DO iq = 1, SIZE(tr)
+!    tr(iq)%keys%<key> = getKey_prv(it, "<key>", tr%keys,  tran0 )   !--- For additional keys
+!  END DO
+END SUBROUTINE setDirectKeys
+!==============================================================================================================================
+
+!==============================================================================================================================
+LOGICAL FUNCTION dispTraSection(message, sname) RESULT(lerr)
+  CHARACTER(LEN=*), INTENT(IN) :: message
+  CHARACTER(LEN=*), INTENT(IN) :: sname
+  INTEGER :: idb, iq, nq
+  INTEGER, ALLOCATABLE :: hadv(:), vadv(:)
+  TYPE(tra), POINTER :: tm(:)
+  lerr = .FALSE.
+  idb = strIdx(dBase(:)%name, sname); IF(idb == 0) RETURN
+  tm => dBase(idb)%trac
+  nq = SIZE(tm)
+  IF(test(getKeyByName_im('hadv', hadv, tm(:)%name, tm(:)%keys),lerr)) RETURN
+  IF(test(getKeyByName_im('vadv', vadv, tm(:)%name, tm(:)%keys),lerr)) RETURN
+  CALL msg(TRIM(message)//':')
+  IF(test(dispTable('iiissis', ['iq        ','hadv      ','vadv      ','short name','parent    ','igen      ','phase     '], &
+                    cat(tm(:)%name,  tm(:)%prnt, tm(:)%phas), cat([(iq, iq=1, nq)],  hadv,  vadv, tm(:)%igen)), lerr)) RETURN
+END FUNCTION dispTraSection
+!==============================================================================================================================
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!== CREATE A SCALAR ALIAS OF THE COMPONENT OF THE TRACERS DESCRIPTOR "t" NAMED "tname" ========================================
+!==============================================================================================================================
+FUNCTION aliasTracer(tname, t) RESULT(out)
+  TYPE(tra),         POINTER    :: out
+  CHARACTER(LEN=*),  INTENT(IN) :: tname
+  TYPE(tra), TARGET, INTENT(IN) :: t(:)
+  INTEGER :: it
+  it = strIdx(t(:)%name, tname)
+  out => NULL(); IF(it /= 0) out => t(it)
+END FUNCTION aliasTracer
+!------------------------------------------------------------------------------------------------------------------------------
+
+
+!==============================================================================================================================
+!=== FROM A LIST OF INDEXES OR NAMES, CREATE A SUBSET OF THE TRACERS DESCRIPTORS LIST "trac" ==================================
+!==============================================================================================================================
+FUNCTION trSubset_Indx(trac,idx) RESULT(out)
+  TYPE(tra), ALLOCATABLE             ::  out(:)
+  TYPE(tra), ALLOCATABLE, INTENT(IN) :: trac(:)
+  INTEGER,                INTENT(IN) ::  idx(:)
+  out = trac(idx)
+  CALL indexUpdate(out)
+END FUNCTION trSubset_Indx
+!------------------------------------------------------------------------------------------------------------------------------
+FUNCTION trSubset_Name(trac,nam) RESULT(out)
+  TYPE(tra), ALLOCATABLE             ::  out(:)
+  TYPE(tra), ALLOCATABLE, INTENT(IN) :: trac(:)
+  CHARACTER(LEN=*),       INTENT(IN) ::  nam(:)
+  out = trac(strIdx(trac(:)%name, nam))
+  CALL indexUpdate(out)
+END FUNCTION trSubset_Name
+!------------------------------------------------------------------------------------------------------------------------------
+
+
+!==============================================================================================================================
+!=== CREATE THE SUBSET OF THE TRACERS DESCRIPTORS LIST "trac" HAVING THE FIRST GENERATION ANCESTOR NAMED "nam" ================
+!==============================================================================================================================
+FUNCTION trSubset_Nam1(trac,nam) RESULT(out)
+  TYPE(tra), ALLOCATABLE             ::  out(:)
+  TYPE(tra), ALLOCATABLE, INTENT(IN) :: trac(:)
+  CHARACTER(LEN=*),       INTENT(IN) ::  nam
+  out = trac(strFind(delPhase(trac(:)%nam1), nam))
+  CALL indexUpdate(out)
+END FUNCTION trSubset_Nam1
+!------------------------------------------------------------------------------------------------------------------------------
+
+
+!==============================================================================================================================
+!=== UPDATE THE INDEXES iparnt(:), idesc=(:) AND igen(:) IN THE TRACERS DESCRIPTOR LIST "tr" (USEFULL FOR SUBSETS) ============
+!==============================================================================================================================
+SUBROUTINE indexUpdate(tr)
+  TYPE(tra), INTENT(INOUT) :: tr(:)
+  INTEGER :: iq, ig, ng, ngen
+  INTEGER, ALLOCATABLE :: ix(:)
+  tr(:)%iprnt = strIdx( tr(:)%name, tr(:)%prnt )                     !--- Parent index
+  ngen = MAXVAL(tr(:)%igen, MASK=.TRUE.)
+  DO iq = 1, SIZE(tr)
+    ng = tr(iq)%igen                                                 !--- Generation of the current tracer
+    ix = idxAncestor(tr, igen = ng); ix = PACK(ix, ix/=0)            !--- Indexes of the tracers with ancestor tr(iq)
+    !--- Childs indexes in growing generation order
+    tr(iq)%idesc = [( PACK(ix, MASK = tr(ix)%igen == ig), ig = ng+1, ngen)]
+    tr(iq)%ndesc =     SUM(  [( COUNT(tr(ix)%igen == ig), ig = ng+1, ngen)] )
+    tr(iq)%nchld =              COUNT(tr(ix)%igen == ng+1)
+  END DO
+END SUBROUTINE indexUpdate
+!------------------------------------------------------------------------------------------------------------------------------
+  
+  
+!==============================================================================================================================
+!=== READ THE ISOTOPES NAMED "iso" FROM THE TRACERS SECTIONS "tr" IN THE FILE "fnam" ; PUT RESULT IN A TRACERS DESCRIPTOR ====
+!===  * SYNTAX IS THE SAME AS IN THE "tracer.def" FILE ; EACH TRACER SECTION CONTAINS ONE LINE EACH OF ITS KNOWN ISOTOPES  ====
+!===  * EACH TRACERS SECTION CAN CONTAIN A "params" VIRTUAL ISOTOPE LINE CONTAINING DEFAULT PARAMETERS FOR THE ISOTOPES    ====
+!===  * IF SOME KEYS ARE FOUND BOTH IN THE "*.def" FILES AND THE "params" SECTION, TEH VALUE FROM "*.def" FILE IS RETAINED ====
+!===  * ON EACH ISOTOPE LINE, A DEFINED KEY CAN BE USED IN THE OTHER KEYS AS A PARAMETER (SIGNLE LEVEL DEPENDENCY !)       ====
+!===  * THE DIFFERENT ISOTOPES SETS (ONE EACH PARENT TRACER) ARE MERGED INTO A SINGLE TRACERS DESCRIPTOR VECTOR            ====
+!===  * 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)
+  CHARACTER(LEN=*),  INTENT(IN)    :: fnam                           !--- Input file name
+  TYPE(iso), TARGET, INTENT(INOUT) :: isot(:)                        !--- Isotopes descriptors (field "prnt" must be defined !)
+  INTEGER :: ik, is, it, idb, nk0, i
+  INTEGER :: nk, ns, nt, ndb, nb0, i0
+  CHARACTER(LEN=256), POINTER     :: k(:), v(:), k0(:), v0(:)
+  CHARACTER(LEN=256), ALLOCATABLE :: vals(:)
+  CHARACTER(LEN=256)     :: val
+  TYPE(kys),    POINTER  ::   ky(:)
+  TYPE(tra),    POINTER  ::    t(:)
+  TYPE(tra), ALLOCATABLE ::   tt(:)
+  TYPE(db),  ALLOCATABLE ::  tdb(:)
+  LOGICAL,   ALLOCATABLE :: liso(:)
+
+  !--- THE INPUT FILE MUST BE PRESENT
+  IF(test(fmsg(testFile(fnam),'Missing isotopes parameters file "'//TRIM(fnam)//'"'),lerr)) RETURN
+
+  !--- READ THE FILE SECTIONS, ONE EACH PARENT TRACER
+  nb0 = SIZE(dBase, DIM=1)+1                                         !--- Next database element index
+  IF(test(readSections(fnam,strStack(isot(:)%prnt,',')),lerr)) RETURN!--- Read sections, one each parent tracer
+  ndb = SIZE(dBase, DIM=1)                                           !--- Current database size
+
+  DO idb = nb0, ndb
+    t => dBase(idb)%trac(:)
+    nt = SIZE(t)                                                     !--- Number of isotopes in the current database section
+
+PRINT*
+PRINT*,'AVANT:'
+DO it=1,SIZE(t); print*,TRIM(t(it)%name)//':  '//strStack([(TRIM(t(it)%keys%key(i))//'='//TRIM(t(it)%keys%val(i)), i=1, SIZE(t(it)%keys%key))]); END DO
+    !--- GET FEW GLOBAL KEYS FROM "def" FILES AND ADD THEM TO THE 'params' SECTION
+    CALL addKeysFromDef(dBase(idb)%trac, 'params')
+
+    !--- SUBSTITUTE THE KEYS DEFINED IN THE 'params' VIRTUAL TRACER ; SUBSTITUTE LOCAL KEYS ; REMOVE 'params' VIRTUAL TRACER
+    CALL subDefault(dBase(idb)%trac, 'params', .TRUE.)
+
+PRINT*
+PRINT*,'AVANT REDUCTION:'
+    t => dBase(idb)%trac(:)
+    DO it=1,SIZE(t); print*,TRIM(t(it)%name)//':  '//strStack([(TRIM(t(it)%keys%key(i))//'='//TRIM(t(it)%keys%val(i)), i=1, SIZE(t(it)%keys%key))]); END DO
+
+    !--- REDUCE THE EXPRESSIONS TO OBTAIN SCALARS
+    DO it=1, nt
+      v => dBase(idb)%trac(it)%keys%val(:)
+      WHERE(reduceExpr(v, vals)) v = vals
+    END DO
+
+PRINT*
+PRINT*,'APRES:'
+    t => dBase(idb)%trac(:)
+    DO it=1,SIZE(t); print*,TRIM(t(it)%name)//':  '//strStack([(TRIM(t(it)%keys%key(i))//'='//TRIM(t(it)%keys%val(i)), i=1, SIZE(t(it)%keys%key))]); END DO
+
+    !--- TRANSFER THE key=val PAIRS TO THE ISOTOPES DESCRIPTOR
+print*
+print*,'isot%prnt = '//strStack(isot%prnt)
+
+    ky => isot(strIdx(isot(:)%prnt, dBase(idb)%name))%keys           !--- Keys of "isot" tracers with parent "dBase(idb)%name"
+print*,'ky%name = '//strStack(ky%name)
+    is=1
+    DO it = 1, nt; IF(it == i0) CYCLE
+print*,'AAAAAA '//strStack(ky%name)
+print*,'       '//TRIM(t(it)%name)
+      is = strIdx(ky(:)%name, t(it)%name)                            !--- Index of the "isot(:)" tracer named "t(it)%name"
+      IF(is == 0) CYCLE                                              !--- Current isotope is not present in "isot" => not needed
+      k => ky(is)%key; k = t(it)%keys%key
+      v => ky(is)%val; v = t(it)%keys%val
+      WHERE(reduceExpr(v, vals)) v = vals
+      DO ik=1, SIZE(k); IF(reduceExpr(v(ik),val)) v(ik) = val; END DO!--- Reduce operations (for substituted variables)
+print*,'(4) '//strStack([(TRIM(k(i))//'='//TRIM(v(i)), i=1, SIZE(k))])
+    END DO
+print*,'(7) i0=',i0
+
+    !--- CHECK FOR MISSING ISOTOPES (NO KEYS ALLOCATED)
+    liso = [( ALLOCATED(ky(is)%key), is=1, SIZE(ky) )]
+print*,'liso=',liso
+    IF(test(checkList(ky(:)%name, &
+      .NOT.liso, 'Check file "'//TRIM(fnam)//'" in section "'//TRIM(dBase(idb)%name)//'"', 'isotopes', 'missing'),lerr)) RETURN
+  END DO
+
+  !--- CLEAN THE DATABASE ENTRIES
+  IF(nb0 == 1) THEN
+    DEALLOCATE(dBase); ALLOCATE(dBase(0))
+  ELSE
+    ALLOCATE(tdb(nb0-1)); tdb(1:nb0-1)=dBase(1:nb0-1); CALL MOVE_ALLOC(FROM=tdb, TO=dBase)
+  END IF
+  lerr = dispIsotopes(isot,'isotopes parameters read from file')
+
+END FUNCTION readIsotopesFile
+!==============================================================================================================================
+
+!==============================================================================================================================
+LOGICAL FUNCTION dispIsotopes(ides, message) RESULT(lerr)
+  TYPE(iso),        INTENT(IN) :: ides(:)                            !--- Isotopes descriptor vector
+  CHARACTER(LEN=*), INTENT(IN) :: message                            !--- Message to display
+  INTEGER :: ik, nk, ip, it, nt
+  CHARACTER(LEN=256) :: prf
+  CHARACTER(LEN=256), ALLOCATABLE :: ttl(:), val(:,:)
+  CALL msg(TRIM(message)//':')
+  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
+    END DO
+    IF(test(fmsg(dispTable(prf, ttl, val, cat([(it,it=1,nt)]), rFmt='(EN8.4)'),'Problem with the table content'), lerr)) RETURN
+    DEALLOCATE(ttl, val)
+  END DO        
+END FUNCTION dispIsotopes
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+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(kys),         INTENT(INOUT) :: ky
+  LOGICAL, OPTIONAL, INTENT(IN)    :: lOverWrite
+  CHARACTER(LEN=256),  ALLOCATABLE :: k(:), v(:)
+  INTEGER :: iky, nky
+  LOGICAL :: lo
+!------------------------------------------------------------------------------------------------------------------------------
+  lo=.FALSE.; IF(PRESENT(lOverWrite)) lo=lOverWrite
+  iky = strIdx(ky%key,key)
+  IF(iky == 0) THEN
+    nky = SIZE(ky%key)
+    IF(nky == 0) THEN
+      ky%key = TRIM(key); ky%val = TRIM(val)
+    ELSE
+      ALLOCATE(k(nky+1)); k(1:nky) = ky%key; k(nky+1) = TRIM(key); CALL MOVE_ALLOC(FROM=k, TO=ky%key)
+      ALLOCATE(v(nky+1)); v(1:nky) = ky%val; v(nky+1) = TRIM(val); CALL MOVE_ALLOC(FROM=v, TO=ky%val)
+    END IF
+  ELSE IF(lo) THEN                                                   !--- Overwriting
+    ky%key(iky) = TRIM(key); ky%val(iky) = TRIM(val)
+  END IF
+END SUBROUTINE addKey_1
+!==============================================================================================================================
+SUBROUTINE addKey_tra(key, val, tr, lOverWrite, tname)
+!------------------------------------------------------------------------------------------------------------------------------
+! Purpose: Add the <key>=<val> pair in all the components of the "tr(itr)%keys" keys descriptor:
+!          * "tname"   specified: for the index "itr" of the tracer named "tname"
+!          * "tname" unspecified: for all the tracers
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=*),           INTENT(IN)    :: key, val
+  TYPE(tra),                  INTENT(INOUT) :: tr(:)
+  LOGICAL,          OPTIONAL, INTENT(IN)    :: lOverWrite
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)    :: tname
+  INTEGER :: itr
+  LOGICAL :: lo
+!------------------------------------------------------------------------------------------------------------------------------
+  lo=.FALSE.; IF(PRESENT(lOverWrite)) lo=lOverWrite
+  IF(PRESENT(tname)) THEN
+    itr = strIdx(tr%name, tname)
+    IF(itr == 0) RETURN
+    CALL addKey_1(key, val, tr(itr)%keys, lo)
+  ELSE
+    DO itr = 1, SIZE(tr); CALL addKey_1(key, val, tr(itr)%keys, lo); END DO
+  END IF
+END SUBROUTINE addKey_tra
+!==============================================================================================================================
+SUBROUTINE addKeysFromDef(tr, tr0)
+  USE ioipsl_getin_p_mod, ONLY : getin_p
+  TYPE(tra), ALLOCATABLE, INTENT(INOUT) :: tr(:)
+  CHARACTER(LEN=*),       INTENT(IN)    :: tr0
+  CHARACTER(LEN=256) :: val
+  INTEGER            :: ik, i0
+  i0 = strIdx(tr%name, tr0)
+  IF(i0 == 0) RETURN
+  DO ik = 1, SIZE(tr(i0)%keys%key)
+    val   =   'zzzz'; CALL  getin_p(tr(i0)%keys%key(ik), val)
+    IF(val /= 'zzzz') CALL addKey_1(tr(i0)%keys%key(ik), val, tr(i0)%keys, .TRUE.)
+  END DO
+END SUBROUTINE addKeysFromDef
+!==============================================================================================================================
+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(tra),        INTENT(INOUT) :: ky(:)
+  CHARACTER(LEN=256), 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) )]
+  k = PACK(ky(itr)%keys%key, MASK=ll); CALL MOVE_ALLOC(FROM=k, TO=ky(itr)%keys%key)
+  v = PACK(ky(itr)%keys%val, MASK=ll); CALL MOVE_ALLOC(FROM=v, TO=ky(itr)%keys%val)
+END SUBROUTINE delKey_1
+!==============================================================================================================================
+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(tra),        INTENT(INOUT) :: ky(:)
+  INTEGER :: iky
+!------------------------------------------------------------------------------------------------------------------------------
+  DO iky = 1, SIZE(ky); CALL delKey_1(iky, keyn, ky); END DO
+END SUBROUTINE delKey
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== 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(tra), OPTIONAL, INTENT(IN) ::  tracers_(:)
+  TYPE(iso), OPTIONAL, INTENT(IN) :: isotopes_(:)
+  IF(PRESENT( tracers_))  tracers =  tracers_
+  IF(PRESENT(isotopes_)) isotopes = isotopes_
+END SUBROUTINE getKey_init
+!==============================================================================================================================
+CHARACTER(LEN=256) FUNCTION fgetKey(itr, keyn, ky, def_val) RESULT(out)
+!------------------------------------------------------------------------------------------------------------------------------
+! Purpose: Internal function ; get a key value in string format (this is the returned argument).
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER,                    INTENT(IN) :: itr
+  CHARACTER(LEN=*),           INTENT(IN) :: keyn
+  TYPE(kys),                  INTENT(IN) :: ky(:)
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: def_val
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: ik
+  ik = 0; IF(itr>0 .AND. itr<=SIZE(ky)) ik = strIdx(ky(itr)%key(:), keyn)
+  out = '';              IF(ik /= 0) out = ky(itr)%val(ik)           !--- Key was found
+  IF(PRESENT(def_val) .AND. ik == 0) out = def_val                   !--- Default value from arguments
+END FUNCTION fgetKey
+!==============================================================================================================================
+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 suffix, then in "isotopes" without.
+  !     * "ky"   specified: try in "ky"      for "tnam" with, then without phase suffix.
+  !    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=256),  INTENT(OUT) :: val
+  CHARACTER(LEN=*),    INTENT(IN)  :: tname
+  TYPE(kys), OPTIONAL, INTENT(IN)  :: ky(:)
+  INTEGER :: is
+  lerr = .FALSE.
+  IF(PRESENT(ky)) THEN
+    val = getKeyByName_prv(keyn,          tname , ky);    IF(val /= '') RETURN !--- "ky" and "tnam"
+    val = getKeyByName_prv(keyn, delPhase(tname), ky)                          !--- "ky" and "tnam" without phase
+  ELSE
+    IF(.NOT.ALLOCATED(tracers))  RETURN
+    val = getKeyByName_prv(keyn, tname, tracers(:)%keys); IF(val /= '') RETURN !--- "tracers" and "tnam"
+    IF(.NOT.ALLOCATED(isotopes)) RETURN
+    IF(SIZE(isotopes) == 0)      RETURN
+    DO is = 1, SIZE(isotopes); IF(strIdx(isotopes(is)%keys(:)%name, tname) /= 0) EXIT; END DO
+    IF(is /= 0) val = getKeyByName_prv(keyn, tname, isotopes(is)%keys(:))      !--- "isotopes" and "tnam" without phase
+  END IF
+
+CONTAINS
+
+FUNCTION getKeyByName_prv(keyn, tname, ky) RESULT(val)
+  CHARACTER(LEN=256)            :: val
+  CHARACTER(LEN=*), INTENT(IN)  :: keyn
+  CHARACTER(LEN=*), INTENT(IN)  :: tname
+  TYPE(kys),        INTENT(IN)  :: ky(:)
+  INTEGER :: itr, iky
+  val = ''; iky = 0
+  itr = strIdx(ky(:)%name, tname);                 IF(itr==0) RETURN           !--- Get the index of the wanted tracer
+  IF(itr /= 0) iky = strIdx(ky(itr)%key(:), keyn); IF(iky==0) RETURN           !--- Wanted key    index
+  val = ky(itr)%val(iky)
+END FUNCTION getKeyByName_prv
+
+END FUNCTION getKeyByName_s1
+!==============================================================================================================================
+LOGICAL FUNCTION getKeyByName_sm(keyn, val, tnam, ky) RESULT(lerr)
+  CHARACTER(LEN=*),                   INTENT(IN)  :: keyn
+  CHARACTER(LEN=256),    ALLOCATABLE, INTENT(OUT) ::  val(:)
+  CHARACTER(LEN=*), TARGET, OPTIONAL, INTENT(IN)  :: tnam(:)
+  TYPE(kys),        TARGET, OPTIONAL, INTENT(IN)  ::   ky(:)
+  CHARACTER(LEN=256), POINTER :: n(:)
+  INTEGER :: iq
+  n => tracers(:)%keys%name; IF(PRESENT(tnam)) n => tnam(:)
+  ALLOCATE(val(SIZE(n)))
+  IF(     PRESENT(ky)) lerr = ANY([(getKeyByName_s1(keyn, val(iq), n(iq), ky), iq=1, SIZE(n))])
+  IF(.NOT.PRESENT(ky)) lerr = ANY([(getKeyByName_s1(keyn, val(iq), n(iq)),     iq=1, SIZE(n))])
+END FUNCTION getKeyByName_sm
+!==============================================================================================================================
+LOGICAL FUNCTION getKeyByName_i1(keyn, val, tnam, ky) RESULT(lerr)
+  CHARACTER(LEN=*),    INTENT(IN)  :: keyn
+  INTEGER,             INTENT(OUT) :: val
+  CHARACTER(LEN=*),    INTENT(IN)  :: tnam
+  TYPE(kys), OPTIONAL, INTENT(IN)  :: ky(:)
+  CHARACTER(LEN=256) :: sval
+  INTEGER :: ierr
+  IF(     PRESENT(ky)) lerr = getKeyByName_s1(keyn, sval, tnam, ky)
+  IF(.NOT.PRESENT(ky)) lerr = getKeyByName_s1(keyn, sval, tnam)
+  IF(test(fmsg(lerr,   'key "'//TRIM(keyn)//'" or tracer "'//TRIM(tnam)//'" is missing'),        lerr)) RETURN
+  READ(sval, *, IOSTAT=ierr) val
+  IF(test(fmsg(ierr/=0,'key "'//TRIM(keyn)//'" of tracer "'//TRIM(tnam)//'" is not an integer'), lerr)) RETURN
+END FUNCTION getKeyByName_i1
+!==============================================================================================================================
+LOGICAL FUNCTION getKeyByName_im(keyn, val, tnam, ky) RESULT(lerr)
+  CHARACTER(LEN=*),                   INTENT(IN)  :: keyn
+  INTEGER,               ALLOCATABLE, INTENT(OUT) ::  val(:)
+  CHARACTER(LEN=*), TARGET, OPTIONAL, INTENT(IN)  :: tnam(:)
+  TYPE(kys),        TARGET, OPTIONAL, INTENT(IN)  ::   ky(:)
+  CHARACTER(LEN=256), POINTER :: n(:)
+  INTEGER :: iq
+  n => tracers(:)%name; IF(PRESENT(tnam)) n => tnam(:)
+  ALLOCATE(val(SIZE(n)))
+  IF(     PRESENT(ky)) lerr = ANY([(getKeyByName_i1(keyn, val(iq), n(iq), ky), iq=1, SIZE(n))])
+  IF(.NOT.PRESENT(ky)) lerr = ANY([(getKeyByName_i1(keyn, val(iq), n(iq)),     iq=1, SIZE(n))])
+END FUNCTION getKeyByName_im
+!==============================================================================================================================
+LOGICAL FUNCTION getKeyByName_r1(keyn, val, tnam, ky) RESULT(lerr)
+  CHARACTER(LEN=*),    INTENT(IN)  :: keyn
+  REAL,                INTENT(OUT) :: val
+  CHARACTER(LEN=*),    INTENT(IN)  :: tnam
+  TYPE(kys), OPTIONAL, INTENT(IN)  :: ky(:)
+  CHARACTER(LEN=256) :: sval
+  INTEGER :: ierr
+  IF(     PRESENT(ky)) lerr = getKeyByName_s1(keyn, sval, tnam, ky)
+  IF(.NOT.PRESENT(ky)) lerr = getKeyByName_s1(keyn, sval, tnam)
+  IF(test(fmsg(lerr,   'key "'//TRIM(keyn)//'" or tracer "'//TRIM(tnam)//'" is missing'),    lerr)) RETURN
+  READ(sval, *, IOSTAT=ierr) val
+  IF(test(fmsg(ierr/=0,'key "'//TRIM(keyn)//'" of tracer "'//TRIM(tnam)//'" is not a real'), lerr)) RETURN
+END FUNCTION getKeyByName_r1
+!==============================================================================================================================
+LOGICAL FUNCTION getKeyByName_rm(keyn, val, tnam, ky) RESULT(lerr)
+  CHARACTER(LEN=*),                   INTENT(IN)  :: keyn
+  REAL,                  ALLOCATABLE, INTENT(OUT) ::  val(:)
+  CHARACTER(LEN=*), TARGET, OPTIONAL, INTENT(IN)  :: tnam(:)
+  TYPE(kys),        TARGET, OPTIONAL, INTENT(IN)  ::   ky(:)
+  CHARACTER(LEN=256), POINTER :: n(:)
+  INTEGER :: iq
+  n => tracers(:)%name;  IF(PRESENT(tnam)) n => tnam(:)
+  ALLOCATE(val(SIZE(n)))
+  IF(     PRESENT(ky)) lerr = ANY([(getKeyByName_r1(keyn, val(iq), n(iq), ky), iq=1, SIZE(n))])
+  IF(.NOT.PRESENT(ky)) lerr = ANY([(getKeyByName_r1(keyn, val(iq), n(iq)),     iq=1, SIZE(n))])
+END FUNCTION getKeyByName_rm
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== REMOVE, IF ANY, OR ADD THE PHASE SUFFIX OF THE TRACER NAMED "s" ==========================================================
+!==============================================================================================================================
+ELEMENTAL CHARACTER(LEN=256) FUNCTION delPhase(s) RESULT(out)
+  CHARACTER(LEN=*), INTENT(IN) :: s
+  INTEGER :: l
+  out = s
+  IF(s == '') RETURN
+  l=LEN_TRIM(s)
+  IF(s(l-1:l-1)=='-' .AND. INDEX(known_phases,s(l:l))/=0) out = s(1:l-2)
+END FUNCTION delPhase
+!------------------------------------------------------------------------------------------------------------------------------
+CHARACTER(LEN=256) FUNCTION addPhase_1(s,pha) RESULT(out)
+  CHARACTER(LEN=*), INTENT(IN) :: s
+  CHARACTER(LEN=1), INTENT(IN) :: pha
+  IF(INDEX(s,'_')==0) THEN; out = TRIM(s)//'-'//pha; RETURN; END IF
+  out = TRIM(strHead(s,'_'))//'-'//pha//TRIM(strTail(s,'_'))
+END FUNCTION addPhase_1
+!------------------------------------------------------------------------------------------------------------------------------
+FUNCTION addPhase_m(s,pha) RESULT(out)
+  CHARACTER(LEN=*),    INTENT(IN) :: s(:)
+  CHARACTER(LEN=1),    INTENT(IN) :: pha
+  CHARACTER(LEN=256), ALLOCATABLE :: out(:)  
+  INTEGER :: k
+  out = [( addPhase_1(s(k), pha), k=1, SIZE(s) )]
+END FUNCTION addPhase_m
+!------------------------------------------------------------------------------------------------------------------------------
+
+
+!==============================================================================================================================
+!=== GET THE NAME(S) OF THE ANCESTOR(S) OF TRACER(S) "tname" AT GENERATION "igen"  IN THE TRACERS DESCRIPTORS LIST "tr" =======
+!==============================================================================================================================
+CHARACTER(LEN=256) FUNCTION ancestor_1(t, tname, igen) RESULT(out)
+  TYPE(tra),         INTENT(IN) :: t(:)
+  CHARACTER(LEN=*),  INTENT(IN) :: tname
+  INTEGER, OPTIONAL, INTENT(IN) :: igen
+  INTEGER :: ig, ix
+  ig = 1; IF(PRESENT(igen)) ig = igen
+  ix = idxAncestor_1(t, tname, ig)
+  out = ''; IF(ix /= 0) out = t(ix)%name
+END FUNCTION ancestor_1
+!------------------------------------------------------------------------------------------------------------------------------
+FUNCTION ancestor_m(t, tname, igen) RESULT(out)
+  CHARACTER(LEN=256), ALLOCATABLE        ::   out(:)
+  TYPE(tra),                  INTENT(IN) ::     t(:)
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: tname(:)
+  INTEGER,          OPTIONAL, INTENT(IN) :: igen
+  INTEGER, ALLOCATABLE :: ix(:)
+  INTEGER :: ig
+  ig = 1; 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" ancestor of "tname"
+  TYPE(tra),         INTENT(IN) :: t(:)
+  CHARACTER(LEN=*),  INTENT(IN) :: tname
+  INTEGER, OPTIONAL, INTENT(IN) :: igen
+  INTEGER :: ig
+  ig = 1; IF(PRESENT(igen)) ig = igen
+  out = strIdx(t(:)%name, tname)
+  IF(out == 0)          RETURN
+  IF(t(out)%igen <= ig) RETURN
+  DO WHILE(t(out)%igen > ig); out = strIdx(t(:)%name, t(out)%prnt); END DO
+END FUNCTION idxAncestor_1
+!------------------------------------------------------------------------------------------------------------------------------
+FUNCTION idxAncestor_m(t, tname, igen) RESULT(out)
+  INTEGER,          ALLOCATABLE          ::   out(:)
+  TYPE(tra),                  INTENT(IN) ::     t(:)
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: tname(:)
+  INTEGER,          OPTIONAL, INTENT(IN) :: igen
+  INTEGER :: ig, ix
+  ig = 1; 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
+!==============================================================================================================================
+
+
+END MODULE readTracFiles_mod
Index: LMDZ6/branches/LMDZ-tracers/libf/misc/strings_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ-tracers/libf/misc/strings_mod.F90	(revision 3852)
+++ LMDZ6/branches/LMDZ-tracers/libf/misc/strings_mod.F90	(revision 3852)
@@ -0,0 +1,1343 @@
+MODULE strings_mod
+
+  IMPLICIT NONE
+
+  PRIVATE
+  PUBLIC :: modname, init_printout, msg, fmsg, get_in, lunout, prt_level
+  PUBLIC :: strLower, strHead, strStack,  strClean,  strIdx,  strCount, strReplace
+  PUBLIC :: strUpper, strTail, strStackm, strReduce, strFind, strParse, cat, find
+  PUBLIC :: dispTable, dispOutliers, dispNameList
+  PUBLIC :: is_numeric, bool2str, int2str, real2str, dble2str
+  PUBLIC :: reduceExpr, str2bool, str2int, str2real, str2dble
+  PUBLIC :: addQuotes, testFile, checkList, removeComment, test
+
+  INTERFACE get_in;     MODULE PROCEDURE getin_s,  getin_i,  getin_r,  getin_l;  END INTERFACE get_in
+  INTERFACE  msg;       MODULE PROCEDURE        msg_1,      msg_l1,       msg_m; END INTERFACE  msg
+  INTERFACE fmsg;       MODULE PROCEDURE       fmsg_1,     fmsg_l1,      fmsg_m; END INTERFACE fmsg
+  INTERFACE strHead;    MODULE PROCEDURE    strHead_1,                strHead_m; END INTERFACE strHead
+  INTERFACE strTail;    MODULE PROCEDURE    strTail_1,                strTail_m; END INTERFACE strTail
+  INTERFACE strStack;   MODULE PROCEDURE   strStack_1,               strStack_2; END INTERFACE strStack
+  INTERFACE strClean;   MODULE PROCEDURE   strClean_1,               strClean_m; END INTERFACE strClean
+  INTERFACE strReduce;  MODULE PROCEDURE  strReduce_1,              strReduce_2; END INTERFACE strReduce
+  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, &
+!                 horzcat_d1,  horzcat_dm,
+                                           horzcat_sm,  horzcat_im,  horzcat_rm; END INTERFACE cat
+  INTERFACE find;       MODULE PROCEDURE      strFind,    find_int,    find_boo; END INTERFACE find
+  INTERFACE dispOutliers; MODULE PROCEDURE dispOutliers_1, dispOutliers_2; END INTERFACE dispOutliers
+  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
+
+  CHARACTER(LEN=256), SAVE :: modname = ''                 !--- Current subroutine name
+  INTEGER,            SAVE :: lunout  = 6                  !--- Printing unit  (default: 6, ie. on screen)
+  INTEGER,            SAVE :: prt_level = 1                !--- Printing level (default: 1, ie. print all)
+
+
+CONTAINS
+
+LOGICAL FUNCTION test(lcond, lout) RESULT(lerr)
+  LOGICAL, INTENT(IN)  :: lcond
+  LOGICAL, INTENT(OUT) :: lout
+  lerr = lcond; lout = lcond
+END FUNCTION test
+
+!==============================================================================================================================
+SUBROUTINE init_printout(lunout_, prt_level_)
+  INTEGER, INTENT(IN) :: lunout_, prt_level_
+  lunout = lunout_
+END SUBROUTINE init_printout
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== Same as getin ; additional last argument: the default value.
+!==============================================================================================================================
+!==============================================================================================================================
+!=== Same as getin ; additional last argument: the default value.
+!==============================================================================================================================
+SUBROUTINE getin_s(nam, val, def)
+USE ioipsl_getin_p_mod, ONLY: getin => getin_p
+  CHARACTER(LEN=*), INTENT(IN)    :: nam
+  CHARACTER(LEN=*), INTENT(INOUT) :: val
+  CHARACTER(LEN=*), INTENT(IN)    :: def
+  val = def; CALL getin(nam, val)
+  WRITE(lunout,*)TRIM(nam)//' = '//TRIM(val)
+END SUBROUTINE getin_s
+!==============================================================================================================================
+SUBROUTINE getin_i(nam, val, def)
+USE ioipsl_getin_p_mod, ONLY: getin => getin_p
+  CHARACTER(LEN=*), INTENT(IN)    :: nam
+  INTEGER,          INTENT(INOUT) :: val
+  INTEGER,          INTENT(IN)    :: def
+  val = def; CALL getin(nam, val)
+  WRITE(lunout,*)TRIM(nam)//' = '//TRIM(int2str(val))
+END SUBROUTINE getin_i
+!==============================================================================================================================
+SUBROUTINE getin_r(nam, val, def)
+USE ioipsl_getin_p_mod, ONLY: getin => getin_p
+  CHARACTER(LEN=*), INTENT(IN)    :: nam
+  REAL,             INTENT(INOUT) :: val
+  REAL,             INTENT(IN)    :: def
+  val = def; CALL getin(nam, val)
+  WRITE(lunout,*)TRIM(nam)//' = '//TRIM(real2str(val))
+END SUBROUTINE getin_r
+!==============================================================================================================================
+SUBROUTINE getin_l(nam, val, def)
+USE ioipsl_getin_p_mod, ONLY: getin => getin_p
+  CHARACTER(LEN=*), INTENT(IN)    :: nam
+  LOGICAL,          INTENT(INOUT) :: val
+  LOGICAL,          INTENT(IN)    :: def
+  val = def; CALL getin(nam, val)
+  WRITE(lunout,*)TRIM(nam)//' = '//TRIM(bool2str(val))
+END SUBROUTINE getin_l
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== Display one or several messages, one each line, starting with the current routine name "modname".
+!==============================================================================================================================
+SUBROUTINE msg_1(str, unit)
+  CHARACTER(LEN=*),  INTENT(IN) :: str
+  INTEGER, OPTIONAL, INTENT(IN) :: unit
+  INTEGER :: unt
+  unt = lunout; IF(PRESENT(unit)) unt = unit
+  WRITE(unt,'(a)') TRIM(modname)//': '//str
+END SUBROUTINE msg_1
+!==============================================================================================================================
+SUBROUTINE msg_l1(ll, str, unit)
+  LOGICAL,           INTENT(IN) :: ll
+  CHARACTER(LEN=*),  INTENT(IN) :: str
+  INTEGER, OPTIONAL, INTENT(IN) :: unit
+  INTEGER :: unt
+  IF(.NOT.ll) RETURN
+  unt = lunout; IF(PRESENT(unit)) unt = unit
+  WRITE(unt,'(a)') TRIM(modname)//': '//str
+END SUBROUTINE msg_l1
+!==============================================================================================================================
+SUBROUTINE msg_m(str, unit, nmax)
+  CHARACTER(LEN=*),  INTENT(IN) :: str(:)
+  INTEGER, OPTIONAL, INTENT(IN) :: unit
+  INTEGER, OPTIONAL, INTENT(IN) :: nmax
+  CHARACTER(LEN=256), ALLOCATABLE :: s(:)
+  INTEGER :: unt, nmx, k
+  unt = lunout ; IF(PRESENT(unit)) unt = unit
+  nmx = 128;     IF(PRESENT(nmax)) nmx = nmax
+  s = strStackm(str, ', ', nmx)
+  DO k=1,SIZE(s); WRITE(unt,'(a)') TRIM(modname)//': '//TRIM(s(k)); END DO
+END SUBROUTINE msg_m
+!==============================================================================================================================
+LOGICAL FUNCTION fmsg_1(str, unit) RESULT(lerr)
+  CHARACTER(LEN=*),  INTENT(IN)  :: str
+  INTEGER, OPTIONAL, INTENT(IN)  :: unit
+  INTEGER :: unt
+  lerr = .TRUE.
+  unt = lunout ; IF(PRESENT(unit)) unt = unit
+  CALL msg_1(str, unt)
+END FUNCTION fmsg_1
+!==============================================================================================================================
+LOGICAL FUNCTION fmsg_l1(li, str, unit) RESULT(lerr)
+  LOGICAL,           INTENT(IN)  :: li
+  CHARACTER(LEN=*),  INTENT(IN)  :: str
+  INTEGER, OPTIONAL, INTENT(IN)  :: unit
+  INTEGER :: unt
+  lerr = li;     IF(.NOT.lerr) RETURN
+  unt = lunout ; IF(PRESENT(unit)) unt = unit
+  CALL msg_l1(lerr, str, unt)
+END FUNCTION fmsg_l1
+!==============================================================================================================================
+LOGICAL FUNCTION fmsg_m(str, unit, nmax) RESULT(lerr)
+  CHARACTER(LEN=*),  INTENT(IN)  :: str(:)
+  INTEGER, OPTIONAL, INTENT(IN)  :: unit
+  INTEGER, OPTIONAL, INTENT(IN)  :: nmax
+  INTEGER :: unt, nmx
+  lerr = .TRUE.
+  unt = lunout ; IF(PRESENT(unit)) unt = unit
+  nmx = 128;     IF(PRESENT(nmax)) nmx = nmax
+  CALL msg_m(str, unt, nmx)
+END FUNCTION fmsg_m
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== Lower/upper case conversion function. ====================================================================================
+!==============================================================================================================================
+ELEMENTAL CHARACTER(LEN=256) FUNCTION strLower(str) RESULT(out)
+  CHARACTER(LEN=*), INTENT(IN) :: str
+  INTEGER :: k
+  out = str
+  DO k=1,LEN_TRIM(str)
+    IF(str(k:k)>='A'.OR.str(k:k)<='Z') out(k:k)=ACHAR(IACHAR(str(k:k))-32)
+  END DO
+END FUNCTION strLower
+!==============================================================================================================================
+ELEMENTAL CHARACTER(LEN=256) FUNCTION strUpper(str) RESULT(out)
+  CHARACTER(LEN=*), INTENT(IN) :: str
+  INTEGER :: k
+  out = str
+  DO k=1,LEN_TRIM(str)
+    IF(str(k:k)>='a'.OR.str(k:k)<='z') out(k:k)=ACHAR(IACHAR(str(k:k))+32)
+  END DO
+END FUNCTION strUpper
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== Extract first or last element of a string "str" delimited by separator "sep" =============================================
+!==============================================================================================================================
+CHARACTER(LEN=256) FUNCTION strHead_1(str,sep) RESULT(out)
+  CHARACTER(LEN=*), INTENT(IN) :: str, sep
+  out = str(1:INDEX(str,sep)-1); IF(out=='') out=str
+END FUNCTION strHead_1
+!==============================================================================================================================
+CHARACTER(LEN=256) FUNCTION strTail_1(str,sep) RESULT(out)
+  CHARACTER(LEN=*), INTENT(IN) :: str, sep
+  out = str(INDEX(str,sep,.TRUE.)+LEN(sep):LEN_TRIM(str)); IF(out=='') out=str
+END FUNCTION strTail_1
+!==============================================================================================================================
+FUNCTION strHead_m(str,sep) RESULT(out)
+  CHARACTER(LEN=256), ALLOCATABLE :: out(:)
+  CHARACTER(LEN=*),    INTENT(IN) :: str(:), sep
+  INTEGER :: k
+  out = [(strHead_1(str(k),sep), k=1, SIZE(str))]
+END FUNCTION strHead_m
+!==============================================================================================================================
+FUNCTION strTail_m(str,sep) RESULT(out)
+  CHARACTER(LEN=256), ALLOCATABLE :: out(:)
+  CHARACTER(LEN=*),    INTENT(IN) :: str(:), sep
+  INTEGER :: k
+  out = [(strTail_1(str(k),sep), k=1, SIZE(str))]
+END FUNCTION strTail_m
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== Concatenates the strings "str(:)" with separator "sep" into a single string using a separator (',' by default). ==========
+!==============================================================================================================================
+FUNCTION strStack_1(str, sep) RESULT(out)
+  CHARACTER(LEN=:),          ALLOCATABLE :: out
+  CHARACTER(LEN=*),           INTENT(IN) :: str(:)
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: sep
+  CHARACTER(LEN=:), ALLOCATABLE :: s
+  INTEGER :: is
+  IF(SIZE(str) == 0) THEN; out = ''; RETURN; END IF
+  ALLOCATE(s, SOURCE=','); IF(PRESENT(sep)) s=sep
+  out = TRIM(str(1)); DO is=2,SIZE(str, DIM=1); out = TRIM(out)//s//' '//TRIM(str(is)); END DO
+END FUNCTION strStack_1
+!==============================================================================================================================
+FUNCTION strStack_2(str1, str2, sep) RESULT(out)
+  CHARACTER(LEN=:),          ALLOCATABLE :: out
+  CHARACTER(LEN=*),           INTENT(IN) :: str1(:), str2(:)
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: sep
+  CHARACTER(LEN=:), ALLOCATABLE :: s
+  INTEGER :: is
+  IF(SIZE(str1) == 0) THEN; out = ''; RETURN; END IF
+  ALLOCATE(s, SOURCE=','); IF(PRESENT(sep)) s=sep
+  out = TRIM(str1(1))//'='//TRIM(str2(1))
+  DO is=2,SIZE(str1); out = TRIM(out)//s//' '//TRIM(str1(is))//'='//TRIM(str2(is)); END DO
+END FUNCTION strStack_2
+!==============================================================================================================================
+!=== Concatenate the strings "str(:)" with separator "sep" into one or several lines of "nmax" characters max (for display) ===
+!==============================================================================================================================
+FUNCTION strStackm(str, sep, nmax) RESULT(out)
+  CHARACTER(LEN=256),        ALLOCATABLE :: out(:)
+  CHARACTER(LEN=*),           INTENT(IN) :: str(:)
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: sep
+  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
+  CHARACTER(LEN=256), ALLOCATABLE :: t(:)
+  CHARACTER(LEN=256) :: sp
+  INTEGER :: is, ns, no, mx, n
+  IF(SIZE(str) == 0) THEN; out = ['']; RETURN; END IF
+  sp =', '; IF(PRESENT(sep )) sp = sep
+  ns = 2  ; IF(PRESENT(sep )) ns = LEN(sep)
+  mx = 256; IF(PRESENT(nmax)) mx = nmax
+  no = 1; out = ['']
+  DO is = 1, SIZE(str)
+    n = LEN_TRIM(str(is)); IF(out(no)/='') n = n+ns+LEN_TRIM(out(no))          !--- Line length after "str(is)" inclusion
+    IF(out(no) == '') THEN
+      out(no) = str(is)                                                        !--- Empty new line: set to "str(is)"
+    ELSE IF(n <= mx) THEN
+      out(no) = TRIM(out(no))//sp(1:ns)//TRIM(str(is))                         !--- Append "str(is)" to the current line
+    ELSE
+      ALLOCATE(t(no+1)); t(1:no) = out; no=no+1; t(no) = str(is)               !--- Full line: "str(si)" put in next line
+      CALL MOVE_ALLOC(FROM=t, TO=out)
+    END IF
+  END DO
+END FUNCTION strStackm
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== String cleaning: replace tabulation by spaces, remove NULL characters and comments. ======================================
+!==============================================================================================================================
+SUBROUTINE strClean_1(str)
+  CHARACTER(LEN=*), INTENT(INOUT) :: str
+  INTEGER :: k, n, m
+  n = LEN(str)
+  DO k = n, 1, -1
+    m = IACHAR(str(k:k))
+    IF(m==9) str(k:k) = ' '                           !--- Replace the tabulations with spaces
+    IF(m==0) str(k:n) = str(k+1:n)//' '               !--- Remove the NULL characters
+  END DO
+  m = INDEX(str,'!')-1; IF(m==-1) m = LEN_TRIM(str)   !--- Remove end of line comment
+  str = ADJUSTL(str(1:m))
+END SUBROUTINE strClean_1
+!==============================================================================================================================
+SUBROUTINE strClean_m(str)
+  CHARACTER(LEN=*), INTENT(INOUT) :: str(:)
+  INTEGER :: k
+  DO k = 1, SIZE(str); CALL strClean_1(str(k)); END DO
+END SUBROUTINE strClean_m
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== strReduce_1(str1)     : Remove duplicated elements of str1.    ===========================================================
+!=== strReduce_2(str1,str2): Append str1 with new elements of str2. ===========================================================
+!==============================================================================================================================
+SUBROUTINE strReduce_1(str1, nb)
+  CHARACTER(LEN=*), ALLOCATABLE, INTENT(INOUT) :: str1(:)
+  INTEGER,          OPTIONAL,    INTENT(OUT)   :: nb
+  CHARACTER(LEN=256), ALLOCATABLE :: s1(:)
+  INTEGER :: k, n, n1
+  CALL MOVE_ALLOC(FROM = str1, TO = s1); CALL strClean(s1)
+  n1 = SIZE(s1, DIM=1)                                     !--- Total nb of  elements in "s1"
+  n  = COUNT( [( ALL(s1(1:k-1)/=s1(k)), k=1, n1 )] )       !--- Nb of unique elements in "s1"
+  ALLOCATE(str1(n)); IF(n==0) RETURN; str1(1) = s1(1)
+  n=1; DO k=2,n1; IF(ANY(s1(1:k-1)==s1(k))) CYCLE; n=n+1; str1(n)=s1(k); END DO
+  IF(PRESENT(nb)) nb = n
+END SUBROUTINE strReduce_1
+!==============================================================================================================================
+SUBROUTINE strReduce_2(str1, str2)
+  CHARACTER(LEN=*),   ALLOCATABLE, INTENT(INOUT) :: str1(:)
+  CHARACTER(LEN=*),                INTENT(IN)    :: str2(:)
+  CHARACTER(LEN=256), ALLOCATABLE :: s1(:), s2(:)
+  INTEGER :: k
+  IF(SIZE(str2)==0) RETURN
+  s2 = str2; CALL strClean(s2)
+  IF(.NOT.ALLOCATED(s2)) RETURN
+  IF(SIZE(s2) == 0) THEN; DEALLOCATE(s2); RETURN; END IF
+  IF(.NOT.ALLOCATED(str1)) THEN
+    str1 = s2
+  ELSE IF(SIZE(str1)==0) THEN
+    str1 = s2
+  ELSE
+    s1 = str1; CALL strClean(s1)
+    str1 = [s1, PACK(s2, MASK= [( ALL(s1(:) /= s2(k)), k=1, SIZE(s2) )] ) ]
+  END IF
+END SUBROUTINE strReduce_2
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== GET THE INDEX OF THE FIRST APPEARANCE IN THE STRING VECTOR "str(:)" OF THE STRING(s) "s" =================================
+!==============================================================================================================================
+INTEGER FUNCTION strIdx_1(str, s) RESULT(out)
+  CHARACTER(LEN=*), INTENT(IN) :: str(:), s
+  DO out = 1, SIZE(str); IF(str(out) == s) EXIT; END DO
+  IF(out == 1+SIZE(str) .OR. SIZE(str)==0) out = 0
+END FUNCTION strIdx_1
+!==============================================================================================================================
+FUNCTION strIdx_m(str, s, n) RESULT(out)
+  CHARACTER(LEN=*),  INTENT(IN)  :: str(:), s(:)
+  INTEGER, OPTIONAL, INTENT(OUT) :: n
+  INTEGER,           ALLOCATABLE :: out(:)
+  INTEGER :: k
+  out = [(strIdx_1(str(:), s(k)), k=1, SIZE(s))]
+  IF(PRESENT(n)) n = COUNT(out(:)/=0)
+END FUNCTION strIdx_m
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== GET THE INDEX LIST OF THE ELEMENTS OF "str(:)" EQUAL TO "s" AND OPTIONALY, ITS LENGTH "n" ================================
+!==============================================================================================================================
+FUNCTION strFind(str, s, n) RESULT(out)
+  CHARACTER(LEN=*),  INTENT(IN)  :: str(:), s
+  INTEGER, OPTIONAL, INTENT(OUT) :: n
+  INTEGER,           ALLOCATABLE :: out(:)
+  INTEGER :: k
+  out = PACK( [(k, k=1, SIZE(str(:), DIM=1))], MASK = str(:) == s )
+  IF(PRESENT(n)) n = SIZE(out(:), DIM=1)
+END FUNCTION strFind
+!==============================================================================================================================
+FUNCTION find_int(i,j,n) RESULT(out)
+  INTEGER,           INTENT(IN)  :: i(:), j
+  INTEGER, OPTIONAL, INTENT(OUT) :: n
+  INTEGER,           ALLOCATABLE :: out(:)
+  INTEGER :: k
+  out = PACK( [(k, k=1, SIZE(i(:), DIM=1))], MASK = i(:) == j )
+  IF(PRESENT(n)) n = SIZE(out(:), DIM=1)
+END FUNCTION find_int
+!==============================================================================================================================
+FUNCTION find_boo(l,n) RESULT(out)
+  LOGICAL,           INTENT(IN)  :: l(:)
+  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"                                 ==========================================
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== GET THE INDEX IN "rawList" OF THE 1ST APPEARANCE OF ONE OF THE "del(:)" SEPARATORS (0 IF NONE OF THEM ARE PRESENT)
+!===  IF lSc == .TRUE.:  * SKIP HEAD SIGNS OR EXPONENTS SIGNS THAT SHOULD NOT BE CONFUSED WITH SEPARATORS
+!===                     * THEN TEST WHETHER THE STRING FROM START TO THE FOUND SEPARATOR IS A CORRECTLY FORMATTED NUMBER
+!==============================================================================================================================
+LOGICAL FUNCTION strIdx_prv(rawList, del, ibeg, idx, idel, lSc) RESULT(lerr)
+  CHARACTER(LEN=*),  INTENT(IN)  :: rawList                          !--- String in which delimiters have to be identified
+  CHARACTER(LEN=*),  INTENT(IN)  :: del(:)                           !--- List of delimiters
+  INTEGER,           INTENT(IN)  :: ibeg                             !--- Start index
+  INTEGER,           INTENT(OUT) :: idx                              !--- Index of the first identified delimiter in "rawList"
+  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, ALLOCATABLE :: ii(:)
+  LOGICAL              :: ll
+  CHARACTER(LEN=256)   :: d
+
+  lerr = .FALSE.
+  idx = strIdx1(rawList, del, ibeg, idel)
+  IF(idx == 0)                  RETURN                               !--- No separator found                      => finished
+  IF(.NOT.PRESENT(lSc))         RETURN; IF(.NOT.lSc ) RETURN         !--- No need to check exceptions for numbers => finished
+  IF(INDEX('+-',del(idel))==0)  RETURN                               !--- No possible sign ambiguity              => finished
+  ll = idx == 1                                                      !--- This is a front sign of a number
+  IF(idx /= 1) ll = INDEX('^*/+-',rawList(idx-1:idx-1)) /= 0         !--- This is a front sign of a number after an operator
+  IF(ll)       idx = strIdx1(rawList, del,idx+1,idel)                !---   => TO THE NEXT DELIMITER
+  IF(idx /= 0) THEN
+    IF(idx /= 1) ll = INDEX('eE' ,rawList(idx-1:idx-1)) /= 0         !--- Exponent sign of a number in scientific notation ?
+    IF(ll)     idx = strIdx1(rawList, del,idx+1,idel)                !--- Identify the next delimiter
+  END IF
+  IF(idx /= 0) lerr = is_numeric(rawList(ibeg:idx              -1))  !--- Check whether the previous word was a number
+  IF(idx == 0) lerr = is_numeric(rawList(ibeg:LEN_TRIM(rawList)-1))  !--- Check whether the previous word was a number
+  IF(idx == 0) idel = 0
+
+CONTAINS
+
+INTEGER FUNCTION strIdx1(str, del, ib, id) RESULT(idx)
+  CHARACTER(LEN=*),  INTENT(IN)  :: str
+  CHARACTER(LEN=*),  INTENT(IN)  :: 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
+END FUNCTION strIdx1
+
+END FUNCTION strIdx_prv
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== Count the number of elements separated by "delimiter" in list "rawList". =================================================
+!==============================================================================================================================
+LOGICAL FUNCTION strCount_11(rawList, delimiter, nb, lSc) RESULT(out)
+  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)
+END FUNCTION strCount_11
+!==============================================================================================================================
+LOGICAL FUNCTION strCount_m1(rawList, delimiter, nb, lSc) RESULT(out)
+  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.
+  ALLOCATE(nb(SIZE(rawList)))
+  DO id = 1, SIZE(rawList)
+    out = out .AND. strCount_1m(rawList(id), [delimiter], nb(id), ll)
+  END DO
+END FUNCTION strCount_m1
+!==============================================================================================================================
+LOGICAL FUNCTION strCount_1m(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
+
+  INTEGER              :: ib, ie, jd, nr
+  LOGICAL              :: ll
+  CHARACTER(LEN=1024)  :: r
+
+  modname = 'strCount'
+  lerr = .FALSE.
+  ll   = .FALSE.; IF(PRESENT(lSc)) ll = lSc
+  r  = TRIM(ADJUSTL(rawList))
+  nr = LEN_TRIM(r); IF(nr == 0) RETURN
+  nb = 1; ib = 1
+  DO
+    lerr = strIdx_prv(r, delimiter, ib, ie, jd, ll)
+    IF(jd == 0) EXIT
+    IF(fmsg(lerr,'"'//TRIM(r(ib:ie-1))//'" is not numeric')) RETURN
+    ib = ie + LEN(delimiter(jd))
+    DO WHILE(r(ib:ib) == ' ' .AND. ib < nr); ib = ib + 1; END DO     !--- Skip spaces before next chain
+    nb = nb + 1
+  END DO
+END FUNCTION strCount_1m
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== Purpose: Parse "delimiter"-separated list "rawList" into the pair keys(:), vals(:).   ====================================
+!===          Corresponding "vals" remains empty if the element does not contain "=" sign. ====================================
+!==============================================================================================================================
+LOGICAL FUNCTION strParse_1(rawList, delimiter, keys, lSc, vals, n) RESULT(lerr)
+  CHARACTER(LEN=*),                          INTENT(IN)  :: rawList, delimiter
+  CHARACTER(LEN=256), ALLOCATABLE,           INTENT(OUT) :: keys(:)
+  LOGICAL,                         OPTIONAL, INTENT(IN)  :: lSc      !--- Take care about numbers in scientific notation
+  CHARACTER(LEN=256), 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=*),                          INTENT(IN)  :: rawList, delimiter(:)
+  CHARACTER(LEN=256),           ALLOCATABLE, INTENT(OUT) :: keys(:)
+  LOGICAL,            OPTIONAL,              INTENT(IN)  :: lSc      !--- Take care about numbers in scientific notation
+  CHARACTER(LEN=256), OPTIONAL, ALLOCATABLE, INTENT(OUT) :: vals(:)
+  INTEGER,            OPTIONAL,              INTENT(OUT) :: n
+  INTEGER,            OPTIONAL, ALLOCATABLE, INTENT(OUT) :: id(:)
+
+  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(strCount_1m(rawList, delimiter, nk, ll), "Couldn't parse list: non-numerical strings were found"),lerr)) RETURN
+
+  !--- FEW ALLOCATIONS
+  ALLOCATE(keys(nk))
+  IF(PRESENT(vals)) ALLOCATE(vals(nk))
+  IF(PRESENT(id))   ALLOCATE(id(nk-1))
+  IF(PRESENT(n)) n = nk
+
+  !--- PARSING
+  r  = TRIM(ADJUSTL(rawList))
+  nr = LEN_TRIM(r); IF(nr == 0) RETURN
+  ib = 1
+  DO ik = 1, nk-1
+    IF(test(fmsg(strIdx_prv(r, delimiter, ib, ie, jd, ll),'Non-numeric values found'),lerr)) RETURN
+    keys(ik) = r(ib:ie-1)
+    IF(PRESENT(vals)) CALL parseKeys(keys(ik), vals(ik))             !--- Parse a <key>=<val> pair
+    IF(PRESENT(id  )) id(ik) = jd                                    !--- Index in "delimiter(:)" of the "ik"th delimiter
+    ib = ie + LEN_TRIM( delimiter(jd) )                              !--- Length of the current delimiter
+    DO WHILE(r(ib:ib) == ' ' .AND. ib < nr); ib = ib + 1; END DO     !--- Skip spaces before next chain
+  END DO
+  keys(nk) = r(ib:nr)
+  IF(PRESENT(vals)) CALL parseKeys(keys(nk), vals(nk))               !--- Parse a <key>=<val> pair
+
+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"
+  val = ADJUSTL(key(ix+1:LEN_TRIM(key)))
+  key = ADJUSTL(key(1:ix-1))
+END SUBROUTINE parseKeys
+
+END FUNCTION strParse_m    
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== String substitution: replace "key" by "val" each time it appears in "str".
+!==============================================================================================================================
+SUBROUTINE strReplace_1(str, key, val, lsurr)
+  CHARACTER(LEN=*),  INTENT(INOUT) :: str        !--- Main string
+  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)
+  i0 = 1
+  DO
+    ns = LEN_TRIM(str)
+    ix = INDEX(str(i0:ns), TRIM(key))            !--- First appearance index of "key" in "s", starting from index "i0"
+    IF(ix == 0) EXIT
+    ix = ix + i0 -1
+    IF(lsur) THEN                                !--- Key must be surrounded by special characters
+      !--- lb=.TRUE.: key is at the very beginning of "str" or located after  a special character
+      lb = ix     ==1;  IF(.NOT.lb) lb = INDEX('+-*/()^', str(ix-1 :ix-1 ))/=0
+      !--- le=.TRUE.: key is at the very end       of "str" or located before a special character
+      le = ix+nk-1==ns; IF(.NOT.le) le = INDEX('+-*/()^', str(ix+nk:ix+nk))/=0
+      IF(.NOT.(lb.AND.le)) THEN; i0 = i0 + nk; CYCLE; END IF
+    END IF
+    str = str(1:ix-1)//TRIM(val)//str(ix+nk:ns)
+  END DO
+END SUBROUTINE strReplace_1
+!==============================================================================================================================
+SUBROUTINE strReplace_m(str, key, val, lsurr)
+  CHARACTER(LEN=*),  INTENT(INOUT) :: str(:)     !--- Main strings vector
+  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
+  INTEGER :: k
+  LOGICAL :: ll
+  ll=.FALSE.; IF(PRESENT(lsurr)) ll=lsurr
+  DO k=1, SIZE(str); CALL strReplace_1(str(k),key,val,ll); END DO
+END SUBROUTINE strReplace_m
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== Contatenate horizontally scalars/vectors of strings/integers/reals into a vector/array ===================================
+!==============================================================================================================================
+FUNCTION horzcat_s1(s0, s1, s2, s3, s4, s5, s6, s7, s8, s9) RESULT(out)
+  CHARACTER(LEN=*),           TARGET, INTENT(IN) :: s0
+  CHARACTER(LEN=*), OPTIONAL, TARGET, INTENT(IN) :: s1, s2, s3, s4, s5, s6, s7, s8, s9
+  CHARACTER(LEN=256), ALLOCATABLE :: out(:)
+  CHARACTER(LEN=256), POINTER     :: s
+  LOGICAL :: lv(10)
+  INTEGER :: iv
+  lv = [   .TRUE.   , PRESENT(s1), PRESENT(s2), PRESENT(s3), PRESENT(s4) , &
+         PRESENT(s5), PRESENT(s6), PRESENT(s7), PRESENT(s8), PRESENT(s9) ]
+  ALLOCATE(out(COUNT(lv)))
+  DO iv=1, COUNT(lv)
+    SELECT CASE(iv)
+      CASE(1); s=> s0; CASE(2); s=> s1; CASE(3); s=> s2; CASE(4); s=> s3; CASE(5); s=> s4
+      CASE(6); s=> s5; CASE(7); s=> s6; CASE(8); s=> s7; CASE(9); s=> s8; CASE(10);s=> s9
+    END SELECT
+    out(iv) = s
+  END DO
+END FUNCTION horzcat_s1
+!==============================================================================================================================
+FUNCTION horzcat_sm(s0, s1, s2, s3, s4, s5, s6, s7, s8, s9) RESULT(out)
+  CHARACTER(LEN=*),           TARGET, DIMENSION(:), INTENT(IN) :: s0
+  CHARACTER(LEN=*), OPTIONAL, TARGET, DIMENSION(:), INTENT(IN) :: s1, s2, s3, s4, s5, s6, s7, s8, s9
+  CHARACTER(LEN=256), ALLOCATABLE :: out(:,:)
+  CHARACTER(LEN=256), POINTER     :: s(:)
+  LOGICAL :: lv(10)
+  INTEGER :: nrow, ncol, iv, n
+  lv = [   .TRUE.   , PRESENT(s1), PRESENT(s2), PRESENT(s3), PRESENT(s4) , &
+         PRESENT(s5), PRESENT(s6), PRESENT(s7), PRESENT(s8), PRESENT(s9) ]
+  nrow = SIZE(s0); ncol=COUNT(lv)
+  ALLOCATE(out(nrow, ncol))
+  DO iv=1, ncol
+    SELECT CASE(iv)
+      CASE(1); s=> s0; CASE(2); s=> s1; CASE(3); s=> s2; CASE(4); s=> s3; CASE(5); s=> s4
+      CASE(6); s=> s5; CASE(7); s=> s6; CASE(8); s=> s7; CASE(9); s=> s8; CASE(10);s=> s9
+    END SELECT
+    n = SIZE(s, DIM=1)
+    IF(n/=nrow) THEN; CALL msg('Can''t concatenate vectors of differing lengths',1); STOP; END IF
+    out(:,iv) = s(:)
+  END DO
+END FUNCTION horzcat_sm
+!==============================================================================================================================
+FUNCTION horzcat_i1(i0, i1, i2, i3, i4, i5, i6, i7, i8, i9) RESULT(out)
+  INTEGER,           TARGET, INTENT(IN) :: i0
+  INTEGER, OPTIONAL, TARGET, INTENT(IN) :: i1, i2, i3, i4, i5, i6, i7, i8, i9
+  INTEGER, ALLOCATABLE :: out(:)
+  INTEGER, POINTER     :: i
+  LOGICAL :: lv(10)
+  INTEGER :: iv
+  lv = [   .TRUE.   , PRESENT(i1), PRESENT(i2), PRESENT(i3), PRESENT(i4) , &
+         PRESENT(i5), PRESENT(i6), PRESENT(i7), PRESENT(i8), PRESENT(i9) ]
+  ALLOCATE(out(COUNT(lv)))
+  DO iv=1, COUNT(lv)
+    SELECT CASE(iv)
+      CASE(1); i=> i0; CASE(2); i=> i1; CASE(3); i=> i2; CASE(4); i=> i3; CASE(5); i=> i4
+      CASE(6); i=> i5; CASE(7); i=> i6; CASE(8); i=> i7; CASE(9); i=> i8; CASE(10);i=> i9
+    END SELECT
+    out(iv) = i
+  END DO
+END FUNCTION horzcat_i1
+!==============================================================================================================================
+FUNCTION horzcat_im(i0, i1, i2, i3, i4, i5, i6, i7, i8, i9) RESULT(out)
+  INTEGER,           TARGET, DIMENSION(:), INTENT(IN) :: i0
+  INTEGER, OPTIONAL, TARGET, DIMENSION(:), INTENT(IN) :: i1, i2, i3, i4, i5, i6, i7, i8, i9
+  INTEGER, ALLOCATABLE :: out(:,:)
+  INTEGER, POINTER     :: i(:)
+  LOGICAL :: lv(10)
+  INTEGER :: nrow, ncol, iv, n
+  lv = [   .TRUE.   , PRESENT(i1), PRESENT(i2), PRESENT(i3), PRESENT(i4) , &
+         PRESENT(i5), PRESENT(i6), PRESENT(i7), PRESENT(i8), PRESENT(i9) ]
+  nrow = SIZE(i0); ncol=COUNT(lv)
+  ALLOCATE(out(nrow, ncol))
+  DO iv=1, ncol
+    SELECT CASE(iv)
+      CASE(1); i=> i0; CASE(2); i=> i1; CASE(3); i=> i2; CASE(4); i=> i3; CASE(5); i=> i4
+      CASE(6); i=> i5; CASE(7); i=> i6; CASE(8); i=> i7; CASE(9); i=> i8; CASE(10);i=> i9
+    END SELECT
+    n = SIZE(i, DIM=1)
+    IF(n/=nrow) THEN; CALL msg('Can''t concatenate vectors of differing lengths',1); STOP; END IF
+    out(:,iv) = i(:)
+  END DO
+END FUNCTION horzcat_im
+!==============================================================================================================================
+FUNCTION horzcat_r1(r0, r1, r2, r3, r4, r5, r6, r7, r8, r9) RESULT(out)
+  REAL,           TARGET, INTENT(IN) :: r0
+  REAL, OPTIONAL, TARGET, INTENT(IN) :: r1, r2, r3, r4, r5, r6, r7, r8, r9
+  REAL, ALLOCATABLE :: out(:)
+  REAL, POINTER     :: r
+  LOGICAL :: lv(10)
+  INTEGER :: iv
+  lv = [   .TRUE.   , PRESENT(r1), PRESENT(r2), PRESENT(r3), PRESENT(r4) , &
+         PRESENT(r5), PRESENT(r6), PRESENT(r7), PRESENT(r8), PRESENT(r9) ]
+  ALLOCATE(out(COUNT(lv)))
+  DO iv=1, COUNT(lv)
+    SELECT CASE(iv)
+      CASE(1); r=> r0; CASE(2); r=> r1; CASE(3); r=> r2; CASE(4); r=> r3; CASE(5); r=> r4
+      CASE(6); r=> r5; CASE(7); r=> r6; CASE(8); r=> r7; CASE(9); r=> r8; CASE(10);r=> r9
+    END SELECT
+    out(iv) = r
+  END DO
+END FUNCTION horzcat_r1
+!==============================================================================================================================
+FUNCTION horzcat_rm(r0, r1, r2, r3, r4, r5, r6, r7, r8, r9) RESULT(out)
+  REAL,           TARGET, DIMENSION(:), INTENT(IN) :: r0
+  REAL, OPTIONAL, TARGET, DIMENSION(:), INTENT(IN) :: r1, r2, r3, r4, r5, r6, r7, r8, r9
+  REAL, ALLOCATABLE :: out(:,:)
+  REAL, POINTER     :: r(:)
+  LOGICAL :: lv(10)
+  INTEGER :: nrow, ncol, iv, n
+  lv = [   .TRUE.   , PRESENT(r1), PRESENT(r2), PRESENT(r3), PRESENT(r4) , &
+         PRESENT(r5), PRESENT(r6), PRESENT(r7), PRESENT(r8), PRESENT(r9) ]
+  nrow = SIZE(r0); ncol=COUNT(lv)
+  ALLOCATE(out(nrow, ncol))
+  DO iv=1, ncol
+    SELECT CASE(iv)
+      CASE(1); r=> r0; CASE(2); r=> r1; CASE(3); r=> r2; CASE(4); r=> r3; CASE(5); r=> r4
+      CASE(6); r=> r5; CASE(7); r=> r6; CASE(8); r=> r7; CASE(9); r=> r8; CASE(10);r=> r9
+    END SELECT
+    n = SIZE(r, DIM=1)
+    IF(n/=nrow) THEN; CALL msg('Can''t concatenate vectors of differing lengths',1); STOP; END IF
+    out(:,iv) = r(:)
+  END DO
+END FUNCTION horzcat_rm
+!==============================================================================================================================
+FUNCTION horzcat_d1(d0, d1, d2, d3, d4, d5, d6, d7, d8, d9) RESULT(out)
+  DOUBLE PRECISION,           TARGET, INTENT(IN) :: d0
+  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)
+  INTEGER :: iv
+  lv = [   .TRUE.   , PRESENT(d1), PRESENT(d2), PRESENT(d3), PRESENT(d4) , &
+         PRESENT(d5), PRESENT(d6), PRESENT(d7), PRESENT(d8), PRESENT(d9) ]
+  ALLOCATE(out(COUNT(lv)))
+  DO iv=1, COUNT(lv)
+    SELECT CASE(iv)
+      CASE(1); d=> d0; CASE(2); d=> d1; CASE(3); d=> d2; CASE(4); d=> d3; CASE(5); d=> d4
+      CASE(6); d=> d5; CASE(7); d=> d6; CASE(8); d=> d7; CASE(9); d=> d8; CASE(10);d=> d9
+    END SELECT
+    out(iv) = d
+  END DO
+END FUNCTION horzcat_d1
+!==============================================================================================================================
+FUNCTION horzcat_dm(d0, d1, d2, d3, d4, d5, d6, d7, d8, d9) RESULT(out)
+  DOUBLE PRECISION,           TARGET, DIMENSION(:), INTENT(IN) :: d0
+  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)
+  INTEGER :: nrow, ncol, iv, n
+  lv = [   .TRUE.   , PRESENT(d1), PRESENT(d2), PRESENT(d3), PRESENT(d4) , &
+         PRESENT(d5), PRESENT(d6), PRESENT(d7), PRESENT(d8), PRESENT(d9) ]
+  nrow = SIZE(d0); ncol=COUNT(lv)
+  ALLOCATE(out(nrow, ncol))
+  DO iv=1, ncol
+    SELECT CASE(iv)
+      CASE(1); d=> d0; CASE(2); d=> d1; CASE(3); d=> d2; CASE(4); d=> d3; CASE(5); d=> d4
+      CASE(6); d=> d5; CASE(7); d=> d6; CASE(8); d=> d7; CASE(9); d=> d8; CASE(10);d=> d9
+    END SELECT
+    n = SIZE(d, DIM=1)
+    IF(n/=nrow) THEN; CALL msg('Can''t concatenate vectors of differing lengths',1); STOP; END IF
+    out(:,iv) = d(:)
+  END DO
+END FUNCTION horzcat_dm
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!--- Display a clean table composed of successive vectors of same length.
+!=== The profile "p" describe in which order to pick up the columns from "s", "i" and "r" for display.
+!==============================================================================================================================
+LOGICAL FUNCTION dispTable(p, titles, s, i, r, rFmt, nmax, unit) RESULT(lerr)
+  CHARACTER(LEN=*),           INTENT(IN)  :: p             !--- DISPLAY MAP: s/i/r
+  CHARACTER(LEN=*),           INTENT(IN)  :: titles(:)     !--- TITLES (ONE EACH COLUMN)
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: s(:,:)        !--- STRINGS
+  INTEGER,          OPTIONAL, INTENT(IN)  :: i(:,:)        !--- INTEGERS
+  REAL,             OPTIONAL, INTENT(IN)  :: r(:,:)        !--- REALS
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: rFmt          !--- Format for reals
+  INTEGER,          OPTIONAL, INTENT(IN)  :: nmax          !--- Display less than "nrow" rows
+  INTEGER,          OPTIONAL, INTENT(IN)  :: unit          !--- Output unit (default: screen)
+
+  CHARACTER(LEN=2048) :: row
+  CHARACTER(LEN=256)  :: rFm, el
+  CHARACTER(LEN=256), ALLOCATABLE :: d(:,:)
+  CHARACTER(LEN=1) :: s1, sp
+  INTEGER :: is, ii, ir, np, nrow, unt, ic
+  INTEGER :: ns, ni, nr, ncol, nmx
+  INTEGER, ALLOCATABLE :: n(:)
+  INTEGER, PARAMETER   :: nm=2                             !--- Space between values & columns
+  LOGICAL :: ls, li, lr
+
+  modname = 'dispTable'
+  rFm = '*';    IF(PRESENT(rFmt)) rFm = rFmt               !--- Specified format for reals
+  unt = lunout; IF(PRESENT(unit)) unt = unit               !--- Specified output unit
+  np = LEN_TRIM(p); ns = 0; ni = 0; nr = 0; ncol = 0
+  ls = PRESENT(s); li = PRESENT(i); lr = PRESENT(r)
+  lerr = .FALSE.; IF(.NOT.ANY([ls,li,lr])) RETURN          !--- Nothing to do
+  sp = '|'                                                 !--- Separator
+
+  !--- CHECK ARGUMENTS COHERENCE
+  lerr = np /= SIZE(titles); IF(fmsg(lerr, 'string "pattern" length and titles list mismatch')) RETURN
+  IF(ls) THEN; ns = SIZE(s, DIM=1); ncol = ncol + SIZE(s, DIM=2)
+    lerr = COUNT([(p(ic:ic)=='s', ic=1, np)]) /= SIZE(s, DIM=2)
+  END IF
+  IF(li) THEN; ni = SIZE(i, DIM=1); ncol = ncol + SIZE(i, DIM=2)
+    lerr = COUNT([(p(ic:ic)=='i', ic=1, np)]) /= SIZE(i, DIM=2)
+  END IF
+  IF(lr) THEN; nr = SIZE(r, DIM=1); ncol = ncol + SIZE(r, DIM=2)
+    lerr = COUNT([(p(ic:ic)=='r', ic=1, np)]) /= SIZE(r, DIM=2)
+  END IF
+
+  IF(fmsg(lerr, 'string "pattern" length and arguments number mismatch')) RETURN
+  lerr = ncol /= SIZE(titles); IF(fmsg(lerr, '"titles" length and arguments number mismatch')) RETURN
+  lerr = ls.AND.li.AND.ns/=ni; IF(fmsg(lerr, 'string and integer arguments lengths mismatch')) RETURN
+  lerr = ls.AND.lr.AND.ns/=nr; IF(fmsg(lerr,    'string and real arguments lengths mismatch')) RETURN
+  lerr = li.AND.lr.AND.ni/=nr; IF(fmsg(lerr,   'integer and real arguments lengths mismatch')) RETURN
+  nrow = MAX(ns,ni,nr)+1
+  nmx = nrow; IF(PRESENT(nmax)) nmx = MIN(nmx,nmax+1)
+
+  !--- Allocate the assembled quantities array
+  ALLOCATE(d(nrow,ncol), n(ncol))
+
+  !--- Assemble the vectors into a strings array in the order indicated by "pattern"
+  is =  1; ii = 1; ir = 1
+  DO ic = 1, ncol
+    d(1,ic) = TRIM(titles(ic))
+    SELECT CASE(p(ic:ic))
+      CASE('s'); d(2:nmx,ic) =          s(:,is)     ; is = is + 1
+      CASE('i'); d(2:nmx,ic) =  int2str(i(:,ii)    ); ii = ii + 1
+      CASE('r'); d(2:nmx,ic) = real2str(r(:,ir),rFm); ir = ir + 1
+    END SELECT
+    n(ic)=0; DO ir=1, nmx; n(ic)=MAX(n(ic), LEN_TRIM(d(ir,ic))); END DO
+  END DO
+  n(:) = n(:) + 2*nm
+
+  !--- Display the strings array as a table
+  DO ir = 1, nmx; row = ''
+    DO ic = 1, ncol; el = d(ir,ic)
+      s1 = sp
+      row = TRIM(row)//REPEAT(' ',nm)//TRIM(el)//REPEAT(' ',n(ic)-LEN_TRIM(el)-nm)//s1
+    END DO
+    nr = LEN_TRIM(row)-1                                             !--- Final separator removed
+    CALL msg(row(1:nr), unt)
+    IF(ir /= 1) CYCLE                                                !--- Titles are underlined
+    row=''; DO ic=1,ncol; row=TRIM(row)//REPEAT('-',n(ic))//'+'; END DO
+    CALL msg(row(1:LEN_TRIM(row)-1), unt)
+  END DO
+
+END FUNCTION dispTable
+!==============================================================================================================================
+
+!==============================================================================================================================
+LOGICAL FUNCTION dispNamelist(unt, p, titles, s, i, r, rFmt, llast) RESULT(lerr)
+  INTEGER,                    INTENT(IN)  :: unt           !--- Output unit
+  CHARACTER(LEN=*),           INTENT(IN)  :: p             !--- DISPLAY MAP: s/i/r
+  CHARACTER(LEN=*),           INTENT(IN)  :: titles(:)     !--- TITLES (ONE EACH COLUMN)
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: s(:,:)        !--- STRINGS
+  INTEGER,          OPTIONAL, INTENT(IN)  :: i(:,:)        !--- INTEGERS
+  REAL,             OPTIONAL, INTENT(IN)  :: r(:,:)        !--- REALS
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: rFmt          !--- Format for reals
+  LOGICAL,          OPTIONAL, INTENT(IN)  :: llast         !--- Last variable: no final ','
+
+  CHARACTER(LEN=256)  :: rFm, el
+  CHARACTER(LEN=256), ALLOCATABLE :: d(:,:)
+  CHARACTER(LEN=:),   ALLOCATABLE :: sp, row
+  INTEGER :: is, ii, ir, nrow, ic
+  INTEGER :: ns, ni, nr, ncol, np
+  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)
+  lerr = .FALSE.; IF(.NOT.ANY([ls,li,lr])) RETURN          !--- Nothing to do
+  la = .FALSE.; IF(PRESENT(llast)) la = llast
+
+  !--- CHECK ARGUMENTS COHERENCE
+  ns = 0; ni = 0; nr = 0; np = LEN_TRIM(p); ncol = 0
+  IF(ls) THEN; ns = SIZE(s, DIM=1); ncol = ncol + SIZE(s, DIM=2)
+    lerr = COUNT([(p(ic:ic)=='s', ic=1, np)]) /= SIZE(s, DIM=2)
+  END IF
+  IF(li) THEN; ni = SIZE(i, DIM=1); ncol = ncol + SIZE(i, DIM=2)
+    lerr = COUNT([(p(ic:ic)=='i', ic=1, np)]) /= SIZE(i, DIM=2)
+  END IF
+  IF(lr) THEN; nr = SIZE(r, DIM=1); ncol = ncol + SIZE(r, DIM=2)
+    lerr = COUNT([(p(ic:ic)=='r', ic=1, np)]) /= SIZE(r, DIM=2)
+  END IF
+  IF(fmsg(lerr, 'string "pattern" length and arguments number mismatch')) RETURN
+  lerr = ncol /= SIZE(titles); IF(fmsg(lerr, '"titles" length and arguments number mismatch')) RETURN
+  lerr = ls.AND.li.AND.ns/=ni; IF(fmsg(lerr, 'string and integer arguments lengths mismatch')) RETURN
+  lerr = ls.AND.lr.AND.ns/=nr; IF(fmsg(lerr,    'string and real arguments lengths mismatch')) RETURN
+  lerr = li.AND.lr.AND.ni/=nr; IF(fmsg(lerr,   'integer and real arguments lengths mismatch')) RETURN
+
+  !--- Allocate the assembled quantities array
+  nrow = MAX(ns,ni,nr)+1
+  ALLOCATE(d(nrow,ncol), n(ncol))
+
+  !--- Assemble the vectors into a strings array in the order indicated by "pattern"
+  is =  1; ii = 1; ir = 1
+  DO ic = 1, ncol
+    d(1,ic) = TRIM(titles(ic))
+    SELECT CASE(p(ic:ic))
+      CASE('s'); d(2:nrow,ic) =          s(:,is)     ; is = is + 1
+      CASE('i'); d(2:nrow,ic) =  int2str(i(:,ii)    ); ii = ii + 1
+      CASE('r'); d(2:nrow,ic) = real2str(r(:,ir),rFm); ir = ir + 1
+    END SELECT
+    n(ic) = 0; DO ir=1, nrow; n(ic)=MAX(n(ic), LEN_TRIM(d(ir,ic))); END DO
+    IF(needQuotes(d(2,ic)) .AND. ic/=1) n(ic) = n(ic) + 2 !--- For quotes, using second line only
+  END DO
+
+  !--- Display the strings array as a table
+  DO ir = 1, nrow
+    row = ''; sp = '   '; IF(TRIM(d(ir,1)) /= '') sp = ' = '
+    DO ic = 1, ncol
+      el = d(ir,ic); IF(ic /= 1) el = addQuotes_1(el)
+      row = row//TRIM(el)//REPEAT(' ',n(ic)-LEN_TRIM(el))//sp
+      sp = '   '; IF(ic /= ncol-1) CYCLE
+      IF(TRIM(d(MIN(ir+1,nrow),1)) /= '' .AND. (ir /= nrow .OR. .NOT.la)) sp = ' , '
+    END DO
+    WRITE(unt,'(a)')TRIM(row)
+  END DO
+
+  !--- End of section
+  IF(la) THEN
+    WRITE(unt,'(a)')'/'
+    WRITE(unt,'(a)')
+  END IF
+
+END FUNCTION dispNameList
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+LOGICAL FUNCTION dispOutliers_1(ll, a, n, err_msg, nam, subn, nmax, unit) RESULT(lerr)
+! Display outliers list in tables
+! If "nam" is supplied, it means the last index is for tracers => one table each tracer for rank > 2.
+  LOGICAL,                    INTENT(IN)  :: ll(:)                   !--- Linearized mask of outliers
+  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)  :: nmax, unit              !--- Maximum number of lines to display (default: all)
+  CHARACTER(LEN=256),         ALLOCATABLE :: ttl(:)
+  LOGICAL,                    ALLOCATABLE :: m(:)
+  INTEGER,                    ALLOCATABLE :: ki(:), kj(:)
+  INTEGER                                 :: i, j, k, rk, rk1, ib, ie, itr, nm, unt, nmx, nv
+  CHARACTER(LEN=256)                      :: mes, sub, fm='(f12.9)', v, s
+  CHARACTER(LEN=256),         ALLOCATABLE :: vnm(:)
+
+  lerr = ANY(ll); IF(.NOT.lerr) RETURN                               !--- No outliers -> finished
+
+  mes = 'outliers found'; IF(PRESENT(err_msg)) mes = err_msg         !--- Error message
+  vnm = ['a'];            IF(PRESENT(nam ))    vnm = nam             !--- Variables names
+  sub = 'dispOutliers';   IF(PRESENT(subn))    sub = subn            !--- Calling subroutine name
+  nmx = SIZE(a);          IF(PRESENT(nmax))    nmx = MIN(nmx,nmax)   !--- Maximum number of lines to print
+  unt = lunout;           IF(PRESENT(unit))    unt = unit            !--- Unit to print messages
+
+  rk = SIZE(n); nv = SIZE(vnm)
+  IF(test(fmsg(nv /= 1 .AND. nv /= n(rk), 'In "'//TRIM(sub)//'": SIZE(nam) /= 1 or =last "n" element'    , unt),lerr)) RETURN
+  IF(test(fmsg(SIZE(a) /= SIZE(ll),       'In "'//TRIM(sub)//'": "ll" and "a" sizes mismatch'            , unt),lerr)) RETURN
+  IF(test(fmsg(SIZE(a) /= PRODUCT(n),     'In "'//TRIM(sub)//'": profile "n" does not match "a" and "ll"', unt),lerr)) RETURN
+
+  WRITE(unt,*)'Outliers detected by '//TRIM(sub)//': '//TRIM(mes)
+
+  !--- SCALAR CASE: single value to display
+  IF(rk==1.AND.n(1)==1) THEN
+    IF(ll(1)) WRITE(unt,'(a," = ",f12.9)')TRIM(nam(1)),a(1); RETURN
+  END IF
+
+  rk1 = rk; IF(nv==1) rk1 = rk-1                                    !--- Rank of each displayed table
+  SELECT CASE(rk1)                                                  !--- Indices list
+    CASE(1,2); ki = [ (i,i=1,n(1)) ]
+    CASE(3);   ki = [((i,i=1,n(1)),j=1,n(2))]; kj = [((j,i=1,n(1)),j=1,n(2))]
+    CASE DEFAULT; WRITE(unt,*)'Sorry: routine "dispOutliers" is limited to rank 3'; RETURN
+  END SELECT
+
+  !--- VECTOR CASE:  table " name | value " (known names)  /  )  /  " i | a(i) " (unknown names)
+  IF(rk==1) THEN
+    ALLOCATE(ttl(2)); ttl(2) = TRIM(vnm(1))//'(i)'; ttl(1) = 'i'
+    IF(nv == 1) lerr = dispTable('sr', ttl,               s=cat(PACK(nam,ll)), r=cat(PACK(a,ll)), rFmt=fm, nmax=nmax)
+    IF(nv /= 1) lerr = dispTable('ir', ['name ','value'], i=cat(PACK(ki,m)),   r=cat(PACK(a,ll)), rFmt=fm, nmax=nmax)
+    CALL msg(lerr,'In '//TRIM(sub)//": can't display outliers table", unt)
+    RETURN
+  END IF
+
+  !--- OTHER CASES: one table for each tracer (last index)
+  ttl = [(ACHAR(k), k = 105, 104+rk)]                                !--- Titles list ('i', 'j', 'k'...)
+  s = strStack( ttl(1:rk-1) )                                        !--- Tracer name dummy indexes: (i, j, k, ...
+
+  DO itr=1,n(rk)
+    nm = PRODUCT(n(1:rk-1))                                          !--- number of elements per tracer
+    ie = itr*nm; ib = ie-nm+1; m=ll(ib:ie)                           !--- section bounds for tracer "itr" ; outlayers mask
+    IF(.NOT.ANY(m)) CYCLE                                            !--- no outlayers for tracer "itr"
+    v = TRIM(vnm(MIN(itr,SIZE(vnm))))//'('//TRIM(s)                  !--- "<name>("
+    IF(nv == 1) ttl(rk) = TRIM(v)//','//int2str(itr)//')'            !--- "<name>(i,j,itr)" (single name)
+    IF(nv /= 1) ttl(rk) = TRIM(v)//')'                               !--- "<nam(itr)>(i,j)" (one name each table/itr index)
+    IF(rk==2) lerr = dispTable('ir',  ttl, i=cat(PACK(ki,m)),            r=cat(PACK(a(ib:ie),m)), rFmt=fm, nmax=nmax)
+    IF(rk==3) lerr = dispTable('iir', ttl, i=cat(PACK(ki,m),PACK(kj,m)), r=cat(PACK(a(ib:ie),m)), rFmt=fm, nmax=nmax)
+    CALL msg(lerr,'In '//TRIM(sub)//": can't display outliers table", unt)
+    IF(lerr) THEN; CALL msg("Can't display outliers table"); RETURN; END IF
+  END DO
+END FUNCTION dispOutliers_1
+!==============================================================================================================================
+LOGICAL FUNCTION dispOutliers_2(ll, a, n, err_msg, nam, subn, nmax, unit) RESULT(lerr)
+! Display outliers list in tables
+! If "nam" is supplied and, it means the last index is for tracers => one table each tracer for rank > 2.
+  LOGICAL,                    INTENT(IN)  :: ll(:)                   !--- Linearized mask of outliers
+  REAL,                       INTENT(IN)  ::  a(:,:)                 !--- Linearized arrays of values stacked along 2nd dim.
+  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)  :: nmax, unit              !--- Maximum number of lines to display (default: all)
+
+  CHARACTER(LEN=256)                      :: mes, sub, fm='(f12.9)', prf
+  CHARACTER(LEN=256),         ALLOCATABLE :: ttl(:), vnm(:)
+  LOGICAL,                    ALLOCATABLE :: m(:)
+  INTEGER,                    ALLOCATABLE :: ki(:), kj(:), kl(:)
+  INTEGER                                 :: i, j, k, rk, ib, ie, itr, nm, nv, unt, nmx, rk1
+  REAL,                       ALLOCATABLE :: val(:,:)
+
+  lerr = ANY(ll); IF(.NOT.lerr) RETURN                                    !--- No outliers -> finished
+  rk = SIZE(n); nv = SIZE(a,2)
+  mes = 'outliers found';         IF(PRESENT(err_msg)) mes = err_msg      !--- Error message
+  vnm = [(ACHAR(k+96),k=1,nv)];   IF(PRESENT(nam ))    vnm = nam          !--- Variables names
+  sub = 'dispOutliers';           IF(PRESENT(subn))    sub = subn         !--- Calling subroutine name
+  nmx = SIZE(a);                  IF(PRESENT(nmax))    nmx = MIN(nmx,nmax)!--- Maximum number of lines to print
+  unt = lunout;                   IF(PRESENT(unit))    unt = unit         !--- Unit to print messages
+  lerr = SIZE(vnm) /= nv;         IF(fmsg(lerr, 'In "dispOutlayers_2": SIZE(nam) /= SIZE(a,2)'  ,unt)) RETURN
+  lerr = SIZE(a,1) /= SIZE(ll);   IF(fmsg(lerr,'In '//TRIM(sub)//': "ll" and "a" sizes mismatch',unt)) RETURN
+  lerr = SIZE(a,1) /= PRODUCT(n); IF(fmsg(lerr,'In '//TRIM(sub)//': profile "n" does not match "a" and "ll"',unt)) RETURN
+
+  SELECT CASE(rk1)                                                   !--- 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)) ]
+    CASE(2); ki = [ ((i,i=1,n(1)),j=1,n(2))];           kj = [ ((j,i=1,n(1)),j=1,n(2))]
+    CASE(3); ki = [(((i,i=1,n(1)),j=1,n(2)),k=1,n(3))]; kj = [(((j,i=1,n(1)),j=1,n(2)),k=1,n(3))]
+             kl = [(((k,i=1,n(1)),j=1,n(2)),k=1,n(3))]
+    CASE DEFAULT; WRITE(unt,*)'Sorry: routine "dispOutliers_2" is limited to rank 3'; RETURN
+  END SELECT
+
+  ttl = [(ACHAR(k), k = 105, 104+rk), vnm]                           !--- Titles list ('i', 'j', 'k'...'var1', 'var2', ...)
+  prf = REPEAT('i',rk)//REPEAT('r',nv)                               !--- Profile
+  ALLOCATE(val(COUNT(ll),nv)); DO k=1, nv; val(:,k) = PACK(a(:,k),ll); END DO
+  IF(rk == 1) lerr = dispTable(prf, ttl, i = cat(PACK(ki,ll)),             r = val, rFmt=fm, nmax=nmax)
+  IF(rk == 2) lerr = dispTable(prf, ttl, i = cat(PACK(ki,ll),PACK(kj,ll)), r = val, rFmt=fm, nmax=nmax)
+  IF(rk == 3) lerr = dispTable(prf, ttl, i = cat(PACK(ki,ll),PACK(kj,ll),PACK(kl,ll)), &
+                     r = val, rFmt=fm, nmax=nmax)
+  CALL msg(lerr,'In '//TRIM(sub)//": can't display outliers table", unt)
+END FUNCTION dispOutliers_2
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== Reduce an algebrical expression (basic operations and parenthesis) to a single number (string format) ====================
+!==============================================================================================================================
+LOGICAL FUNCTION reduceExpr_1(str, val) RESULT(out)
+  CHARACTER(LEN=*),   INTENT(IN)  :: str
+  CHARACTER(LEN=256), INTENT(OUT) :: val
+  CHARACTER(LEN=1024)              :: s, vv
+  CHARACTER(LEN=256) :: v
+!  CHARACTER(LEN=:),    ALLOCATABLE :: v
+  CHARACTER(LEN=1024), ALLOCATABLE :: vl(:)
+  INTEGER,             ALLOCATABLE :: ip(:)
+  INTEGER :: nl, nn, i, j, im, ix
+  LOGICAL :: ll
+  modname = 'reduceExpr_1'
+  s = str
+
+  !--- Check wether the parenthesis are correctly formed
+  ll = strCount(s,'(',nl)
+  ll = strCount(s,')',nn)
+  out = nl == nn
+  IF(fmsg(.NOT.out, 'Mismatching number of opening and closing parenthesis: '//TRIM(s))) RETURN
+  nl = 2*nl-1
+
+  !--- Build vectors ; vl: parenthesis-separated elements ; ip: parenthesis types (1: opening, 2: closing)
+  ALLOCATE(ip(nl-1),vl(nl))
+  j = 1; im = 1
+  DO i = 1, LEN_TRIM(str)
+    ix = INDEX('()', str(i:i))
+    IF(ix == 0) CYCLE
+    ip(j) = ix
+    vl(j) = str(im:i-1)
+    j = j + 1; im = i + 1
+  END DO
+  vl(j) = str(im:LEN_TRIM(str))
+
+  !--- Search for each opening/closing parenthesis pair
+  DO WHILE(nl > 1)
+    i = 1; DO WHILE(ip(i) /= 1 .OR. ip(i+1) /= 2); i = i + 1; END DO !IF(i > SIZE(ip)+1) EXIT;END DO
+    out = reduceExpr_basic(vl(i+1), v); IF(.NOT. out) RETURN
+    v = TRIM(vl(i))//TRIM(v); IF(i+2<=nl) v=TRIM(v)//TRIM(vl(i+2))
+    vv = v//REPEAT(' ',768)
+    IF(i == 1) THEN;         ip =  ip(3:nl-1);                vl = [            vv, vl(4  :nl)]
+    ELSE IF(i == nl-1) THEN; ip =  ip(1:nl-2);                vl = [vl(1:nl-3), vv            ]
+    ELSE;                    ip = [ip(1: i-1), ip(i+2:nl-1)]; vl = [vl(1: i-1), vv, vl(i+3:nl)]; END IF
+    nl = SIZE(vl)
+  END DO
+  out = reduceExpr_basic(vl(1), val)
+END FUNCTION reduceExpr_1
+
+
+!==============================================================================================================================
+!=== Reduce a simple algebrical expression (basic operations, no parenthesis) to a single number (string format) ==============
+!==============================================================================================================================
+LOGICAL FUNCTION reduceExpr_basic(str, val) RESULT(out)
+  CHARACTER(LEN=*),   INTENT(IN)  :: str
+  CHARACTER(LEN=256), INTENT(OUT) :: val
+  DOUBLE PRECISION,   ALLOCATABLE :: vl(:)
+  INTEGER,            ALLOCATABLE :: id(:)
+  CHARACTER(LEN=256), ALLOCATABLE :: ky(:)
+  CHARACTER(LEN=1),   ALLOCATABLE :: op(:)
+
+  CHARACTER(LEN=1024) :: s
+  DOUBLE PRECISION :: v, vm, vp
+  INTEGER      :: i, ni, io
+  LOGICAL :: ll
+
+  modname = 'reduceExpr_basic'
+  out = .TRUE.
+  print*,TRIM(str),is_numeric(str)
+!  val =''
+print*,"coincoin 1",len(str),len_trim(str)
+  IF(is_numeric(str)) THEN; val = TRIM(str); RETURN; END IF
+print*,"coincoin 2"
+  op = ['^','/','*','+','-']                                                   !--- List of recognized operations
+  s = str
+  ll = strParse_m(s, op, ky, .TRUE., id = id)                                  !--- Parse the values
+  vl = str2dble(ky)                                                            !--- Conversion to doubles
+  out = ALL(vl < HUGE(1.d0))
+  IF(fmsg(.NOT.out,'Some values are non-numeric in: '//TRIM(s))) RETURN        !--- Non-numerical values found
+  DO io = 1, SIZE(op)                                                          !--- Loop on known operators (order matters !)
+    DO i = SIZE(id), 1, -1                                                     !--- Loop on found operators
+      ni = SIZE(id)
+      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(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
+      IF(i == ni) THEN; id =  id(1:ni-1);     ELSE; id = [id(1:i-1),    id(i+1:ni  )]; END IF
+    END DO
+  END DO
+  val = dble2str(vl(1))
+
+END FUNCTION reduceExpr_basic
+!==============================================================================================================================
+
+!==============================================================================================================================
+FUNCTION reduceExpr_m(str, val) RESULT(out)
+  CHARACTER(LEN=*),    INTENT(IN)  :: str(:)
+  CHARACTER(LEN=256),  INTENT(OUT) :: val(SIZE(str))
+  LOGICAL                          :: out(SIZE(str))
+  INTEGER :: i
+  out = [(reduceExpr_1(str(i), val(i)), i=1, SIZE(str))]
+END FUNCTION reduceExpr_m
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== Check whether a string is a number or not ================================================================================
+!==============================================================================================================================
+ELEMENTAL LOGICAL FUNCTION is_numeric(str) RESULT(out)
+  CHARACTER(LEN=*), INTENT(IN) :: str
+  REAL    :: x
+  INTEGER :: e
+  CHARACTER(LEN=12) :: fmt
+  WRITE(fmt,'("(f",i0,".0)")') LEN_TRIM(str)
+  READ(str,fmt,IOSTAT=e) x
+  out = e==0
+END FUNCTION is_numeric
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== Convert a string into a logical/integer integer or an integer/real into a string =========================================
+!==============================================================================================================================
+ELEMENTAL LOGICAL FUNCTION str2bool(str) RESULT(out)
+  CHARACTER(LEN=*), INTENT(IN) :: str
+  INTEGER :: ierr
+  READ(str,*,IOSTAT=ierr) out
+  IF(ierr==0) RETURN
+  out = ANY(['t     ','true  ','.true.','y     ','yes   ']==strLower(str))
+END FUNCTION str2bool
+!==============================================================================================================================
+ELEMENTAL INTEGER FUNCTION str2int(str) RESULT(out)
+  CHARACTER(LEN=*), INTENT(IN) :: str
+  INTEGER :: ierr
+  READ(str,*,IOSTAT=ierr) out
+  IF(ierr/=0) out = -HUGE(1)
+END FUNCTION str2int
+!==============================================================================================================================
+ELEMENTAL REAL FUNCTION str2real(str) RESULT(out)
+  CHARACTER(LEN=*), INTENT(IN) :: str
+  INTEGER :: ierr
+  READ(str,*,IOSTAT=ierr) out
+  IF(ierr/=0) out = -HUGE(1.)
+END FUNCTION str2real
+!==============================================================================================================================
+ELEMENTAL DOUBLE PRECISION FUNCTION str2dble(str) RESULT(out)
+  CHARACTER(LEN=*), INTENT(IN) :: str
+  INTEGER :: ierr
+  READ(str,*,IOSTAT=ierr) out
+  IF(ierr/=0) out = -HUGE(1.d0)
+END FUNCTION str2dble
+!==============================================================================================================================
+ELEMENTAL CHARACTER(LEN=256) FUNCTION bool2str(b) RESULT(out)
+  LOGICAL, INTENT(IN) :: b
+  WRITE(out,*)b
+  out = ADJUSTL(out)
+END FUNCTION bool2str
+!==============================================================================================================================
+ELEMENTAL CHARACTER(LEN=256) FUNCTION int2str(i) RESULT(out)
+  INTEGER, INTENT(IN) :: i
+  WRITE(out,*)i
+  out = ADJUSTL(out)
+END FUNCTION int2str
+!==============================================================================================================================
+ELEMENTAL CHARACTER(LEN=256) FUNCTION real2str(r,fmt) RESULT(out)
+  REAL,                       INTENT(IN) :: r
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: fmt
+  IF(     PRESENT(fmt)) WRITE(out,fmt)r
+  IF(.NOT.PRESENT(fmt)) WRITE(out, * )r
+  out = ADJUSTL(out)
+END FUNCTION real2str
+!==============================================================================================================================
+ELEMENTAL CHARACTER(LEN=256) FUNCTION dble2str(d,fmt) RESULT(out)
+  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
+  out = ADJUSTL(out)
+END FUNCTION dble2str
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+FUNCTION addQuotes_1(s) RESULT(out)
+  CHARACTER(LEN=*), INTENT(IN)  :: s
+  CHARACTER(LEN=:), ALLOCATABLE :: out
+  IF(needQuotes(s)) THEN; out = "'"//TRIM(s)//"'"; ELSE; out = s; END IF
+END FUNCTION addQuotes_1
+!==============================================================================================================================
+FUNCTION addQuotes_m(s) RESULT(out)
+  CHARACTER(LEN=*), INTENT(IN)  :: s(:)
+  CHARACTER(LEN=:), ALLOCATABLE :: out(:)
+  INTEGER :: k, n
+  n = MAXVAL(LEN_TRIM(s), MASK=.TRUE.)
+  ALLOCATE(CHARACTER(LEN=n) :: out(SIZE(s)))
+  DO k=1,SIZE(s)
+    IF(needQuotes(s(k))) THEN; out(k) = "'"//TRIM(s(k))//"'"; ELSE; out(k) = s(k); END IF
+  END DO
+END FUNCTION addQuotes_m
+!==============================================================================================================================
+ELEMENTAL LOGICAL FUNCTION needQuotes(s) RESULT(out)
+  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
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== DISPLAY "<message>: the following <items> are <reason>" FOLLOWED BY THE LIST OF <str> FOR WHICH <lerr>==T. ===============
+!==============================================================================================================================
+LOGICAL FUNCTION checkList(str, lerr, message, items, reason, nmax) RESULT(out)
+! Purpose: Messages in case a list contains wrong elements (indicated by lerr boolean vector).
+! Note:    Return value "out" is .TRUE. if there are errors (ie at least one element of "lerr" is TRUE).
+  CHARACTER(LEN=*),   INTENT(IN)  :: str(:)
+  LOGICAL,            INTENT(IN)  :: lerr(:)
+  CHARACTER(LEN=*),   INTENT(IN)  :: message, items, reason
+  INTEGER,  OPTIONAL, INTENT(IN)  :: nmax
+  CHARACTER(LEN=256), ALLOCATABLE :: s(:)
+  INTEGER :: i, nmx
+  nmx = 256; IF(PRESENT(nmax)) nmx=nmax
+  out = ANY(lerr); IF(.NOT.out) RETURN
+  CALL msg(TRIM(message)//': the following '//TRIM(items)//' are '//TRIM(reason)//':')
+  s = strStackm(PACK(str, MASK=lerr), ', ',nmx)
+  DO i=1,SIZE(s,DIM=1); CALL msg(s(i)); END DO
+END FUNCTION checkList
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== Remove comment in line "str", ie all the characters from the first "#" sign found in "str". ==============================
+!==============================================================================================================================
+SUBROUTINE removeComment(str)
+  CHARACTER(LEN=*), INTENT(INOUT) :: str
+  INTEGER :: ix
+  ix = INDEX(str,'# '); IF(ix /= 0) str = str(1:ix-1)//REPEAT(' ',LEN(str)-ix+1)
+END SUBROUTINE removeComment
+!==============================================================================================================================
+
+
+
+END MODULE strings_mod
Index: LMDZ6/branches/LMDZ-tracers/libf/misc/trac_types_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ-tracers/libf/misc/trac_types_mod.F90	(revision 3852)
+++ LMDZ6/branches/LMDZ-tracers/libf/misc/trac_types_mod.F90	(revision 3852)
@@ -0,0 +1,52 @@
+MODULE trac_types_mod
+
+!=== TRACERS DESCRIPTOR DERIVED TYPE AND ASSOCIATED ROUTINES INTERFACES =======================================================
+  PRIVATE
+  PUBLIC :: tra, iso, db, kys
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE kys                                                           !=== TYPE FOR A SET OF KEYS ASSOCIATED TO AN ELEMENT
+    CHARACTER(LEN=256)              :: name                          !--- Tracer name
+    CHARACTER(LEN=256), ALLOCATABLE :: key(:)                        !--- Keys string list
+    CHARACTER(LEN=256), ALLOCATABLE :: val(:)                        !--- Corresponding values string list
+  END TYPE kys
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE tra                                                           !=== TYPE FOR SINGLE TRACER
+    CHARACTER(LEN=256)   :: name = ''                                !--- Name
+    CHARACTER(LEN=256)   :: nam1 = ''                                !--- Generation 1 ancestor name
+    CHARACTER(LEN=256)   :: prnt = ''                                !--- Parent name
+    CHARACTER(LEN=256)   :: lnam = ''                                !--- Long name (with adv. scheme)
+    CHARACTER(LEN=256)   :: type = 'tracer'                          !--- Type (so far: 'tracer'/'tag')
+    CHARACTER(LEN=256)   :: phas = 'g'                               !--- Phase ('g'as/'l'iquid/'s'olid)
+    CHARACTER(LEN=256)   :: comp                                     !--- Coma-separated list of components (Ex: lmdz,inca)
+    INTEGER              :: iadv = 10                                !--- Advection scheme used
+    INTEGER              :: igen = 1                                 !--- Generation number (>=1)
+    INTEGER              :: itr  = 0                                 !--- Index in tr_seri (0: not in physics)
+    INTEGER              :: iprnt = 0                                !--- Parent index
+    INTEGER, ALLOCATABLE :: idesc(:)                                 !--- Descendants index (in growing generation order)
+    INTEGER              :: ndesc = 0                                !--- Number of descendants (all generations)
+    INTEGER              :: nchld = 0                                !--- Number of childs    (first generation)
+    INTEGER              :: iso_igr = 0                              !--- Isotopes group index in isotopes(:)
+    INTEGER              :: iso_num = 0                              !--- Isotope  name  index in isotopes(iso_igr)%trac(:)
+    INTEGER              :: iso_zon = 0                              !--- Isotope  zone  index in isotopes(iso_igr)%zone(:)
+    INTEGER              :: iso_pha = 0                              !--- Isotope  phase index in isotopes(iso_igr)%phas
+    TYPE(kys)            :: keys                                     !--- <key>=<val> pairs vector
+  END TYPE tra
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE db                                                            !=== TYPE FOR TRACERS SECTION
+    CHARACTER(LEN=256)     :: name                                   !--- Section name
+    TYPE(tra), ALLOCATABLE :: trac(:)                                !--- Tracers descriptors
+  END TYPE db
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE iso                                                           !=== TYPE FOR ISOTOPES FAMILY DESCENDING ON TRACER "prnt"
+    CHARACTER(LEN=256)              :: prnt                          !--- Isotopes family name (parent tracer name ; ex: H2O)
+    LOGICAL                         :: check=.FALSE.                 !--- Triggering of the checking routines
+    TYPE(kys),          ALLOCATABLE :: keys(:)                       !--- Isotopes keys/values pairs list (length: niso)
+    CHARACTER(LEN=256), ALLOCATABLE :: trac(:)                       !--- Isotopes + tagging tracers list (length: nitr)
+    CHARACTER(LEN=256), ALLOCATABLE :: zone(:)                       !--- Geographic tagging zones names list
+    CHARACTER(LEN=256)              :: phas = 'g'                    !--- Phases list: [g][l][s]
+    INTEGER                         :: niso=0, nzon=0, nitr=0, npha=0!--- Number of isotopes, zones, total isotopes and phases
+    INTEGER,            ALLOCATABLE :: iTraPha(:,:)   ! (iqiso)      !--- Idx in "trac(1:niso)" = f(name(1:nitr)),phas)
+    INTEGER,            ALLOCATABLE :: iZonIso(:,:)   ! (index_trac) !--- Idx in "trac(1:nitr)" = f(zone, name(1:niso))
+  END TYPE iso
+
+END MODULE trac_types_mod
