Index: LMDZ6/branches/LMDZ-tracers/libf/misc/readTracFiles_mod.f90
===================================================================
--- LMDZ6/branches/LMDZ-tracers/libf/misc/readTracFiles_mod.f90	(revision 3852)
+++ LMDZ6/branches/LMDZ-tracers/libf/misc/readTracFiles_mod.f90	(revision 3891)
@@ -1,6 +1,6 @@
 MODULE readTracFiles_mod
 
-  USE strings_mod,    ONLY: msg, testFile,  strFind, strStack, strCount,   strHead, removeComment, dispTable, fmsg, &
-                            cat, checkList, strIdx,  strParse, strReplace, strTail,  reduceExpr, modname, find, test
+  USE strings_mod, ONLY: msg, testFile,  strFind, strStack, strReduce,  strHead, strCount,   find, dispTable, fmsg, &
+          removeComment, cat, checkList, strIdx,  strParse, strReplace, strTail, reduceExpr, test, modname, get_in
   USE trac_types_mod, ONLY : tra, iso, db, kys
 
@@ -9,4 +9,5 @@
   PRIVATE
 
+  PUBLIC :: initIsotopes
   PUBLIC :: readTracersFiles, aliasTracer, tracersSubset, indexUpdate     !--- TOOLS ASSOCIATED TO TRACERS  DESCRIPTORS
   PUBLIC :: readIsotopesFile                                              !--- TOOLS ASSOCIATED TO ISOTOPES DESCRIPTORS
@@ -82,5 +83,5 @@
 !------------------------------------------------------------------------------------------------------------------------------
   lerr = .FALSE.
-  modname = 'readTracersFiles'
+!  modname = 'readTracersFiles'
   IF(.NOT.ALLOCATED(dBase)) ALLOCATE(dBase(0))
 
@@ -225,5 +226,5 @@
   CHARACTER(LEN=256), ALLOCATABLE :: sec(:)
   INTEGER,            ALLOCATABLE ::  ix(:)
-  INTEGER :: n0, idb, ndb
+  INTEGER :: n0, idb, ndb, i, j
   LOGICAL :: ll
 !------------------------------------------------------------------------------------------------------------------------------
@@ -272,5 +273,5 @@
       ll = strParse(str,' ', keys = s, vals = v, n = n)              !--- Parse <key>=<val> pairs
       tt = dBase(ndb)%trac(:)
-      tmp%name = s(1); tmp%keys = kys(s(1), s(2:n), v(2:n))
+      tmp%name = s(1); tmp%comp=secn; tmp%keys = kys(s(1), s(2:n), v(2:n))
       dBase(ndb)%trac = [tt(:), tmp]
       DEALLOCATE(tt)
@@ -294,14 +295,14 @@
   TYPE(tra), ALLOCATABLE, TARGET, INTENT(INOUT) :: t(:)
   CHARACTER(LEN=*),               INTENT(IN)    :: defName
-  INTEGER :: i0, it, k
-  TYPE(kys), POINTER     :: k0
+  INTEGER :: jd, it, k
+  TYPE(kys), POINTER :: ky
   TYPE(tra), ALLOCATABLE :: tt(:)
-  i0 = strIdx(t(:)%name, defName)
-  IF(i0 == 0) RETURN
-  k0 => t(i0)%keys
-  DO k = 1, SIZE(k0%key)                                             !--- Loop on the keys of the tracer named "defName"
-    CALL addKey_tra(TRIM(k0%key(k)), TRIM(k0%val(k)), t)             !--- Add key to all the tracers (no overwriting)
-  END DO
-  tt = [t(1:i0-1),t(i0+1:SIZE(t))]; CALL MOVE_ALLOC(FROM=tt, TO=t)   !--- Remove the virtual tracer named "defName"
+  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)
+  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
 !==============================================================================================================================
@@ -338,6 +339,4 @@
 
 
-
-
 !==============================================================================================================================
 LOGICAL FUNCTION expandSection(tr, sname, fname) RESULT(lerr)
@@ -451,5 +450,5 @@
   CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: fname                    !--- File name
   CHARACTER(LEN=256) :: mesg
-  CHARACTER(LEN=256) :: bp(SIZE(tr, DIM=1)), pha
+  CHARACTER(LEN=256) :: bp(SIZE(tr, DIM=1)), pha                     !--- Bad phases list, phases of current tracer
   CHARACTER(LEN=1) :: p
   INTEGER :: ip, np, iq, nq
@@ -814,17 +813,20 @@
   
 !==============================================================================================================================
-!=== READ THE ISOTOPES NAMED "iso" FROM THE TRACERS SECTIONS "tr" IN THE FILE "fnam" ; PUT RESULT IN A TRACERS DESCRIPTOR ====
-!===  * SYNTAX IS THE SAME AS IN THE "tracer.def" FILE ; EACH TRACER SECTION CONTAINS ONE LINE EACH OF ITS KNOWN ISOTOPES  ====
-!===  * EACH TRACERS SECTION CAN CONTAIN A "params" VIRTUAL ISOTOPE LINE CONTAINING DEFAULT PARAMETERS FOR THE ISOTOPES    ====
-!===  * IF SOME KEYS ARE FOUND BOTH IN THE "*.def" FILES AND THE "params" SECTION, TEH VALUE FROM "*.def" FILE IS RETAINED ====
-!===  * ON EACH ISOTOPE LINE, A DEFINED KEY CAN BE USED IN THE OTHER KEYS AS A PARAMETER (SIGNLE LEVEL DEPENDENCY !)       ====
-!===  * THE DIFFERENT ISOTOPES SETS (ONE EACH PARENT TRACER) ARE MERGED INTO A SINGLE TRACERS DESCRIPTOR VECTOR            ====
-!===  * THE ROUTINE GIVES AN ERROR IF A REQUIRED ISOTOPE IS NOT AVAILABLE IN THE DATABASE STORED IN "fnam"                 ====
-!==============================================================================================================================
-
+!=== READ FILE "fnam" TO APPEND THE "dBase" TRACERS DATABASE WITH AS MUCH SECTIONS AS PARENTS NAMES IN "isot(:)%prnt":     ====
+!===  * Each section dBase(i)%name contains the isotopes "dBase(i)%trac(:)" descending on "dBase(i)%name"="iso(i)%prnt"    ====
+!===  * 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):                             ====
+!===      prnt,   nzon, zone(:),   niso, keys(:)%name,   nitr, trac(:),   npha, phas,  iTraPha(:,:),  iZonPhi(:,:)         ====
+!===  * 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(iso), TARGET, INTENT(INOUT) :: isot(:)                        !--- Isotopes descriptors (field "prnt" must be defined !)
-  INTEGER :: ik, is, it, idb, nk0, i
+  INTEGER :: ik, is, it, idb, nk0, i, iis
   INTEGER :: nk, ns, nt, ndb, nb0, i0
   CHARACTER(LEN=256), POINTER     :: k(:), v(:), k0(:), v0(:)
@@ -832,6 +834,5 @@
   CHARACTER(LEN=256)     :: val
   TYPE(kys),    POINTER  ::   ky(:)
-  TYPE(tra),    POINTER  ::    t(:)
-  TYPE(tra), ALLOCATABLE ::   tt(:)
+  TYPE(tra),    POINTER  ::   tt(:), t
   TYPE(db),  ALLOCATABLE ::  tdb(:)
   LOGICAL,   ALLOCATABLE :: liso(:)
@@ -844,12 +845,7 @@
   IF(test(readSections(fnam,strStack(isot(:)%prnt,',')),lerr)) RETURN!--- Read sections, one each parent tracer
   ndb = SIZE(dBase, DIM=1)                                           !--- Current database size
-
   DO idb = nb0, ndb
-    t => dBase(idb)%trac(:)
-    nt = SIZE(t)                                                     !--- Number of isotopes in the current database section
-
-PRINT*
-PRINT*,'AVANT:'
-DO it=1,SIZE(t); print*,TRIM(t(it)%name)//':  '//strStack([(TRIM(t(it)%keys%key(i))//'='//TRIM(t(it)%keys%val(i)), i=1, SIZE(t(it)%keys%key))]); END DO
+   iis = idb-nb0+1
+
     !--- GET FEW GLOBAL KEYS FROM "def" FILES AND ADD THEM TO THE 'params' SECTION
     CALL addKeysFromDef(dBase(idb)%trac, 'params')
@@ -858,45 +854,20 @@
     CALL subDefault(dBase(idb)%trac, 'params', .TRUE.)
 
-PRINT*
-PRINT*,'AVANT REDUCTION:'
-    t => dBase(idb)%trac(:)
-    DO it=1,SIZE(t); print*,TRIM(t(it)%name)//':  '//strStack([(TRIM(t(it)%keys%key(i))//'='//TRIM(t(it)%keys%val(i)), i=1, SIZE(t(it)%keys%key))]); END DO
-
-    !--- REDUCE THE EXPRESSIONS TO OBTAIN SCALARS
-    DO it=1, nt
-      v => dBase(idb)%trac(it)%keys%val(:)
-      WHERE(reduceExpr(v, vals)) v = vals
+    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)
+      is = strIdx(isot(iis)%keys(:)%name, dBase(idb)%trac(it)%name)  !--- Index of the "isot(iis)%keys(:)%name" tracer named "t%name"
+      IF(is == 0) CYCLE
+      t => dBase(idb)%trac(it)
+      liso = reduceExpr(t%keys%val, vals)                            !--- Reduce expressions (for substituted variables)
+      isot(iis)%keys(is)%key = PACK(t%keys%key, MASK=liso)
+      isot(iis)%keys(is)%val = PACK(  vals,     MASK=liso)
     END DO
 
-PRINT*
-PRINT*,'APRES:'
-    t => dBase(idb)%trac(:)
-    DO it=1,SIZE(t); print*,TRIM(t(it)%name)//':  '//strStack([(TRIM(t(it)%keys%key(i))//'='//TRIM(t(it)%keys%val(i)), i=1, SIZE(t(it)%keys%key))]); END DO
-
-    !--- TRANSFER THE key=val PAIRS TO THE ISOTOPES DESCRIPTOR
-print*
-print*,'isot%prnt = '//strStack(isot%prnt)
-
-    ky => isot(strIdx(isot(:)%prnt, dBase(idb)%name))%keys           !--- Keys of "isot" tracers with parent "dBase(idb)%name"
-print*,'ky%name = '//strStack(ky%name)
-    is=1
-    DO it = 1, nt; IF(it == i0) CYCLE
-print*,'AAAAAA '//strStack(ky%name)
-print*,'       '//TRIM(t(it)%name)
-      is = strIdx(ky(:)%name, t(it)%name)                            !--- Index of the "isot(:)" tracer named "t(it)%name"
-      IF(is == 0) CYCLE                                              !--- Current isotope is not present in "isot" => not needed
-      k => ky(is)%key; k = t(it)%keys%key
-      v => ky(is)%val; v = t(it)%keys%val
-      WHERE(reduceExpr(v, vals)) v = vals
-      DO ik=1, SIZE(k); IF(reduceExpr(v(ik),val)) v(ik) = val; END DO!--- Reduce operations (for substituted variables)
-print*,'(4) '//strStack([(TRIM(k(i))//'='//TRIM(v(i)), i=1, SIZE(k))])
-    END DO
-print*,'(7) i0=',i0
-
     !--- CHECK FOR MISSING ISOTOPES (NO KEYS ALLOCATED)
-    liso = [( ALLOCATED(ky(is)%key), is=1, SIZE(ky) )]
-print*,'liso=',liso
-    IF(test(checkList(ky(:)%name, &
-      .NOT.liso, 'Check file "'//TRIM(fnam)//'" in section "'//TRIM(dBase(idb)%name)//'"', 'isotopes', 'missing'),lerr)) RETURN
+    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
 
@@ -907,8 +878,100 @@
     ALLOCATE(tdb(nb0-1)); tdb(1:nb0-1)=dBase(1:nb0-1); CALL MOVE_ALLOC(FROM=tdb, TO=dBase)
   END IF
-  lerr = dispIsotopes(isot,'isotopes parameters read from file')
+  lerr = dispIsotopes(isot, 'Isotopes parameters read from file')
 
 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(tra), ALLOCATABLE, TARGET, INTENT(INOUT) :: trac(:)
+  TYPE(iso), ALLOCATABLE, TARGET, INTENT(INOUT) :: isot(:)
+  CHARACTER(LEN=256), ALLOCATABLE :: p(:), str(:)                    !--- Temporary storage
+  CHARACTER(LEN=256) :: iname
+  CHARACTER(LEN=1)   :: ph                                           !--- Phase
+  INTEGER :: nbIso, ic, ip, iq, it, iz
+  LOGICAL, ALLOCATABLE :: ll(:)                                      !--- Mask
+  TYPE(tra), POINTER   ::  t(:), t1
+  TYPE(iso), POINTER   ::  s
+
+  t => trac
+
+  p = PACK(delPhase(t%prnt), MASK = t%type=='tracer' .AND. t%igen==2)!--- Parents of 2nd generation 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(:)%prnt = p
+  DO ic = 1, SIZE(p)                                                 !--- Loop on isotopes classes
+    s => isot(ic)
+    iname = s%prnt                                                   !--- 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(:)%prnt) == iname .AND. t(:)%phas == 'g'
+    str = PACK(delPhase(t(:)%name), MASK = ll)                       !--- Effectively found isotopes of "iname"
+    s%niso = SIZE(str)                                               !--- Number of "effectively found isotopes of "iname"
+    ALLOCATE(s%keys(s%niso))
+    FORALL(it = 1:s%niso) s%keys(it)%name = str(it)
+
+    !=== Geographic tagging tracers descending on tracer "iname": mask, names, number
+    ll = t(:)%type=='tag'    .AND. delPhase(t(:)%nam1) == iname .AND. t(:)%igen == 3
+    s%zone = PACK(strTail(t(:)%name,'_'), MASK = ll)                 !--- Tagging zones names  for isotopes category "iname"
+    CALL strReduce(s%zone)
+    s%nzon = SIZE(s%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)
+    s%nitr = s%niso + SIZE(str)                                      !--- Number of isotopes + their geographic tracers [ntraciso]
+    ALLOCATE(s%trac(s%nitr))
+    FORALL(it = 1:s%niso) s%trac(it) = s%keys(it)%name
+    FORALL(it = s%niso+1:s%nitr) s%trac(it) = str(it-s%niso)
+
+    !=== Phases for tracer "iname"
+    s%phas = ''
+    DO ip = 1, nphases; ph = known_phases(ip:ip); IF(strIdx(t%name,addPhase(iname, ph)) /= 0) s%phas = TRIM(s%phas)//ph; END DO
+    s%npha = LEN_TRIM(s%phas)                                        !--- 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%nam1) /= iname) CYCLE                            !--- Only deal with tracers descending on "iname"
+      t1%iso_igr = ic                                                 !--- Isotopes family       idx in list "isotopes(:)%prnt"
+      t1%iso_num = strIdx(s%trac, delPhase(strHead(t1%name,'_')))     !--- Current isotope       idx in effective isotopes list
+      t1%iso_zon = strIdx(s%zone,          strTail(t1%name,'_') )     !--- Current isotope zone  idx in effective zones    list
+      t1%iso_pha =  INDEX(s%phas,TRIM(t1%phas))                       !--- Current isotope phase idx in effective phases   list
+      IF(t1%igen /= 3) t1%iso_zon = 0                                 !--- Skip possible generation 2 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
+    s%iTraPha = RESHAPE( [( (strIdx(t(:)%name,  addPhase(s%trac(it),s%phas(ip:ip))),     it=1, s%nitr), ip=1, s%npha)], &
+                         [s%nitr, s%npha] )
+
+    !=== Table used to get ix (index in tagging tracers isotopes list, size nitr) from the zone and isotope indexes
+    s%iZonIso = RESHAPE( [( (strIdx(s%trac(:), TRIM(s%trac(it))//'_'//TRIM(s%zone(iz))), iz=1, s%nzon), it=1, s%niso)], &
+                         [s%nzon, s%niso] )
+  END DO
+  
+  !=== Indexes, in dynamical tracers list, of the tracers transmitted to phytrac (nqtottr non-vanishing elements)
+  ll = delPhase(t%name)/='H2O' .AND. t%iso_num ==0              !--- Mask of tracers passed to the physics
+  t(:)%itr = UNPACK([(iq,iq=1,COUNT(ll))], ll, [(0, iq=1, SIZE(t))])
+
+  !=== 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)) CALL abort_gcm(modname,'Problem when reading isotopes parameters',1)
+
+END SUBROUTINE initIsotopes
+!==============================================================================================================================
+
 
 !==============================================================================================================================
@@ -955,49 +1018,37 @@
   IF(iky == 0) THEN
     nky = SIZE(ky%key)
-    IF(nky == 0) THEN
-      ky%key = TRIM(key); ky%val = TRIM(val)
-    ELSE
-      ALLOCATE(k(nky+1)); k(1:nky) = ky%key; k(nky+1) = TRIM(key); CALL MOVE_ALLOC(FROM=k, TO=ky%key)
-      ALLOCATE(v(nky+1)); v(1:nky) = ky%val; v(nky+1) = TRIM(val); CALL MOVE_ALLOC(FROM=v, TO=ky%val)
-    END IF
+    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) = TRIM(key); ky%val(iky) = TRIM(val)
+    ky%key(iky) = key; ky%val(iky) = val
   END IF
 END SUBROUTINE addKey_1
 !==============================================================================================================================
-SUBROUTINE addKey_tra(key, val, tr, lOverWrite, tname)
-!------------------------------------------------------------------------------------------------------------------------------
-! Purpose: Add the <key>=<val> pair in all the components of the "tr(itr)%keys" keys descriptor:
-!          * "tname"   specified: for the index "itr" of the tracer named "tname"
-!          * "tname" unspecified: for all the tracers
-!------------------------------------------------------------------------------------------------------------------------------
-  CHARACTER(LEN=*),           INTENT(IN)    :: key, val
-  TYPE(tra),                  INTENT(INOUT) :: tr(:)
-  LOGICAL,          OPTIONAL, INTENT(IN)    :: lOverWrite
-  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)    :: tname
+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(kys),         INTENT(INOUT) :: ky(:)
+  LOGICAL, OPTIONAL, INTENT(IN)    :: lOverWrite
   INTEGER :: itr
   LOGICAL :: lo
 !------------------------------------------------------------------------------------------------------------------------------
   lo=.FALSE.; IF(PRESENT(lOverWrite)) lo=lOverWrite
-  IF(PRESENT(tname)) THEN
-    itr = strIdx(tr%name, tname)
-    IF(itr == 0) RETURN
-    CALL addKey_1(key, val, tr(itr)%keys, lo)
-  ELSE
-    DO itr = 1, SIZE(tr); CALL addKey_1(key, val, tr(itr)%keys, lo); END DO
-  END IF
-END SUBROUTINE addKey_tra
-!==============================================================================================================================
-SUBROUTINE addKeysFromDef(tr, tr0)
-  USE ioipsl_getin_p_mod, ONLY : getin_p
-  TYPE(tra), ALLOCATABLE, INTENT(INOUT) :: tr(:)
+  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(tra), ALLOCATABLE, INTENT(INOUT) :: t(:)
   CHARACTER(LEN=*),       INTENT(IN)    :: tr0
   CHARACTER(LEN=256) :: val
-  INTEGER            :: ik, i0
-  i0 = strIdx(tr%name, tr0)
-  IF(i0 == 0) RETURN
-  DO ik = 1, SIZE(tr(i0)%keys%key)
-    val   =   'zzzz'; CALL  getin_p(tr(i0)%keys%key(ik), val)
-    IF(val /= 'zzzz') CALL addKey_1(tr(i0)%keys%key(ik), val, tr(i0)%keys, .TRUE.)
+  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, 'zzzz')
+    IF(val /= 'zzzz') CALL addKey_1(t(jd)%keys%key(ik), val, t(jd)%keys, .TRUE.)
   END DO
 END SUBROUTINE addKeysFromDef
@@ -1064,6 +1115,6 @@
 LOGICAL FUNCTION getKeyByName_s1(keyn, val, tname, ky) RESULT(lerr)
   !--- Purpose: Get the value of the key named "keyn" for the tracer named "tnam".
-  !     * "ky" unspecified: try in "tracers" for "tnam" with phase suffix, then in "isotopes" without.
-  !     * "ky"   specified: try in "ky"      for "tnam" with, then without phase suffix.
+  !     * "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
@@ -1074,6 +1125,6 @@
   lerr = .FALSE.
   IF(PRESENT(ky)) THEN
-    val = getKeyByName_prv(keyn,          tname , ky);    IF(val /= '') RETURN !--- "ky" and "tnam"
-    val = getKeyByName_prv(keyn, delPhase(tname), ky)                          !--- "ky" and "tnam" without phase
+    val = getKeyByName_prv(keyn, tname , ky);    IF(val /= '') RETURN          !--- "ky" and "tnam"
+    val = getKeyByName_prv(keyn, delPhase(strHead(tname,'_')), ky)             !--- "ky" and "tnam" without phase
   ELSE
     IF(.NOT.ALLOCATED(tracers))  RETURN
@@ -1081,5 +1132,5 @@
     IF(.NOT.ALLOCATED(isotopes)) RETURN
     IF(SIZE(isotopes) == 0)      RETURN
-    DO is = 1, SIZE(isotopes); IF(strIdx(isotopes(is)%keys(:)%name, tname) /= 0) EXIT; END DO
+    DO is = 1, SIZE(isotopes); IF(strIdx(isotopes(is)%keys(:)%name, delPhase(strHead(tname,'_'))) /= 0) EXIT; END DO
     IF(is /= 0) val = getKeyByName_prv(keyn, tname, isotopes(is)%keys(:))      !--- "isotopes" and "tnam" without phase
   END IF
@@ -1175,9 +1226,13 @@
 ELEMENTAL CHARACTER(LEN=256) FUNCTION delPhase(s) RESULT(out)
   CHARACTER(LEN=*), INTENT(IN) :: s
-  INTEGER :: l
+  INTEGER :: l, i
   out = s
   IF(s == '') RETURN
-  l=LEN_TRIM(s)
-  IF(s(l-1:l-1)=='-' .AND. INDEX(known_phases,s(l:l))/=0) out = s(1:l-2)
+  i = INDEX(s, '_'); l = LEN_TRIM(s)
+  IF(i == 0) THEN
+    IF(s(l-1:l-1)=='-' .AND. INDEX(known_phases,s(l:l)) /= 0) out = s(1:l-2)
+  ELSE; i=i-1
+    IF(s(i-1:i-1)=='-' .AND. INDEX(known_phases,s(i:i)) /= 0) out = s(1:i-2)//s(i+1:l)
+  END IF
 END FUNCTION delPhase
 !------------------------------------------------------------------------------------------------------------------------------
@@ -1185,6 +1240,10 @@
   CHARACTER(LEN=*), INTENT(IN) :: s
   CHARACTER(LEN=1), INTENT(IN) :: pha
-  IF(INDEX(s,'_')==0) THEN; out = TRIM(s)//'-'//pha; RETURN; END IF
-  out = TRIM(strHead(s,'_'))//'-'//pha//TRIM(strTail(s,'_'))
+  INTEGER :: l, i
+  out = s
+  IF(s == '') RETURN
+  i = INDEX(s, '_'); l = LEN_TRIM(s)
+  IF(i == 0) out =  TRIM(s)//'-'//pha
+  IF(i /= 0) out = s(1:i-1)//'-'//pha//'_'//s(i+1:l)
 END FUNCTION addPhase_1
 !------------------------------------------------------------------------------------------------------------------------------
