MODULE readTracFiles_mod_ico

  USE strings_mod_ico,    ONLY: msg, testFile,  strFind, strStack, strReduce,  strHead, strCount, find, fmsg, reduceExpr, &
             removeComment, cat, checkList, str2int, strParse, strReplace, strTail, strIdx, maxlen, test, dispTable, get_in
  USE trac_types_mod_ico, ONLY: trac_type, isot_type, keys_type

  IMPLICIT NONE

  PRIVATE

  PUBLIC :: initIsotopes, maxlen, trac_type, isot_type, keys_type
  PUBLIC :: readTracersFiles, indexUpdate, setGeneration             !--- TOOLS ASSOCIATED TO TRACERS  DESCRIPTORS
  PUBLIC :: readIsotopesFile                                         !--- TOOLS ASSOCIATED TO ISOTOPES DESCRIPTORS
  PUBLIC :: getKey_init, getKey, fGetKey, setDirectKeys              !--- GET/SET KEYS FROM/TO tracers & isotopes

  PUBLIC :: addPhase, new2oldName,  getPhase, &                      !--- FUNCTIONS RELATED TO THE PHASES
            delPhase, old2newName, getiPhase, &                      !--- + ASSOCIATED VARIABLES
            known_phases, old_phases, phases_sep, phases_names, nphases

  PUBLIC :: oldH2OIso, newH2OIso                                     !--- NEEDED FOR BACKWARD COMPATIBILITY (OLD traceur.def)

  PUBLIC :: tran0, idxAncestor, ancestor                             !--- GENERATION 0 TRACER + TOOLS FOR GENERATIONS
!------------------------------------------------------------------------------------------------------------------------------
  TYPE :: dataBase_type                                              !=== TYPE FOR TRACERS SECTION
    CHARACTER(LEN=maxlen)  :: name                                   !--- Section name
    TYPE(trac_type), ALLOCATABLE :: trac(:)                          !--- Tracers descriptors
  END TYPE dataBase_type
!------------------------------------------------------------------------------------------------------------------------------
  INTERFACE getKey
    MODULE PROCEDURE getKeyByName_s1, getKeyByName_i1, getKeyByName_r1, getKeyByName_sm, getKeyByName_im, getKeyByName_rm
  END INTERFACE getKey
!------------------------------------------------------------------------------------------------------------------------------
  INTERFACE fGetKey;       MODULE PROCEDURE fgetKeyByIndex_s1, fgetKeyByName_s1; END INTERFACE fGetKey
  INTERFACE tracersSubset; MODULE PROCEDURE trSubset_Indx, trSubset_Name, trSubset_gen0Name; END INTERFACE tracersSubset
  INTERFACE idxAncestor;   MODULE PROCEDURE idxAncestor_1, idxAncestor_m; END INTERFACE idxAncestor
  INTERFACE    ancestor;   MODULE PROCEDURE    ancestor_1,    ancestor_m; END INTERFACE    ancestor
  INTERFACE    addPhase;   MODULE PROCEDURE   addPhase_s1,   addPhase_sm,   addPhase_i1,   addPhase_im; END INTERFACE addPhase
  INTERFACE old2newName;   MODULE PROCEDURE old2newName_1, old2newName_m; END INTERFACE old2newName
  INTERFACE new2oldName;   MODULE PROCEDURE new2oldName_1, new2oldName_m; END INTERFACE new2oldName
!------------------------------------------------------------------------------------------------------------------------------

  !=== 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   = 'vlirb'          !--- Old phases for water (no separator)
  CHARACTER(LEN=maxlen), PARAMETER :: known_phases = 'glsrb'          !--- 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  ', 'blosnow']
  CHARACTER(LEN=1), SAVE :: phases_sep  =  '_'                       !--- Phase separator
  LOGICAL,          SAVE :: tracs_merge = .TRUE.                     !--- Merge/stack tracers lists
  LOGICAL,          SAVE :: lSortByGen  = .TRUE.                     !--- Sort by growing generation

  !--- KEPT JUST TO MANAGE OLD WATER ISOTOPES NAMES
  !--- Apart from that context, on limitaion on isotopes names (as long as they have a corresponding line in isotopes_params.def)
  CHARACTER(LEN=maxlen), SAVE :: oldH2OIso(5) = ['eau',     'HDO',     'O18',     'O17',     'HTO'    ]
  CHARACTER(LEN=maxlen), SAVE :: newH2OIso(5) = ['H2[16]O', 'H[2]HO ', 'H2[18]O', 'H2[17]O', 'H[3]HO ']


  !=== LOCAL COPIES OF THE TRACERS AND ISOTOPES DESCRIPTORS, USED BY getKey (INITIALIZED IN getKey_init)
  TYPE(trac_type), ALLOCATABLE, TARGET, SAVE ::  tracers(:)
  TYPE(isot_type), ALLOCATABLE, TARGET, SAVE :: isotopes(:)

  CHARACTER(LEN=maxlen) :: modname

CONTAINS

!==============================================================================================================================
!==============================================================================================================================
!=== READ ONE OR SEVERAL TRACER FILES AND FILL A "tr" TRACERS DESCRIPTOR DERIVED TYPE.
!=== THE RETURN VALUE fType DEPENDS ON WHAT IS FOUND:
!===  0: NO ADEQUATE FILE FOUND ; DEFAULT VALUES MUST BE USED
!===  1: AN "OLD STYLE" TRACERS FILE "traceur.def":
!===    First line: <nb tracers>     Other lines: <hadv> <vadv> <tracer name> [<parent name>]
!===  2: A  "NEW STYLE" TRACERS FILE  "tracer.def" WITH SEVERAL SECTIONS.
!===  3: SEVERAL  "  "  TRACERS FILES "tracer_<component>.def" WITH A SINGLE SECTION IN EACH.
!=== REMARKS:
!===  * EACH SECTION BEGINS WITH A "&<section name> LINE
!===  * DEFAULT VALUES FOR ALL THE SECTIONS OF THE FILE ARE DEFINED IN THE SPECIAL SECTION "&default"
!===  * EACH SECTION LINE HAS THE STRUCTURE:  <name(s)>  <key1>=<value1> <key2>=<value2> ...
!===  * SO FAR, THE DEFINED KEYS ARE: parent, phases, hadv, vadv, type
!===  * <name> AND <parent> CAN BE LISTS OF COMA-SEPARATED TRACERS ; THE ROUTINE EXPAND THESE FACTORIZATIONS.
!=== FUNCTION RETURN VALUE "lerr" IS FALSE IN CASE SOMETHING WENT WRONG.
!=== ABOUT THE KEYS:
!     * The "keys" component (of type 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, fType, tracs) RESULT(lerr)
!------------------------------------------------------------------------------------------------------------------------------
  CHARACTER(LEN=*),             INTENT(IN)  :: type_trac              !--- List of components used
  INTEGER,                      INTENT(OUT) :: fType                  !--- Type of input file found
  TYPE(trac_type), ALLOCATABLE, INTENT(OUT) :: tracs(:)
  CHARACTER(LEN=maxlen),  ALLOCATABLE :: s(:), sections(:), trac_files(:)
  CHARACTER(LEN=maxlen) :: str, fname, mesg
  INTEGER               :: is, nsec, ierr, it, ntrac, ns, ip, ix
  LOGICAL, ALLOCATABLE  :: ll(:), lGen3(:)
!------------------------------------------------------------------------------------------------------------------------------
  lerr = .FALSE.
  modname = 'readTracersFiles'
  IF(.NOT.ALLOCATED(dBase)) ALLOCATE(dBase(0))

  !--- Required sections + corresponding files names (new style single section case)
  IF(test(strParse(type_trac, '|', sections), lerr)) RETURN          !--- Parse "type_trac" list

  nsec = SIZE(sections, DIM=1)
  ALLOCATE(trac_files(nsec)); DO is=1, nsec; trac_files(is) = 'tracer_'//TRIM(sections(is))//'.def'; END DO

  !--- LOOK AT AVAILABLE FILES
  ll = .NOT.testFile(trac_files)
  fType = 0
  IF(.NOT.testFile('traceur.def') .AND. nsec==1) fType = 1           !--- OLD STYLE FILE
  IF(.NOT.testFile('tracer.def'))                fType = 2           !--- NEW STYLE ; SINGLE  FILE, SEVERAL SECTIONS
  IF(ALL(ll))                                    fType = 3           !--- NEW STYLE ; SEVERAL FILES, SINGLE SECTION USED
  IF(ANY(ll) .AND. fType/=3) THEN                                    !--- MISSING FILES
    IF(test(checkList(trac_files, .NOT.ll, 'Failed reading tracers description', 'files', 'missing'), lerr)) RETURN
  END IF

  !--- CHECK WHETHER type_trac AND FILE TYPE ARE COMPATIBLE
  IF(test(fmsg('No multiple sections for the old format "traceur.def"', ll = SIZE(sections)>1 .AND. fType==1), lerr)) RETURN

  !--- TELLS WHAT WAS IS ABOUT TO BE USED
  IF (fmsg('No adequate tracers description file(s) found ; default values will be used',          modname, fType==0)) RETURN
  CALL msg('Trying to read old-style tracers description file "traceur.def"',                      modname, fType==1)
  CALL msg('Trying to read the new style multi-sections tracers description file "tracer.def"',    modname, fType==2)
  CALL msg('Trying to read the new style single section tracers description files "tracer_*.def"', modname, fType==3)

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  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>]
      ALLOCATE(tracs(ntrac))
      DO it=1,ntrac                                                  !=== READ RAW DATA: loop on the line/tracer number
        READ(90,'(a)',IOSTAT=ierr) str
        IF(test(fmsg('Invalid format for "' //TRIM(fname)//'"', modname, ierr>0), lerr)) RETURN
        IF(test(fmsg('Not enough lines in "'//TRIM(fname)//'"', modname, ierr<0), lerr)) RETURN
        ll = strParse(str, ' ', s, n=ns)
        CALL msg('This file is for air tracers only',           modname, ns == 3 .AND. it == 1)
        CALL msg('This files specifies the transporting fluid', modname, ns == 4 .AND. it == 1)
        tracs(it)%name   = old2newName(s(3), ip)                     !--- Set %name:   name   of the tracer
        tracs(it)%parent = tran0                                     !--- Default transporting fluid name
        IF(ns == 4) tracs(it)%parent = old2newName(s(4))             !--- Set %parent: parent of the tracer
        tracs(it)%phase = known_phases(ip:ip)                        !--- Set %phase:  tracer phase (default: "g"azeous)
        tracs(it)%component = TRIM(type_trac)                        !--- Set %component: model component name
        tracs(it)%keys%key = ['hadv', 'vadv']                        !--- Set %keys%key
        tracs(it)%keys%val = s(1:2)                                  !--- Set %keys%val
      END DO
      CLOSE(90)
      CALL setGeneration(tracs)                                      !--- Set %iGeneration and %gen0Name
      WHERE(tracs%iGeneration == 2) tracs%type = 'tag'               !--- Set %type:        'tracer' or 'tag'
      IF(test(checkTracers(tracs, fname, fname), lerr)) RETURN       !--- Detect orphans and check phases
      IF(test(checkUnique (tracs, fname, fname), lerr)) RETURN       !--- Detect repeated tracers
      CALL sortTracers  (tracs)                                      !--- Sort the tracers
      tracs(:)%keys%name = tracs(:)%name                             !--- Copy tracers names in keys components
    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    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; 
    tracs = dBase(1)%trac
  ELSE IF(tracs_merge) THEN
    CALL msg('The multiple required sections will be MERGED.',    modname)
    IF(test(mergeTracers(dBase, tracs), lerr)) RETURN
  ELSE
    CALL msg('The multiple required sections will be CUMULATED.', modname)
    IF(test(cumulTracers(dBase, tracs), lerr)) RETURN
  END IF
  WHERE(tracs%gen0Name(1:3) /= 'H2O') tracs%isInPhysics=.TRUE.       !--- Set %isInPhysics: passed to physics
  CALL setDirectKeys(tracs)                                          !--- Set %iqParent, %iqDescen, %nqDescen, %nqChilds
END FUNCTION readTracersFiles
!==============================================================================================================================

!==============================================================================================================================
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
  LOGICAL,  ALLOCATABLE :: lTg(:)                                    !--- Tagging tracers mask
  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
    CALL setGeneration   (dBase(idb)%trac)                           !---                 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, i, j
  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
  CHARACTER(LEN=maxlen) :: secn
  INTEGER               :: ierr, n
!------------------------------------------------------------------------------------------------------------------------------
  IF(.NOT.ALLOCATED(dBase)) ALLOCATE(dBase(0))
  OPEN(90, FILE=fnam, FORM='formatted', STATUS='old')
  DO; READ(90,'(a)', IOSTAT=ierr)str
    IF(ierr    /= 0 ) EXIT                                           !--- Finished: error or end of file
    IF(str(1:1)=='#') CYCLE                                          !--- Skip comments lines
    CALL removeComment(str)                                          !--- Skip comments at the end of a line
    IF(str     == '') CYCLE                                          !--- Skip empty line (probably at the end of the file)
    IF(str(1:1)=='&') THEN                                           !=== SECTION HEADER LINE
      ndb  = SIZE(dBase)                                             !--- Number of sections so far
      secn = str(2:LEN_TRIM(str))//' '                               !--- Current section name
      IF(ANY(dBase(:)%name == secn)) CYCLE                           !--- Already known section
      IF(secn(1:7) == 'version')     CYCLE                           !--- Skip the "version" special section
      ndb = ndb + 1                                                  !--- Extend database
      ALLOCATE(tdb(ndb))
      tdb(1:ndb-1)  = dBase
      tdb(ndb)%name = secn
      ALLOCATE(tdb(ndb)%trac(0))
      CALL MOVE_ALLOC(FROM=tdb, TO=dBase)
    ELSE                                                             !=== TRACER LINE
      ll = strParse(str,' ', keys = s, vals = v, n = n)              !--- Parse <key>=<val> pairs
      tt = dBase(ndb)%trac(:)
      tmp%name = s(1); tmp%keys = 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)                   !--- Add key to all the tracers (no overwriting)
    DO it = 1, SIZE(t); CALL addKey_1(ky%key(k), ky%val(k), t(it)%keys); END DO
  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, i
  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

    !--- 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
  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  CALL delKey(['parent','type  '], tr)

  ALLOCATE(ttr(nq))
  iq = 1
  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  DO it = 1, nt                                                      !=== EXPAND TRACERS AND PARENTS NAMES LISTS
  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    ll = strParse(tr(it)%name,   ',', ta, n=ntr)                     !--- Number of tracers
    ll = strParse(tr(it)%parent, ',', pa, n=npr)                     !--- Number of parents
    DO ipr=1,npr                                                     !--- Loop on parents list elts
      DO itr=1,ntr                                                   !--- Loop on tracers list elts
        i = iq+itr-1+(ipr-1)*ntr
        ttr(i)%name   = TRIM(ta(itr))
        ttr(i)%parent = TRIM(pa(ipr))
        ttr(i)%keys%name = ta(itr)
        ttr(i)%keys%key  = tr(it)%keys%key
        ttr(i)%keys%val  = tr(it)%keys%val
!        ttr(i)%keys   = keys_type(ta(itr), tr(it)%keys%key, tr(it)%keys%val)
      END DO
    END DO
    ttr(iq:iq+ntr*npr-1)%type      = tr(it)%type                     !--- Duplicating type
    ttr(iq:iq+ntr*npr-1)%component = tr(it)%component                !--- Duplicating type
    iq = iq + ntr*npr
  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  END DO
  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  DEALLOCATE(ta,pa)
  CALL MOVE_ALLOC(FROM=ttr, TO=tr)

END FUNCTION expandSection
!==============================================================================================================================

!==============================================================================================================================
SUBROUTINE setGeneration(tr)
!------------------------------------------------------------------------------------------------------------------------------
! Purpose: Determine, for each tracer of "tr(:)":
!   * %iGeneration: the generation number 
!   * %gen0Name:    the generation 0 ancestor name
!------------------------------------------------------------------------------------------------------------------------------
  TYPE(trac_type),     INTENT(INOUT) :: tr(:)                        !--- Tracer derived type vector
  INTEGER                            :: iq, nq, ig
  LOGICAL,               ALLOCATABLE :: lg(:)
  CHARACTER(LEN=maxlen), ALLOCATABLE :: prn(:)
!------------------------------------------------------------------------------------------------------------------------------
  tr(:)%iGeneration = -1                                             !--- error if -1
  nq = SIZE(tr, DIM=1)                                               !--- Number of tracers lines
  lg = tr(:)%parent == tran0                                         !--- Flag for generation 0 tracers
  WHERE(lg) tr(:)%iGeneration = 0                                    !--- Generation 0 tracers

  !=== Determine generation for each tracer
  ig=-1; prn = [tran0]
  DO                                                                 !--- Update current generation flag
    IF(ig/=-1) prn = PACK( tr(:)%name, MASK=tr(:)%iGeneration == ig)
    lg(:) = [(ANY(prn(:) == tr(iq)%parent), iq=1, nq)]               !--- Current generation tracers flag
    IF( ALL( .NOT. lg ) ) EXIT                                       !--- Empty current generation
    ig = ig+1; WHERE(lg) tr(:)%iGeneration = ig
  END DO
  tr(:)%gen0Name = ancestor(tr)                                      !--- First generation ancestor name

END SUBROUTINE setGeneration
!==============================================================================================================================

!==============================================================================================================================
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, trn
  CHARACTER(LEN=1) :: p
  INTEGER :: ip, np, iq, jq, nq, it, nt, nc, i, n
  LOGICAL :: lTg, lEx
!------------------------------------------------------------------------------------------------------------------------------
  nq = SIZE(tr, DIM=1)
  nt = 0
  DO iq = 1, nq                                                      !--- GET THE NUMBER OF TRACERS
    IF(tr(iq)%iGeneration /= 0) CYCLE                                !--- Only deal with generation 0 tracers
    nc = COUNT(tr(:)%gen0Name==tr(iq)%name .AND. tr%iGeneration/=0)  !--- Number of childs of tr(iq)
    tr(iq)%phase = fgetKey(iq, 'phases', tr(:)%keys)                 !--- Phases list      of tr(iq)
    np = LEN_TRIM(tr(iq)%phase)                                      !--- Number of phases of tr(iq)
    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
    lTg = tr(iq)%type=='tag'                                         !--- Current tracer is a tag
    i0 = strFind(tr(:)%name, TRIM(tr(iq)%gen0Name), n)               !--- Indexes of first generation ancestor copies
    np = SUM([( LEN_TRIM(tr(i0(i))%phase),i=1,n )], 1)               !--- Number of phases for current tracer tr(iq)
    lEx = np>1                                                       !--- Phase suffix only required if phases number is > 1
    IF(lTg) lEx = lEx .AND. tr(iq)%iGeneration>0                     !--- No phase suffix for generation 0 tags
    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)
        trn = TRIM(tr(iq)%name); nam = trn                           !--- Tracer name (regular case)
        IF(lTg) nam = TRIM(tr(iq)%parent)                            !--- Parent name (tagging case)
        IF(lEx) nam = addPhase(nam, p )                              !--- Phase extension needed
        IF(lTg) nam = TRIM(nam)//'_'//TRIM(trn)                      !--- <parent>_<name> for tags
        ttr(it) = tr(iq)                                             !--- Same <key>=<val> pairs
        ttr(it)%name      = TRIM(nam)                                !--- Name with possibly phase suffix
        ttr(it)%keys%name = TRIM(nam)                                !--- Name inside the keys decriptor
        ttr(it)%phase     = p                                        !--- Single phase entry
        IF(lEx .AND. tr(iq)%iGeneration>0) THEN
          ttr(it)%parent   = addPhase(ttr(it)%parent,   p)
          ttr(it)%gen0Name = addPhase(ttr(it)%gen0Name, p)
        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 childs sorted by growing generation, one after the other.
!   TO BE ADDED IF NECESSARY: HIGHER MOMENTS AT THE END
!------------------------------------------------------------------------------------------------------------------------------
  TYPE(trac_type), ALLOCATABLE, INTENT(INOUT) :: tr(:)         !--- Tracer derived type vector
  INTEGER,         ALLOCATABLE :: iy(:), iz(:)
  INTEGER :: ig, ng, iq, jq, ip, nq, n, ix(SIZE(tr)), k
!------------------------------------------------------------------------------------------------------------------------------
  nq = SIZE(tr)
  DO ip = nphases, 1, -1
    iq = strIdx(tr(:)%name, addPhase('H2O', ip))
    IF(iq == 0) CYCLE
    tr = [tr(iq), tr(1:iq-1), tr(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)" childs 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     :: t1(:), t2(:)
  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
  !----------------------------------------------------------------------------------------------------------------------------
    t1 => sections(is)%trac(:)
    !--------------------------------------------------------------------------------------------------------------------------
    DO iq=1, nt(is)                                                  !=== LOOP ON TRACERS
    !--------------------------------------------------------------------------------------------------------------------------
      tnam = TRIM(t1(iq)%name)                                       !--- Original name
      IF(COUNT(t1%name == tnam) == 1) CYCLE                          !--- Current tracer is not duplicated: finished
      tnam_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, %nqChilds
  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(:)
  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([(fgetKey(iq, 'hadv',  tm(:)%keys, '10'), iq=1, nq)])
  vadv = str2int([(fgetKey(iq, 'vadv',  tm(:)%keys, '10'), iq=1, nq)])
  phas =         [(fgetKey(iq, 'phases',tm(:)%keys, 'g' ), iq=1, nq)]
  CALL msg(TRIM(message)//':', modname)
  IF(tm(1)%parent == '') THEN
    IF(test(dispTable('iiiss',   ['iq    ','hadv  ','vadv  ','name  ','phase '], cat(tm%name, phas), cat([(iq, iq=1, nq)], &
                                            hadv,    vadv),                 sub=modname), lerr)) RETURN
  ELSE
    IF(test(dispTable('iiissis', ['iq    ','hadv  ','vadv  ','name  ','parent','igen  ','phase '], cat(tm%name, tm%parent, &
           tm%phase), cat([(iq, iq=1, nq)], hadv,    vadv, tm%iGeneration), 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, ng, igen, ngen
  INTEGER, ALLOCATABLE :: ix(:)
  tr(:)%iqParent = strIdx( tr(:)%name, tr(:)%parent )                !--- Parent index
  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))
    ix = idxAncestor(tr, igen=ig)                                    !--- Ancestor of generation "ng" for each tr
    DO igen = ig+1, ngen
      tr(iq)%iqDescen = [tr(iq)%iqDescen, find(ix==iq .AND. tr%iGeneration==igen)]
      tr(iq)%nqDescen = SIZE(tr(iq)%iqDescen)
      IF(igen == ig+1) tr(iq)%nqChilds=tr(iq)%nqDescen
    END DO
  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 (initIsotopes):                             ====
!===      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(fnam, isot) RESULT(lerr)
  CHARACTER(LEN=*),        INTENT(IN)    :: fnam                     !--- Input file name
  TYPE(isot_type), TARGET, INTENT(INOUT) :: isot(:)                  !--- Isotopes descriptors (field "prnt" must be defined !)
  INTEGER :: ik, is, it, idb, nk0, i, iis
  INTEGER :: nk, ns, nt, ndb, nb0, i0
  CHARACTER(LEN=maxlen), POINTER     :: k(:), v(:), k0(:), v0(:)
  CHARACTER(LEN=maxlen), ALLOCATABLE :: vals(:)
  CHARACTER(LEN=maxlen)              :: val, modname
  TYPE(keys_type),           POINTER ::   ky(:)
  TYPE(trac_type),           POINTER ::   tt(:), t
  TYPE(dataBase_type),   ALLOCATABLE ::  tdb(:)
  LOGICAL,               ALLOCATABLE :: liso(:)
  modname = 'readIsotopesFile'

  !--- THE INPUT FILE MUST BE PRESENT
  IF(test(fmsg('Missing isotopes parameters file "'//TRIM(fnam)//'"', modname, testFile(fnam)),lerr)) RETURN

  !--- 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
      liso = reduceExpr(t%keys%val, vals)                            !--- Reduce expressions (for substituted variables)
      IF(test(ANY(liso), lerr)) RETURN                               !--- Some non-numerical elements were found
      isot(iis)%keys(is)%key = PACK(t%keys%key, MASK=.NOT.liso)
      isot(iis)%keys(is)%val = PACK(  vals,     MASK=.NOT.liso)
    END DO

    !--- CHECK FOR MISSING ISOTOPES (NO KEYS ALLOCATED)
    liso = [( ALLOCATED(isot(iis)%keys(is)%key), is=1, SIZE(isot(iis)%keys) )]
    IF(test(checkList(isot(iis)%keys(:)%name, .NOT.liso, &
      'Check file "'//TRIM(fnam)//'" in section "'//TRIM(dBase(idb)%name)//'"', 'isotopes', 'missing'),lerr)) RETURN
  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(isot, 'Isotopes parameters read from file', modname)

END FUNCTION readIsotopesFile
!==============================================================================================================================

!==============================================================================================================================
!=== IF ISOTOPES (2ND GENERATION TRACERS) ARE DETECTED:                                                                     ===
!===    * COMPUTE MOST OF THE RELATED QUANTITIES ("isot" COMPONENTS).                                                       ===
!===    * COMPUTE FEW ISOTOPES-DEDICATED "trac" COMPONENTS                                                                  ===
!===    * CALL readIsotopesFile 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.  ===
!==============================================================================================================================
SUBROUTINE initIsotopes(trac, isot)
  TYPE(trac_type), ALLOCATABLE, TARGET, INTENT(INOUT) :: trac(:)
  TYPE(isot_type), ALLOCATABLE, TARGET, INTENT(INOUT) :: isot(:)
  CHARACTER(LEN=maxlen), ALLOCATABLE :: p(:), str(:)                 !--- Temporary storage
  CHARACTER(LEN=maxlen) :: iname
  CHARACTER(LEN=1)   :: ph                                           !--- Phase
  INTEGER :: nbIso, ic, ip, iq, it, iz
  LOGICAL, ALLOCATABLE :: ll(:)                                      !--- Mask
  TYPE(trac_type), POINTER   ::  t(:), t1
  TYPE(isot_type), POINTER   ::  i

  t => trac

  p = PACK(delPhase(t%parent), MASK = t%type=='tracer' .AND. t%iGeneration==1) !--- Parents of generation 1 isotopes
  CALL strReduce(p, nbIso)
  ALLOCATE(isot(nbIso))

  IF(nbIso==0) RETURN                                                !=== NO ISOTOPES: FINISHED

  !--- ISOTOPES RELATED VARIABLES ; NULL OR EMPTY IF NO ISOTOPES
  isot(:)%parent = p
  DO ic = 1, SIZE(p)                                                 !--- Loop on isotopes classes
    i => isot(ic)
    iname = i%parent                                                 !--- Current isotopes class name (parent tracer name)

    !=== Isotopes childs of tracer "iname": mask, names, number (same for each phase of "iname")
    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 childs 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 => trac(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 "isotopes_params.def" FILE
  !    DONE HERE, AND NOT ONLY IN "infotrac_phy", BECAUSE SOME PHYSICAL PARAMS ARE NEEDED FOR RESTARTS (tnat AND alpha_ideal)
  IF(readIsotopesFile('isotopes_params.def',isot)) THEN 
     STOP 'Problem when reading isotopes parameters in initIsotopes' 
  ENDIF

END SUBROUTINE initIsotopes
!==============================================================================================================================


!==============================================================================================================================
LOGICAL FUNCTION dispIsotopes(ides, message, modname) RESULT(lerr)
  TYPE(isot_type),  INTENT(IN) :: ides(:)                            !--- Isotopes descriptor vector
  CHARACTER(LEN=*), INTENT(IN) :: message                            !--- Message to display
  CHARACTER(LEN=*), INTENT(IN) :: modname                            !--- Calling subroutine name
  INTEGER :: ik, nk, ip, it, nt
  CHARACTER(LEN=maxlen) :: prf
  CHARACTER(LEN=maxlen), ALLOCATABLE :: ttl(:), val(:,:)
  CALL msg(TRIM(message)//':', modname)
  DO ip = 1, SIZE(ides)                                              !--- Loop on parents tracers
    nk = SIZE(ides(ip)%keys(1)%key)                                  !--- Same keys for each isotope
    nt = SIZE(ides(ip)%keys)                                         !--- Number of isotopes
    prf = 'i'//REPEAT('s',nk+1)                                      !--- Profile for table printing
    ALLOCATE(ttl(nk+2), val(nt,nk+1))
    ttl(1:2) = ['iq  ','name']; ttl(3:nk+2) = ides(ip)%keys(1)%key(:)!--- Titles line with keys names
    val(:,1) = ides(ip)%keys(:)%name                                 !--- Values table 1st column: isotopes names  
    DO ik = 1, nk
      DO it = 1, nt
        val(it,ik+1) = ides(ip)%keys(it)%val(ik)                     !--- Other columns: keys values
      END DO
    END DO
    IF(test(fmsg('Problem with the table content', modname, dispTable(prf, ttl, val, cat([(it,it=1,nt)]), rFmt='(EN8.4)',     &
       sub=modname)), lerr)) RETURN
    DEALLOCATE(ttl, val)
  END DO        
END FUNCTION dispIsotopes
!==============================================================================================================================


!==============================================================================================================================
SUBROUTINE addKey_1(key, val, ky, lOverWrite)
!------------------------------------------------------------------------------------------------------------------------------
! Purpose: Add the <key>=<val> pair in the "ky" keys descriptor.
!------------------------------------------------------------------------------------------------------------------------------
  CHARACTER(LEN=*),  INTENT(IN)    :: key, val
  TYPE(keys_type),   INTENT(INOUT) :: ky
  LOGICAL, OPTIONAL, INTENT(IN)    :: lOverWrite
  CHARACTER(LEN=maxlen), ALLOCATABLE :: k(:), v(:)
  INTEGER :: iky, nky
  LOGICAL :: lo
!------------------------------------------------------------------------------------------------------------------------------
  lo=.FALSE.; IF(PRESENT(lOverWrite)) lo=lOverWrite
  iky = strIdx(ky%key,key)
  IF(iky == 0) THEN
    nky = SIZE(ky%key)
    IF(nky == 0) THEN; ky%key = [key]; ky%val = [val]; ELSE; ky%key = [ky%key, key]; ky%val = [ky%val, val]; END IF
  ELSE IF(lo) THEN                                                   !--- Overwriting
    ky%key(iky) = key; ky%val(iky) = val
  END IF
END SUBROUTINE addKey_1
!==============================================================================================================================
SUBROUTINE addKey_m(key, val, ky, lOverWrite)
!------------------------------------------------------------------------------------------------------------------------------
! Purpose: Add the <key>=<val> pair in all the components of the "ky" keys descriptor.
!------------------------------------------------------------------------------------------------------------------------------
  CHARACTER(LEN=*),  INTENT(IN)    :: key, val
  TYPE(keys_type),   INTENT(INOUT) :: ky(:)
  LOGICAL, OPTIONAL, INTENT(IN)    :: lOverWrite
  INTEGER :: itr
  LOGICAL :: lo
!------------------------------------------------------------------------------------------------------------------------------
  lo=.FALSE.; IF(PRESENT(lOverWrite)) lo=lOverWrite
  DO itr = 1, SIZE(ky); CALL addKey_1(key, val, ky(itr), lo); END DO
END SUBROUTINE addKey_m
!==============================================================================================================================
SUBROUTINE addKeysFromDef(t, tr0)
!------------------------------------------------------------------------------------------------------------------------------
! Purpose: The values of the keys of the tracer named "tr0" are overwritten by the values found in the *.def files, if any.
!------------------------------------------------------------------------------------------------------------------------------
  TYPE(trac_type), ALLOCATABLE, INTENT(INOUT) :: t(:)
  CHARACTER(LEN=*),             INTENT(IN)    :: tr0
  CHARACTER(LEN=maxlen) :: val
  INTEGER               :: ik, jd
  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
!==============================================================================================================================
SUBROUTINE delKey_1(itr, keyn, ky)
!------------------------------------------------------------------------------------------------------------------------------
! Purpose: Internal routine.
!   Remove <key>=<val> pairs in the "itr"th component of the "ky" keys descriptor.
!------------------------------------------------------------------------------------------------------------------------------
  INTEGER,          INTENT(IN)    :: itr
  CHARACTER(LEN=*), INTENT(IN)    :: keyn(:)
  TYPE(trac_type),  INTENT(INOUT) :: ky(:)
  CHARACTER(LEN=maxlen), ALLOCATABLE :: k(:), v(:)
  LOGICAL,               ALLOCATABLE :: ll(:)
  INTEGER :: iky
!------------------------------------------------------------------------------------------------------------------------------
  IF(itr<=0 .OR. itr>SIZE(ky, DIM=1)) RETURN                          !--- Index is out of range
  ll = [( ALL(keyn/=ky(itr)%keys%key(iky)), iky=1, SIZE(ky(itr)%keys%key) )]
  k = PACK(ky(itr)%keys%key, MASK=ll); CALL MOVE_ALLOC(FROM=k, TO=ky(itr)%keys%key)
  v = PACK(ky(itr)%keys%val, MASK=ll); CALL MOVE_ALLOC(FROM=v, TO=ky(itr)%keys%val)
END SUBROUTINE delKey_1
!==============================================================================================================================
SUBROUTINE delKey(keyn, ky)
!------------------------------------------------------------------------------------------------------------------------------
! Purpose: Internal routine.
!   Remove <key>=<val> pairs in all the components of the "t" tracers descriptor.
!------------------------------------------------------------------------------------------------------------------------------
  CHARACTER(LEN=*), INTENT(IN)    :: keyn(:)
  TYPE(trac_type),  INTENT(INOUT) :: ky(:)
  INTEGER :: iky
!------------------------------------------------------------------------------------------------------------------------------
  DO iky = 1, SIZE(ky); CALL delKey_1(iky, keyn, ky); END DO
END SUBROUTINE delKey
!==============================================================================================================================


!==============================================================================================================================
!=== PUBLIC ROUTINES: GET A KEY FROM A <key>=<val> LIST ; VECTORS, TRACER AND DATABASE VERSIONS ===============================
!=== BEWARE !!! IF THE "ky" ARGUMENT IS NOT PRESENT, THEN THE VARIABLES "tracers" AND "isotopes" ARE USED. ====================
!===     THEY ARE LOCAL TO THIS MODULE, SO MUST MUST BE INITIALIZED FIRST USING the "getKey_init" ROUTINE  ====================
!==============================================================================================================================
SUBROUTINE getKey_init(tracers_, isotopes_)
  TYPE(trac_type), OPTIONAL, INTENT(IN) ::  tracers_(:)
  TYPE(isot_type), OPTIONAL, INTENT(IN) :: isotopes_(:)
  IF(PRESENT( tracers_))  tracers =  tracers_
  IF(PRESENT(isotopes_)) isotopes = isotopes_
END SUBROUTINE getKey_init
!==============================================================================================================================
CHARACTER(LEN=maxlen) FUNCTION fgetKeyByIndex_s1(itr, keyn, ky, def_val) RESULT(val)
!------------------------------------------------------------------------------------------------------------------------------
! Purpose: Internal function ; get a string formatted key value (returned argument) from its name and the tracer index.
!------------------------------------------------------------------------------------------------------------------------------
  INTEGER,                    INTENT(IN) :: itr
  CHARACTER(LEN=*),           INTENT(IN) :: keyn
  TYPE(keys_type),            INTENT(IN) :: ky(:)
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: def_val
!------------------------------------------------------------------------------------------------------------------------------
  INTEGER :: iky
  iky = 0;  IF(itr >  0 .AND. itr <= SIZE(ky)) iky = strIdx(ky(itr)%key(:), keyn)
  val = ''; IF(iky /= 0) val = ky(itr)%val(iky)                      !--- Key was found
  IF(PRESENT(def_val) .AND. iky == 0) val = def_val                  !--- Default value from arguments
END FUNCTION fgetKeyByIndex_s1
!==============================================================================================================================
CHARACTER(LEN=maxlen) FUNCTION fgetKeyByName_s1(tname, keyn, ky, def_val, lerr) RESULT(val)
!------------------------------------------------------------------------------------------------------------------------------
! Purpose: Internal function ; get a string formatted key value (returned argument) from its name and the tracer name.
!------------------------------------------------------------------------------------------------------------------------------
  CHARACTER(LEN=*),           INTENT(IN)  :: tname, keyn
  TYPE(keys_type),            INTENT(IN)  :: ky(:)
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: def_val
  LOGICAL,          OPTIONAL, INTENT(OUT) :: lerr
!------------------------------------------------------------------------------------------------------------------------------
  INTEGER :: iky, itr
  val = ''; iky = 0
  itr = strIdx(ky(:)%name, tname)                                    !--- Get the index of the wanted tracer
  IF(PRESENT(lerr)) lerr = itr==0; IF(itr == 0) RETURN
  IF(itr >  0 .AND. itr <= SIZE(ky)) iky = strIdx(ky(itr)%key(:), keyn)
  IF(iky /= 0) val = ky(itr)%val(iky)                                !--- Key was found
  IF(PRESENT(def_val) .AND. iky == 0) val = def_val                  !--- Default value from arguments
END FUNCTION fgetKeyByName_s1
!==============================================================================================================================
LOGICAL FUNCTION getKeyByName_s1(keyn, val, tname, ky) RESULT(lerr)
  !--- Purpose: Get the value of the key named "keyn" for the tracer named "tnam".
  !     * "ky" unspecified: try in "tracers" for "tnam" with phase and tagging suffixes, then in "isotopes" without.
  !     * "ky"   specified: try in "ky"      for "tnam" with phase and tagging suffixes, then without.
  !    The returned error code is always .FALSE.: an empty string is returned when the key hasn't been found.
  CHARACTER(LEN=*),          INTENT(IN)  :: keyn
  CHARACTER(LEN=maxlen),     INTENT(OUT) :: val
  CHARACTER(LEN=*),          INTENT(IN)  :: tname
  TYPE(keys_type), OPTIONAL, INTENT(IN)  :: ky(:)
  CHARACTER(LEN=maxlen) :: tnam
  INTEGER, ALLOCATABLE  :: is(:)
  INTEGER :: i, itr
  tnam = delPhase(strHead(tname,'_',.FALSE.))                        !--- Remove tag and phase
  IF(PRESENT(ky)) THEN
    val = fgetKeyByName_s1(tname, keyn, ky, lerr=lerr)               !--- "ky" and "tname"
    IF(val /= '' .OR. lerr)      RETURN
    val = fgetKeyByName_s1(tnam,  keyn, ky, lerr=lerr)               !--- "ky" and "tnam"
  ELSE
    IF(.NOT.ALLOCATED(tracers))  RETURN
    val = fgetKeyByName_s1(tname, keyn, tracers(:)%keys, lerr=lerr)  !--- "tracers" and "tname"
    IF(val /= ''.AND..NOT.lerr)  RETURN
    IF(.NOT.ALLOCATED(isotopes)) RETURN
    IF(SIZE(isotopes) == 0)      RETURN
    !--- Search the "is" isotopes class index of the isotope named "tnam"
    is = find([(ANY(isotopes(i)%keys(:)%name == tnam), i=1, SIZE(isotopes))])
    IF(test(SIZE(is) == 0,lerr)) RETURN
    val = fgetKeyByName_s1(tname, keyn, isotopes(is(1))%keys(:),lerr=lerr)!--- "isotopes" and "tnam"
  END IF
END FUNCTION getKeyByName_s1
!==============================================================================================================================
LOGICAL FUNCTION getKeyByName_sm(keyn, val, tname, ky) RESULT(lerr)
  CHARACTER(LEN=*),                   INTENT(IN)  :: keyn
  CHARACTER(LEN=maxlen), ALLOCATABLE, INTENT(OUT) ::   val(:)
  CHARACTER(LEN=*), TARGET, OPTIONAL, INTENT(IN)  :: tname(:)
  TYPE(keys_type),  TARGET, OPTIONAL, INTENT(IN)  ::    ky(:)
  TYPE(keys_type),           POINTER :: k(:)
  CHARACTER(LEN=maxlen), ALLOCATABLE :: n(:)
  INTEGER :: iq, nq
  IF(test(.NOT.(PRESENT(tname).OR.PRESENT(ky)), lerr)) RETURN
  IF(PRESENT(ky   )) nq = SIZE(ky%name)
  IF(PRESENT(tname)) nq = SIZE(  tname)
  ALLOCATE(val(nq))
  IF(PRESENT(tname)) THEN
    IF(     PRESENT(ky)) lerr = ANY([(getKeyByName_s1(keyn, val(iq), tname(iq),   ky), iq=1, nq)])
    IF(.NOT.PRESENT(ky)) lerr = ANY([(getKeyByName_s1(keyn, val(iq), tname(iq)      ), iq=1, nq)])
  ELSE;                  lerr = ANY([(getKeyByName_s1(keyn, val(iq), ky(iq)%name, ky), iq=1, nq)])
  END IF
END FUNCTION getKeyByName_sm
!==============================================================================================================================
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
  IF(     PRESENT(ky)) lerr = getKeyByName_s1(keyn, sval, tname, ky)
  IF(.NOT.PRESENT(ky)) lerr = getKeyByName_s1(keyn, sval, tname)
  IF(test(fmsg('key "'//TRIM(keyn)//'" or tracer "'//TRIM(tname)//'" is missing',        modname, lerr), lerr)) RETURN
  READ(sval, *, IOSTAT=ierr) val
  IF(test(fmsg('key "'//TRIM(keyn)//'" of tracer "'//TRIM(tname)//'" is not an integer', modname, lerr), lerr)) RETURN
END FUNCTION getKeyByName_i1
!==============================================================================================================================
LOGICAL FUNCTION getKeyByName_im(keyn, val, tname, ky) RESULT(lerr)
  CHARACTER(LEN=*),                   INTENT(IN)  :: keyn
  INTEGER,               ALLOCATABLE, INTENT(OUT) ::   val(:)
  CHARACTER(LEN=*), TARGET, OPTIONAL, INTENT(IN)  :: tname(:)
  TYPE(keys_type),  TARGET, OPTIONAL, INTENT(IN)  ::    ky(:)
  TYPE(keys_type),           POINTER :: k(:)
  CHARACTER(LEN=maxlen), ALLOCATABLE :: n(:)
  INTEGER :: iq, nq
  IF(test(.NOT.(PRESENT(tname).OR.PRESENT(ky)), lerr)) RETURN
  IF(PRESENT(ky   )) nq = SIZE(ky%name)
  IF(PRESENT(tname)) nq = SIZE(  tname)
  ALLOCATE(val(nq))
  IF(PRESENT(tname)) THEN
    IF(     PRESENT(ky)) lerr = ANY([(getKeyByName_i1(keyn, val(iq), tname(iq),   ky), iq=1, nq)])
    IF(.NOT.PRESENT(ky)) lerr = ANY([(getKeyByName_i1(keyn, val(iq), tname(iq)      ), iq=1, nq)])
  ELSE;                  lerr = ANY([(getKeyByName_i1(keyn, val(iq), ky(iq)%name, ky), iq=1, nq)])
  END IF
END FUNCTION getKeyByName_im
!==============================================================================================================================
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
  IF(     PRESENT(ky)) lerr = getKeyByName_s1(keyn, sval, tname, ky)
  IF(.NOT.PRESENT(ky)) lerr = getKeyByName_s1(keyn, sval, tname)
  IF(test(fmsg('key "'//TRIM(keyn)//'" or tracer "'//TRIM(tname)//'" is missing',    modname, lerr), lerr)) RETURN
  READ(sval, *, IOSTAT=ierr) val
  IF(test(fmsg('key "'//TRIM(keyn)//'" of tracer "'//TRIM(tname)//'" is not a real', modname, lerr), lerr)) RETURN
END FUNCTION getKeyByName_r1
!==============================================================================================================================
LOGICAL FUNCTION getKeyByName_rm(keyn, val, tname, ky) RESULT(lerr)
  CHARACTER(LEN=*),                   INTENT(IN)  :: keyn
  REAL,                  ALLOCATABLE, INTENT(OUT) ::   val(:)
  CHARACTER(LEN=*), TARGET, OPTIONAL, INTENT(IN)  :: tname(:)
  TYPE(keys_type),  TARGET, OPTIONAL, INTENT(IN)  ::    ky(:)
  TYPE(keys_type),           POINTER :: k(:)
  CHARACTER(LEN=maxlen), ALLOCATABLE :: n(:)
  INTEGER :: iq, nq
  IF(test(.NOT.(PRESENT(tname).OR.PRESENT(ky)), lerr)) RETURN
  IF(PRESENT(ky   )) nq = SIZE(ky%name)
  IF(PRESENT(tname)) nq = SIZE(  tname)
  ALLOCATE(val(nq))
  IF(PRESENT(tname)) THEN
    IF(     PRESENT(ky)) lerr = ANY([(getKeyByName_r1(keyn, val(iq), tname(iq),   ky), iq=1, nq)])
    IF(.NOT.PRESENT(ky)) lerr = ANY([(getKeyByName_r1(keyn, val(iq), tname(iq)      ), iq=1, nq)])
  ELSE;                  lerr = ANY([(getKeyByName_r1(keyn, val(iq), ky(iq)%name, ky), iq=1, nq)])
  END IF
END FUNCTION getKeyByName_rm
!==============================================================================================================================


!==============================================================================================================================
!=== 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 :: l, i, ix
  CHARACTER(LEN=maxlen) :: sh, st
  out = s
  IF(s == '') RETURN                                                           !--- Empty string: nothing to do

  !--- Special case: old phases for water, no phases separator
  i = INDEX(s,'_'); sh = s; IF(i/=0) sh=s(1:i-1); st='H2O'; IF(i/=0) st='H2O_'//s(i+1:LEN_TRIM(s))
  IF(ANY([('H2O'//old_phases(ix:ix), ix=1, nphases)] == sh)) THEN; out=st; RETURN; END IF

  !--- Index of found phase in "known_phases"
  ix = MAXLOC( [(i, i=1,nphases)], MASK=[( INDEX(s, phases_sep//known_phases(i:i))/=0, i=1, nphases)], DIM=1 )
  IF(ix == 0) RETURN                                                           !--- No phase pattern found
  i = INDEX(s, phases_sep//known_phases(ix:ix))                                !--- Index of <sep><pha> pattern in "str"
  l = LEN_TRIM(s)
  IF(i == l-1) THEN                                                            !--- <var><sep><pha>       => return <var>
    out = s(1:l-2)
  ELSE IF(s(i+2:i+2) == '_') THEN                                              !--- <var><sep><pha>_<tag> => return <var>_<tag>
    out = s(1:i-1)//s(i+2:l)
  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) RETURN                                                           !--- Null 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
!------------------------------------------------------------------------------------------------------------------------------


!------------------------------------------------------------------------------------------------------------------------------
CHARACTER(LEN=maxlen) FUNCTION old2newName_1(oldName, iPhase) RESULT(newName)
  !--- Convert an old style name into a new one.
  !    Only usable with old style "traceur.def" files, in which only water isotopes are allowed.
  !    In these files, H2O descendants names are: H2O<phase>[_<isotope>][_<tag>], with:
  !    phase = v, l or i ; isotope = eau, HDO, O18, O17 or HTO.
  CHARACTER(LEN=*),  INTENT(IN)  :: oldName
  INTEGER, OPTIONAL, INTENT(OUT) :: iPhase
  CHARACTER(LEN=maxlen), ALLOCATABLE :: tmp(:)
  INTEGER :: ix, ip, it, nt
  LOGICAL :: lerr, lH2O
  newName = oldName
  IF(PRESENT(iPhase)) iPhase = 1                                               !--- Default: gaseous phase
  lH2O=.FALSE.
  IF(LEN_TRIM(oldName) > 3) THEN
    lH2O = oldName(1:3)=='H2O' .AND. INDEX(old_phases,oldName(4:4))/=0         !--- H2O<phase>*,  with phase=="v", "l", "i" or "r"
    IF(LEN_TRIM(oldName) > 4) lH2O = lH2O .AND. oldName(5:5) == '_'            !--- H2O<phase>_*, with phase=="v", "l", "i" or "r"
  END IF
  IF(.NOT.lH2O) RETURN
  IF(LEN_TRIM(oldName)>3) THEN; IF(INDEX(old_Phases,oldName(4:4))==0) RETURN; END IF 
  lerr = strParse(oldName, '_', tmp, n=nt)
  ip = strIdx([('H2O'//old_phases(ip:ip), ip=1, nphases)], tmp(1))             !--- Phase index (/=0 if any)
  IF(PRESENT(iPhase)) iPhase = ip
  newName = addPhase('H2O', ip)                                                !--- Water
  IF(nt == 1) RETURN                                                           !--- Water: finished
  ix = strIdx(oldH2OIso, tmp(2))                                               !--- Index in the known isotopes list
  IF(ix == 0) newName = addPhase(tmp(2),        ip)                            !--- Not an isotope
  IF(ix /= 0) newName = addPhase(newH2OIso(ix), ip)                            !--- Isotope
  IF(nt == 3) newName = TRIM(newName)//'_'//TRIM(tmp(3))                       !--- Tagging tracer
END FUNCTION old2newName_1
!------------------------------------------------------------------------------------------------------------------------------
FUNCTION old2newName_m(oldName, iPhase) RESULT(newName)
  CHARACTER(LEN=*),  INTENT(IN)  :: oldName(:)
  INTEGER, OPTIONAL, INTENT(OUT) :: iPhase
  CHARACTER(LEN=maxlen)          :: newName(SIZE(oldName))
  INTEGER :: i
  newName = [(old2newName_1(oldName(i), iPhase), i=1, SIZE(oldName))]
END FUNCTION old2newName_m
!------------------------------------------------------------------------------------------------------------------------------

!------------------------------------------------------------------------------------------------------------------------------
CHARACTER(LEN=maxlen) FUNCTION new2oldName_1(newName, iPhase) RESULT(oldName)
  !--- Convert a new style name into an old one.
  !    Only convertable names are water descendants names H2O_<phase>, <isotope>_<phase>, <isotope>_<phase>_<tag>, with:
  !    phase = g, l or s ; isotope = H2[16]O, H[2]O, H2<[18]O, H2[17]O or H[3]O.
  CHARACTER(LEN=*),  INTENT(IN)    :: newName
  INTEGER, OPTIONAL, INTENT(OUT)   :: iPhase
  INTEGER :: ix, ip, it, nt
  LOGICAL :: lH2O
  CHARACTER(LEN=maxlen) :: tag
  ix = strIdx([(addPhase('H2O',ip), ip=1, nphases)], newName)                  !--- Phase index for H2O_<phase>
  IF(ix /= 0) THEN; oldName = 'H2O'//old_phases(ix:ix); RETURN; END IF         !--- H2O_<phase> case
  ix = strIdx(newH2OIso, strHead(newName, phases_sep, .TRUE.))                 !--- Isotope index
  IF(ix == 0) THEN; oldName = newName;                  RETURN; END IF         !--- Not a water descendant
  ip = getiPhase(newName)                                                      !--- Phase index
  oldName = TRIM(oldH2OIso(ix))//old_phases(ip:ip)                             !--- <isotope>_<phase>
  tag = strTail(delPhase(newName), TRIM(newH2OIso(ix)))                        !--- Get "_<tag>" if any
  IF(tag /= delPhase(newName) .AND. tag /= '') oldName = TRIM(oldName)//tag    !--- Tagging tracer
END FUNCTION new2oldName_1
!------------------------------------------------------------------------------------------------------------------------------
FUNCTION new2oldName_m(newName, iPhase) RESULT(oldName)
  CHARACTER(LEN=*),  INTENT(IN)  :: newName(:)
  INTEGER, OPTIONAL, INTENT(OUT) :: iPhase
  CHARACTER(LEN=maxlen)          :: oldName(SIZE(newName))
  INTEGER :: i
  oldName = [(new2oldName_1(newName(i), iPhase), i=1, SIZE(newName))]
END FUNCTION new2oldName_m
!------------------------------------------------------------------------------------------------------------------------------


!==============================================================================================================================
!=== GET THE NAME(S) OF THE ANCESTOR(S) OF TRACER(S) "tname" AT GENERATION "igen"  IN THE TRACERS DESCRIPTORS LIST "tr" =======
!==============================================================================================================================
CHARACTER(LEN=maxlen) FUNCTION ancestor_1(t, tname, igen) RESULT(out)
  TYPE(trac_type),   INTENT(IN) :: t(:)
  CHARACTER(LEN=*),  INTENT(IN) :: tname
  INTEGER, OPTIONAL, INTENT(IN) :: igen
  INTEGER :: ig, ix
  ig = 0; IF(PRESENT(igen)) ig = igen
  ix = idxAncestor_1(t, tname, ig)
  out = ''; IF(ix /= 0) out = t(ix)%name
END FUNCTION ancestor_1
!------------------------------------------------------------------------------------------------------------------------------
FUNCTION ancestor_m(t, tname, igen) RESULT(out)
  CHARACTER(LEN=maxlen), ALLOCATABLE     ::   out(:)
  TYPE(trac_type),            INTENT(IN) ::     t(:)
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: tname(:)
  INTEGER,          OPTIONAL, INTENT(IN) :: igen
  INTEGER, ALLOCATABLE :: ix(:)
  INTEGER :: ig
  ig = 0; IF(PRESENT(igen)) ig = igen
  IF(     PRESENT(tname)) ix = idxAncestor_m(t, tname,     ig)
  IF(.NOT.PRESENT(tname)) ix = idxAncestor_m(t, t(:)%name, ig)
  ALLOCATE(out(SIZE(ix))); out(:) = ''
  WHERE(ix /= 0) out = t(ix)%name
END FUNCTION ancestor_m
!==============================================================================================================================


!==============================================================================================================================
!=== GET THE INDEX(ES) OF THE ANCESTOR(S) OF TRACER(S) "tname" AT GENERATION "igen"  IN THE TRACERS DESCRIPTORS LIST "tr" =====
!==============================================================================================================================
INTEGER FUNCTION idxAncestor_1(t, tname, igen) RESULT(out)
! Return the name of the generation "igen" (>=0) ancestor of "tname"
  TYPE(trac_type),   INTENT(IN) :: t(:)
  CHARACTER(LEN=*),  INTENT(IN) :: tname
  INTEGER, OPTIONAL, INTENT(IN) :: igen
  INTEGER :: ig
  ig = 0; IF(PRESENT(igen)) ig = igen
  out = strIdx(t(:)%name, tname)
  IF(out == 0)                 RETURN            !--- Tracer not found
  IF(t(out)%iGeneration <= ig) RETURN            !--- Tracer has a lower generation number than asked generation 'igen"
  DO WHILE(t(out)%iGeneration > ig); out = strIdx(t(:)%name, t(out)%parent); END DO
END FUNCTION idxAncestor_1
!------------------------------------------------------------------------------------------------------------------------------
FUNCTION idxAncestor_m(t, tname, igen) RESULT(out)
  INTEGER,          ALLOCATABLE          ::   out(:)
  TYPE(trac_type),            INTENT(IN) ::     t(:)
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: tname(:)
  INTEGER,          OPTIONAL, INTENT(IN) :: igen
  INTEGER :: ig, ix
  ig = 0; IF(PRESENT(igen)) ig = igen
  IF(     PRESENT(tname)) out = [(idxAncestor_1(t, tname(ix),  ig), ix=1, SIZE(tname))]
  IF(.NOT.PRESENT(tname)) out = [(idxAncestor_1(t, t(ix)%name, ig), ix=1, SIZE(t))]
END FUNCTION idxAncestor_m
!==============================================================================================================================


END MODULE readTracFiles_mod_ico
