Index: LMDZ6/branches/Ocean_skin/libf/misc/readTracFiles_mod.f90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/misc/readTracFiles_mod.f90	(revision 4368)
+++ LMDZ6/branches/Ocean_skin/libf/misc/readTracFiles_mod.f90	(revision 4368)
@@ -0,0 +1,1987 @@
+MODULE readTracFiles_mod
+
+  USE strings_mod,    ONLY: msg, find, get_in, str2int, dispTable, testFile, strReduce,  strFind, strStack, strHead,  &
+       test, removeComment, cat, fmsg, maxlen, int2str, checkList, strParse, strReplace, strTail, strCount, strIdx, reduceExpr
+
+  IMPLICIT NONE
+
+  PRIVATE
+
+  PUBLIC :: maxlen                                              !--- PARAMETER FOR CASUAL STRING LENGTH
+  PUBLIC :: tracers                                             !--- TRACERS  DESCRIPTION DATABASE
+  PUBLIC :: trac_type, setGeneration, indexUpdate               !--- TRACERS  DESCRIPTION ASSOCIATED TOOLS
+  PUBLIC :: testTracersFiles, readTracersFiles                  !--- TRACERS FILES READING ROUTINES
+  PUBLIC :: getKey, fGetKey, fGetKeys, addKey, setDirectKeys    !--- TOOLS TO GET/SET KEYS FROM/TO  tracers & isotopes
+  PUBLIC :: getKeysDBase,    setKeysDBase                       !--- TOOLS TO GET/SET THE DATABASE (tracers & isotopes)
+
+  PUBLIC :: addPhase, getiPhase,  old_phases, phases_sep, &     !--- FUNCTIONS RELATED TO THE PHASES
+   nphases, delPhase, getPhase, known_phases, phases_names      !--- + ASSOCIATED VARIABLES
+
+  PUBLIC :: oldH2OIso, newH2OIso, old2newH2O, new2oldH2O        !--- H2O ISOTOPES BACKWARD COMPATIBILITY (OLD traceur.def)
+  PUBLIC :: oldHNO3,   newHNO3                                  !--- HNO3 REPRO   BACKWARD COMPATIBILITY (OLD start.nc)
+
+  PUBLIC :: tran0, idxAncestor, ancestor                        !--- GENERATION 0 TRACER + TOOLS FOR GENERATIONS
+
+  !=== FOR ISOTOPES: GENERAL
+  PUBLIC :: isot_type, readIsotopesFile, isoSelect              !--- ISOTOPES DESCRIPTION TYPE + READING ROUTINE
+  PUBLIC :: ixIso, nbIso                                        !--- INDEX OF SELECTED ISOTOPES CLASS + NUMBER OF CLASSES
+
+  !=== FOR ISOTOPES: H2O FAMILY ONLY
+  PUBLIC :: iH2O
+
+  !=== FOR ISOTOPES: DEPENDING ON THE SELECTED ISOTOPES CLASS
+  PUBLIC :: isotope, isoKeys                                    !--- SELECTED ISOTOPES DATABASE + ASSOCIATED KEYS
+  PUBLIC :: isoName, isoZone, isoPhas                           !--- ISOTOPES AND TAGGING ZONES NAMES AND PHASES
+  PUBLIC :: niso,    nzone,   nphas,   ntiso                    !---  " " NUMBERS + ISOTOPES AND TAGGING TRACERS NUMBERS
+  PUBLIC :: itZonIso                                            !--- Idx IN isoName(1:niso) = f(tagging idx, isotope idx)
+  PUBLIC :: iqIsoPha                                            !--- Idx IN qx(1:nqtot)     = f(isotope idx,   phase idx)
+  PUBLIC :: isoCheck                                            !--- FLAG TO RUN ISOTOPES CHECKING ROUTINES
+
+  PUBLIC :: maxTableWidth
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE :: keys_type                                        !=== TYPE FOR A SET OF KEYS ASSOCIATED TO AN ELEMENT
+    CHARACTER(LEN=maxlen)              :: name             !--- Tracer name
+    CHARACTER(LEN=maxlen), ALLOCATABLE :: key(:)           !--- Keys string list
+    CHARACTER(LEN=maxlen), ALLOCATABLE :: val(:)           !--- Corresponding values string list
+  END TYPE keys_type
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE :: trac_type                                        !=== TYPE FOR A SINGLE TRACER NAMED "name"
+    CHARACTER(LEN=maxlen) :: name        = ''              !--- Name of the tracer
+    CHARACTER(LEN=maxlen) :: gen0Name    = ''              !--- First generation ancestor name
+    CHARACTER(LEN=maxlen) :: parent      = ''              !--- Parent name
+    CHARACTER(LEN=maxlen) :: longName    = ''              !--- Long name (with advection scheme suffix)
+    CHARACTER(LEN=maxlen) :: type        = 'tracer'        !--- Type  (so far: 'tracer' / 'tag')
+    CHARACTER(LEN=maxlen) :: phase       = 'g'             !--- Phase ('g'as / 'l'iquid / 's'olid)
+    CHARACTER(LEN=maxlen) :: component   = ''              !--- Coma-separated list of components (Ex: lmdz,inca)
+    INTEGER               :: iGeneration = -1              !--- Generation number (>=0)
+    INTEGER               :: iqParent    = 0               !--- Parent index
+    INTEGER,  ALLOCATABLE :: iqDescen(:)                   !--- Descendants index (in growing generation order)
+    INTEGER               :: nqDescen    = 0               !--- Number of descendants (all generations)
+    INTEGER               :: nqChildren  = 0               !--- Number of children  (first generation)
+    TYPE(keys_type)       :: keys                          !--- <key>=<val> pairs vector
+    INTEGER               :: iadv        = 10              !--- Advection scheme used
+    LOGICAL               :: isAdvected  = .FALSE.         !--- "true" tracers: iadv > 0.   COUNT(isAdvected )=nqtrue
+    LOGICAL               :: isInPhysics = .TRUE.          !--- "true" tracers: in tr_seri. COUNT(isInPhysics)=nqtottr
+    INTEGER               :: iso_iGroup  = 0               !--- Isotopes group index in isotopes(:)
+    INTEGER               :: iso_iName   = 0               !--- Isotope  name  index in isotopes(iso_iGroup)%trac(:)
+    INTEGER               :: iso_iZone   = 0               !--- Isotope  zone  index in isotopes(iso_iGroup)%zone(:)
+    INTEGER               :: iso_iPhase  = 0               !--- Isotope  phase index in isotopes(iso_iGroup)%phase
+  END TYPE trac_type
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE :: isot_type                                        !=== TYPE FOR AN ISOTOPES FAMILY DESCENDING ON TRACER "parent"
+    CHARACTER(LEN=maxlen)              :: parent           !--- Isotopes family name (parent tracer name ; ex: H2O)
+    LOGICAL                            :: check=.FALSE.    !--- Triggering of the checking routines
+    TYPE(keys_type),       ALLOCATABLE :: keys(:)          !--- Isotopes keys/values pairs list     (length: niso)
+    CHARACTER(LEN=maxlen), ALLOCATABLE :: trac(:)          !--- Isotopes + tagging tracers list     (length: ntiso)
+    CHARACTER(LEN=maxlen), ALLOCATABLE :: zone(:)          !--- Geographic tagging zones names list (length: nzone)
+    CHARACTER(LEN=maxlen)              :: phase = 'g'      !--- Phases list: [g][l][s]              (length: nphas)
+    INTEGER                            :: niso  = 0        !--- Number of isotopes, excluding tagging tracers
+    INTEGER                            :: nzone = 0        !--- Number of geographic tagging zones
+    INTEGER                            :: ntiso = 0        !--- Number of isotopes, including tagging tracers
+    INTEGER                            :: nphas = 0        !--- Number phases
+    INTEGER,               ALLOCATABLE :: iqIsoPha(:,:)    !--- Idx in "tracers(1:nqtot)" = f(name(1:ntiso)),phas)
+                                                           !---        "iqIsoPha" former name: "iqiso"
+    INTEGER,               ALLOCATABLE :: itZonIso(:,:)    !--- Idx in "trac(1:ntiso)" = f(zone, name(1:niso))
+                                                           !---        "itZonIso" former name: "index_trac"
+  END TYPE isot_type
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE :: dataBase_type                                         !=== TYPE FOR TRACERS SECTION
+    CHARACTER(LEN=maxlen)  :: name                              !--- Section name
+    TYPE(trac_type), ALLOCATABLE :: trac(:)                     !--- Tracers descriptors
+  END TYPE dataBase_type
+!------------------------------------------------------------------------------------------------------------------------------
+  INTERFACE getKey
+    MODULE PROCEDURE getKeyByName_s1, getKeyByName_s1m, getKeyByName_sm, getKey_sm, &
+                     getKeyByName_i1, getKeyByName_i1m, getKeyByName_im, getKey_im, &
+                     getKeyByName_r1, getKeyByName_r1m, getKeyByName_rm, getKey_rm, &
+                     getKeyByName_l1, getKeyByName_l1m, getKeyByName_lm, getKey_lm
+  END INTERFACE getKey
+!------------------------------------------------------------------------------------------------------------------------------
+  INTERFACE    isoSelect;  MODULE PROCEDURE  isoSelectByIndex,  isoSelectByName; END INTERFACE isoSelect
+  INTERFACE  old2newH2O;   MODULE PROCEDURE  old2newH2O_1,  old2newH2O_m;        END INTERFACE old2newH2O
+  INTERFACE  new2oldH2O;   MODULE PROCEDURE  new2oldH2O_1,  new2oldH2O_m;        END INTERFACE new2oldH2O
+  INTERFACE fGetKey;       MODULE PROCEDURE fgetKeyIdx_s1, fgetKeyNam_s1;        END INTERFACE fGetKey
+  INTERFACE tracersSubset; MODULE PROCEDURE trSubset_Indx, trSubset_Name, trSubset_gen0Name; END INTERFACE tracersSubset
+  INTERFACE idxAncestor;   MODULE PROCEDURE idxAncestor_1, idxAncestor_m, idxAncestor_mt;    END INTERFACE idxAncestor
+  INTERFACE    ancestor;   MODULE PROCEDURE    ancestor_1,    ancestor_m,    ancestor_mt;    END INTERFACE    ancestor
+  INTERFACE      addKey;   MODULE PROCEDURE      addKey_1; END INTERFACE addKey!,      addKey_m,     addKey_mm;     END INTERFACE addKey
+  INTERFACE    addPhase;   MODULE PROCEDURE   addPhase_s1,   addPhase_sm,   addPhase_i1,   addPhase_im; END INTERFACE addPhase
+!------------------------------------------------------------------------------------------------------------------------------
+
+  !=== MAIN DATABASE: files sections descriptors
+  TYPE(dataBase_type), SAVE, ALLOCATABLE, TARGET :: dBase(:)
+
+  !--- SOME PARAMETERS THAT ARE NOT LIKELY TO CHANGE OFTEN
+  CHARACTER(LEN=maxlen), SAVE      :: tran0        = 'air'      !--- Default transporting fluid
+  CHARACTER(LEN=maxlen), PARAMETER :: old_phases   = 'vlir'     !--- Old phases for water (no separator)
+  CHARACTER(LEN=maxlen), PARAMETER :: known_phases = 'glsr'     !--- Known phases initials
+  INTEGER, PARAMETER :: nphases = LEN_TRIM(known_phases)        !--- Number of phases
+  CHARACTER(LEN=maxlen), SAVE      :: phases_names(nphases) &   !--- Known phases names
+                                             = ['gaseous', 'liquid ', 'solid  ', 'cloud  ']
+  CHARACTER(LEN=1), SAVE :: phases_sep  =  '_'                  !--- Phase separator
+  LOGICAL,          SAVE :: tracs_merge = .TRUE.                !--- Merge/stack tracers lists
+  LOGICAL,          SAVE :: lSortByGen  = .TRUE.                !--- Sort by growing generation
+  CHARACTER(LEN=maxlen), SAVE :: isoFile = 'isotopes_params.def'!--- Name of the isotopes parameters file
+
+  !--- CORRESPONDANCE BETWEEN OLD AND NEW WATER NAMES
+  CHARACTER(LEN=maxlen), SAVE :: oldH2OIso(5) = ['eau',     'HDO',     'O18',     'O17',     'HTO'    ]
+  CHARACTER(LEN=maxlen), SAVE :: newH2OIso(5) = ['H2[16]O', 'H[2]HO ', 'H2[18]O', 'H2[17]O', 'H[3]HO ']
+
+  !--- CORRESPONDANCE BETWEEN OLD AND NEW HNO3 RELATED SPECIES NAMES
+  CHARACTER(LEN=maxlen), SAVE ::   oldHNO3(2) = ['HNO3_g ', 'HNO3   ']
+  CHARACTER(LEN=maxlen), SAVE ::   newHNO3(2) = ['HNO3   ', 'HNO3tot']
+
+  !=== TRACERS AND ISOTOPES DESCRIPTORS, USED BY getKey
+  TYPE(trac_type), ALLOCATABLE, TARGET, SAVE ::  tracers(:)
+  TYPE(isot_type), ALLOCATABLE, TARGET, SAVE :: isotopes(:)
+
+  !=== ALIASES OF VARIABLES FROM SELECTED ISOTOPES FAMILY EMBEDDED IN "isotope" (isotopes(ixIso))
+  TYPE(isot_type),         SAVE, POINTER :: isotope             !--- CURRENTLY SELECTED ISOTOPES FAMILY DESCRIPTOR
+  INTEGER,                 SAVE          :: ixIso, iH2O         !--- Index of the selected isotopes family and H2O family
+  INTEGER,                 SAVE          :: nbIso               !--- Number of isotopes classes
+  LOGICAL,                 SAVE          :: isoCheck            !--- Flag to trigger the checking routines
+  TYPE(keys_type),         SAVE, POINTER :: isoKeys(:)          !--- ONE SET OF KEYS FOR EACH ISOTOPE (LISTED IN isoName)
+  CHARACTER(LEN=maxlen),   SAVE, POINTER :: isoName(:),   &     !--- ISOTOPES NAMES FOR THE CURRENTLY SELECTED FAMILY
+                                            isoZone(:),   &     !--- TAGGING ZONES  FOR THE CURRENTLY SELECTED FAMILY
+                                            isoPhas             !--- USED PHASES    FOR THE CURRENTLY SELECTED FAMILY
+  INTEGER,                 SAVE          ::  niso, nzone, &     !--- NUMBER OF ISOTOPES, TAGGING ZONES AND PHASES
+                                            nphas, ntiso        !--- NUMBER OF PHASES AND ISOTOPES + ISOTOPIC TAGGING TRACERS
+  INTEGER,                 SAVE, POINTER ::itZonIso(:,:), &     !--- INDEX IN "isoTrac" AS f(tagging zone idx,  isotope idx)
+                                           iqIsoPha(:,:)        !--- INDEX IN "qx"      AS f(isotopic tracer idx, phase idx)
+
+  INTEGER,    PARAMETER :: maxTableWidth = 192                  !--- Maximum width of a table displayed with "dispTable"
+  CHARACTER(LEN=maxlen) :: modname
+
+CONTAINS
+
+!==============================================================================================================================
+!==============================================================================================================================
+!=== READ ONE OR SEVERAL TRACER FILES AND FILL A "tr" TRACERS DESCRIPTOR DERIVED TYPE.
+!=== THE RETURNED VALUE fType DEPENDS ON WHAT IS FOUND:
+!===  0: NO ADEQUATE FILE FOUND ; DEFAULT VALUES MUST BE USED
+!===  1: AN "OLD STYLE" TRACERS FILE "traceur.def":
+!===    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 keys_type) 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 (longName, iadv, isAdvected for now).
+!     * 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: iqParent, iqDescen).
+!     * If you need to convert a %key/%val pair into a direct-access key, add the corresponding line in "setDirectKeys".
+!==============================================================================================================================
+LOGICAL FUNCTION readTracersFiles(type_trac, fTyp, lRepr) RESULT(lerr)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=*),  INTENT(IN)  :: type_trac                        !--- List of components used
+  INTEGER, OPTIONAL, INTENT(OUT) :: fTyp                             !--- Type of input file found
+  LOGICAL, OPTIONAL, INTENT(IN)  :: lRepr                            !--- Activate the HNNO3 exceptions for REPROBUS
+  CHARACTER(LEN=maxlen),  ALLOCATABLE :: s(:), sections(:), trac_files(:)
+  CHARACTER(LEN=maxlen) :: str, fname, tname, pname, cname
+  INTEGER               :: nsec, ierr, it, ntrac, ns, ip, ix, fType
+  LOGICAL :: lRep
+  TYPE(keys_type), POINTER :: k
+!------------------------------------------------------------------------------------------------------------------------------
+  lerr = .FALSE.
+  modname = 'readTracersFiles'
+  IF(.NOT.ALLOCATED(dBase)) ALLOCATE(dBase(0))
+  lRep=.FALSE.; IF(PRESENT(lRepr)) lRep = lRepr
+
+  !--- Required sections + corresponding files names (new style single section case) for tests
+  IF(test(testTracersFiles(modname, type_trac, fType, .TRUE., trac_files, sections), lerr)) RETURN
+  IF(PRESENT(fTyp)) fTyp = fType
+  nsec = SIZE(sections)
+
+  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+  SELECT CASE(fType)                         !--- Set %name, %genOName, %parent, %type, %phase, %component, %iGeneration, %keys
+  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+    CASE(1)                                                          !=== 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('Invalid format for "'//TRIM(fname)//'"', modname, ierr /= 0), lerr)) RETURN
+
+      !--- READ THE REMAINING LINES: <hadv> <vadv> <tracer> [<transporting fluid>]
+      IF(ALLOCATED(tracers)) DEALLOCATE(tracers)
+      ALLOCATE(tracers(ntrac))
+      DO it=1,ntrac                                                  !=== READ RAW DATA: loop on the line/tracer number
+        READ(90,'(a)',IOSTAT=ierr) str
+        IF(test(fmsg('Invalid format for "' //TRIM(fname)//'"', modname, ierr>0), lerr)) RETURN
+        IF(test(fmsg('Not enough lines in "'//TRIM(fname)//'"', modname, ierr<0), lerr)) RETURN
+        lerr = strParse(str, ' ', s, ns)
+        CALL msg('This file is for air tracers only',           modname, ns == 3 .AND. it == 1)
+        CALL msg('This files specifies the transporting fluid', modname, ns == 4 .AND. it == 1)
+        k => tracers(it)%keys
+
+        !=== NAME OF THE TRACER
+        tname = old2newH2O(s(3), ip)
+        ix = strIdx(oldHNO3, s(3))
+        IF(ix /= 0 .AND. lRep) tname = newHNO3(ix)                   !--- Exception for HNO3 (REPROBUS ONLY)
+        tracers(it)%name = tname                                     !--- Set %name
+        CALL addKey_1('name', tname, k)                              !--- Set the name of the tracer
+        tracers(it)%keys%name = tname                                !--- Copy tracers names in keys components
+
+        !=== NAME OF THE COMPONENT
+        cname = type_trac                                            !--- Name of the model component
+        IF(ANY([(addPhase('H2O', ip), ip = 1, nphases)] == tname)) cname = 'lmdz'
+        tracers(it)%component = cname                                !--- Set %component
+        CALL addKey_1('component', cname, k)                         !--- Set the name of the model component
+
+        !=== NAME OF THE PARENT
+        pname = tran0                                                !--- Default name: default transporting fluid (air)
+        IF(ns == 4) THEN
+          pname = old2newH2O(s(4))
+          ix = strIdx(oldHNO3, s(4))
+          IF(ix /= 0 .AND. lRep) pname = newHNO3(ix)                 !--- Exception for HNO3 (REPROBUS ONLY)
+        END IF
+        tracers(it)%parent = pname                                   !--- Set %parent
+        CALL addKey_1('parent', pname, k)
+
+        !=== PHASE AND ADVECTION SCHEMES NUMBERS
+        tracers(it)%phase = known_phases(ip:ip)                      !--- Set %phase:  tracer phase (default: "g"azeous)
+        CALL addKey_1('phase', known_phases(ip:ip), k)               !--- Set the phase  of the tracer (default: "g"azeous)
+        CALL addKey_1('hadv', s(1),  k)                              !--- Set the horizontal advection schemes number
+        CALL addKey_1('vadv', s(2),  k)                              !--- Set the vertical   advection schemes number
+      END DO
+      CLOSE(90)
+      IF(test(setGeneration(tracers), lerr)) RETURN                  !--- Set %iGeneration and %gen0Name
+      WHERE(tracers%iGeneration == 2) tracers(:)%type = 'tag'        !--- Set %type:        'tracer' or 'tag'
+      DO it=1,ntrac
+        CALL addKey_1('type', tracers(it)%type, tracers(it)%keys)    !--- Set the type of tracer
+      END DO
+      IF(test(checkTracers(tracers, fname, fname), lerr)) RETURN     !--- Detect orphans and check phases
+      IF(test(checkUnique (tracers, fname, fname), lerr)) RETURN     !--- Detect repeated tracers
+      CALL sortTracers    (tracers)                                  !--- Sort the tracers
+    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+    CASE(2); IF(test(feedDBase(["tracer.def"], [type_trac], modname), lerr)) RETURN !=== SINGLE   FILE, MULTIPLE SECTIONS
+    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+    CASE(3); IF(test(feedDBase(  trac_files  ,  sections,   modname), lerr)) RETURN !=== MULTIPLE FILES, SINGLE  SECTION
+  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+  END SELECT
+  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+  IF(ALL([2,3] /= fType)) RETURN
+
+  IF(nsec  == 1) THEN; 
+    tracers = dBase(1)%trac
+  ELSE IF(tracs_merge) THEN
+    CALL msg('The multiple required sections will be MERGED.',    modname)
+    IF(test(mergeTracers(dBase, tracers), lerr)) RETURN
+  ELSE
+    CALL msg('The multiple required sections will be CUMULATED.', modname)
+    IF(test(cumulTracers(dBase, tracers), lerr)) RETURN
+  END IF
+  CALL setDirectKeys(tracers)                                        !--- Set %iqParent, %iqDescen, %nqDescen, %nqChildren
+END FUNCTION readTracersFiles
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+LOGICAL FUNCTION testTracersFiles(modname, type_trac, fType, lDisp, tracf, sects) RESULT(lerr)
+  CHARACTER(LEN=*),                             INTENT(IN)  :: modname, type_trac
+  INTEGER,                                      INTENT(OUT) :: fType
+  LOGICAL,                                      INTENT(IN)  :: lDisp
+  CHARACTER(LEN=maxlen), ALLOCATABLE, OPTIONAL, INTENT(OUT) :: tracf(:), sects(:)
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: trac_files(:), sections(:)
+  LOGICAL, ALLOCATABLE :: ll(:)
+  INTEGER :: is, nsec
+
+  !--- PARSE "type_trac" LIST AND DETERMINE THE TRACERS FILES NAMES (FOR CASE 3: MULTIPLE FILES, SINNGLE SECTION PER FILE)
+  IF(test(strParse(type_trac, '|', sections,  n=nsec), lerr)) RETURN !--- Parse "type_trac" list
+  IF(PRESENT(sects)) sects = sections
+  ALLOCATE(trac_files(nsec)); DO is=1, nsec; trac_files(is) = 'tracer_'//TRIM(sections(is))//'.def'; END DO
+  IF(PRESENT(tracf)) tracf = trac_files
+
+  nsec = SIZE(trac_files, DIM=1)
+  ll = .NOT.testFile(trac_files)
+  fType = 0
+  IF(.NOT.testFile('traceur.def') .AND. nsec==1) fType = 1           !--- OLD STYLE FILE
+  IF(.NOT.testFile('tracer.def'))                fType = 2           !--- NEW STYLE ; SINGLE  FILE, SEVERAL SECTIONS
+  IF(ALL(ll))                                    fType = 3           !--- NEW STYLE ; SEVERAL FILES, SINGLE SECTION USED
+  IF(.NOT.test(lDisp, lerr)) RETURN                                  !--- NO CHECKING/DISPLAY NEEDED: JUST GET type_trac,fType
+  IF(ANY(ll) .AND. fType/=3) THEN                                    !--- MISSING FILES
+    IF(test(checkList(trac_files, .NOT.ll, 'Failed reading tracers description', 'files', 'missing'), lerr)) RETURN
+  END IF
+
+  !--- CHECK WHETHER type_trac AND FILE TYPE ARE COMPATIBLE
+  IF(test(fmsg('No multiple sections for the old format "traceur.def"', ll = nsec>1 .AND. fType==1), lerr)) RETURN
+
+  !--- TELLS WHAT WAS IS ABOUT TO BE USED
+  CALL msg('Trying to read old-style tracers description file "traceur.def"',                      modname, fType==1)
+  CALL msg('Trying to read the new style multi-sections tracers description file "tracer.def"',    modname, fType==2)
+  CALL msg('Trying to read the new style single section tracers description files "tracer_*.def"', modname, fType==3)
+END FUNCTION testTracersFiles
+!==============================================================================================================================
+
+!==============================================================================================================================
+LOGICAL FUNCTION feedDBase(fnames, snames, modname) RESULT(lerr)
+! Purpose: Read the sections "snames(is)" (pipe-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(:)                         !--- Pipe-deparated list of sections (one list each file)
+  CHARACTER(LEN=*), INTENT(IN)  :: modname                           !--- Calling routine name
+  INTEGER,  ALLOCATABLE :: ndb(:)                                    !--- Number of sections for each file
+  INTEGER,  ALLOCATABLE :: ixf(:)                                    !--- File index for each section of the expanded list
+  CHARACTER(LEN=maxlen) :: 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)                                               !--- Set %name, %keys
+    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                    !--- FILE AND SECTION NAMES
+    lerr = ANY([(dispTraSection('RAW CONTENT OF SECTION "'//TRIM(snm)//'"', snm, modname), idb=1, SIZE(dBase))])
+    IF(test(expandSection(dBase(idb)%trac, snm, fnm), lerr)) RETURN  !--- EXPAND NAMES ;  set %parent, %type, %component
+    IF(test(setGeneration(dBase(idb)%trac),           lerr)) RETURN  !---                 set %iGeneration,   %genOName
+    IF(test(checkTracers (dBase(idb)%trac, snm, fnm), lerr)) RETURN  !--- CHECK ORPHANS AND PHASES
+    IF(test(checkUnique  (dBase(idb)%trac, snm, fnm), lerr)) RETURN  !--- CHECK TRACERS UNIQUENESS
+    CALL expandPhases    (dBase(idb)%trac)                           !--- EXPAND PHASES ; set %phase
+    CALL sortTracers     (dBase(idb)%trac)                           !--- SORT TRACERS
+    lerr = ANY([(dispTraSection('EXPANDED CONTENT OF SECTION "'//TRIM(snm)//'"', snm, modname), idb=1, SIZE(dBase))])
+  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+  END DO
+  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+END FUNCTION feedDBase
+!------------------------------------------------------------------------------------------------------------------------------
+
+!------------------------------------------------------------------------------------------------------------------------------
+LOGICAL FUNCTION readSections(fnam,snam,defName) RESULT(lerr)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=*),           INTENT(IN) :: fnam                     !--- File name
+  CHARACTER(LEN=*),           INTENT(IN) :: snam                     !--- Pipe-separated sections list
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: defName                  !--- Special section (default values) name
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE(dataBase_type),   ALLOCATABLE :: tdb(:)
+  CHARACTER(LEN=maxlen), 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=maxlen), ALLOCATABLE ::  s(:), v(:)
+  TYPE(trac_type),       ALLOCATABLE :: tt(:)
+  TYPE(trac_type)       :: tmp
+  CHARACTER(LEN=1024)   :: str, str2
+  CHARACTER(LEN=maxlen) :: secn
+  INTEGER               :: ierr, n
+!------------------------------------------------------------------------------------------------------------------------------
+  IF(.NOT.ALLOCATED(dBase)) ALLOCATE(dBase(0))
+  OPEN(90, FILE=fnam, FORM='formatted', STATUS='old')
+  DO; str=''
+    DO
+      READ(90,'(a)', IOSTAT=ierr)str2                                !--- Read a full line
+      str=TRIM(str)//' '//TRIM(str2)                                 !--- Append "str" with the current line
+      n=LEN_TRIM(str); IF(n == 0) EXIT                               !--- Empty line (probably end of file)
+      IF(IACHAR(str(n:n))  /= 92) EXIT                               !--- No "\" continuing line symbol found => end of line
+      str = str(1:n-1)                                               !--- Remove the "\" continuing line symbol
+    END DO
+    str = ADJUSTL(str)                                               !--- Remove the front space
+    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,' ', s, n, v)                                !--- Parse <key>=<val> pairs
+      tt = dBase(ndb)%trac(:)
+      tmp%name = s(1); tmp%keys = keys_type(s(1), s(2:n), v(2:n))    !--- Set %name and %keys
+      dBase(ndb)%trac = [tt(:), tmp]
+      DEALLOCATE(tt)
+!      dBase(ndb)%trac = [dBase(ndb)%trac(:), tra(name=s(1), keys=keys_type(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(trac_type), ALLOCATABLE, TARGET, INTENT(INOUT) :: t(:)
+  CHARACTER(LEN=*),                     INTENT(IN)    :: defName
+  INTEGER :: jd, it, k
+  TYPE(keys_type), POINTER :: ky
+  TYPE(trac_type), ALLOCATABLE :: tt(:)
+  jd = strIdx(t(:)%name, defName)
+  IF(jd == 0) RETURN
+  ky => t(jd)%keys
+  DO k = 1, SIZE(ky%key)                                             !--- Loop on the keys of the tracer named "defName"
+!   CALL addKey_m(ky%key(k), ky%val(k), t(:)%keys, .FALSE.)           !--- Add key to all the tracers (no overwriting)
+    DO it = 1, SIZE(t); CALL addKey_1(ky%key(k), ky%val(k), t(it)%keys, .FALSE.); END DO
+  END DO
+  tt = [t(1:jd-1),t(jd+1:SIZE(t))]; CALL MOVE_ALLOC(FROM=tt, TO=t)   !--- Remove the virtual tracer named "defName"
+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(trac_type), ALLOCATABLE, TARGET, INTENT(INOUT) :: t(:)
+  CHARACTER(LEN=*),                     INTENT(IN)    :: defName
+  LOGICAL,                              INTENT(IN)    :: lSubLocal
+  INTEGER :: i0, it, ik
+  TYPE(keys_type), POINTER     :: k0, ky
+  TYPE(trac_type), 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(trac_type), ALLOCATABLE, INTENT(INOUT) :: tr(:)                 !--- Tracer derived type vector
+  CHARACTER(LEN=*),             INTENT(IN)    :: sname
+  CHARACTER(LEN=*), OPTIONAL,   INTENT(IN)    :: fname
+  TYPE(trac_type),       ALLOCATABLE :: ttr(:)
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: ta(:), pa(:)
+  CHARACTER(LEN=maxlen) :: msg1, modname
+  INTEGER :: it, nt, iq, nq, itr, ntr, ipr, npr
+  LOGICAL :: ll
+  modname = 'expandSection'
+  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)%parent    = fgetKey(it, 'parent', tr(:)%keys,  tran0  )
+    tr(it)%type      = fgetKey(it, 'type'  , tr(:)%keys, 'tracer')
+    tr(it)%component = sname
+!   CALL addKey_m('component', sname, tr(:)%keys)
+    DO iq=1,SIZE(tr); CALL addKey_1('component', sname, tr(iq)%keys); END DO
+
+    !--- Determine the number of tracers and parents ; coherence checking
+    ll = strCount(tr(it)%name,   ',', ntr)
+    ll = strCount(tr(it)%parent, ',', 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', modname); RETURN
+    END IF
+    nq = nq + ntr*npr                 
+  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+  END DO
+  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+  ALLOCATE(ttr(nq))
+  iq = 1
+  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+  DO it = 1, nt                                                      !=== EXPAND TRACERS AND PARENTS NAMES LISTS
+  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+    ll = strParse(tr(it)%name,   ',', ta, ntr)                       !--- Number of tracers
+    ll = strParse(tr(it)%parent, ',', pa, npr)                       !--- Number of parents
+    DO ipr=1,npr                                                     !--- Loop on parents list elts
+      DO itr=1,ntr                                                   !--- Loop on tracers list elts
+        ttr(iq)%keys%key  = tr(it)%keys%key
+        ttr(iq)%keys%val  = tr(it)%keys%val
+        ttr(iq)%keys%name = ta(itr)
+        ttr(iq)%name      = TRIM(ta(itr));    CALL addKey_1('name',      ta(itr),          ttr(iq)%keys)
+        ttr(iq)%parent    = TRIM(pa(ipr));    CALL addKey_1('parent',    pa(ipr),          ttr(iq)%keys)
+        ttr(iq)%type      = tr(it)%type;      CALL addKey_1('type',      tr(it)%type,      ttr(iq)%keys)
+        ttr(iq)%component = tr(it)%component; CALL addKey_1('component', tr(it)%component, ttr(iq)%keys)
+        iq = iq+1
+      END DO
+    END DO
+  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+  END DO
+  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+  DEALLOCATE(ta,pa)
+  CALL MOVE_ALLOC(FROM=ttr, TO=tr)
+
+END FUNCTION expandSection
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+LOGICAL FUNCTION setGeneration(tr) RESULT(lerr)
+!------------------------------------------------------------------------------------------------------------------------------
+! Purpose: Determine, for each tracer of "tr(:)":
+!   * %iGeneration: the generation number 
+!   * %gen0Name:    the generation 0 ancestor name
+!          Check also for orphan tracers (tracers not descending on "tran0").
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE(trac_type),     INTENT(INOUT) :: tr(:)                        !--- Tracer derived type vector
+  INTEGER                            :: iq, jq, ig
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: parent(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen) :: modname
+  modname = 'setGeneration'
+  IF(test(fmsg('missing "parent" attribute', modname, getKey('parent', parent, ky=tr(:)%keys)), lerr)) RETURN
+  DO iq = 1, SIZE(tr)
+    jq = iq; ig = 0
+    DO WHILE(parent(jq) /= tran0)
+      jq = strIdx(tr(:)%name, parent(jq))
+      IF(test(fmsg('Orphan tracer "'//TRIM(tr(iq)%name)//'"', modname, jq == 0), lerr)) RETURN
+      ig = ig + 1
+    END DO
+    tr(iq)%gen0Name = tr(jq)%name; CALL addKey_1('gen0Name',    tr(iq)%gen0Name,   tr(iq)%keys)
+    tr(iq)%iGeneration = ig;       CALL addKey_1('iGeneration', TRIM(int2str(ig)), tr(iq)%keys)
+  END DO
+END FUNCTION setGeneration
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+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(trac_type),            INTENT(IN) :: tr(:)                    !--- Tracer derived type vector
+  CHARACTER(LEN=*),           INTENT(IN) :: sname                    !--- Section name
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: fname                    !--- File name
+  CHARACTER(LEN=maxlen) :: mesg
+  CHARACTER(LEN=maxlen) :: bp(SIZE(tr, DIM=1)), pha                  !--- Bad phases list, phases of current tracer
+  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%iGeneration==-1, mesg, 'tracers', 'orphan'), lerr)) RETURN
+
+  !=== CHECK PHASES
+  DO iq=1,nq; IF(tr(iq)%iGeneration/=0) CYCLE                        !--- Generation O 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%iGeneration==0 .AND. bp/='', mesg, 'tracers phases', 'unknown')
+END FUNCTION checkTracers
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+LOGICAL FUNCTION checkUnique(tr, sname, fname) RESULT(lerr)
+!------------------------------------------------------------------------------------------------------------------------------
+! Purpose: Make sure that tracers are not repeated.
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE(trac_type),            INTENT(IN) :: tr(:)                    !--- Tracer derived type vector
+  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=maxlen) :: 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(tr(iq)%type == 'tag') CYCLE                         !--- Tags can be repeated
+    tnam = TRIM(tr(iq)%name)
+    ll = tr(:)%name==TRIM(tnam)                                      !--- Mask for current tracer name
+    IF(COUNT(ll)==1 ) CYCLE                                          !--- Tracer is not repeated
+    IF(tr(iq)%iGeneration>0) THEN
+      tdup(iq) = tnam                                                !--- gen>0: 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". Phases are not repeated for a tracer, thanks to "checkUnique".
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE(trac_type), ALLOCATABLE, INTENT(INOUT) :: tr(:)               !--- Tracer derived type vector
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE(trac_type), ALLOCATABLE :: ttr(:)
+  INTEGER,   ALLOCATABLE ::  i0(:)
+  CHARACTER(LEN=maxlen)  :: nam, pha, tname
+  CHARACTER(LEN=1) :: p
+  INTEGER :: ip, np, iq, jq, nq, it, nt, nc, i, n
+  LOGICAL :: lTag, lExt
+!------------------------------------------------------------------------------------------------------------------------------
+  nq = SIZE(tr, DIM=1)
+  nt = 0
+  DO iq = 1, nq                                                      !--- GET THE NUMBER OF TRACERS
+    IF(tr(iq)%iGeneration /= 0) CYCLE                                !--- Only deal with generation 0 tracers
+    nc = COUNT(tr(:)%gen0Name==tr(iq)%name .AND. tr%iGeneration/=0)  !--- Number of children of tr(iq)
+    tr(iq)%phase = fgetKey(iq, 'phases', tr(:)%keys)                 !--- Phases list        of tr(iq)
+    np = LEN_TRIM(tr(iq)%phase)                                      !--- Number of phases   of tr(iq)
+    nt = nt + (1+nc) * np                                            !--- Number of tracers after expansion
+  END DO
+  ALLOCATE(ttr(nt))                                                  !--- Version  of "tr" after phases expansion
+  it = 1                                                             !--- Current "ttr(:)" index
+  DO iq = 1, nq                                                      !--- Loop on "tr(:)" indexes
+    lTag = tr(iq)%type=='tag'                                        !--- Current tracer is a tag
+    i0 = strFind(tr(:)%name, TRIM(tr(iq)%gen0Name), n)               !--- Indexes of first generation ancestor copies
+    np = SUM([( LEN_TRIM(tr(i0(i))%phase),i=1,n )], 1)               !--- Number of phases for current tracer tr(iq)
+    lExt = np>1                                                      !--- Phase suffix only required if phases number is > 1
+    IF(lTag) lExt = lExt .AND. tr(iq)%iGeneration>0                  !--- No phase suffix for generation 0 tags
+    DO i=1,n                                                         !=== LOOP ON GENERATION 0 ANCESTORS
+      jq = i0(i)                                                     !--- tr(jq): ith tracer with same gen 0 ancestor as tr(iq)
+      IF(tr(iq)%iGeneration==0) jq=iq                                !--- Generation 0: count the current tracer phases only
+      pha = tr(jq)%phase                                             !--- Phases list for tr(jq)
+      DO ip = 1, LEN_TRIM(pha)                                       !=== LOOP ON PHASES LISTS
+        p = pha(ip:ip)
+        tname = TRIM(tr(iq)%name); nam = tname                       !--- Tracer name (regular case)
+        IF(lTag) nam = TRIM(tr(iq)%parent)                           !--- Parent name (tagging case)
+        IF(lExt) nam = addPhase(nam, p )                             !--- Phase extension needed
+        IF(lTag) nam = TRIM(nam)//'_'//TRIM(tname)                   !--- <parent>_<name> for tags
+        ttr(it) = tr(iq)                                             !--- Same <key>=<val> pairs
+        ttr(it)%name      = TRIM(nam)                                !--- Name with possibly phase suffix
+        ttr(it)%keys%name = TRIM(nam)                                !--- Name inside the keys decriptor
+        ttr(it)%phase     = p                                        !--- Single phase entry
+        CALL addKey_1('name', nam, ttr(it)%keys)
+        CALL addKey_1('phase', p,  ttr(it)%keys)
+        IF(lExt .AND. tr(iq)%iGeneration>0) THEN
+          ttr(it)%parent   = addPhase(tr(iq)%parent,   p)
+          ttr(it)%gen0Name = addPhase(tr(iq)%gen0Name, p)
+          CALL addKey_1('parent',   ttr(it)%parent,   ttr(it)%keys)
+          CALL addKey_1('gen0Name', ttr(it)%gen0Name, ttr(it)%keys)
+        END IF
+        it = it+1
+      END DO
+      IF(tr(iq)%iGeneration==0) EXIT                                 !--- Break phase loop for gen 0
+    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:
+!  * Put water at the beginning of the vector, in the "known_phases" order.
+!  * lGrowGen == T: in ascending generations numbers.
+!  * lGrowGen == F: tracer + its children sorted by growing generation, one after the other.
+!   TO BE ADDED IF NECESSARY: HIGHER MOMENTS AT THE END
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE(trac_type), ALLOCATABLE, INTENT(INOUT) :: tr(:)         !--- Tracer derived type vector
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE(trac_type), ALLOCATABLE        :: tr2(:)
+  INTEGER,         ALLOCATABLE        :: iy(:), iz(:)
+  INTEGER                             :: ig, ng, iq, jq, ip, nq, n, ix(SIZE(tr)), k
+!  tr2 is introduced in order to cope with a bug in gfortran 4.8.5 compiler
+!------------------------------------------------------------------------------------------------------------------------------
+  nq = SIZE(tr)
+  DO ip = nphases, 1, -1
+    iq = strIdx(tr(:)%name, addPhase('H2O', ip))
+    IF(iq == 0) CYCLE
+    tr2 = tr(:)
+    tr = [tr2(iq), tr2(1:iq-1), tr2(iq+1:nq)]
+  END DO
+  IF(lSortByGen) THEN
+    iq = 1
+    ng = MAXVAL(tr(:)%iGeneration, MASK=.TRUE., DIM=1)               !--- Number of generations
+    DO ig = 0, ng                                                    !--- Loop on generations
+      iy = PACK([(k, k=1, nq)], MASK=tr(:)%iGeneration==ig)          !--- Generation ig tracers indexes
+      n = SIZE(iy)
+      ix(iq:iq+n-1) = iy                                             !--- Stack growing generations idxs
+      iq = iq + n
+    END DO
+  ELSE
+    iq = 1
+    DO jq = 1, nq                                                    !--- Loop on generation 0 tracers
+      IF(tr(jq)%iGeneration /= 0) CYCLE                              !--- Skip generations /= 0
+      ix(iq) = jq                                                    !--- Generation 0 ancestor index first
+      iq = iq + 1                                                    !--- Next "iq" for next generations tracers
+      iy = strFind(tr(:)%gen0Name, TRIM(tr(jq)%name))                !--- Indexes of "tr(jq)" children in "tr(:)"
+      ng = MAXVAL(tr(iy)%iGeneration, MASK=.TRUE., DIM=1)            !--- Number of generations of the "tr(jq)" family
+      DO ig = 1, ng                                                  !--- Loop   on generations of the "tr(jq)" family
+        iz = find(tr(iy)%iGeneration, 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(dataBase_type),  TARGET, INTENT(IN)  :: sections(:)
+  TYPE(trac_type), ALLOCATABLE, INTENT(OUT) ::       tr(:)
+  TYPE(trac_type), POINTER ::   t1(:),   t2(:)
+  INTEGER,     ALLOCATABLE :: ixct(:), ixck(:)
+  INTEGER :: is, k1, k2, nk2, i1, i2, nt2
+  CHARACTER(LEN=maxlen) :: s1, v1, v2, tnam, knam, modname
+  modname = 'mergeTracers'
+  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)//'":', modname)
+    CALL msg(t1(PACK(ixct, MASK = ixct/=0))%name, modname, nmax=128) !--- Display duplicates (the 128 first at most)
+    !--------------------------------------------------------------------------------------------------------------------------
+    DO i2=1,nt2; tnam = TRIM(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('Parent name'//TRIM(s1), modname, t1(i1)%parent      /= t2(i2)%parent),      lerr)) RETURN
+      IF(test(fmsg('Type'       //TRIM(s1), modname, t1(i1)%type        /= t2(i2)%type),        lerr)) RETURN
+      IF(test(fmsg('Generation' //TRIM(s1), modname, t1(i1)%iGeneration /= t2(i2)%iGeneration), 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)%component = TRIM(tr(i1)%component)//','//TRIM(tr(i2)%component)
+
+      !--- 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), modname)
+      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), modname)
+      END DO
+      !------------------------------------------------------------------------------------------------------------------------
+    END DO
+    !--------------------------------------------------------------------------------------------------------------------------
+  END DO
+  CALL sortTracers(tr)
+
+END FUNCTION mergeTracers
+!==============================================================================================================================
+
+!==============================================================================================================================
+LOGICAL FUNCTION cumulTracers(sections, tr) RESULT(lerr)
+  TYPE(dataBase_type),  TARGET, INTENT(IN)  :: sections(:)
+  TYPE(trac_type), ALLOCATABLE, INTENT(OUT) ::       tr(:)
+  TYPE(trac_type), POINTER     :: t(:)
+  INTEGER,   ALLOCATABLE :: nt(:)
+  CHARACTER(LEN=maxlen)  :: 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
+  !----------------------------------------------------------------------------------------------------------------------------
+    t => sections(is)%trac(:)
+    !--------------------------------------------------------------------------------------------------------------------------
+    DO iq=1, nt(is)                                                  !=== LOOP ON TRACERS
+    !--------------------------------------------------------------------------------------------------------------------------
+      tnam = TRIM(t(iq)%name)                                        !--- Original name
+      IF(COUNT(t%name == tnam) == 1) CYCLE                           !--- Current tracer is not duplicated: finished
+      tnam_new = TRIM(tnam)//'_'//TRIM(sections(is)%name)            !--- Same with section extension
+      nq = SUM(nt(1:is-1))                                           !--- Number of tracers in previous sections
+      ns = nt(is)                                                    !--- Number of tracers in the current section
+      tr(iq + nq)%name = TRIM(tnam_new)                              !--- Modify tracer name
+      WHERE(tr(1+nq:ns+nq)%parent==tnam) tr(1+nq:ns+nq)%parent=tnam_new  !--- Modify parent name
+    !--------------------------------------------------------------------------------------------------------------------------
+    END DO
+  !----------------------------------------------------------------------------------------------------------------------------
+  END DO
+  !----------------------------------------------------------------------------------------------------------------------------
+  CALL sortTracers(tr)
+END FUNCTION cumulTracers
+!==============================================================================================================================
+
+!==============================================================================================================================
+SUBROUTINE setDirectKeys(tr)
+  TYPE(trac_type), INTENT(INOUT) :: tr(:)
+
+  !--- Update %iqParent, %iqDescen, %nqDescen, %nqChildren
+  CALL indexUpdate(tr)
+
+  !--- Extract some direct-access keys
+!  DO iq = 1, SIZE(tr)
+!    tr(iq)%keys%<key> = getKey_prv(it, "<key>", tr%keys, <default_value> )
+!  END DO
+END SUBROUTINE setDirectKeys
+!==============================================================================================================================
+
+!==============================================================================================================================
+LOGICAL FUNCTION dispTraSection(message, sname, modname) RESULT(lerr)
+  CHARACTER(LEN=*), INTENT(IN) :: message, sname, modname
+  INTEGER :: idb, iq, nq
+  INTEGER, ALLOCATABLE :: hadv(:), vadv(:)
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: phas(:), prnt(:)
+  TYPE(trac_type), POINTER :: tm(:)
+  lerr = .FALSE.
+  idb = strIdx(dBase(:)%name, sname); IF(idb == 0) RETURN
+  tm => dBase(idb)%trac
+  nq = SIZE(tm)
+  !--- BEWARE ! Can't use the "getKeyByName" functions yet.
+  !             Names must first include the phases for tracers defined on multiple lines.
+  hadv = str2int(fgetKeys('hadv',  tm(:)%keys, '10'))
+  vadv = str2int(fgetKeys('vadv',  tm(:)%keys, '10'))
+  prnt =         fgetKeys('parent',tm(:)%keys,  '' )
+  IF(getKey('phases', phas, ky=tm(:)%keys)) phas = fGetKeys('phase', tm(:)%keys, 'g')
+  CALL msg(TRIM(message)//':', modname)
+  IF(ALL(prnt == 'air')) THEN
+    IF(test(dispTable('iiiss',   ['iq    ','hadv  ','vadv  ','name  ','phase '],                   cat(tm%name,       phas),  &
+                 cat([(iq, iq=1, nq)], hadv, vadv),                 nColMax=maxTableWidth, nHead=2, sub=modname), lerr)) RETURN
+  ELSE IF(ALL(tm%iGeneration == -1)) THEN
+    IF(test(dispTable('iiisss', ['iq    ','hadv  ','vadv  ','name  ','parent','phase '],           cat(tm%name, prnt, phas),  &
+                 cat([(iq, iq=1, nq)], hadv, vadv),                 nColMax=maxTableWidth, nHead=2, sub=modname), lerr)) RETURN
+  ELSE
+    IF(test(dispTable('iiissis', ['iq    ','hadv  ','vadv  ','name  ','parent','igen  ','phase '], cat(tm%name, prnt, phas),  &
+                 cat([(iq, iq=1, nq)], hadv, vadv, tm%iGeneration), nColMax=maxTableWidth, nHead=2, sub=modname), lerr)) RETURN
+  END IF
+END FUNCTION dispTraSection
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!== CREATE A SCALAR ALIAS OF THE COMPONENT OF THE TRACERS DESCRIPTOR "t" NAMED "tname" ========================================
+!==============================================================================================================================
+FUNCTION aliasTracer(tname, t) RESULT(out)
+  TYPE(trac_type),         POINTER    :: out
+  CHARACTER(LEN=*),        INTENT(IN) :: tname
+  TYPE(trac_type), 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(trac_type), ALLOCATABLE             ::  out(:)
+  TYPE(trac_type), 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(trac_type), ALLOCATABLE             ::  out(:)
+  TYPE(trac_type), 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_gen0Name(trac,nam) RESULT(out)
+  TYPE(trac_type), ALLOCATABLE             ::  out(:)
+  TYPE(trac_type), ALLOCATABLE, INTENT(IN) :: trac(:)
+  CHARACTER(LEN=*),             INTENT(IN) ::  nam
+  out = trac(strFind(delPhase(trac(:)%gen0Name), nam))
+  CALL indexUpdate(out)
+END FUNCTION trSubset_gen0Name
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== UPDATE THE INDEXES iqParent, iqDescend AND iGeneration IN THE TRACERS DESCRIPTOR LIST "tr" (USEFULL FOR SUBSETS) =========
+!==============================================================================================================================
+SUBROUTINE indexUpdate(tr)
+  TYPE(trac_type), INTENT(INOUT) :: tr(:)
+  INTEGER :: iq, ig, igen, ngen, ix(SIZE(tr))
+  tr(:)%iqParent = strIdx( tr(:)%name, tr(:)%parent )                !--- Parent index
+  DO iq = 1, SIZE(tr); CALL addKey_1('iqParent', int2str(tr(iq)%iqParent), tr(iq)%keys); END DO
+  ngen = MAXVAL(tr(:)%iGeneration, MASK=.TRUE.)
+  DO iq = 1, SIZE(tr)
+    ig = tr(iq)%iGeneration
+    IF(ALLOCATED(tr(iq)%iqDescen)) DEALLOCATE(tr(iq)%iqDescen)
+    ALLOCATE(tr(iq)%iqDescen(0))
+    CALL idxAncestor(tr, ix, ig)                                     !--- Ancestor of generation "ng" for each tr
+    DO igen = ig+1, ngen
+      tr(iq)%iqDescen = [tr(iq)%iqDescen, find(ix==iq .AND. tr%iGeneration==igen)]
+      tr(iq)%nqDescen = SIZE(tr(iq)%iqDescen)
+      IF(igen == ig+1) THEN
+        tr(iq)%nqChildren = tr(iq)%nqDescen
+        CALL addKey_1('nqChildren', int2str(tr(iq)%nqChildren), tr(iq)%keys)
+      END IF
+    END DO
+    CALL addKey_1('iqDescen', strStack(int2str(tr(iq)%iqDescen)), tr(iq)%keys)
+    CALL addKey_1('nqDescen',          int2str(tr(iq)%nqDescen) , tr(iq)%keys)
+  END DO
+END SUBROUTINE indexUpdate
+!==============================================================================================================================
+  
+  
+!==============================================================================================================================
+!=== READ FILE "fnam" TO APPEND THE "dBase" TRACERS DATABASE WITH AS MUCH SECTIONS AS PARENTS NAMES IN "isot(:)%parent":   ====
+!===  * Each section dBase(i)%name contains the isotopes "dBase(i)%trac(:)" descending on "dBase(i)%name"="iso(i)%parent"  ====
+!===  * For each isotopes class, the <key>=<val> vector of each tracer is moved into the isotopes descriptor "isot"        ====
+!=== NOTES:                                                                                                                ====
+!===  * Most of the "isot" components have been defined in the calling routine (readIsotopes):                             ====
+!===      parent,  nzone, zone(:),  niso, keys(:)%name,  ntiso, trac(:),  nphas, phas,  iqIsoPha(:,:),  itZonPhi(:,:)      ====
+!===  * Same syntax for isotopes file and "tracer.def": a tracers section contains one line for each of its isotopes       ====
+!===  * Each tracers section can contain a "params" virtual isotope line of isotopes parameters default values             ====
+!===  * In case keys are found both in the "params" section and the "*.def" file, the later value is retained              ====
+!===  * On each isotope line, defined keys can be used for other keys defintions (single level depth substitution)         ====
+!===  * The routine gives an error if a required isotope is not available in the database stored in "fnam"                 ====
+!==============================================================================================================================
+LOGICAL FUNCTION readIsotopesFile_prv(fnam, isot) RESULT(lerr)
+  CHARACTER(LEN=*),        INTENT(IN)    :: fnam                     !--- Input file name
+  TYPE(isot_type), TARGET, INTENT(INOUT) :: isot(:)                  !--- Isotopes descriptors (field %parent must be defined!)
+  INTEGER :: is, iis, it, idb, ndb, nb0
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: vals(:)
+  CHARACTER(LEN=maxlen)              :: modname
+  TYPE(trac_type),           POINTER ::   tt(:), t
+  TYPE(dataBase_type),   ALLOCATABLE ::  tdb(:)
+  modname = 'readIsotopesFile'
+
+  !--- THE INPUT FILE MUST BE PRESENT
+  IF(test(fmsg('Missing isotopes parameters file "'//TRIM(fnam)//'"', modname, testFile(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(:)%parent,'|')),lerr)) RETURN !--- Read sections, one each parent tracer
+  ndb = SIZE(dBase, DIM=1)                                           !--- Current database size
+  DO idb = nb0, ndb
+    iis = idb-nb0+1
+
+    !--- 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.)
+
+    tt => dBase(idb)%trac
+
+    !--- REDUCE THE EXPRESSIONS TO OBTAIN SCALARS AND TRANSFER THEM TO THE "isot" ISOTOPES DESCRIPTORS VECTOR
+    DO it = 1, SIZE(dBase(idb)%trac)
+      t => dBase(idb)%trac(it)
+      is = strIdx(isot(iis)%keys(:)%name, t%name)                    !--- Index in "isot(iis)%keys(:)%name" of isotope "t%name"
+      IF(is == 0) CYCLE
+      IF(test(ANY(reduceExpr(t%keys%val, vals)), lerr)) RETURN       !--- Reduce expressions ; detect non-numerical elements
+      isot(iis)%keys(is)%key = t%keys%key
+      isot(iis)%keys(is)%val = vals
+    END DO
+
+    !--- CHECK FOR MISSING ISOTOPES (NO KEYS ALLOCATED)
+    IF(test(checkList(isot(iis)%keys(:)%name, .NOT.[( ALLOCATED(isot(iis)%keys(is)%key), is=1, SIZE(isot(iis)%keys) )], &
+      'Check file "'//TRIM(fnam)//'" in section "'//TRIM(dBase(idb)%name)//'"', 'isotopes', 'missing'), lerr)) RETURN
+  END DO
+
+  !--- 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
+
+  !--- GET THE isoCheck ENTRY FROM THE *.DEF FILES (MIGHT BE CHANGED TO A CLASS-DEPENDANT KEYWORD)
+  CALL get_in('ok_iso_verif', isot(strIdx(isot%parent, 'H2O'))%check, .FALSE.)
+
+  lerr = dispIsotopes()
+
+CONTAINS
+
+!------------------------------------------------------------------------------------------------------------------------------
+LOGICAL FUNCTION dispIsotopes() RESULT(lerr)
+  INTEGER :: ik, nk, ip, it, nt
+  CHARACTER(LEN=maxlen) :: prf
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: ttl(:), val(:,:)
+  CALL msg('Isotopes parameters read from file "'//TRIM(fnam)//'":', modname)
+  DO ip = 1, SIZE(isot)                                              !--- Loop on parents tracers
+    nk = SIZE(isot(ip)%keys(1)%key)                                  !--- Same keys for each isotope
+    nt = SIZE(isot(ip)%keys)                                         !--- Number of isotopes
+    prf = 'i'//REPEAT('s',nk+1)                                      !--- Profile for table printing
+    ALLOCATE(ttl(nk+2), val(nt,nk+1))
+    ttl(1:2) = ['it  ','name']; ttl(3:nk+2) = isot(ip)%keys(1)%key(:)!--- Titles line with keys names
+    val(:,1) = isot(ip)%keys(:)%name                                 !--- Values table 1st column: isotopes names  
+    DO ik = 1, nk
+      DO it = 1, nt
+        val(it,ik+1) = isot(ip)%keys(it)%val(ik)                     !--- Other columns: keys values
+      END DO
+    END DO
+    IF(test(fmsg('Problem with the table content', modname, dispTable(prf, ttl, val, &
+            cat([(it,it=1,nt)]), rFmt='(EN8.4)', nColMax=maxTableWidth, nHead=2, sub=modname)), lerr)) RETURN
+    DEALLOCATE(ttl, val)
+  END DO        
+END FUNCTION dispIsotopes
+!------------------------------------------------------------------------------------------------------------------------------
+
+END FUNCTION readIsotopesFile_prv
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== IF ISOTOPES (2ND GENERATION TRACERS) ARE DETECTED:                                                                     ===
+!===    * COMPUTE MOST OF THE RELATED QUANTITIES ("isot" COMPONENTS).                                                       ===
+!===    * COMPUTE FEW ISOTOPES-DEDICATED "trac" COMPONENTS                                                                  ===
+!===    * CALL readIsotopesFile_prv TO GET PHYSICAL QUANTITIES (<key>=<val> PAIRS)                                          ===
+!===      NOTE: THIS IS DONE HERE (IN A ROUTINE CALLED BY THE DYNAMIC), BECAUSE THE DYNAMIC NEEDS FEW PHYSICAL PARAMETERS.  ===
+!==============================================================================================================================
+LOGICAL FUNCTION readIsotopesFile(iNames) RESULT(lerr)
+  CHARACTER(LEN=maxlen), TARGET, OPTIONAL, INTENT(IN)  :: iNames(:)
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: p(:), str(:)                 !--- Temporary storage
+  CHARACTER(LEN=maxlen) :: iName, modname
+  CHARACTER(LEN=1)   :: ph                                           !--- Phase
+  INTEGER :: ic, ip, iq, it, iz
+  LOGICAL, ALLOCATABLE :: ll(:)                                      !--- Mask
+  TYPE(trac_type), POINTER   ::  t(:), t1
+  TYPE(isot_type), POINTER   ::  i
+  lerr = .FALSE.
+  modname = 'readIsotopesFile'
+
+  t => tracers
+
+  !--- GET FROM "tracers" THE FULL LIST OF AVAILABLE ISOTOPES CLASSES
+  p = PACK(delPhase(t%parent), MASK = t%type=='tracer' .AND. t%iGeneration==1)
+  CALL strReduce(p, nbIso)
+
+  !--- CHECK WHETHER NEEDED ISOTOPES CLASSES "iNames" ARE AVAILABLE OR NOT
+  IF(PRESENT(iNames)) THEN
+    DO it = 1, SIZE(iNames)
+      IF(test(fmsg('No isotopes class "'//TRIM(iNames(it))//'" found among tracers', modname, ALL(p /= iNames(it))), lerr)) RETURN
+    END DO
+    p = iNames; nbIso = SIZE(p)
+  END IF
+  IF(ALLOCATED(isotopes)) DEALLOCATE(isotopes)
+  ALLOCATE(isotopes(nbIso))
+
+  IF(nbIso==0) RETURN                                                !=== NO ISOTOPES: FINISHED
+
+  !--- ISOTOPES RELATED VARIABLES ; NULL OR EMPTY IF NO ISOTOPES
+  isotopes(:)%parent = p
+  DO ic = 1, SIZE(p)                                                 !--- Loop on isotopes classes
+    i => isotopes(ic)
+    iname = i%parent                                                 !--- Current isotopes class name (parent tracer name)
+
+    !=== Isotopes children of tracer "iname": mask, names, number (same for each phase of "iname")
+    ll = t(:)%type=='tracer' .AND. delPhase(t(:)%parent) == iname .AND. t(:)%phase == 'g'
+    str = PACK(delPhase(t(:)%name), MASK = ll)                       !--- Effectively found isotopes of "iname"
+    i%niso = SIZE(str)                                               !--- Number of "effectively found isotopes of "iname"
+    ALLOCATE(i%keys(i%niso))
+    FORALL(it = 1:i%niso) i%keys(it)%name = str(it)
+
+    !=== Geographic tagging tracers descending on tracer "iname": mask, names, number
+    ll = t(:)%type=='tag'    .AND. delPhase(t(:)%gen0Name) == iname .AND. t(:)%iGeneration == 2
+    i%zone = PACK(strTail(t(:)%name,'_'), MASK = ll)                 !--- Tagging zones names  for isotopes category "iname"
+    CALL strReduce(i%zone)
+    i%nzone = SIZE(i%zone)                                           !--- Tagging zones number for isotopes category "iname"
+
+    !=== Geographic tracers of the isotopes children of tracer "iname" (same for each phase of "iname")
+    !    NOTE: One might like to create a similar variable for 2nd generation tagging tracers (tagging the gen1 tracers)
+    str = PACK(delPhase(t(:)%name), MASK=ll)
+    CALL strReduce(str)
+    i%ntiso = i%niso + SIZE(str)                                     !--- Number of isotopes + their geographic tracers [ntiso]
+    ALLOCATE(i%trac(i%ntiso))
+    FORALL(it = 1:i%niso) i%trac(it) = i%keys(it)%name
+    FORALL(it = i%niso+1:i%ntiso) i%trac(it) = str(it-i%niso)
+
+    !=== Phases for tracer "iname"
+    i%phase = ''
+    DO ip = 1, nphases; ph = known_phases(ip:ip); IF(strIdx(t%name,addPhase(iname, ph)) /= 0) i%phase = TRIM(i%phase)//ph; END DO
+    i%nphas = LEN_TRIM(i%phase)                                       !--- Equal to "nqo" for water
+
+    !=== Tables giving the index in a table of effectively found items for each dynamical tracer (1<=iq<=nqtot)
+    DO iq = 1, SIZE(t)
+      t1 => tracers(iq)
+      IF(delPhase(t1%gen0Name)/=iname .OR. t1%iGeneration==0) CYCLE  !--- Only deal with tracers descending on "iname"
+      t1%iso_iGroup = ic                                             !--- Isotopes family       idx in list "isotopes(:)%parent"
+      t1%iso_iName  = strIdx(i%trac, delPhase(strHead(t1%name,'_'))) !--- Current isotope       idx in effective isotopes list
+      t1%iso_iZone  = strIdx(i%zone,          strTail(t1%name,'_') ) !--- Current isotope zone  idx in effective zones    list
+      t1%iso_iPhase =  INDEX(i%phase,TRIM(t1%phase))                 !--- Current isotope phase idx in effective phases   list
+      IF(t1%iGeneration /= 2) t1%iso_iZone = 0                       !--- Skip possible generation 1 tagging tracers
+    END DO
+
+    !=== Table used to get iq (index in dyn array, size nqtot) from the isotope and phase indexes ; the full isotopes list
+    !    (including tagging tracers) is sorted this way:  iso1, iso2, ..., iso1_zone1, iso2_zone1, ..., iso1_zoneN, iso2_zoneN
+    i%iqIsoPha = RESHAPE( [( (strIdx(t%name,  addPhase(i%trac(it),i%phase(ip:ip))),       it=1, i%ntiso), ip=1, i%nphas)], &
+                         [i%ntiso, i%nphas] )
+    !=== Table used to get ix (index in tagging tracers isotopes list, size ntiso) from the zone and isotope indexes
+    i%itZonIso = RESHAPE( [( (strIdx(i%trac(:), TRIM(i%trac(it))//'_'//TRIM(i%zone(iz))), iz=1, i%nzone), it=1, i%niso )], &
+                         [i%nzone, i%niso] )
+  END DO
+
+  !=== READ PHYSICAL PARAMETERS FROM isoFile FILE
+  IF(test(readIsotopesFile_prv(isoFile, isotopes), lerr)) RETURN
+
+  !=== CHECK CONSISTENCY
+  IF(test(testIsotopes(), lerr)) RETURN
+
+  !=== SELECT FIRST ISOTOPES CLASS OR, IF POSSIBLE, WATER CLASS
+  IF(.NOT.test(isoSelect(1, .TRUE.), lerr)) THEN; IF(isotope%parent == 'H2O') iH2O = ixIso; END IF
+
+CONTAINS
+
+!------------------------------------------------------------------------------------------------------------------------------
+LOGICAL FUNCTION testIsotopes() RESULT(lerr)     !--- MAKE SURE MEMBERS OF AN ISOTOPES FAMILY ARE PRESENT IN THE SAME PHASES
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: ix, it, ip, np, iz, nz
+  TYPE(isot_type), POINTER :: i
+  DO ix = 1, nbIso
+    i => isotopes(ix)
+    !--- Check whether each isotope and tagging isotopic tracer is present in the same number of phases
+    DO it = 1, i%ntiso
+      np = SUM([(COUNT(tracers(:)%name == addPhase(i%trac(it), i%phase(ip:ip))), ip=1, i%nphas)])
+      IF(test(fmsg(TRIM(int2str(np))//' phases instead of '//TRIM(int2str(i%nphas))//' for '//TRIM(i%trac(it)), &
+        modname, np /= i%nphas), lerr)) RETURN
+    END DO
+    DO it = 1, i%niso
+      nz = SUM([(COUNT(i%trac == TRIM(i%trac(it))//'_'//i%zone(iz)), iz=1, i%nzone)])
+      IF(test(fmsg(TRIM(int2str(nz))//' tagging zones instead of '//TRIM(int2str(i%nzone))//' for '//TRIM(i%trac(it)), &
+        modname, nz /= i%nzone), lerr)) RETURN
+    END DO
+  END DO
+END FUNCTION testIsotopes
+!------------------------------------------------------------------------------------------------------------------------------
+
+END FUNCTION readIsotopesFile
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== THE ROUTINE isoSelect IS USED TO SWITCH FROM AN ISOTOPE FAMILY TO ANOTHER: ISOTOPES DEPENDENT PARAMETERS ARE UPDATED
+!     Single generic "isoSelect" routine, using the predefined index of the parent (fast version) or its name (first call).
+!==============================================================================================================================
+LOGICAL FUNCTION isoSelectByName(iName, lVerbose) RESULT(lerr)
+   IMPLICIT NONE
+   CHARACTER(LEN=*),  INTENT(IN) :: iName
+   LOGICAL, OPTIONAL, INTENT(IN) :: lVerbose
+   INTEGER :: iIso
+   LOGICAL :: lV
+   lV = .FALSE.; IF(PRESENT(lVerbose)) lV = lVerbose
+   iIso = strIdx(isotopes(:)%parent, iName)
+   IF(test(iIso == 0, lerr)) THEN
+      niso = 0; ntiso = 0; nzone = 0; nphas = 0; isoCheck=.FALSE.
+      CALL msg('no isotope family named "'//TRIM(iName)//'"', ll=lV)
+      RETURN
+   END IF
+   lerr = isoSelectByIndex(iIso, lV)
+END FUNCTION isoSelectByName
+!==============================================================================================================================
+LOGICAL FUNCTION isoSelectByIndex(iIso, lVerbose) RESULT(lerr)
+   IMPLICIT NONE
+   INTEGER,           INTENT(IN) :: iIso
+   LOGICAL, OPTIONAL, INTENT(IN) :: lVerbose
+   LOGICAL :: lV
+   lv = .FALSE.; IF(PRESENT(lVerbose)) lv = lVerbose
+   lerr = .FALSE.
+   IF(iIso == ixIso) RETURN                                          !--- Nothing to do if the index is already OK
+   lerr = iIso<=0 .OR. iIso>SIZE(isotopes)
+   CALL msg('Inconsistent isotopes family index '//TRIM(int2str(iIso))//': should be > 0 and <= '&
+          //TRIM(int2str(SIZE(isotopes)))//'"', ll = lerr .AND. lV)
+   IF(lerr) RETURN
+   ixIso = iIso                                                      !--- Update currently selected family index
+   isotope  => isotopes(ixIso)                                       !--- Select corresponding component
+   isoKeys  => isotope%keys;     niso     = isotope%niso
+   isoName  => isotope%trac;     ntiso    = isotope%ntiso
+   isoZone  => isotope%zone;     nzone    = isotope%nzone
+   isoPhas  => isotope%phase;    nphas    = isotope%nphas
+   itZonIso => isotope%itZonIso; isoCheck = isotope%check
+   iqIsoPha => isotope%iqIsoPha
+END FUNCTION isoSelectByIndex
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== ADD THE <key>=<val> PAIR TO THE "ky[(:)]" KEY[S] DESCRIPTOR[S] OR THE <key>=<val(:)> PAIRS TO THE "ky(:)" KEYS DESCRIPTORS
+!==============================================================================================================================
+SUBROUTINE addKey_1(key, val, ky, lOverWrite)
+  CHARACTER(LEN=*),  INTENT(IN)    :: key, val
+  TYPE(keys_type),   INTENT(INOUT) :: ky
+  LOGICAL, OPTIONAL, INTENT(IN)    :: lOverWrite
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: k(:), v(:)
+  INTEGER :: iky, nky
+  LOGICAL :: lo
+  lo=.TRUE.; IF(PRESENT(lOverWrite)) lo=lOverWrite
+  IF(.NOT.ALLOCATED(ky%key)) THEN
+    ALLOCATE(ky%key(1)); ky%key(1)=key
+    ALLOCATE(ky%val(1)); ky%val(1)=val
+    RETURN
+  END IF
+  iky = 0; IF(ALLOCATED(ky%key)) iky = strIdx(ky%key,key)
+  IF(iky == 0) THEN
+    nky = 0; IF(ALLOCATED(ky%key)) nky = SIZE(ky%key)
+    ALLOCATE(k(nky+1)); k(1:nky) = ky%key; k(nky+1) = key; ky%key = k
+    ALLOCATE(v(nky+1)); v(1:nky) = ky%val; v(nky+1) = val; ky%val = v
+  ELSE IF(lo) THEN
+    ky%key(iky) = key; ky%val(iky) = val
+  END IF
+END SUBROUTINE addKey_1
+!==============================================================================================================================
+SUBROUTINE addKey_m(key, val, ky, lOverWrite)
+  CHARACTER(LEN=*),  INTENT(IN)    :: key, val
+  TYPE(keys_type),   INTENT(INOUT) :: ky(:)
+  LOGICAL, OPTIONAL, INTENT(IN)    :: lOverWrite
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: itr
+  DO itr = 1, SIZE(ky)
+    CALL addKey_1(key, val, ky(itr), lOverWrite)
+  END DO
+END SUBROUTINE addKey_m
+!==============================================================================================================================
+SUBROUTINE addKey_mm(key, val, ky, lOverWrite)
+  CHARACTER(LEN=*),  INTENT(IN)    :: key, val(:)
+  TYPE(keys_type),   INTENT(INOUT) :: ky(:)
+  LOGICAL, OPTIONAL, INTENT(IN)    :: lOverWrite
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: itr
+  DO itr = 1, SIZE(ky); CALL addKey_1(key, val(itr), ky(itr), lOverWrite); END DO
+END SUBROUTINE addKey_mm
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== OVERWRITE THE KEYS OF THE TRACER NAMED "tr0" WITH THE VALUES FOUND IN THE *.def FILES, IF ANY. ===========================
+!==============================================================================================================================
+SUBROUTINE addKeysFromDef(t, tr0)
+  TYPE(trac_type), ALLOCATABLE, INTENT(INOUT) :: t(:)
+  CHARACTER(LEN=*),             INTENT(IN)    :: tr0
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen) :: val
+  INTEGER               :: ik, jd
+  jd = strIdx(t%name, tr0)
+  IF(jd == 0) RETURN
+  DO ik = 1, SIZE(t(jd)%keys%key)
+    CALL get_in(t(jd)%keys%key(ik), val, '*none*')
+    IF(val /= '*none*') CALL addKey_1(t(jd)%keys%key(ik), val, t(jd)%keys, .TRUE.)
+  END DO
+END SUBROUTINE addKeysFromDef
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== REMOVE THE KEYS NAMED "keyn(:)" FROM EITHER THE "itr"th OR ALL THE KEYS DESCRIPTORS OF "ky(:)" ===========================
+!==============================================================================================================================
+SUBROUTINE delKey_1(itr, keyn, ky)
+  INTEGER,          INTENT(IN)    :: itr
+  CHARACTER(LEN=*), INTENT(IN)    :: keyn(:)
+  TYPE(trac_type),  INTENT(INOUT) :: ky(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: k(:), v(:)
+  LOGICAL,               ALLOCATABLE :: ll(:)
+  INTEGER :: iky
+  IF(itr<=0 .OR. itr>SIZE(ky, DIM=1)) RETURN                          !--- Index is out of range
+  ll = [( ALL(keyn/=ky(itr)%keys%key(iky)), iky=1, SIZE(ky(itr)%keys%key) )]
+  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)
+  CHARACTER(LEN=*), INTENT(IN)    :: keyn(:)
+  TYPE(trac_type),  INTENT(INOUT) :: ky(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: iky
+  DO iky = 1, SIZE(ky); CALL delKey_1(iky, keyn, ky); END DO
+END SUBROUTINE delKey
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!================ GET THE VALUE OF A KEY FROM A "keys_type" DERIVED TYPE ; THE RESULT IS THE RETURNED VALUE ===================
+!==============================================================================================================================
+CHARACTER(LEN=maxlen) FUNCTION fgetKeyIdx_s1(itr, keyn, ky, def_val, lerr) RESULT(val)
+  INTEGER,                    INTENT(IN)  :: itr
+  CHARACTER(LEN=*),           INTENT(IN)  :: keyn
+  TYPE(keys_type),            INTENT(IN)  :: ky(:)
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: def_val
+  LOGICAL,          OPTIONAL, INTENT(OUT) :: lerr
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: iky
+  LOGICAL :: ler
+  iky = 0; val = ''
+  IF(.NOT.test(itr <= 0 .OR. itr > SIZE(ky), ler)) iky = strIdx(ky(itr)%key(:), keyn)    !--- Correct index
+  IF(.NOT.test(iky == 0, ler))                     val = ky(itr)%val(iky)                !--- Found key
+  IF(iky == 0) THEN
+    IF(.NOT.test(.NOT.PRESENT(def_val), ler))      val = def_val                         !--- Default value
+  END IF
+  IF(PRESENT(lerr)) lerr = ler
+END FUNCTION fgetKeyIdx_s1
+!==============================================================================================================================
+CHARACTER(LEN=maxlen) FUNCTION fgetKeyNam_s1(tname, keyn, ky, def_val, lerr) RESULT(val)
+  CHARACTER(LEN=*),           INTENT(IN)  :: tname, keyn
+  TYPE(keys_type),            INTENT(IN)  :: ky(:)
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: def_val
+  LOGICAL,          OPTIONAL, INTENT(OUT) :: lerr
+!------------------------------------------------------------------------------------------------------------------------------
+  val = fgetKeyIdx_s1(strIdx(ky(:)%name, tname), keyn, ky, def_val, lerr)
+END FUNCTION fgetKeyNam_s1
+!==============================================================================================================================
+FUNCTION fgetKeys(keyn, ky, def_val, lerr) RESULT(val)
+CHARACTER(LEN=maxlen),        ALLOCATABLE :: val(:)
+  CHARACTER(LEN=*),           INTENT(IN)  :: keyn
+  TYPE(keys_type),            INTENT(IN)  :: ky(:)
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: def_val
+  LOGICAL,          OPTIONAL, INTENT(OUT) :: lerr
+!------------------------------------------------------------------------------------------------------------------------------
+  LOGICAL :: ler(SIZE(ky))
+  INTEGER :: it
+  val = [(fgetKeyIdx_s1(it, keyn, ky, def_val, ler(it)), it = 1, SIZE(ky))]
+  IF(PRESENT(lerr)) lerr = ANY(ler)
+END FUNCTION fgetKeys
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!========== GET THE VALUE OF A KEY FROM A "keys_type" DERIVED TYPE ; THE RETURNED VALUE IS THE ERROR CODE        ==============
+!==========  The key "keyn" is searched in: 1)           "ky(:)%name" (if given)                                 ==============
+!==========                                 2)      "tracers(:)%name"                                            ==============
+!==========                                 3) "isotope%keys(:)%name"                                            ==============
+!==========  for the tracer[s] "tname[(:)]" (if given) or all the available tracers from the used set otherwise. ==============
+!==========  The type of the returned value(s) can be string, integer or real, scalar or vector                  ==============
+!==============================================================================================================================
+LOGICAL FUNCTION getKeyByName_s1(keyn, val, tname, ky) RESULT(lerr)
+  CHARACTER(LEN=*),          INTENT(IN)  :: keyn
+  CHARACTER(LEN=maxlen),     INTENT(OUT) :: val
+  CHARACTER(LEN=*),          INTENT(IN)  :: tname
+  TYPE(keys_type), OPTIONAL, INTENT(IN)  :: ky(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen) :: tnam
+  tnam = strHead(delPhase(tname),'_',.FALSE.)                                            !--- Remove tag and phase
+  IF(PRESENT(ky)) THEN                                                                   !=== KEY FROM "ky"
+               val = fgetKeyNam_s1(tname, keyn, ky,           lerr=lerr)                 !--- "ky" and "tname"
+    IF( lerr ) val = fgetKeyNam_s1(tnam,  keyn, ky,           lerr=lerr)                 !--- "ky" and "tnam"
+  ELSE
+    IF(         .NOT.test(.NOT.ALLOCATED(tracers ), lerr)) lerr = SIZE(tracers ) == 0    !=== KEY FROM "tracers"
+    IF(.NOT.lerr) THEN
+               val = fgetKeyNam_s1(tname, keyn, tracers%keys, lerr=lerr)                 !--- "ky" and "tname"
+      IF(lerr) val = fgetKeyNam_s1(tnam,  keyn, tracers%keys, lerr=lerr)                 !--- "ky" and "tnam"
+    END IF
+    IF(lerr.AND..NOT.test(.NOT.ALLOCATED(isotopes), lerr)) lerr = SIZE(isotopes) == 0    !=== KEY FROM "isotope"
+    IF(.NOT.lerr) THEN
+               val = fgetKeyNam_s1(tname, keyn, isotope%keys, lerr=lerr)                 !--- "ky" and "tname"
+      IF(lerr) val = fgetKeyNam_s1(tnam,  keyn, isotope%keys, lerr=lerr)                 !--- "ky" and "tnam"
+    END IF
+  END IF
+END FUNCTION getKeyByName_s1
+!==============================================================================================================================
+LOGICAL FUNCTION getKeyByName_s1m(keyn, val, tname, ky) RESULT(lerr)
+  CHARACTER(LEN=*),                   INTENT(IN)  :: keyn
+  CHARACTER(LEN=maxlen), ALLOCATABLE, INTENT(OUT) :: val(:)
+  CHARACTER(LEN=*),                   INTENT(IN)  :: tname
+  TYPE(keys_type),          OPTIONAL, INTENT(IN)  :: ky(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen) :: sval
+  lerr = getKeyByName_s1(keyn, sval, tname, ky)
+  IF(test(fmsg('missing key "'//TRIM(keyn)//'" for tracer "'//TRIM(tname)//'"', modname, lerr), lerr)) RETURN
+  lerr = strParse(sval, ',', val)
+END FUNCTION getKeyByName_s1m
+!==============================================================================================================================
+LOGICAL FUNCTION getKeyByName_sm(keyn, val, tname, ky, nam) RESULT(lerr)
+  CHARACTER(LEN=*),                             INTENT(IN)  :: keyn
+  CHARACTER(LEN=maxlen),           ALLOCATABLE, INTENT(OUT) :: val(:)
+  CHARACTER(LEN=*),                             INTENT(IN)  :: tname(:)
+  TYPE(keys_type),       OPTIONAL, TARGET,      INTENT(IN)  :: ky(:)
+  CHARACTER(LEN=maxlen), OPTIONAL, ALLOCATABLE, INTENT(OUT) :: nam(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE(keys_type), POINTER ::  keys(:)
+  LOGICAL :: lk, lt, li
+  INTEGER :: iq, nq
+
+  !--- DETERMINE THE DATABASE TO BE USED (ky, tracers or isotope)
+  lk = PRESENT(ky)
+  lt = .NOT.lk .AND. ALLOCATED(tracers);  IF(lt) lt = SIZE(tracers)  /= 0; IF(lt) lt = ANY(tracers(1)%keys%key(:) == keyn)
+  li = .NOT.lt .AND. ALLOCATED(isotopes); IF(li) li = SIZE(isotopes) /= 0; IF(li) li = ANY(isotope%keys(1)%key(:) == keyn)
+
+  !--- LINK "keys" TO THE RIGHT DATABASE
+  IF(test(.NOT.ANY([lk,lt,li]), lerr)) RETURN
+  IF(lk) keys => ky(:)
+  IF(lt) keys => tracers(:)%keys
+  IF(li) keys => isotope%keys(:)
+
+  !--- GET THE DATA
+  nq = SIZE(tname)
+  ALLOCATE(val(nq))
+  lerr = ANY([(getKeyByName_s1(keyn, val(iq), tname(iq), keys(:)), iq=1, nq)])
+  IF(PRESENT(nam)) nam = tname(:)
+
+END FUNCTION getKeyByName_sm
+!==============================================================================================================================
+LOGICAL FUNCTION getKey_sm(keyn, val, ky, nam) RESULT(lerr)
+  CHARACTER(LEN=*),                             INTENT(IN)  :: keyn
+  CHARACTER(LEN=maxlen),           ALLOCATABLE, INTENT(OUT) :: val(:)
+  TYPE(keys_type),       OPTIONAL,              INTENT(IN)  :: ky(:)
+  CHARACTER(LEN=maxlen), OPTIONAL, ALLOCATABLE, INTENT(OUT) :: nam(:)
+!------------------------------------------------------------------------------------------------------------------------------
+! Note: if names are repeated, getKeyByName_sm can't be used ; this routine, using indexes, must be used instead.
+  IF(PRESENT(ky)) THEN                                                                   !=== KEY FROM "ky"
+    val = fgetKeys(keyn, ky, lerr=lerr)
+    IF(PRESENT(nam)) nam = ky(:)%name
+  ELSE
+    IF(         .NOT.test(.NOT.ALLOCATED(tracers ), lerr)) lerr = SIZE(tracers ) == 0    !=== KEY FROM "tracers"
+    IF(.NOT.lerr) val = fgetKeys(keyn, tracers%keys, lerr=lerr)
+    IF(.NOT.lerr.AND.PRESENT(nam)) nam = tracers(:)%keys%name
+    IF(lerr.AND..NOT.test(.NOT.ALLOCATED(isotopes), lerr)) lerr = SIZE(isotopes) == 0    !=== KEY FROM "isotope"
+    IF(.NOT.lerr) val = fgetKeys(keyn, isotope%keys, lerr=lerr)
+    IF(.NOT.lerr.AND.PRESENT(nam)) nam = isotope%keys(:)%name
+  END IF
+END FUNCTION getKey_sm
+!==============================================================================================================================
+LOGICAL FUNCTION getKeyByName_i1(keyn, val, tname, ky) RESULT(lerr)
+  CHARACTER(LEN=*),          INTENT(IN)  :: keyn
+  INTEGER,                   INTENT(OUT) :: val
+  CHARACTER(LEN=*),          INTENT(IN)  :: tname
+  TYPE(keys_type), OPTIONAL, INTENT(IN)  :: ky(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen) :: sval
+  INTEGER :: ierr
+  lerr = getKeyByName_s1(keyn, sval, tname, ky)
+  IF(test(fmsg('key "'//TRIM(keyn)//'" or tracer "'//TRIM(tname)//'" is missing', modname, lerr), lerr)) RETURN
+  READ(sval, *, IOSTAT=ierr) val
+  IF(test(fmsg('key "'//TRIM(keyn)//'" of tracer "'//TRIM(tname)//'" is not an integer', modname, ierr/=0), lerr)) RETURN
+END FUNCTION getKeyByName_i1
+!==============================================================================================================================
+LOGICAL FUNCTION getKeyByName_i1m(keyn, val, tname, ky) RESULT(lerr)
+  CHARACTER(LEN=*),          INTENT(IN)  :: keyn
+  INTEGER,      ALLOCATABLE, INTENT(OUT) :: val(:)
+  CHARACTER(LEN=*),          INTENT(IN)  :: tname
+  TYPE(keys_type), OPTIONAL, INTENT(IN)  ::  ky(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: sval(:)
+  INTEGER :: ierr, iq, nq
+  IF(test(getKeyByName_s1m(keyn, sval, tname, ky), lerr)) RETURN
+  nq = SIZE(sval); ALLOCATE(val(nq))
+  lerr = .FALSE.; DO iq=1, nq; READ(sval(iq), *, IOSTAT=ierr) val(iq); lerr = lerr .OR. ierr /= 0; END DO
+  IF(test(fmsg('key "'//TRIM(keyn)//'" of tracer "'//TRIM(tname)//'" is not an integer', modname, lerr), lerr)) RETURN
+END FUNCTION getKeyByName_i1m
+!==============================================================================================================================
+LOGICAL FUNCTION getKeyByName_im(keyn, val, tname, ky, nam) RESULT(lerr)
+  CHARACTER(LEN=*),                             INTENT(IN)  ::  keyn
+  INTEGER,                         ALLOCATABLE, INTENT(OUT) ::   val(:)
+  CHARACTER(LEN=*),                             INTENT(IN)  :: tname(:)
+  TYPE(keys_type),       OPTIONAL,              INTENT(IN)  ::    ky(:)
+  CHARACTER(LEN=maxlen), OPTIONAL, ALLOCATABLE, INTENT(OUT) ::   nam(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: sval(:), names(:)
+  INTEGER :: ierr, iq, nq
+  IF(test(getKeyByName_sm(keyn, sval, tname, ky, names), lerr)) RETURN
+  nq = SIZE(sval); ALLOCATE(val(nq))
+  DO iq = 1, nq                                                      !--- CONVERT THE KEYS TO INTEGERS
+    READ(sval(iq), *, IOSTAT=ierr) val(iq)
+    IF(test(fmsg('key "'//TRIM(keyn)//'" of tracer "'//TRIM(names(iq))//'" is not an integer', modname, ierr/=0), lerr)) RETURN
+  END DO
+  IF(PRESENT(nam)) nam = names(:)
+END FUNCTION getKeyByName_im
+!==============================================================================================================================
+LOGICAL FUNCTION getKey_im(keyn, val, ky, nam) RESULT(lerr)
+  CHARACTER(LEN=*),                             INTENT(IN)  :: keyn
+  INTEGER,                         ALLOCATABLE, INTENT(OUT) ::  val(:)
+  TYPE(keys_type),       OPTIONAL,              INTENT(IN)  ::   ky(:)
+  CHARACTER(LEN=maxlen), OPTIONAL, ALLOCATABLE, INTENT(OUT) ::  nam(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: sval(:), names(:)
+  INTEGER :: ierr, iq, nq
+  IF(test(getKey_sm(keyn, sval, ky, names), lerr)) RETURN
+  nq = SIZE(sval); ALLOCATE(val(nq))
+  DO iq = 1, nq
+    READ(sval(iq), *, IOSTAT=ierr) val(iq)
+    IF(test(fmsg('key "'//TRIM(keyn)//'" of tracer "'//TRIM(names(iq))//'" is not an integer', modname, ierr/=0), lerr)) RETURN
+  END DO
+  IF(PRESENT(nam)) nam = names
+END FUNCTION getKey_im
+!==============================================================================================================================
+LOGICAL FUNCTION getKeyByName_r1(keyn, val, tname, ky) RESULT(lerr)
+  CHARACTER(LEN=*),          INTENT(IN)  :: keyn
+  REAL,                      INTENT(OUT) :: val
+  CHARACTER(LEN=*),          INTENT(IN)  :: tname
+  TYPE(keys_type), OPTIONAL, INTENT(IN)  :: ky(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen) :: sval
+  INTEGER :: ierr
+  lerr = getKeyByName_s1(keyn, sval, tname, ky)
+  IF(test(fmsg('key "'//TRIM(keyn)//'" or tracer "'//TRIM(tname)//'" is missing',    modname, lerr), lerr)) RETURN
+  READ(sval, *, IOSTAT=ierr) val
+  IF(test(fmsg('key "'//TRIM(keyn)//'" of tracer "'//TRIM(tname)//'" is not a real', modname, ierr/=0), lerr)) RETURN
+END FUNCTION getKeyByName_r1
+!==============================================================================================================================
+LOGICAL FUNCTION getKeyByName_r1m(keyn, val, tname, ky) RESULT(lerr)
+  CHARACTER(LEN=*),           INTENT(IN)  :: keyn
+  REAL,          ALLOCATABLE, INTENT(OUT) :: val(:)
+  CHARACTER(LEN=*),           INTENT(IN)  :: tname
+  TYPE(keys_type),  OPTIONAL, INTENT(IN)  ::  ky(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: sval(:)
+  INTEGER :: ierr, iq, nq
+  IF(test(getKeyByName_s1m(keyn, sval, tname, ky), lerr)) RETURN
+  nq = SIZE(sval); ALLOCATE(val(nq))
+  lerr = .FALSE.; DO iq=1, nq; READ(sval(iq), *, IOSTAT=ierr) val(iq); lerr = lerr .OR. ierr /= 0; END DO
+  IF(test(fmsg('key "'//TRIM(keyn)//'" of tracer "'//TRIM(tname)//'" is not a vector of reals', modname, lerr), lerr)) RETURN
+END FUNCTION getKeyByName_r1m
+!==============================================================================================================================
+LOGICAL FUNCTION getKeyByName_rm(keyn, val, tname, ky, nam) RESULT(lerr)
+  CHARACTER(LEN=*),                             INTENT(IN)  ::  keyn
+  REAL,                            ALLOCATABLE, INTENT(OUT) ::   val(:)
+  CHARACTER(LEN=*),                             INTENT(IN)  :: tname(:)
+  TYPE(keys_type),       OPTIONAL,              INTENT(IN)  ::    ky(:)
+  CHARACTER(LEN=maxlen), OPTIONAL, ALLOCATABLE, INTENT(OUT) ::   nam(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: sval(:), names(:)
+  INTEGER :: ierr, iq, nq
+  IF(test(getKeyByName_sm(keyn, sval, tname, ky, names), lerr)) RETURN
+  nq = SIZE(sval); ALLOCATE(val(nq))
+  DO iq = 1, nq                                                      !--- CONVERT THE KEYS TO INTEGERS
+    READ(sval(iq), *, IOSTAT=ierr) val(iq)
+    IF(test(fmsg('key "'//TRIM(keyn)//'" of tracer "'//TRIM(names(iq))//'" is not a real', modname, ierr/=0), lerr)) RETURN
+  END DO
+  IF(PRESENT(nam)) nam = names
+END FUNCTION getKeyByName_rm
+!==============================================================================================================================
+LOGICAL FUNCTION getKey_rm(keyn, val, ky, nam) RESULT(lerr)
+  CHARACTER(LEN=*),                             INTENT(IN)  :: keyn
+  REAL,                            ALLOCATABLE, INTENT(OUT) ::  val(:)
+  TYPE(keys_type),       OPTIONAL,              INTENT(IN)  ::   ky(:)
+  CHARACTER(LEN=maxlen), OPTIONAL, ALLOCATABLE, INTENT(OUT) ::  nam(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: sval(:), names(:)
+  INTEGER :: ierr, iq, nq
+  IF(test(getKey_sm(keyn, sval, ky, names), lerr)) RETURN
+  nq = SIZE(sval); ALLOCATE(val(nq))
+  DO iq = 1, nq                                                      !--- CONVERT THE KEYS TO INTEGERS
+    READ(sval(iq), *, IOSTAT=ierr) val(iq)
+    IF(test(fmsg('key "'//TRIM(keyn)//'" of tracer "'//TRIM(names(iq))//'" is not a real', modname, ierr/=0), lerr)) RETURN
+  END DO
+  IF(PRESENT(nam)) nam = names
+END FUNCTION getKey_rm
+!==============================================================================================================================
+LOGICAL FUNCTION getKeyByName_l1(keyn, val, tname, ky) RESULT(lerr)
+  USE strings_mod, ONLY: str2bool
+  CHARACTER(LEN=*),          INTENT(IN)  :: keyn
+  LOGICAL,                   INTENT(OUT) :: val
+  CHARACTER(LEN=*),          INTENT(IN)  :: tname
+  TYPE(keys_type), OPTIONAL, INTENT(IN)  :: ky(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen) :: sval
+  lerr = getKeyByName_s1(keyn, sval, tname, ky)
+  IF(test(fmsg('key "'//TRIM(keyn)//'" or tracer "'//TRIM(tname)//'" is missing', modname, lerr), lerr)) RETURN
+  val = str2bool(sval)
+END FUNCTION getKeyByName_l1
+!==============================================================================================================================
+LOGICAL FUNCTION getKeyByName_l1m(keyn, val, tname, ky) RESULT(lerr)
+  USE strings_mod, ONLY: str2bool
+  CHARACTER(LEN=*),           INTENT(IN)  :: keyn
+  LOGICAL,       ALLOCATABLE, INTENT(OUT) :: val(:)
+  CHARACTER(LEN=*),           INTENT(IN)  :: tname
+  TYPE(keys_type),  OPTIONAL, INTENT(IN)  ::  ky(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: sval(:)
+  INTEGER :: iq, nq
+  IF(test(getKeyByName_s1m(keyn, sval, tname, ky), lerr)) RETURN
+  nq = SIZE(sval); ALLOCATE(val(nq))
+  lerr = .FALSE.; DO iq=1, nq; val(iq)=str2bool(sval(iq)); END DO
+END FUNCTION getKeyByName_l1m
+!==============================================================================================================================
+LOGICAL FUNCTION getKeyByName_lm(keyn, val, tname, ky, nam) RESULT(lerr)
+  USE strings_mod, ONLY: str2bool
+  CHARACTER(LEN=*),                             INTENT(IN)  ::  keyn
+  LOGICAL,                         ALLOCATABLE, INTENT(OUT) ::   val(:)
+  CHARACTER(LEN=*),                             INTENT(IN)  :: tname(:)
+  TYPE(keys_type),       OPTIONAL,              INTENT(IN)  ::    ky(:)
+  CHARACTER(LEN=maxlen), OPTIONAL, ALLOCATABLE, INTENT(OUT) ::   nam(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: sval(:)
+  INTEGER :: iq, nq
+  IF(test(getKeyByName_sm(keyn, sval, tname, ky, nam), lerr)) RETURN
+  nq = SIZE(sval); ALLOCATE(val(nq))
+  lerr = .FALSE.; DO iq=1, nq; val(iq)=str2bool(sval(iq)); END DO
+END FUNCTION getKeyByName_lm
+!==============================================================================================================================
+LOGICAL FUNCTION getKey_lm(keyn, val, ky, nam) RESULT(lerr)
+  USE strings_mod, ONLY: str2bool
+  CHARACTER(LEN=*),                             INTENT(IN)  :: keyn
+  LOGICAL,                         ALLOCATABLE, INTENT(OUT) ::  val(:)
+  TYPE(keys_type),       OPTIONAL,              INTENT(IN)  ::   ky(:)
+  CHARACTER(LEN=maxlen), OPTIONAL, ALLOCATABLE, INTENT(OUT) ::  nam(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: sval(:)
+  INTEGER :: iq, nq
+  IF(test(getKey_sm(keyn, sval, ky, nam), lerr)) RETURN
+  nq = SIZE(sval); ALLOCATE(val(nq))
+  lerr = .FALSE.; DO iq=1, nq; val(iq)=str2bool(sval(iq)); END DO
+END FUNCTION getKey_lm
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== ROUTINES TO GET OR PUT TRACERS AND ISOTOPES DATABASES, SINCE tracers AND isotopes ARE PRIVATE VARIABLES ==================
+!==============================================================================================================================
+SUBROUTINE setKeysDBase(tracers_, isotopes_, isotope_)
+  TYPE(trac_type), OPTIONAL, INTENT(IN) ::  tracers_(:)
+  TYPE(isot_type), OPTIONAL, INTENT(IN) :: isotopes_(:)
+  TYPE(isot_type), OPTIONAL, INTENT(IN) :: isotope_
+!------------------------------------------------------------------------------------------------------------------------------
+  TYPE(isot_type), ALLOCATABLE :: iso(:)
+  INTEGER :: ix, nbIso
+  IF(PRESENT( tracers_)) THEN;  tracers =  tracers_; ELSE; ALLOCATE( tracers(0)); END IF
+  IF(PRESENT(isotopes_)) THEN; isotopes = isotopes_; ELSE; ALLOCATE(isotopes(0)); END IF
+  IF(PRESENT(isotope_ )) THEN
+    ix = strIdx(isotopes(:)%parent, isotope_%parent)
+    IF(ix /= 0) THEN
+      isotopes(ix) = isotope_
+    ELSE
+      nbIso = SIZE(isotopes); ALLOCATE(iso(nbIso+1)); iso(1:nbIso) = isotopes; iso(nbIso+1) = isotope_
+      CALL MOVE_ALLOC(FROM=iso, TO=isotopes)
+    END IF
+  END IF
+END SUBROUTINE setKeysDBase
+!==============================================================================================================================
+SUBROUTINE getKeysDBase(tracers_, isotopes_, isotope_)
+  TYPE(trac_type), OPTIONAL, ALLOCATABLE, INTENT(OUT) ::  tracers_(:)
+  TYPE(isot_type), OPTIONAL, ALLOCATABLE, INTENT(OUT) :: isotopes_(:)
+  TYPE(isot_type), OPTIONAL,              INTENT(OUT) :: isotope_
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: ix
+  IF(PRESENT( tracers_)) THEN;  tracers_ =  tracers; ELSE; ALLOCATE( tracers_(0)); END IF
+  IF(PRESENT(isotopes_)) THEN; isotopes_ = isotopes; ELSE; ALLOCATE(isotopes_(0)); END IF
+  IF(PRESENT(isotope_ )) THEN; ix = strIdx(isotopes(:)%parent, isotope%parent); IF(ix /= 0) isotope_=isotopes(ix); END IF
+END SUBROUTINE getKeysDBase
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== REMOVE, IF ANY, OR ADD THE PHASE SUFFIX OF THE TRACER NAMED "s" ==========================================================
+!==============================================================================================================================
+ELEMENTAL CHARACTER(LEN=maxlen) FUNCTION delPhase(s) RESULT(out)
+  CHARACTER(LEN=*), INTENT(IN) :: s
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: ix, ip, ns
+  out = s; ns = LEN_TRIM(s)
+  IF(s == '')               RETURN                                             !--- Empty string: nothing to do
+  IF(s(1:3)=='H2O' .AND. INDEX(old_phases,s(4:4))/=0 .AND. (ns==4 .OR. s(5:5)=='_')) THEN
+    out='H2O'//s(5:ns)                                                         !--- H2O<phase>[_<iso>][_<tag>]
+  ELSE IF(s(ns-1:ns-1)==phases_sep .AND. INDEX(known_phases,s(ns:ns))/=0) THEN
+    out = s(1:ns-2); RETURN                                                    !--- <var><phase_sep><phase>
+  ELSE; DO ip = 1, nphases; ix = INDEX(s, phases_sep//known_phases(ip:ip)//'_'); IF(ix /= 0) EXIT; END DO
+    IF(ip /= nphases+1) out = s(1:ix-1)//s(ix+2:ns)                            !--- <var><phase_sep><phase>_<tag>
+  END IF
+END FUNCTION delPhase
+!==============================================================================================================================
+CHARACTER(LEN=maxlen) FUNCTION addPhase_s1(s,pha) RESULT(out)
+  CHARACTER(LEN=*),           INTENT(IN) :: s
+  CHARACTER(LEN=1),           INTENT(IN) :: pha
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: l, i
+  out = s
+  IF(s == '') RETURN                                                           !--- Empty string: nothing to do
+  i = INDEX(s, '_')                                                            !--- /=0 for <var>_<tag> tracers names
+  l = LEN_TRIM(s)
+  IF(i == 0) out =  TRIM(s)//phases_sep//pha                                   !--- <var>       => return <var><sep><pha>
+  IF(i /= 0) out = s(1:i-1)//phases_sep//pha//'_'//s(i+1:l)                    !--- <var>_<tag> => return <var><sep><pha>_<tag>
+END FUNCTION addPhase_s1
+!==============================================================================================================================
+FUNCTION addPhase_sm(s,pha) RESULT(out)
+  CHARACTER(LEN=*),           INTENT(IN) :: s(:)
+  CHARACTER(LEN=1),           INTENT(IN) :: pha
+  CHARACTER(LEN=maxlen),     ALLOCATABLE :: out(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: k
+  out = [( addPhase_s1(s(k), pha), k=1, SIZE(s) )]
+END FUNCTION addPhase_sm
+!==============================================================================================================================
+CHARACTER(LEN=maxlen) FUNCTION addPhase_i1(s,ipha,phases) RESULT(out)
+  CHARACTER(LEN=*),           INTENT(IN) :: s
+  INTEGER,                    INTENT(IN) :: ipha
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: phases
+!------------------------------------------------------------------------------------------------------------------------------
+  out = s
+  IF(s == '') RETURN                                                           !--- Empty string: nothing to do
+  IF(ipha == 0 .OR. ipha > nphases) RETURN                                     !--- Absurd index: no phase to add
+  IF(     PRESENT(phases)) out = addPhase_s1(s,       phases(ipha:ipha))
+  IF(.NOT.PRESENT(phases)) out = addPhase_s1(s, known_phases(ipha:ipha))
+END FUNCTION addPhase_i1
+!==============================================================================================================================
+FUNCTION addPhase_im(s,ipha,phases) RESULT(out)
+  CHARACTER(LEN=*),           INTENT(IN) :: s(:)
+  INTEGER,                    INTENT(IN) :: ipha
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: phases
+  CHARACTER(LEN=maxlen),     ALLOCATABLE :: out(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: k
+  IF(     PRESENT(phases)) out = [( addPhase_i1(s(k), ipha,       phases), k=1, SIZE(s) )]
+  IF(.NOT.PRESENT(phases)) out = [( addPhase_i1(s(k), ipha, known_phases), k=1, SIZE(s) )]
+END FUNCTION addPhase_im
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== GET PHASE INDEX IN THE POSSIBLE PHASES LIST OR IN A SPECIFIED LIST ("phases") ============================================
+!==============================================================================================================================
+INTEGER FUNCTION getiPhase(tname, phases) RESULT(iPhase)
+  CHARACTER(LEN=*),           INTENT(IN)  :: tname
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: phases
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen) :: phase
+  IF(     PRESENT(phases)) phase = getPhase(tname,       phases, iPhase)
+  IF(.NOT.PRESENT(phases)) phase = getPhase(tname, known_phases, iPhase)
+END FUNCTION getiPhase
+!==============================================================================================================================
+CHARACTER(LEN=1) FUNCTION getPhase(tname, phases, iPhase) RESULT(phase)
+  CHARACTER(LEN=*),           INTENT(IN)  :: tname
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: phases
+  INTEGER,          OPTIONAL, INTENT(OUT) :: iPhase
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: ip
+  phase = TRIM(strHead(strTail(tname, phases_sep, .TRUE.), '_', .TRUE.))
+  IF(     PRESENT(phases)) ip = INDEX(      phases, phase)
+  IF(.NOT.PRESENT(phases)) ip = INDEX(known_phases, phase)
+  IF(ip == 0) phase = 'g'
+  IF(PRESENT(iPhase)) iPhase = ip
+END FUNCTION getPhase
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!============ CONVERT WATER-DERIVED NAMES FROM FORMER TO CURRENT CONVENTION ; OTHER NAMES ARE LEFT UNTOUCHED ==================
+!======= NAMES STRUCTURE: H2O[<phase>][_<isotope>][_<tag>] (<phase> from "old_phases", <isotope> from "oldH2OIso") ============
+!==============================================================================================================================
+CHARACTER(LEN=maxlen) FUNCTION old2newH2O_1(oldName, iPhase) RESULT(newName)
+  CHARACTER(LEN=*),  INTENT(IN)  :: oldName
+  INTEGER, OPTIONAL, INTENT(OUT) :: iPhase
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: tmp(:)
+  INTEGER :: ix, ip, nt
+  LOGICAL :: lerr
+  newName = oldName
+  IF(PRESENT(iPhase)) iPhase = 1                                               !--- Default: gaseous phase
+  lerr = strParse(oldName, '_', tmp, nt)                                       !--- Parsing: 1 up to 3 elements.
+  ip = strIdx( [('H2O'//old_phases(ix:ix), ix=1, nphases)], tmp(1) )           !--- Phase index
+  IF(ip /= 0 .AND. PRESENT(iPhase)) iPhase = ip                                !--- Returning phase index
+  IF(ip == 0 .AND. tmp(1) /= 'H2O')   RETURN                                   !--- Not an old-style water-related species
+  IF(nt == 1) THEN
+    newName = addPhase('H2O',ip)                                               !=== WATER WITH OR WITHOUT PHASE
+  ELSE
+    ix = strIdx(oldH2OIso, tmp(2))                                             !--- Index in the known isotopes list
+    IF(ix /= 0) tmp(2) = newH2OIso(ix)                                         !--- Move to new isotope name
+    IF(ip /= 0) tmp(2) = addPhase(tmp(2), ip)                                  !--- Add phase to isotope name
+    newName = TRIM(strStack(tmp(2:nt), '_'))                                   !=== WATER ISOTOPE OR TAGGING TRACER
+  END IF
+END FUNCTION old2newH2O_1
+!==============================================================================================================================
+FUNCTION old2newH2O_m(oldName) RESULT(newName)
+  CHARACTER(LEN=*), INTENT(IN) :: oldName(:)
+  CHARACTER(LEN=maxlen)        :: newName(SIZE(oldName))
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: i
+  newName = [(old2newH2O_1(oldName(i)), i=1, SIZE(oldName))]
+END FUNCTION old2newH2O_m
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!============ CONVERT WATER-DERIVED NAMES FROM CURRENT TO FORMER CONVENTION ; OTHER NAMES ARE LEFT UNTOUCHED ==================
+!==== NAMES STRUCTURE: <var>[<phase_sep><phase>][_<tag>] (<phase> from "known_phases", <var> = 'H2O' or from "newH2OIso") =====
+!==============================================================================================================================
+CHARACTER(LEN=maxlen) FUNCTION new2oldH2O_1(newName, iPhase) RESULT(oldName)
+  CHARACTER(LEN=*),  INTENT(IN)  :: newName
+  INTEGER, OPTIONAL, INTENT(OUT) :: iPhase
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: ix, ip
+  CHARACTER(LEN=maxlen) :: var
+  oldName = newName
+  IF(PRESENT(iPhase)) iPhase = 1                                               !--- Default: gaseous phase
+  ip = getiPhase(newName)                                                      !--- Phase index
+  IF(ip /= 0 .AND. PRESENT(iPhase)) iPhase = ip                                !--- Returning phase index
+  ix = strIdx(newH2OIso, strHead(newName, '_'))                                !--- Index in the known H2O isotopes list
+  IF(ix /= 0) THEN
+    oldName = 'H2O'//'_'//TRIM(oldH2OIso(ix))                                  !=== WATER ISOTOPE WITHOUT PHASE
+    IF(newH2OIso(ix)/=newName) oldName=TRIM(oldName)//'_'//strTail(newName,'_')!=== WATER ISOTOPIC TAGGING TRACER WITHOUT PHASE
+  END IF
+  IF(ix /= 0 .OR. ip == 0)           RETURN
+  oldName = 'H2O'//old_phases(ip:ip)
+  IF(newName == addPhase('H2O', ip)) RETURN                                    !=== WATER WITH PHASE
+  var = TRIM(strHead(newName, phases_sep, .TRUE.))                             !--- Head variable name   (no phase)
+  ix = strIdx(newH2OIso, var)                                                  !--- Index in the known H2O isotopes list
+  IF(ix == 0)                        RETURN                                    !=== H2O[vli]_<var> (<var> /= H2O isotope)
+  oldName = TRIM(oldName)//'_'//TRIM(oldH2OIso(ix))                            !=== WATER ISOTOPE WITH PHASE
+  var = addPhase(var, ip)                                                      !--- Head variable with phase
+  IF(newName /= var) oldName = TRIM(oldName)//strTail(newName, TRIM(var))      !=== WATER ISOTOPIC TAGGING TRACER
+END FUNCTION new2oldH2O_1
+!==============================================================================================================================
+FUNCTION new2oldH2O_m(newName) RESULT(oldName)
+  CHARACTER(LEN=*), INTENT(IN) :: newName(:)
+  CHARACTER(LEN=maxlen)        :: oldName(SIZE(newName))
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: i
+  oldName = [(new2oldH2O_1(newName(i)), i=1, SIZE(newName))]
+END FUNCTION new2oldH2O_m
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== GET THE NAME(S) OF THE ANCESTOR(S) OF TRACER(S) "tname" AT GENERATION "igen"  IN THE TRACERS DESCRIPTORS LIST "tr" =======
+!==============================================================================================================================
+SUBROUTINE ancestor_1(t, out, tname, igen)
+  TYPE(trac_type),       INTENT(IN)  :: t(:)
+  CHARACTER(LEN=maxlen), INTENT(OUT) :: out
+  CHARACTER(LEN=*),      INTENT(IN)  :: tname
+  INTEGER,     OPTIONAL, INTENT(IN)  :: igen
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: ix
+  CALL idxAncestor_1(t, ix, tname, igen)
+  out = ''; IF(ix /= 0) out = t(ix)%name
+END SUBROUTINE ancestor_1
+!==============================================================================================================================
+SUBROUTINE ancestor_mt(t, out, tname, igen)
+  TYPE(trac_type),       INTENT(IN)  :: t(:)
+  CHARACTER(LEN=*),      INTENT(IN)  :: tname(:)
+  CHARACTER(LEN=maxlen), INTENT(OUT) :: out(SIZE(tname))
+  INTEGER,     OPTIONAL, INTENT(IN)  :: igen
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: ix(SIZE(tname))
+  CALL idxAncestor_mt(t, ix, tname, igen)
+  out(:) = ''; WHERE(ix /= 0) out = t(ix)%name
+END SUBROUTINE ancestor_mt
+!==============================================================================================================================
+SUBROUTINE ancestor_m(t, out, igen)
+  TYPE(trac_type),       INTENT(IN)  :: t(:)
+  CHARACTER(LEN=maxlen), INTENT(OUT) :: out(SIZE(t))
+  INTEGER,     OPTIONAL, INTENT(IN)  :: igen
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: ix(SIZE(t))
+  CALL idxAncestor_m(t, ix, igen)
+  out(:) = ''; WHERE(ix /= 0) out = t(ix)%name
+END SUBROUTINE ancestor_m
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== GET THE INDEX(ES) OF THE GENERATION "igen" ANCESTOR(S) OF "tname" (t%name IF UNSPECIFIED) IN THE "t" LIST ================
+!==============================================================================================================================
+SUBROUTINE idxAncestor_1(t, idx, tname, igen)
+  TYPE(trac_type),   INTENT(IN)  :: t(:)
+  INTEGER,           INTENT(OUT) :: idx
+  CHARACTER(LEN=*),  INTENT(IN)  :: tname
+  INTEGER, OPTIONAL, INTENT(IN)  :: igen
+  INTEGER :: ig
+  ig = 0; IF(PRESENT(igen)) ig = igen
+  idx = strIdx(t(:)%name, tname)
+  IF(idx == 0)                 RETURN            !--- Tracer not found
+  IF(t(idx)%iGeneration <= ig) RETURN            !--- Tracer has a lower generation number than asked generation 'igen"
+  DO WHILE(t(idx)%iGeneration > ig); idx = strIdx(t(:)%name, t(idx)%parent); END DO
+END SUBROUTINE idxAncestor_1
+!------------------------------------------------------------------------------------------------------------------------------
+SUBROUTINE idxAncestor_mt(t, idx, tname, igen)
+  TYPE(trac_type),   INTENT(IN)  :: t(:)
+  CHARACTER(LEN=*),  INTENT(IN)  :: tname(:)
+  INTEGER,           INTENT(OUT) :: idx(SIZE(tname))
+  INTEGER, OPTIONAL, INTENT(IN)  :: igen
+  INTEGER :: ix
+  DO ix = 1, SIZE(tname); CALL idxAncestor_1(t, idx(ix), tname(ix), igen); END DO
+END SUBROUTINE idxAncestor_mt
+!------------------------------------------------------------------------------------------------------------------------------
+SUBROUTINE idxAncestor_m(t, idx, igen)
+  TYPE(trac_type),   INTENT(IN)  :: t(:)
+  INTEGER,           INTENT(OUT) :: idx(SIZE(t))
+  INTEGER, OPTIONAL, INTENT(IN)  :: igen
+  INTEGER :: ix
+  DO ix = 1, SIZE(t); CALL idxAncestor_1(t, idx(ix), t(ix)%name, igen); END DO
+END SUBROUTINE idxAncestor_m
+!==============================================================================================================================
+
+
+END MODULE readTracFiles_mod
Index: LMDZ6/branches/Ocean_skin/libf/misc/strings_mod.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/misc/strings_mod.F90	(revision 4368)
+++ LMDZ6/branches/Ocean_skin/libf/misc/strings_mod.F90	(revision 4368)
@@ -0,0 +1,1507 @@
+MODULE strings_mod
+
+  IMPLICIT NONE
+
+  PRIVATE
+  PUBLIC :: maxlen, init_printout, msg, fmsg, get_in, lunout, prt_level
+  PUBLIC :: strLower, strHead, strStack,  strCount, strReduce,  strClean, strIdx
+  PUBLIC :: strUpper, strTail, strStackm, strParse, strReplace, strFind, find, cat
+  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_m; END INTERFACE  msg
+  INTERFACE fmsg;       MODULE PROCEDURE       fmsg_1,                   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 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 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
+
+  INTEGER, PARAMETER :: maxlen    = 256                    !--- Standard maximum length for strings
+  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_
+  prt_level = prt_level_
+END SUBROUTINE init_printout
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== Same as getin ; additional last argument: the default value.
+!==============================================================================================================================
+SUBROUTINE getin_s(nam, val, def)
+USE ioipsl_getincom, ONLY: getin
+  CHARACTER(LEN=*), INTENT(IN)    :: nam
+  CHARACTER(LEN=*), INTENT(INOUT) :: val
+  CHARACTER(LEN=*), INTENT(IN)    :: def
+  val = def; CALL getin(nam, val)
+  IF(val/=def) WRITE(lunout,*)TRIM(nam)//' = '//TRIM(val)
+END SUBROUTINE getin_s
+!==============================================================================================================================
+SUBROUTINE getin_i(nam, val, def)
+USE ioipsl_getincom, ONLY: getin
+  CHARACTER(LEN=*), INTENT(IN)    :: nam
+  INTEGER,          INTENT(INOUT) :: val
+  INTEGER,          INTENT(IN)    :: def
+  val = def; CALL getin(nam, val)
+  IF(val/=def) WRITE(lunout,*)TRIM(nam)//' = '//TRIM(int2str(val))
+END SUBROUTINE getin_i
+!==============================================================================================================================
+SUBROUTINE getin_r(nam, val, def)
+USE ioipsl_getincom, ONLY: getin
+  CHARACTER(LEN=*), INTENT(IN)    :: nam
+  REAL,             INTENT(INOUT) :: val
+  REAL,             INTENT(IN)    :: def
+  val = def; CALL getin(nam, val)
+  IF(val/=def) WRITE(lunout,*)TRIM(nam)//' = '//TRIM(real2str(val))
+END SUBROUTINE getin_r
+!==============================================================================================================================
+SUBROUTINE getin_l(nam, val, def)
+USE ioipsl_getincom, ONLY: getin
+  CHARACTER(LEN=*), INTENT(IN)    :: nam
+  LOGICAL,          INTENT(INOUT) :: val
+  LOGICAL,          INTENT(IN)    :: def
+  val = def; CALL getin(nam, val)
+  IF(val.NEQV.def) 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, modname, ll, unit)
+  !--- Display a simple message "str". Optional parameters:
+  !    * "modname": module name, displayed in front of the message (with ": " separator) if present.
+  !    * "ll":      message trigger ; message is displayed only if ll==.TRUE.
+  !    * "unit":    write unit (by default: "lunout")
+  CHARACTER(LEN=*),           INTENT(IN) :: str
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: modname
+  LOGICAL,          OPTIONAL, INTENT(IN) :: ll
+  INTEGER,          OPTIONAL, INTENT(IN) :: unit
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen) :: subn
+  INTEGER :: unt
+  subn = '';    IF(PRESENT(modname)) subn = modname
+  IF(PRESENT(ll)) THEN; IF(.NOT.ll) RETURN; END IF
+  unt = lunout; IF(PRESENT(unit)) unt = unit
+  IF(subn == '') WRITE(unt,'(a)') str                                          !--- Simple message
+  IF(subn /= '') WRITE(unt,'(a)') TRIM(subn)//': '//str                        !--- Routine name provided
+END SUBROUTINE msg_1
+!==============================================================================================================================
+SUBROUTINE msg_m(str, modname, ll, unit, nmax)
+  !--- Same as msg_1 with multiple strings that are stacked (separator: coma) on up to "nmax" full lines.
+  CHARACTER(LEN=*),           INTENT(IN) :: str(:)
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: modname
+  LOGICAL,          OPTIONAL, INTENT(IN) :: ll
+  INTEGER,          OPTIONAL, INTENT(IN) :: unit
+  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: s(:)
+  CHARACTER(LEN=maxlen) :: subn
+  INTEGER :: unt, nmx, k
+  LOGICAL :: l
+  subn = '';    IF(PRESENT(modname)) subn = modname
+  l   = .TRUE.; IF(PRESENT(ll))     l = ll
+  unt = lunout; IF(PRESENT(unit)) unt = unit
+  nmx = 128;    IF(PRESENT(nmax)) nmx = nmax
+  s = strStackm(str, ', ', nmx)
+  DO k=1,SIZE(s); CALL msg_1(s(k), subn,  l,   unt); END DO
+END SUBROUTINE msg_m
+!==============================================================================================================================
+LOGICAL FUNCTION fmsg_1(str, modname, ll, unit) RESULT(l)
+  CHARACTER(LEN=*),           INTENT(IN) :: str
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: modname
+  LOGICAL,          OPTIONAL, INTENT(IN) :: ll
+  INTEGER,          OPTIONAL, INTENT(IN) :: unit
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen) :: subn
+  INTEGER :: unt
+  subn = '';    IF(PRESENT(modname)) subn = modname
+  l   = .TRUE.; IF(PRESENT(ll))     l = ll
+  unt = lunout; IF(PRESENT(unit)) unt = unit
+  CALL msg_1(str, subn, l, unt)
+END FUNCTION fmsg_1
+!==============================================================================================================================
+LOGICAL FUNCTION fmsg_m(str, modname, ll, unit, nmax) RESULT(l)
+  CHARACTER(LEN=*),           INTENT(IN)  :: str(:)
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: modname
+  LOGICAL,          OPTIONAL, INTENT(IN) :: ll
+  INTEGER,          OPTIONAL, INTENT(IN) :: unit
+  INTEGER,          OPTIONAL, INTENT(IN)  :: nmax
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen) :: subn
+  INTEGER :: unt, nmx
+  subn = '';    IF(PRESENT(modname)) subn = modname
+  l   = .TRUE.; IF(PRESENT(ll))     l = ll
+  unt = lunout; IF(PRESENT(unit)) unt = unit
+  nmx = 128;    IF(PRESENT(nmax)) nmx = nmax
+  CALL msg_m(str, subn, l, unt, nmx)
+END FUNCTION fmsg_m
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== Lower/upper case conversion function. ====================================================================================
+!==============================================================================================================================
+ELEMENTAL CHARACTER(LEN=maxlen) 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' .AND. str(k:k)<='Z') out(k:k)=ACHAR(IACHAR(str(k:k))+32)
+  END DO
+END FUNCTION strLower
+!==============================================================================================================================
+ELEMENTAL CHARACTER(LEN=maxlen) 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' .AND. str(k:k)<='z') out(k:k)=ACHAR(IACHAR(str(k:k))-32)
+  END DO
+END FUNCTION strUpper
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== Extract the substring in front of the last (first if lFirst==TRUE) occurrence of separator "sep" in "str" ================
+!=== Examples for str='a_b_c' and sep='_' and the bash command with the same effect:                           ================
+!===    * strHead(..,.FALSE.) = 'a'           ${str%%$sep*}                                                    ================
+!===    * strHead(..,.TRUE.)  = 'a_b'         ${str%$sep*}                                                     ================
+!==============================================================================================================================
+CHARACTER(LEN=maxlen) FUNCTION strHead_1(str,sep,lFirst) RESULT(out)
+  CHARACTER(LEN=*),           INTENT(IN) :: str
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: sep
+  LOGICAL,          OPTIONAL, INTENT(IN) :: lFirst
+!------------------------------------------------------------------------------------------------------------------------------
+  LOGICAL :: lf
+  lf = .FALSE.; IF(PRESENT(lFirst)) lf=lFirst
+  IF(PRESENT(sep)) THEN
+    out = str(1:INDEX(str,sep,.NOT.lf)-1)
+  ELSE
+    out = str(1:INDEX(str,'/',.NOT.lf)-1)
+  END IF
+  IF(out == '') out = str
+END FUNCTION strHead_1
+!==============================================================================================================================
+FUNCTION strHead_m(str,sep,lFirst) RESULT(out)
+  CHARACTER(LEN=maxlen),     ALLOCATABLE :: out(:)
+  CHARACTER(LEN=*),           INTENT(IN) :: str(:)
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: sep
+  LOGICAL,          OPTIONAL, INTENT(IN) :: lFirst
+!------------------------------------------------------------------------------------------------------------------------------
+  LOGICAL :: lf
+  INTEGER :: k
+  lf = .FALSE.; IF(PRESENT(lFirst)) lf=lFirst
+  IF(PRESENT(sep)) THEN
+    out = [(strHead_1(str(k), sep,   lf), k=1, SIZE(str))]
+  ELSE
+    out = [(strHead_1(str(k), lFirst=lf), k=1, SIZE(str))]
+  END IF
+END FUNCTION strHead_m
+!==============================================================================================================================
+!=== Extract the substring following the last (first if lFirst==TRUE) occurrence of separator "sep" in "str"   ================
+!=== Examples for str='a_b_c' and sep='_' and the bash command with the same effect:                           ================
+!===    * strTail(..,.FALSE.) = 'c'           ${str#*$sep}                                                     ================
+!===    * strTail(..,.TRUE.)  = 'b_c'         ${str##*$sep}                                                    ================
+!==============================================================================================================================
+CHARACTER(LEN=maxlen) FUNCTION strTail_1(str,sep,lFirst) RESULT(out)
+  CHARACTER(LEN=*),           INTENT(IN) :: str
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: sep
+  LOGICAL,          OPTIONAL, INTENT(IN) :: lFirst
+!------------------------------------------------------------------------------------------------------------------------------
+  LOGICAL :: lf
+  lf = .FALSE.; IF(PRESENT(lFirst)) lf=lFirst
+  IF(PRESENT(sep)) THEN
+    out = str(INDEX(str,sep,.NOT.lf)+LEN(sep):LEN_TRIM(str))
+  ELSE
+    out = str(INDEX(str,'/',.NOT.lf)+1:LEN_TRIM(str))
+  END IF
+  IF(out == '') out = str
+END FUNCTION strTail_1
+!==============================================================================================================================
+FUNCTION strTail_m(str,sep,lFirst) RESULT(out)
+  CHARACTER(LEN=maxlen),     ALLOCATABLE :: out(:)
+  CHARACTER(LEN=*),           INTENT(IN) :: str(:)
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: sep
+  LOGICAL,          OPTIONAL, INTENT(IN) :: lFirst
+!------------------------------------------------------------------------------------------------------------------------------
+  LOGICAL :: lf
+  INTEGER :: k
+  lf = .FALSE.; IF(PRESENT(lFirst)) lf=lFirst
+  IF(PRESENT(sep)) THEN
+    out = [(strTail_1(str(k), sep,   lf), k=1, SIZE(str))]
+  ELSE
+    out = [(strTail_1(str(k), lFirst=lf), k=1, SIZE(str))]
+  END IF
+END FUNCTION strTail_m
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== Concatenates the strings "str(:)" with separator "sep" into a single string using a separator (',' by default). ==========
+!==============================================================================================================================
+FUNCTION strStack(str, sep, mask) RESULT(out)
+  CHARACTER(LEN=:),          ALLOCATABLE :: out
+  CHARACTER(LEN=*),           INTENT(IN) :: str(:)
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: sep
+  LOGICAL,          OPTIONAL, INTENT(IN) :: mask(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=:), ALLOCATABLE :: s
+  INTEGER :: is, i0
+  IF(SIZE(str) == 0) THEN; out = ''; RETURN; END IF
+  ALLOCATE(s, SOURCE=', '); IF(PRESENT(sep)) s=sep
+  IF(PRESENT(mask)) THEN
+    IF(ALL(.NOT.mask)) THEN; out = ''; RETURN; END IF
+    i0 = 0; DO WHILE(.NOT.mask(i0+1)); i0 = i0+1; END DO
+    out = str(i0); DO is=i0+1,SIZE(str, DIM=1); IF(.NOT.mask(is)) CYCLE; out = TRIM(out)//s//TRIM(str(is)); END DO
+  ELSE
+    out = str(1); DO is=2,SIZE(str, DIM=1); out = TRIM(out)//s//TRIM(str(is)); END DO
+  END IF
+END FUNCTION strStack
+!==============================================================================================================================
+!=== 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=maxlen),     ALLOCATABLE :: out(:)
+  CHARACTER(LEN=*),           INTENT(IN) :: str(:)
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: sep
+  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: t(:)
+  CHARACTER(LEN=maxlen) :: 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(str, nb)
+  CHARACTER(LEN=*), ALLOCATABLE, INTENT(INOUT) :: str(:)
+  INTEGER,          OPTIONAL,    INTENT(OUT)   :: nb
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: s1(:)
+  INTEGER :: k, n, n1
+  IF(PRESENT(nb)) nb = 0
+  CALL MOVE_ALLOC(FROM=str, 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(str(n))
+  IF(n==0) RETURN
+  str(1) = s1(1)
+  n=1; DO k=2,n1; IF(ANY(s1(1:k-1)==s1(k))) CYCLE; n=n+1; str(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=maxlen), 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[(:)]" ============================
+!=== OPTIONALY: GET THE NUMBER OF FOUND ELEMENTS "n". NB: UNFOUND => INDEX=0                       ============================
+!==============================================================================================================================
+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 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 :: idx0                                                    !--- Used to display an identified non-numeric string
+  lerr = .FALSE.
+  idx = strIdx1(rawList, del, ibeg, idel)                            !--- idx/=0: del(idel) is at position "idx" in "rawList" 
+  IF(.NOT.PRESENT(lSc))               RETURN                         !--- No need to check exceptions for numbers => finished
+  IF(.NOT.        lSc )               RETURN                         !--- No need to check exceptions for numbers => finished
+
+  !=== No delimiter found: the whole string must be a valid number
+  IF(idx == 0) THEN                                                  !--- No element of "del" in "rawList"
+    lerr = .NOT.is_numeric(rawList(ibeg:LEN_TRIM(rawList)))          !--- String must be a number
+    IF(lerr) idx = LEN_TRIM(rawList); RETURN                         !--- Set idx so that rawList(ibeg:idx-1) = whole string
+  END IF
+
+  IF(test(idx == 1 .AND. INDEX('+-',del(idel)) /= 0, lerr)) RETURN   !--- The front delimiter is different from +/-: error
+  IF(     idx /= 1 .AND. is_numeric(rawList(ibeg:idx-1)))   RETURN   !--- The input string head is a valid number
+
+  !=== The string part in front of the 1st delimiter is not a valid number: search for next delimiter index "idx"
+  idx0 = idx ; idx = strIdx1(rawList, del, idx+1, idel)              !--- Keep start index because idx is recycled
+  IF(idx == 0) THEN
+    lerr = .NOT.is_numeric(rawList(ibeg:LEN_TRIM(rawList)))          !--- No other delimiter: whole string must be a valid numb
+    IF(lerr) idx = idx0; RETURN
+  END IF
+  lerr = .NOT.is_numeric(rawList(ibeg:idx-1))
+
+CONTAINS
+
+!------------------------------------------------------------------------------------------------------------------------------
+INTEGER FUNCTION strIdx1(str, del, ib, id) RESULT(i)
+!--- Get the index of the first appereance of one of the delimiters "del(:)" in "str" starting from position "ib".
+!--- "id" is the index in "del(:)" of the first delimiter found.
+  IMPLICIT NONE
+  CHARACTER(LEN=*),  INTENT(IN)  :: str, del(:)
+  INTEGER,           INTENT(IN)  :: ib
+  INTEGER,           INTENT(OUT) :: id
+!------------------------------------------------------------------------------------------------------------------------------
+  DO i = ib, LEN_TRIM(str); id = strIdx(del, str(i:i)); IF(id /= 0) EXIT; END DO
+  IF(i > LEN_TRIM(str)) THEN; i = 0; id = 0; END IF
+END FUNCTION strIdx1
+
+END FUNCTION strIdx_prv
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+!=== Count the number of elements separated by "delimiter" in list "rawList". =================================================
+!==============================================================================================================================
+LOGICAL FUNCTION strCount_11(rawList, delimiter, nb, lSc) RESULT(lerr)
+  CHARACTER(LEN=*),  INTENT(IN)  :: rawList
+  CHARACTER(LEN=*),  INTENT(IN)  :: delimiter
+  INTEGER,           INTENT(OUT) :: nb
+  LOGICAL, OPTIONAL, INTENT(IN)  :: lSc
+!------------------------------------------------------------------------------------------------------------------------------
+  LOGICAL :: ll
+  ll = .FALSE.; IF(PRESENT(lSc)) ll = lSc
+  lerr = strCount_1m(rawList, [delimiter], nb, ll)
+END FUNCTION strCount_11
+!==============================================================================================================================
+LOGICAL FUNCTION strCount_m1(rawList, delimiter, nb, lSc) RESULT(lerr)
+  CHARACTER(LEN=*),     INTENT(IN)  :: rawList(:)
+  CHARACTER(LEN=*),     INTENT(IN)  :: delimiter
+  INTEGER, ALLOCATABLE, INTENT(OUT) :: nb(:)
+  LOGICAL,    OPTIONAL, INTENT(IN)  :: lSc
+!------------------------------------------------------------------------------------------------------------------------------
+  LOGICAL :: ll
+  INTEGER :: id
+  ll  = .FALSE.; IF(PRESENT(lSc)) ll = lSc .AND. INDEX('+-', delimiter) /= 0
+  lerr = .TRUE.
+  ALLOCATE(nb(SIZE(rawList)))
+  DO id = 1, SIZE(rawList)
+    lerr = lerr .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
+  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(fmsg('"'//TRIM(r(ib:ie-1))//'" is not numeric', ll=lerr)) RETURN
+    IF(ie == 0 .OR. jd == 0) EXIT
+    ib = ie + LEN(delimiter(jd))
+    DO WHILE(r(ib:ib) == ' ' .AND. ib < nr); ib = ib + 1; END DO     !--- Skip spaces before next chain
+    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(rawList, delimiter, keys, n, vals) RESULT(lerr)
+  CHARACTER(LEN=*),                             INTENT(IN)  :: rawList, delimiter
+  CHARACTER(LEN=maxlen), ALLOCATABLE,           INTENT(OUT) :: keys(:)
+  INTEGER,                            OPTIONAL, INTENT(OUT) :: n
+  CHARACTER(LEN=maxlen), ALLOCATABLE, OPTIONAL, INTENT(OUT) :: vals(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=1024) :: r
+  INTEGER :: nr, nk
+  lerr = .FALSE.
+  r  = TRIM(ADJUSTL(rawList))
+  nr = LEN_TRIM(r); IF(nr == 0) THEN; keys = ['']; RETURN; END IF
+  CALL strParse_prv(nk)                                              !--- COUNT THE ELEMENTS
+  ALLOCATE(keys(nk))
+  IF(PRESENT(vals)) THEN
+    ALLOCATE(vals(nk)); CALL strParse_prv(nk, keys, vals)            !--- PARSE THE KEYS
+  ELSE
+    CALL strParse_prv(nk, keys)                                      !--- PARSE THE KEYS
+  END IF
+  IF(PRESENT(n)) n = nk
+
+CONTAINS
+
+!------------------------------------------------------------------------------------------------------------------------------
+SUBROUTINE strParse_prv(nkeys, keys, vals)
+!--- * Get the number of elements after parsing ("nkeys" only is present)
+!--- * Parse the <key>=<val> pairs and store result in "keys" and "vals" (already allocated)
+  IMPLICIT NONE
+  INTEGER,                         INTENT(OUT) :: nkeys
+  CHARACTER(LEN=maxlen), OPTIONAL, INTENT(OUT) :: keys(:)
+  CHARACTER(LEN=maxlen), OPTIONAL, INTENT(OUT) :: vals(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: ib, ie
+  nkeys = 1; ib = 1
+  DO
+    ie = INDEX(rawList(ib:nr), delimiter)+ib-1                       !--- Determine the next separator start index
+    IF(ie == ib-1) EXIT
+    IF(PRESENT(keys)) keys(nkeys) = r(ib:ie-1)                       !--- Get the ikth key
+    IF(PRESENT(vals)) CALL parseKeys(keys(nkeys), vals(nkeys))       !--- Parse the ikth <key>=<val> pair
+    ib = ie + LEN(delimiter)
+    DO WHILE(r(ib:ib) == ' ' .AND. ib < nr); ib = ib + 1; END DO     !--- Skip spaces before next chain
+    nkeys = nkeys+1
+  END DO
+  IF(PRESENT(keys)) keys(nkeys) = r(ib:nr)                           !--- Get the last key
+  IF(PRESENT(vals)) CALL parseKeys(keys(nkeys), vals(nkeys))         !--- Parse the last <key>=<val> pair
+END SUBROUTINE strParse_prv
+
+!------------------------------------------------------------------------------------------------------------------------------
+SUBROUTINE parseKeys(key, val)
+  CHARACTER(LEN=*), INTENT(INOUT) :: key
+  CHARACTER(LEN=*), INTENT(OUT)   :: val
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: ix
+  ix = INDEX(key, '='); IF(ix == 0) RETURN                           !--- First "=" index in "key"
+  val = ADJUSTL(key(ix+1:LEN_TRIM(key)))
+  key = ADJUSTL(key(1:ix-1))
+END SUBROUTINE parseKeys
+
+END FUNCTION strParse
+!==============================================================================================================================
+LOGICAL FUNCTION strParse_m(rawList, delimiter, keys, n, vals, lSc, id) RESULT(lerr)
+  CHARACTER(LEN=*),                             INTENT(IN)  :: rawList, delimiter(:)
+  CHARACTER(LEN=maxlen),           ALLOCATABLE, INTENT(OUT) :: keys(:)  !--- Parsed keys vector
+  INTEGER,               OPTIONAL,              INTENT(OUT) :: n        !--- Length of the parsed vector
+  CHARACTER(LEN=maxlen), OPTIONAL, ALLOCATABLE, INTENT(OUT) :: vals(:)  !--- Values for <name>=<value> keys
+  LOGICAL,               OPTIONAL,              INTENT(IN)  :: lSc      !--- Take care about numbers in scientific notation
+  INTEGER,               OPTIONAL, ALLOCATABLE, INTENT(OUT) :: id(:)    !--- Indexes of the separators in "delimiter(:)" vector
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=1024) :: r
+  INTEGER :: nr, ik, nk, ib, ie, jd
+  LOGICAL :: ll
+  ll = .FALSE.; IF(PRESENT(lSc)) ll = lSc
+  IF(test(fmsg("Couldn't parse list: non-numerical strings were found", ll=strCount_1m(rawList, delimiter, nk, ll)),lerr)) RETURN
+
+  !--- 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('Non-numeric values found', ll=strIdx_prv(r, delimiter, ib, ie, jd, ll)),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
+!------------------------------------------------------------------------------------------------------------------------------
+  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=maxlen), ALLOCATABLE :: out(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen), 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=maxlen), ALLOCATABLE :: out(:,:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen), 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"); 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"); 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"); 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"); 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.
+!===  * nRowMax lines are displayed (default: all lines)
+!===  * nColMax characters (default: as long as needed) are displayed at most on a line. If the effective total length is
+!===    higher, several partial tables are displayed ; the nHead (default: 1) first columns are included in each sub-table.
+!==============================================================================================================================
+LOGICAL FUNCTION dispTable(p, titles, s, i, r, rFmt, nRowMax, nColMax, nHead, unit, sub) 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)  :: nRowMax       !--- Display at most "nRowMax" rows
+  INTEGER,          OPTIONAL, INTENT(IN)  :: nColMax       !--- Display at most "nColMax" characters each line
+  INTEGER,          OPTIONAL, INTENT(IN)  :: nHead         !--- Head columns repeated for multiple tables display
+  INTEGER,          OPTIONAL, INTENT(IN)  :: unit          !--- Output unit (default: screen)
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: sub           !--- Subroutine name
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=2048) :: row
+  CHARACTER(LEN=maxlen)  :: rFm, el, subn
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: d(:,:)
+  CHARACTER(LEN=1) :: s1, sp
+  INTEGER :: is, ii, ir, it, k, nmx,  unt, ic, np
+  INTEGER :: ns, ni, nr, nt, l, ncol, nHd, ib, l0
+  INTEGER, ALLOCATABLE :: n(:), ncmx(:)
+  INTEGER, PARAMETER   :: nm=1                             !--- Space between values & columns
+  LOGICAL :: ls, li, lr
+  subn = '';    IF(PRESENT(sub)) subn = sub
+  rFm = '*';    IF(PRESENT(rFmt)) rFm = rFmt               !--- Specified format for reals
+  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('display map "p" length and titles list mismatch', subn, lerr)) RETURN
+  IF(ls) THEN
+    ns = SIZE(s, 1); ncol = ncol + SIZE(s, 2); lerr = COUNT([(p(ic:ic)=='s', ic=1, np)]) /= SIZE(s, 2)
+  END IF
+  IF(li) THEN
+    ni = SIZE(i, 1); ncol = ncol + SIZE(i, 2); lerr = COUNT([(p(ic:ic)=='i', ic=1, np)]) /= SIZE(i, 2)
+  END IF
+  IF(lr) THEN
+    nr = SIZE(r, 1); ncol = ncol + SIZE(r, 2); lerr = COUNT([(p(ic:ic)=='r', ic=1, np)]) /= SIZE(r, 2)
+  END IF
+  IF(fmsg('display map "p" length and arguments number mismatch', subn, lerr)) RETURN
+  lerr = ncol /= SIZE(titles); IF(fmsg('"titles" length and arguments number mismatch', subn, lerr)) RETURN
+  lerr = ls.AND.li.AND.ns/=ni; IF(fmsg('string and integer arguments lengths mismatch', subn, lerr)) RETURN
+  lerr = ls.AND.lr.AND.ns/=nr; IF(fmsg(   'string and real arguments lengths mismatch', subn, lerr)) RETURN
+  lerr = li.AND.lr.AND.ni/=nr; IF(fmsg(  'integer and real arguments lengths mismatch', subn, lerr)) RETURN
+  nmx = MAX(ns,ni,nr)+1; IF(PRESENT(nRowMax)) nmx = MIN(nmx,nRowMax+1)
+
+  !--- Allocate the assembled quantities array
+  ALLOCATE(d(nmx,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
+  END DO
+  CALL cleanZeros(d)
+  DO ic = 1, ncol
+    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
+
+  !--- Build the vector of max column index in case the rows are too long (table must be displayed in multiple parts)
+  nHd = 1; IF(PRESENT(nHead)) nHd = nHead
+  IF(.NOT.PRESENT(nColMax)) THEN
+    nt = 1; ncmx = [ncol]
+  ELSE
+    nt = 1; l0 = SUM(n(1:nHd)+1)+1
+    IF(PRESENT(sub)) l0=l0+LEN_TRIM(subn)+1
+    !--- Count the number of table parts
+    l = l0; DO ic = nHd+1, ncol; l = l+n(ic)+1; IF(l>=nColMax) THEN; nt = nt+1; l = l0+n(ic)+1; END IF; END DO
+    !--- Get the index of the last column for each table part
+    ALLOCATE(ncmx(nt)); k = 1
+    l = l0; DO ic = nHd+1, ncol; l = l+n(ic)+1; IF(l>=nColMax) THEN; ncmx(k) = ic-1; l = l0+n(ic)+1; k = k+1; END IF; END DO
+    ncmx(nt) = ncol
+  END IF
+      
+  !--- Display the strings array as a table
+  DO it = 1, nt
+    DO ir = 1, nmx; row = ''
+      DO ic = 1, nHd; el = d(ir,ic)
+        s1 = sp
+        row = TRIM(row)//REPEAT(' ',nm)//TRIM(el)//REPEAT(' ',n(ic)-LEN_TRIM(el)-nm)//s1
+      END DO
+      ib = nHd+1; IF(it>1) ib = ncmx(it-1)+1
+      DO ic = ib, ncmx(it); 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), subn, unit=unt)
+      IF(ir /= 1) CYCLE                                              !--- Titles only are underlined
+      row=''; DO ic=1,nHd; row=TRIM(row)//REPEAT('-',n(ic))//'+'; END DO
+      DO ic = ib,ncmx(it); row=TRIM(row)//REPEAT('-',n(ic))//'+'; END DO
+      CALL msg(row(1:LEN_TRIM(row)-1), subn, unit=unt)
+    END DO
+    CALL msg('', subn, unit=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=maxlen)  :: rFm, el
+  CHARACTER(LEN=maxlen), 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
+  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('display map "p" length and arguments number mismatch', ll=lerr)) RETURN
+  lerr = ncol /= SIZE(titles); IF(fmsg('"titles" length and arguments number mismatch', ll=lerr)) RETURN
+  lerr = ls.AND.li.AND.ns/=ni; IF(fmsg('string and integer arguments lengths mismatch', ll=lerr)) RETURN
+  lerr = ls.AND.lr.AND.ns/=nr; IF(fmsg(   'string and real arguments lengths mismatch', ll=lerr)) RETURN
+  lerr = li.AND.lr.AND.ni/=nr; IF(fmsg(  'integer and real arguments lengths mismatch', ll=lerr)) 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
+  END DO
+  CALL cleanZeros(d)
+  DO ic = 1, ncol
+    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, nRowmax, nColMax, nHead, 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)  :: nRowMax                 !--- Maximum number of lines to display    (default: all)
+  INTEGER,          OPTIONAL, INTENT(IN)  :: nColMax                 !--- Maximum number of characters per line (default: 2048)
+  INTEGER,          OPTIONAL, INTENT(IN)  :: nHead                   !--- Number of front columns to duplicate  (default: 1)
+  INTEGER,          OPTIONAL, INTENT(IN)  :: unit                    !--- Output unit                           (def: lunout)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen),      ALLOCATABLE :: ttl(:)
+  LOGICAL,                    ALLOCATABLE :: m(:)
+  INTEGER,                    ALLOCATABLE :: ki(:), kj(:)
+  INTEGER                                 :: i, j, k, rk, rk1, ib, ie, itr, nm, unt, nRmx, nCmx, nHd, nv
+  CHARACTER(LEN=maxlen)                   :: mes, sub, fm='(f12.9)', v, s
+  CHARACTER(LEN=maxlen),      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
+  nRmx= SIZE(a);          IF(PRESENT(nRowMax)) nRmx=MIN(nRmx,nRowMax)!-- Maximum number of lines to print
+  nCmx= 2048;             IF(PRESENT(nColMax)) nCmx=MIN(nCmx,nColMax)!-- Maximum number of characters each line
+  nHd = 1;                IF(PRESENT(nHead))   nHd = nHead           !--- Number of front columns to duplicate
+  unt = lunout;           IF(PRESENT(unit))    unt = unit            !--- Unit to print messages
+
+  rk = SIZE(n); nv = SIZE(vnm)
+  IF(test(fmsg('SIZE(nam) /= 1 and /= last "n" element', sub, nv /= 1 .AND. nv /= n(rk), unt),lerr)) RETURN
+  IF(test(fmsg('ll" and "a" sizes mismatch',             sub, SIZE(a) /= SIZE(ll),       unt),lerr)) RETURN
+  IF(test(fmsg('profile "n" does not match "a" and "ll', sub, SIZE(a) /= PRODUCT(n),     unt),lerr)) RETURN
+  CALL msg(mes, sub, unit=unt)
+
+  !--- 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, nRowMax=nRmx, nColMax=nCmx, nHead=nHd, unit=unt, sub=sub)
+    IF(nv /= 1) lerr = dispTable('ir', ['name ','value'], i=cat(PACK(ki,m)),   r=cat(PACK(a,ll)), &
+                                 rFmt=fm, nRowMax=nRmx, nColMax=nCmx, nHead=nHd, unit=unt, sub=sub)
+    CALL msg("can't display outliers table", sub, lerr, 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, nRowMax=nRmx, nColMax=nCmx, nHead=nHd, unit=unt, sub=sub)
+    IF(rk==3) lerr = dispTable('iir', ttl, i=cat(PACK(ki,m),PACK(kj,m)), r=cat(PACK(a(ib:ie),m)), &
+                                 rFmt=fm, nRowMax=nRmx, nColMax=nCmx, nHead=nHd, unit=unt, sub=sub)
+    CALL msg("can't display outliers table", sub, lerr, unt)
+    IF(lerr) RETURN
+  END DO
+END FUNCTION dispOutliers_1
+!==============================================================================================================================
+LOGICAL FUNCTION dispOutliers_2(ll, a, n, err_msg, nam, subn, nRowMax, nColMax, nHead, 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)  :: nRowMax                 !--- Maximum number of lines to display    (default: all)
+  INTEGER,          OPTIONAL, INTENT(IN)  :: nColMax                 !--- Maximum number of characters per line (default: 2048)
+  INTEGER,          OPTIONAL, INTENT(IN)  :: nHead                   !--- Number of front columns to duplicate  (default: 1)
+  INTEGER,          OPTIONAL, INTENT(IN)  :: unit                    !--- Output unit                           (def: lunout)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen)                   :: mes, sub, fm='(f12.9)', prf
+  CHARACTER(LEN=maxlen),      ALLOCATABLE :: ttl(:), vnm(:)
+  INTEGER,                    ALLOCATABLE :: ki(:), kj(:), kl(:)
+  INTEGER                                 :: i, j, k, rk, nv, unt, nRmx, nCmx, nHd
+  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
+  nRmx= SIZE(a);          IF(PRESENT(nRowMax)) nRmx=MIN(nRmx,nRowMax)!-- Maximum number of lines to print
+  nCmx= 2048;             IF(PRESENT(nColMax)) nCmx=MIN(nCmx,nColMax)!-- Maximum number of characters each line
+  nHd = 1;                IF(PRESENT(nHead))   nHd = nHead           !--- Number of front columns to duplicate
+  unt = lunout;                  IF(PRESENT(unit))    unt = unit     !--- Unit to print messages
+  lerr= SIZE(vnm) /= nv;         IF(fmsg('SIZE(nam) /= SIZE(a,2)',                  sub, lerr, unt)) RETURN
+  lerr= SIZE(a,1) /= SIZE(ll);   IF(fmsg('"ll" and "a" sizes mismatch',             sub, lerr, unt)) RETURN
+  lerr= SIZE(a,1) /= PRODUCT(n); IF(fmsg('profile "n" does not match "a" and "ll"', sub, lerr, unt)) RETURN
+
+  SELECT CASE(rk)                                                   !--- Indices list
+    CASE(0); IF(ll(1)) THEN; WRITE(unt,'(a,", ",a," = ",2f12.9)')TRIM(vnm(1)),TRIM(vnm(2)),a(1,1),a(1,2); RETURN; END IF
+    CASE(1); ki = [  (i,i=1,n(1)) ]
+    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, nRowMax=nRmx, nColMax=nCmx, nHead=nHd, unit=unt, sub=sub)
+  IF(rk == 2) lerr = dispTable(prf, ttl, i = cat(PACK(ki,ll),PACK(kj,ll)),             r = val, &
+                                 rFmt=fm, nRowMax=nRmx, nColMax=nCmx, nHead=nHd, unit=unt, sub=sub)
+  IF(rk == 3) lerr = dispTable(prf, ttl, i = cat(PACK(ki,ll),PACK(kj,ll),PACK(kl,ll)), r = val, &
+                                 rFmt=fm, nRowMax=nRmx, nColMax=nCmx, nHead=nHd, unit=unt, sub=sub)
+  CALL msg("can't display outliers table", sub, lerr, 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(lerr)
+  CHARACTER(LEN=*),      INTENT(IN)  :: str
+  CHARACTER(LEN=maxlen), INTENT(OUT) :: val
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=maxlen)              :: v
+  CHARACTER(LEN=1024)                :: s, vv
+  CHARACTER(LEN=1024), ALLOCATABLE :: vl(:)
+  INTEGER,             ALLOCATABLE :: ip(:)
+  INTEGER :: nl, nn, i, j, im, ix
+  LOGICAL :: ll
+  s = str
+
+  !--- Check wether the parenthesis are correctly formed
+  ll = strCount(s,'(',nl)
+  ll = strCount(s,')',nn)
+  lerr = nl /= nn
+  IF(fmsg('Mismatching number of opening and closing parenthesis: '//TRIM(s), ll=lerr)) 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
+    IF(test(reduceExpr_basic(vl(i+1), v), lerr)) 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
+  lerr = 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(lerr)
+  CHARACTER(LEN=*),      INTENT(IN)  :: str
+  CHARACTER(LEN=*),      INTENT(OUT) :: val
+  DOUBLE PRECISION,      ALLOCATABLE :: vl(:)
+  INTEGER,               ALLOCATABLE :: id(:)
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: ky(:)
+  CHARACTER(LEN=1),      ALLOCATABLE :: op(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  CHARACTER(LEN=1024) :: s
+  DOUBLE PRECISION :: v, vm, vp
+  INTEGER      :: i, ni, io
+  lerr = .FALSE.
+  IF(is_numeric(str)) THEN; val=TRIM(str); RETURN; END IF
+  op = ['^','/','*','+','-']                                                   !--- List of recognized operations
+  s = str
+  IF(test(strParse_m(s, op, ky, lSc=.TRUE., id = id), lerr)) RETURN            !--- Parse the values
+  vl = str2dble(ky)                                                            !--- Conversion to doubles
+  lerr = ANY(vl >= HUGE(1.d0))
+  IF(fmsg('Some values are non-numeric in: '//TRIM(s), ll=lerr)) 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(lerr)
+  LOGICAL,               ALLOCATABLE              :: lerr(:)
+  CHARACTER(LEN=*),                   INTENT(IN)  :: str(:)
+  CHARACTER(LEN=maxlen), ALLOCATABLE, INTENT(OUT) :: val(:)
+!------------------------------------------------------------------------------------------------------------------------------
+  INTEGER :: i
+  ALLOCATE(lerr(SIZE(str)),val(SIZE(str)))
+  lerr(:) = [(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
+  IF(TRIM(str) == '') THEN; out = .FALSE.; RETURN; END IF
+  WRITE(fmt,'("(f",i0,".0)")') LEN_TRIM(str)
+  READ(str,fmt,IOSTAT=e) x
+  out = e==0 .AND. INDEX('Ee',str(LEN_TRIM(str):LEN_TRIM(str)))==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=maxlen) FUNCTION bool2str(b) RESULT(out)
+  LOGICAL, INTENT(IN) :: b
+  WRITE(out,*)b
+  out = ADJUSTL(out)
+END FUNCTION bool2str
+!==============================================================================================================================
+ELEMENTAL CHARACTER(LEN=maxlen) FUNCTION int2str(i, nDigits) RESULT(out)
+  INTEGER,           INTENT(IN) :: i
+  INTEGER, OPTIONAL, INTENT(IN) :: nDigits
+!------------------------------------------------------------------------------------------------------------------------------
+  WRITE(out,*)i
+  out = ADJUSTL(out)
+  IF(.NOT.PRESENT(nDigits)) RETURN
+  IF(nDigits > LEN_TRIM(out)) out = REPEAT('0', nDigits - LEN_TRIM(out))//TRIM(out)
+END FUNCTION int2str
+!==============================================================================================================================
+ELEMENTAL CHARACTER(LEN=maxlen) 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=maxlen) 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
+!==============================================================================================================================
+ELEMENTAL SUBROUTINE cleanZeros(s)
+  CHARACTER(LEN=*), INTENT(INOUT) :: s
+  INTEGER :: ls, ix, i
+  IF(is_numeric(s)) THEN
+    ls = LEN_TRIM(s)
+    ix = MAX(INDEX(s,'E'),INDEX(s,'e'),INDEX(s,'D'),INDEX(s,'d'))
+    IF(ix == 0) THEN
+      DO ix = ls,1,-1; IF(s(ix:ix) /= '0') EXIT; END DO; s=s(1:ix+1)
+    ELSE IF(INDEX(s,'.')/=0) THEN
+      i = ix-1; DO WHILE(s(i:i) == '0'); i = i-1; END DO; s=s(1:i)//s(ix:ls)
+    END IF
+  END IF
+END SUBROUTINE cleanZeros
+!==============================================================================================================================
+
+
+!==============================================================================================================================
+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=maxlen), 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/Ocean_skin/libf/misc/wxios.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/misc/wxios.F90	(revision 3617)
+++ LMDZ6/branches/Ocean_skin/libf/misc/wxios.F90	(revision 4368)
@@ -140,4 +140,5 @@
 !            CALL wxios_context_init()
 !        ENDIF
+         WRITE(*,*)'END of WXIOS_INIT', g_comm , g_ctx_name 
 
     END SUBROUTINE wxios_init
