Index: LMDZ6/trunk/libf/dyn3d/dynetat0.F90
===================================================================
--- LMDZ6/trunk/libf/dyn3d/dynetat0.F90	(revision 4120)
+++ LMDZ6/trunk/libf/dyn3d/dynetat0.F90	(revision 4120)
@@ -0,0 +1,225 @@
+SUBROUTINE dynetat0(fichnom,vcov,ucov,teta,q,masse,ps,phis,time)
+!
+!-------------------------------------------------------------------------------
+! Authors: P. Le Van , L.Fairhead
+!-------------------------------------------------------------------------------
+! Purpose: Initial state reading.
+!-------------------------------------------------------------------------------
+  USE infotrac,    ONLY: nqtot, tracers, niso, iqiso, iso_indnum, iso_num, tnat, alpha_ideal, ok_isotopes, iH2O
+  USE strings_mod, ONLY: maxlen, msg, strStack, real2str
+  USE netcdf,      ONLY: NF90_OPEN,  NF90_NOWRITE, NF90_INQ_VARID, &
+                         NF90_CLOSE, NF90_GET_VAR, NF90_NoErr
+  USE readTracFiles_mod, ONLY: new2oldName
+  USE control_mod, ONLY: planet_type
+  USE assert_eq_m, ONLY: assert_eq
+  USE comvert_mod, ONLY: pa,preff
+  USE comconst_mod, ONLY: cpp, daysec, dtvr, g, im, jm, kappa, lllm, omeg, rad
+  USE logic_mod, ONLY: fxyhypb, ysinus
+  USE serre_mod, ONLY: clon, clat, grossismx, grossismy
+  USE temps_mod, ONLY: annee_ref, day_ini, day_ref, itau_dyn, start_time
+  USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
+
+  IMPLICIT NONE
+  include "dimensions.h"
+  include "paramet.h"
+  include "comgeom2.h"
+  include "description.h"
+  include "iniprint.h"
+!===============================================================================
+! Arguments:
+  CHARACTER(LEN=*), INTENT(IN) :: fichnom          !--- FILE NAME
+  REAL, INTENT(OUT) ::  vcov(iip1,jjm, llm)        !--- V COVARIANT WIND
+  REAL, INTENT(OUT) ::  ucov(iip1,jjp1,llm)        !--- U COVARIANT WIND
+  REAL, INTENT(OUT) ::  teta(iip1,jjp1,llm)        !--- POTENTIAL TEMP.
+  REAL, INTENT(OUT) ::     q(iip1,jjp1,llm,nqtot)  !--- TRACERS
+  REAL, INTENT(OUT) :: masse(iip1,jjp1,llm)        !--- MASS PER CELL
+  REAL, INTENT(OUT) ::    ps(iip1,jjp1)            !--- GROUND PRESSURE
+  REAL, INTENT(OUT) ::  phis(iip1,jjp1)            !--- GEOPOTENTIAL
+!===============================================================================
+! Local variables:
+  CHARACTER(LEN=maxlen) :: mesg, var, modname, oldVar
+  INTEGER, PARAMETER :: length=100
+  INTEGER :: iq, fID, vID, idecal, iqParent, iName, iZone, iPhase
+  REAL    :: time, tab_cntrl(length)               !--- RUN PARAMS TABLE
+!-------------------------------------------------------------------------------
+  modname="dynetat0"
+
+!--- Initial state file opening
+  var=fichnom
+  CALL err(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
+  CALL get_var1("controle",tab_cntrl)
+
+!!! AS: idecal is a hack to be able to read planeto starts...
+!!!     .... while keeping everything OK for LMDZ EARTH
+  IF(planet_type=="generic") THEN
+    CALL msg('NOTE NOTE NOTE : Planeto-like start files', modname)
+    idecal = 4
+    annee_ref  = 2000
+  ELSE
+    CALL msg('NOTE NOTE NOTE : Earth-like start files', modname)
+    idecal = 5
+    annee_ref  = tab_cntrl(5)
+  END IF
+  im         = tab_cntrl(1)
+  jm         = tab_cntrl(2)
+  lllm       = tab_cntrl(3)
+  day_ref    = tab_cntrl(4)
+  rad        = tab_cntrl(idecal+1)
+  omeg       = tab_cntrl(idecal+2)
+  g          = tab_cntrl(idecal+3)
+  cpp        = tab_cntrl(idecal+4)
+  kappa      = tab_cntrl(idecal+5)
+  daysec     = tab_cntrl(idecal+6)
+  dtvr       = tab_cntrl(idecal+7)
+  etot0      = tab_cntrl(idecal+8)
+  ptot0      = tab_cntrl(idecal+9)
+  ztot0      = tab_cntrl(idecal+10)
+  stot0      = tab_cntrl(idecal+11)
+  ang0       = tab_cntrl(idecal+12)
+  pa         = tab_cntrl(idecal+13)
+  preff      = tab_cntrl(idecal+14)
+!
+  clon       = tab_cntrl(idecal+15)
+  clat       = tab_cntrl(idecal+16)
+  grossismx  = tab_cntrl(idecal+17)
+  grossismy  = tab_cntrl(idecal+18)
+!
+  IF ( tab_cntrl(idecal+19)==1. )  THEN
+    fxyhypb  = .TRUE.
+!   dzoomx   = tab_cntrl(25)
+!   dzoomy   = tab_cntrl(26)
+!   taux     = tab_cntrl(28)
+!   tauy     = tab_cntrl(29)
+  ELSE
+    fxyhypb = .FALSE.
+    ysinus  = tab_cntrl(idecal+22)==1.
+  END IF
+
+  day_ini    = tab_cntrl(30)
+  itau_dyn   = tab_cntrl(31)
+  start_time = tab_cntrl(32)
+
+!-------------------------------------------------------------------------------
+  CALL msg('rad, omeg, g, cpp, kappa = '//TRIM(strStack(real2str([rad,omeg,g,cpp,kappa]))), modname)
+  CALL check_dim(im,iim,'im','im')
+  CALL check_dim(jm,jjm,'jm','jm')
+  CALL check_dim(lllm,llm,'lm','lllm')
+  CALL get_var1("rlonu",rlonu)
+  CALL get_var1("rlatu",rlatu)
+  CALL get_var1("rlonv",rlonv)
+  CALL get_var1("rlatv",rlatv)
+  CALL get_var2("cu"   ,cu)
+  CALL get_var2("cv"   ,cv)
+  CALL get_var2("aire" ,aire)
+  var="temps"
+  IF(NF90_INQ_VARID(fID,var,vID)/=NF90_NoErr) THEN
+    CALL msg('missing field <temps> ; trying with <Time>', modname)
+    var="Time"
+    CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
+  END IF
+  CALL err(NF90_GET_VAR(fID,vID,time),"get",var)
+  CALL get_var2("phisinit",phis)
+  CALL get_var3("ucov",ucov)
+  CALL get_var3("vcov",vcov)
+  CALL get_var3("teta",teta)
+  CALL get_var3("masse",masse)
+  CALL get_var2("ps",ps)
+
+!--- Tracers
+  DO iq=1,nqtot
+    var = tracers(iq)%name
+    oldVar = new2oldName(var)
+    !--------------------------------------------------------------------------------------------------------------------------
+    IF(NF90_INQ_VARID(fID, var, vID) == NF90_NoErr) THEN                                 !=== REGULAR CASE
+      CALL err(NF90_GET_VAR(fID,vID,q(:,:,:,iq)),"get",var)
+    !--------------------------------------------------------------------------------------------------------------------------
+    ELSE IF(NF90_INQ_VARID(fID, oldVar, vID) == NF90_NoErr) THEN                         !=== OLD NAME
+      CALL msg('Tracer <'//TRIM(var)//'> is missing => initialized to <'//TRIM(oldVar)//'>', modname)
+      CALL err(NF90_GET_VAR(fID,vID,q(:,:,:,iq)),"get",oldVar)
+    !--------------------------------------------------------------------------------------------------------------------------
+#ifdef INCA
+    ELSE IF(NF90_INQ_VARID(fID, 'OX',   vID) == NF90_NoErr .AND. var == 'O3') THEN       !=== INCA: OX INSTEAD OF O3
+      CALL msg('Tracer <O3> is missing => initialized to <OX>', modname)
+      CALL err(NF90_GET_VAR(fID,vID,q(:,:,:,iq)),"get",'OX')
+    !--------------------------------------------------------------------------------------------------------------------------
+#endif
+    ELSE IF(tracers(iq)%iso_iGroup == iH2O .AND. niso > 0) THEN                          !=== WATER ISOTOPES
+!     iName    = tracers(iq)%iso_iName  ! (next commit)
+      iName    = iso_num(iq)
+      iPhase   = tracers(iq)%iso_iPhase
+      iqParent = tracers(iq)%iqParent
+      IF(tracers(iq)%iso_iZone == 0) THEN
+         CALL msg('Tracer <'//TRIM(var)//'> is missing => initialized with a simplified Rayleigh distillation law.', modname)
+         q(:,:,:,iq) = q(:,:,:,iqParent)*tnat(iName)*(q(:,:,:,iqParent)/30.e-3)**(alpha_ideal(iName)-1.)
+      ELSE
+         CALL msg('Tracer <'//TRIM(var)//'> is missing => initialized to its parent isotope concentration.', modname)
+         q(:,:,:,iq) = q(:,:,:,iqiso(iso_indnum(iq),iPhase))
+      END IF
+    !--------------------------------------------------------------------------------------------------------------------------
+    ELSE                                                                                 !=== MISSING: SET TO 0
+      CALL msg('Tracer <'//TRIM(var)//'> is missing => initialized to zero', modname)
+      q(:,:,:,iq)=0.
+    !--------------------------------------------------------------------------------------------------------------------------
+    END IF
+  END DO
+
+  CALL err(NF90_CLOSE(fID),"close",fichnom)
+  day_ini=day_ini+INT(time)
+  time=time-INT(time)
+
+
+  CONTAINS
+
+
+SUBROUTINE check_dim(n1,n2,str1,str2)
+  INTEGER,          INTENT(IN) :: n1, n2
+  CHARACTER(LEN=*), INTENT(IN) :: str1, str2
+  CHARACTER(LEN=maxlen) :: s1, s2
+  IF(n1/=n2) THEN
+    s1='value of '//TRIM(str1)//' ='
+    s2=' read in starting file differs from parametrized '//TRIM(str2)//' ='
+    WRITE(mesg,'(10x,a,i4,2x,a,i4)')TRIM(ADJUSTL(s1)),n1,TRIM(ADJUSTL(s2)),n2
+    CALL ABORT_gcm(TRIM(modname),TRIM(mesg),1)
+  END IF
+END SUBROUTINE check_dim
+
+
+SUBROUTINE get_var1(var,v)
+  CHARACTER(LEN=*), INTENT(IN)  :: var
+  REAL,             INTENT(OUT) :: v(:)
+  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
+  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
+END SUBROUTINE get_var1
+
+
+SUBROUTINE get_var2(var,v)
+  CHARACTER(LEN=*), INTENT(IN)  :: var
+  REAL,             INTENT(OUT) :: v(:,:)
+  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
+  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
+END SUBROUTINE get_var2
+
+
+SUBROUTINE get_var3(var,v)
+  CHARACTER(LEN=*), INTENT(IN)  :: var
+  REAL,             INTENT(OUT) :: v(:,:,:)
+  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
+  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
+END SUBROUTINE get_var3
+
+
+SUBROUTINE err(ierr,typ,nam)
+  INTEGER,          INTENT(IN) :: ierr   !--- NetCDF ERROR CODE
+  CHARACTER(LEN=*), INTENT(IN) :: typ    !--- TYPE OF OPERATION
+  CHARACTER(LEN=*), INTENT(IN) :: nam    !--- FIELD/FILE NAME
+  IF(ierr==NF90_NoERR) RETURN
+  SELECT CASE(typ)
+    CASE('inq');   mesg="Field <"//TRIM(nam)//"> is missing"
+    CASE('get');   mesg="Reading failed for <"//TRIM(nam)//">"
+    CASE('open');  mesg="File opening failed for <"//TRIM(nam)//">"
+    CASE('close'); mesg="File closing failed for <"//TRIM(nam)//">"
+  END SELECT
+  CALL ABORT_gcm(TRIM(modname),TRIM(mesg),ierr)
+END SUBROUTINE err
+
+END SUBROUTINE dynetat0
Index: LMDZ6/trunk/libf/dyn3d/dynetat0.f90
===================================================================
--- LMDZ6/trunk/libf/dyn3d/dynetat0.f90	(revision 4119)
+++ 	(revision )
@@ -1,208 +1,0 @@
-SUBROUTINE dynetat0(fichnom,vcov,ucov,teta,q,masse,ps,phis,time)
-!
-!-------------------------------------------------------------------------------
-! Authors: P. Le Van , L.Fairhead
-!-------------------------------------------------------------------------------
-! Purpose: Initial state reading.
-!-------------------------------------------------------------------------------
-  USE infotrac,    ONLY: nqtot, tracers, iqiso, iso_indnum, tnat, alpha_ideal, &
-                         ok_isotopes
-  USE strings_mod, ONLY: maxlen
-  USE netcdf,      ONLY: NF90_OPEN,  NF90_NOWRITE, NF90_INQ_VARID, NF90_NoErr, &
-                         NF90_CLOSE, NF90_GET_VAR
-  USE control_mod, ONLY: planet_type
-  USE assert_eq_m, ONLY: assert_eq
-  USE comvert_mod, ONLY: pa,preff
-  USE comconst_mod, ONLY: cpp, daysec, dtvr, g, im, jm, kappa, lllm, omeg, rad
-  USE logic_mod, ONLY: fxyhypb, ysinus
-  USE serre_mod, ONLY: clon, clat, grossismx, grossismy
-  USE temps_mod, ONLY: annee_ref, day_ini, day_ref, itau_dyn, start_time
-  USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
-
-  IMPLICIT NONE
-  include "dimensions.h"
-  include "paramet.h"
-  include "comgeom2.h"
-  include "description.h"
-  include "iniprint.h"
-!===============================================================================
-! Arguments:
-  CHARACTER(LEN=*), INTENT(IN) :: fichnom          !--- FILE NAME
-  REAL, INTENT(OUT) ::  vcov(iip1,jjm, llm)        !--- V COVARIANT WIND
-  REAL, INTENT(OUT) ::  ucov(iip1,jjp1,llm)        !--- U COVARIANT WIND
-  REAL, INTENT(OUT) ::  teta(iip1,jjp1,llm)        !--- POTENTIAL TEMP.
-  REAL, INTENT(OUT) ::     q(iip1,jjp1,llm,nqtot)  !--- TRACERS
-  REAL, INTENT(OUT) :: masse(iip1,jjp1,llm)        !--- MASS PER CELL
-  REAL, INTENT(OUT) ::    ps(iip1,jjp1)            !--- GROUND PRESSURE
-  REAL, INTENT(OUT) ::  phis(iip1,jjp1)            !--- GEOPOTENTIAL
-!===============================================================================
-! Local variables:
-  CHARACTER(LEN=maxlen) :: msg, var, modname
-  INTEGER, PARAMETER :: length=100
-  INTEGER :: iq, fID, vID, idecal, iqParent, iName, iZone, iPhase
-  REAL    :: time, tab_cntrl(length)               !--- RUN PARAMS TABLE
-!-------------------------------------------------------------------------------
-  modname="dynetat0"
-
-!--- Initial state file opening
-  var=fichnom
-  CALL err(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
-  CALL get_var1("controle",tab_cntrl)
-
-!!! AS: idecal is a hack to be able to read planeto starts...
-!!!     .... while keeping everything OK for LMDZ EARTH
-  IF(planet_type=="generic") THEN
-    WRITE(lunout,*)'NOTE NOTE NOTE : Planeto-like start files'
-    idecal = 4
-    annee_ref  = 2000
-  ELSE
-    WRITE(lunout,*)'NOTE NOTE NOTE : Earth-like start files'
-    idecal = 5
-    annee_ref  = tab_cntrl(5)
-  END IF
-  im         = tab_cntrl(1)
-  jm         = tab_cntrl(2)
-  lllm       = tab_cntrl(3)
-  day_ref    = tab_cntrl(4)
-  rad        = tab_cntrl(idecal+1)
-  omeg       = tab_cntrl(idecal+2)
-  g          = tab_cntrl(idecal+3)
-  cpp        = tab_cntrl(idecal+4)
-  kappa      = tab_cntrl(idecal+5)
-  daysec     = tab_cntrl(idecal+6)
-  dtvr       = tab_cntrl(idecal+7)
-  etot0      = tab_cntrl(idecal+8)
-  ptot0      = tab_cntrl(idecal+9)
-  ztot0      = tab_cntrl(idecal+10)
-  stot0      = tab_cntrl(idecal+11)
-  ang0       = tab_cntrl(idecal+12)
-  pa         = tab_cntrl(idecal+13)
-  preff      = tab_cntrl(idecal+14)
-!
-  clon       = tab_cntrl(idecal+15)
-  clat       = tab_cntrl(idecal+16)
-  grossismx  = tab_cntrl(idecal+17)
-  grossismy  = tab_cntrl(idecal+18)
-!
-  IF ( tab_cntrl(idecal+19)==1. )  THEN
-    fxyhypb  = .TRUE.
-!   dzoomx   = tab_cntrl(25)
-!   dzoomy   = tab_cntrl(26)
-!   taux     = tab_cntrl(28)
-!   tauy     = tab_cntrl(29)
-  ELSE
-    fxyhypb = .FALSE.
-    ysinus  = tab_cntrl(idecal+22)==1.
-  END IF
-
-  day_ini    = tab_cntrl(30)
-  itau_dyn   = tab_cntrl(31)
-  start_time = tab_cntrl(32)
-
-!-------------------------------------------------------------------------------
-  WRITE(lunout,*)TRIM(modname)//': rad,omeg,g,cpp,kappa',rad,omeg,g,cpp,kappa
-  CALL check_dim(im,iim,'im','im')
-  CALL check_dim(jm,jjm,'jm','jm')
-  CALL check_dim(lllm,llm,'lm','lllm')
-  CALL get_var1("rlonu",rlonu)
-  CALL get_var1("rlatu",rlatu)
-  CALL get_var1("rlonv",rlonv)
-  CALL get_var1("rlatv",rlatv)
-  CALL get_var2("cu"   ,cu)
-  CALL get_var2("cv"   ,cv)
-  CALL get_var2("aire" ,aire)
-  var="temps"
-  IF(NF90_INQ_VARID(fID,var,vID)/=NF90_NoErr) THEN
-    WRITE(lunout,*)TRIM(modname)//": missing field <temps>"
-    WRITE(lunout,*)TRIM(modname)//": trying with <Time>"; var="Time"
-    CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
-  END IF
-  CALL err(NF90_GET_VAR(fID,vID,time),"get",var)
-  CALL get_var2("phisinit",phis)
-  CALL get_var3("ucov",ucov)
-  CALL get_var3("vcov",vcov)
-  CALL get_var3("teta",teta)
-  CALL get_var3("masse",masse)
-  CALL get_var2("ps",ps)
-
-!--- Tracers
-  DO iq=1,nqtot
-    var=TRIM(tracers(iq)%name)
-    IF(NF90_INQ_VARID(fID,var,vID)==NF90_NoErr) THEN
-      CALL err(NF90_GET_VAR(fID,vID,q(:,:,:,iq)),"get",var); CYCLE
-    END IF
-    WRITE(lunout,*)TRIM(modname)//": Tracer <"//TRIM(var)//"> is missing"
-    WRITE(lunout,*)"         It is hence initialized to zero"
-    q(:,:,:,iq)=0.
-   !--- CRisi: for isotops, theoretical initialization using very simplified
-   !           Rayleigh distillation law.
-    iName = tracers(iq)%iso_iName
-    IF(.NOT.ok_isotopes .OR. iName<=0) CYCLE
-    iZone = tracers(iq)%iso_iZone
-    iPhase= tracers(iq)%iso_iPhase
-    iqParent = tracers(iq)%iqParent
-    IF(iZone==0) q(:,:,:,iq) = q(:,:,:,iqParent)*tnat(iName)    &
-                             *(q(:,:,:,iqParent)/30.e-3)**(alpha_ideal(iName)-1.)
-    IF(iZone==1) q(:,:,:,iq) = q(:,:,:,iqiso(iso_indnum(iq),iPhase))
-  END DO
-
-  CALL err(NF90_CLOSE(fID),"close",fichnom)
-  day_ini=day_ini+INT(time)
-  time=time-INT(time)
-
-
-  CONTAINS
-
-
-SUBROUTINE check_dim(n1,n2,str1,str2)
-  INTEGER,          INTENT(IN) :: n1, n2
-  CHARACTER(LEN=*), INTENT(IN) :: str1, str2
-  CHARACTER(LEN=maxlen) :: s1, s2
-  IF(n1/=n2) THEN
-    s1='value of '//TRIM(str1)//' ='
-    s2=' read in starting file differs from parametrized '//TRIM(str2)//' ='
-    WRITE(msg,'(10x,a,i4,2x,a,i4)')TRIM(ADJUSTL(s1)),n1,TRIM(ADJUSTL(s2)),n2
-    CALL ABORT_gcm(TRIM(modname),TRIM(msg),1)
-  END IF
-END SUBROUTINE check_dim
-
-
-SUBROUTINE get_var1(var,v)
-  CHARACTER(LEN=*), INTENT(IN)  :: var
-  REAL,             INTENT(OUT) :: v(:)
-  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
-  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
-END SUBROUTINE get_var1
-
-
-SUBROUTINE get_var2(var,v)
-  CHARACTER(LEN=*), INTENT(IN)  :: var
-  REAL,             INTENT(OUT) :: v(:,:)
-  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
-  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
-END SUBROUTINE get_var2
-
-
-SUBROUTINE get_var3(var,v)
-  CHARACTER(LEN=*), INTENT(IN)  :: var
-  REAL,             INTENT(OUT) :: v(:,:,:)
-  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
-  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
-END SUBROUTINE get_var3
-
-
-SUBROUTINE err(ierr,typ,nam)
-  INTEGER,          INTENT(IN) :: ierr   !--- NetCDF ERROR CODE
-  CHARACTER(LEN=*), INTENT(IN) :: typ    !--- TYPE OF OPERATION
-  CHARACTER(LEN=*), INTENT(IN) :: nam    !--- FIELD/FILE NAME
-  IF(ierr==NF90_NoERR) RETURN
-  SELECT CASE(typ)
-    CASE('inq');   msg="Field <"//TRIM(nam)//"> is missing"
-    CASE('get');   msg="Reading failed for <"//TRIM(nam)//">"
-    CASE('open');  msg="File opening failed for <"//TRIM(nam)//">"
-    CASE('close'); msg="File closing failed for <"//TRIM(nam)//">"
-  END SELECT
-  CALL ABORT_gcm(TRIM(modname),TRIM(msg),1)
-END SUBROUTINE err
-
-END SUBROUTINE dynetat0
Index: LMDZ6/trunk/libf/dyn3d/iniacademic.F90
===================================================================
--- LMDZ6/trunk/libf/dyn3d/iniacademic.F90	(revision 4119)
+++ LMDZ6/trunk/libf/dyn3d/iniacademic.F90	(revision 4120)
@@ -6,5 +6,5 @@
   USE filtreg_mod, ONLY: inifilr
   USE infotrac,    ONLY: nqtot, niso_possibles, ok_isotopes, ok_iso_verif, tnat, alpha_ideal, &
-                         iqiso, tracers, iso_indnum
+                         iqiso, tracers, iso_indnum, iso_num
   USE control_mod, ONLY: day_step,planet_type
   use exner_hyb_m, only: exner_hyb
@@ -22,4 +22,5 @@
   USE temps_mod, ONLY: annee_ref, day_ini, day_ref
   USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
+  USE readTracFiles_mod, ONLY: addPhase
 
   !   Author:    Frederic Hourdin      original: 15/01/93
@@ -62,5 +63,5 @@
   real tetastrat ! potential temperature in the stratosphere, in K
   real tetajl(jjp1,llm)
-  INTEGER i,j,l,lsup,ij, iq, iName, iZone, iPhase, iqParent
+  INTEGER i,j,l,lsup,ij, iq, iName, iPhase, iqParent
 
   REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
@@ -276,19 +277,19 @@
            do iq=1,nqtot
               q(:,:,iq)=0.
-!              IF(tracers(iq)%name == 'H2O'//phases_sep//'g') q(:,:,iq)=1.e-10
-!              IF(tracers(iq)%name == 'H2O'//phases_sep//'l') q(:,:,iq)=1.e-15
-              IF(tracers(iq)%name == 'H2Ov') q(:,:,iq)=1.e-10
-              IF(tracers(iq)%name == 'H2Ol') q(:,:,iq)=1.e-15
+              IF(tracers(iq)%name == addPhase('H2O', 'g')) q(:,:,iq)=1.e-10
+              IF(tracers(iq)%name == addPhase('H2O', 'l')) q(:,:,iq)=1.e-15
 
               ! CRisi: init des isotopes
               ! distill de Rayleigh très simplifiée
-              iName = tracers(iq)%iso_iName
+!             iName    = tracers(iq)%iso_iName  ! (next commit)
+              iName    = iso_num(iq)
               if (.NOT.ok_isotopes .OR. iName <= 0) CYCLE
-              iZone    = tracers(iq)%iso_iZone
               iPhase   = tracers(iq)%iso_iPhase
               iqParent = tracers(iq)%iqParent
-              if (iZone == 0) q(:,:,iq) = q(:,:,iqParent)*tnat(iName) &
-                                        *(q(:,:,iqParent)/30.e-3)**(alpha_ideal(iName)-1)
-              if (iZone == 1) q(:,:,iq) = q(:,:,iqiso(iso_indnum(iq),iPhase))
+              IF(tracers(iq)%iso_iZone == 0) THEN
+                 q(:,:,iq) = q(:,:,iqParent)*tnat(iName)*(q(:,:,iqParent)/30.e-3)**(alpha_ideal(iName)-1.)
+              ELSE
+                 q(:,:,iq) = q(:,:,iqiso(iso_indnum(iq),iPhase))
+              END IF
            enddo
         else
Index: LMDZ6/trunk/libf/dyn3d/leapfrog.F
===================================================================
--- LMDZ6/trunk/libf/dyn3d/leapfrog.F	(revision 4119)
+++ LMDZ6/trunk/libf/dyn3d/leapfrog.F	(revision 4120)
@@ -451,5 +451,5 @@
 c+jld
 
-c  Diagnostique de conservation de l'énergie : initialisation
+c  Diagnostique de conservation de l'energie : initialisation
          IF (ip_ebil_dyn.ge.1 ) THEN 
           ztit='bil dyn'
@@ -498,5 +498,5 @@
        
 c
-c  Diagnostique de conservation de l'énergie : difference
+c  Diagnostique de conservation de l'energie : difference
          IF (ip_ebil_dyn.ge.1 ) THEN 
           ztit='bil phys'
Index: LMDZ6/trunk/libf/dyn3d_common/infotrac.F90
===================================================================
--- LMDZ6/trunk/libf/dyn3d_common/infotrac.F90	(revision 4119)
+++ LMDZ6/trunk/libf/dyn3d_common/infotrac.F90	(revision 4120)
@@ -3,11 +3,9 @@
 MODULE infotrac
 
-   USE       strings_mod, ONLY: msg, find, strIdx,  strFind, strParse, dispTable, int2str,  reduceExpr,  &
+   USE       strings_mod, ONLY: msg, find, strIdx,  strFind, strParse, dispTable, int2str,  reduceExpr, &
                           cat, fmsg, test, strTail, strHead, strStack, strReduce, bool2str, maxlen, testFile
-   USE readTracFiles_mod, ONLY: trac_type, readTracersFiles, addPhase,  phases_sep,  nphases, ancestor,  &
-                                isot_type, readIsotopesFile, delPhase,   old_phases, getKey_init, tran0, &
-                                keys_type, initIsotopes,  indexUpdate, known_phases, getKey, setGeneration, &
-                                new2oldPhase
-
+   USE readTracFiles_mod, ONLY: trac_type, readTracersFiles, addPhase, indexUpdate,  nphases, ancestor,  &
+                                isot_type, old2newName,      delPhase,               getKey_init, tran0, &
+                                keys_type, initIsotopes,     getPhase, known_phases, getKey, setGeneration
    IMPLICIT NONE
 
@@ -23,19 +21,21 @@
    PUBLIC :: isotopes,  nbIso                              !--- Derived type, full isotopes families database + nb of families
    PUBLIC :: isoSelect, ixIso                              !--- Isotopes family selection tool + selected family index
-   PUBLIC :: min_qParent, min_qMass, min_ratio             !--- Min. values for various isotopic quantities
    !=== FOR ISOTOPES: Specific to water
    PUBLIC :: iH2O, tnat, alpha_ideal                       !--- H2O isotopes index, natural abundance, fractionning coeff.
+   PUBLIC :: min_qParent, min_qMass, min_ratio             !--- Min. values for various isotopic quantities
    !=== FOR ISOTOPES: Depending on the selected isotopes family
    PUBLIC :: isotope, isoKeys                              !--- Selected isotopes database + associated keys (cf. getKey)
    PUBLIC :: isoName, isoZone, isoPhas                     !--- Isotopes and tagging zones names, phases
    PUBLIC :: niso,    nzone,   nphas,   ntiso              !---  " " numbers + isotopes & tagging tracers number
-   PUBLIC :: iZonIso, iTraPha                              !--- 2D index tables to get "iq" index
+   PUBLIC :: itZonIso, index_trac                          !--- idx "it" (in "isoName(1:niso)") = function(tagging idx, isotope idx)
+   PUBLIC :: iqTraPha, iqiso                               !--- idx "iq" (in "qx") = function(isotope idx, phase idx) + aliases
    PUBLIC :: isoCheck                                      !--- Run isotopes checking routines
    !=== FOR BOTH TRACERS AND ISOTOPES
    PUBLIC :: getKey                                        !--- Get a key from "tracers" or "isotope"
 
-   PUBLIC :: ntraciso, ntraceurs_zone, iqiso
-   PUBLIC :: ok_isotopes, ok_iso_verif, ok_isotrac, ok_init_iso, use_iso
-   PUBLIC :: index_trac, iso_indnum, indnum_fn_num, niso_possibles
+   !=== OLD QUANTITIES OR ALIASES FOR OLDER NAMES (TO BE REMOVED SOON)
+   PUBLIC :: ntraciso, ntraceurs_zone
+   PUBLIC :: ok_isotopes, ok_iso_verif, ok_isotrac, use_iso
+   PUBLIC :: iso_num, iso_indnum, indnum_fn_num, niso_possibles
    PUBLIC :: qperemin, masseqmin, ratiomin
 
@@ -93,19 +93,15 @@
 !=== DERIVED TYPE EMBEDDING MOST OF THE ISOTOPES-RELATED QUANTITIES (LENGTH: nbIso, NUMBER OF ISOTOPES FAMILIES)
 !    Each entry is accessible using "%" sign.
-!  |-----------------+--------------------------------------------------+----------------+-----------------+
-!  |  entry | length | Meaning                                          | Former name    | Possible values |
-!  |-----------------+--------------------------------------------------+----------------+-----------------+
-!  | parent          | Parent tracer (isotopes family name)             |                |                 |
-!  | keys   | niso   | Isotopes keys/values pairs list + number         |                |                 |
-!  | trac   | ntiso  | Isotopes + tagging tracers list + number         |                |                 |
-!  | zone   | nzone  | Geographic tagging zones   list + number         |                |                 |
-!  | phase  | nphas  | Phases                     list + number         |                | [g][l][s], 1:3  |
-!  | niso            | Number of isotopes, excluding tagging tracers    |                |                 |
-!  | ntiso           | Number of isotopes, including tagging tracers    | ntraciso       |                 |
-!  | nzone           | Number of geographic tagging zones               | ntraceurs_zone |                 |
-!  | nphas           | Number of phases                                 |                |                 |
-!  | iTraPha         | Index in "trac(1:niso)" = f(name(1:ntiso)),phas) | iqiso          | 1:niso          |
-!  | iZonIso         | Index in "trac(1:ntiso)" = f(zone, name(1:niso)) | index_trac     | 1:nzone         |
-!  |-----------------+--------------------------------------------------+----------------+-----------------+
+!  |-----------------+--------------------------------------------------+--------------------+-----------------+
+!  |  entry | length | Meaning                                          |    Former name     | Possible values |
+!  |-----------------+--------------------------------------------------+--------------------+-----------------+
+!  | parent          | Parent tracer (isotopes family name)             |                    |                 |
+!  | keys   | niso   | Isotopes keys/values pairs list + number         |                    |                 |
+!  | trac   | ntiso  | Isotopes + tagging tracers list + number         | / | ntraciso       |                 |
+!  | zone   | nzone  | Geographic tagging zones   list + number         | / | ntraceurs_zone |                 |
+!  | phase  | nphas  | Phases                     list + number         |                    | [g][l][s], 1:3  |
+!  | iqTraPha        | Index in "qx"           = f(name(1:ntiso)),phas) | iqiso              | 1:nqtot         |
+!  | itZonIso        | Index in "trac(1:ntiso)"= f(zone, name(1:niso))  | index_trac         | 1:ntiso         |
+!  +-----------------+--------------------------------------------------+--------------------+-----------------+
 
    REAL, PARAMETER :: min_qParent = 1.e-30, min_qMass = 1.e-18, min_ratio = 1.e-16 ! MVals et CRisi
@@ -127,33 +123,32 @@
    TYPE(isot_type),         SAVE, POINTER :: isotope            !--- CURRENTLY SELECTED ISOTOPES FAMILY DESCRIPTOR
    INTEGER,                 SAVE          :: ixIso, iH2O        !--- Index of the selected isotopes family and H2O family
-   LOGICAL,                 SAVE          :: isoCheck           !--- Flag to trigger the checking routines
+   LOGICAL,                 SAVE, POINTER :: isoCheck           !--- Flag to trigger the checking routines
    TYPE(keys_type),         SAVE, POINTER :: isoKeys(:)         !--- ONE SET OF KEYS FOR EACH ISOTOPE (LISTED IN isoName)
    CHARACTER(LEN=maxlen),   SAVE, POINTER :: isoName(:),   &    !--- ISOTOPES NAMES FOR THE CURRENTLY SELECTED FAMILY
                                              isoZone(:),   &    !--- TAGGING ZONES  FOR THE CURRENTLY SELECTED FAMILY
                                              isoPhas            !--- USED PHASES    FOR THE CURRENTLY SELECTED FAMILY
-   INTEGER, TARGET,         SAVE          ::  niso, nzone, &    !--- NUMBER OF ISOTOPES, TAGGING ZONES AND PHASES
-                                             nphas, ntiso       !--- NUMBER OF PHASES AND ISOTOPES + ISOTOPIC TAGGING TRACERS
-   INTEGER,                 SAVE, POINTER :: iZonIso(:,:)       !--- INDEX IN "isoTrac" AS f(tagging zone, isotope)
-   INTEGER,                 SAVE, POINTER :: iTraPha(:,:)       !--- INDEX IN "isoTrac" AS f(isotopic tracer, phase)
-   INTEGER, ALLOCATABLE,    SAVE ::  index_trac(:,:) ! numero ixt en fn izone, indnum entre 1 et niso
-   INTEGER, ALLOCATABLE,    SAVE ::  iqiso(:,:)      ! donne indice iq en fn de (ixt,phase) 
-
-   !--- Aliases for older names
-   INTEGER, POINTER, SAVE :: ntraciso, ntraceurs_zone
-   REAL,             SAVE :: qperemin, masseqmin, ratiomin
-
-! CRisi: cas particulier des isotopes
-   INTEGER, PARAMETER :: niso_possibles = 5
-   LOGICAL, SAVE      :: ok_isotopes, ok_iso_verif, ok_isotrac, ok_init_iso
+   INTEGER,                 SAVE, POINTER ::  niso, nzone, &    !--- NUMBER OF ISOTOPES, TAGGING ZONES AND PHASES
+                                             nphas, ntiso, &    !--- NUMBER OF PHASES AND ISOTOPES + ISOTOPIC TAGGING TRACERS
+                                            itZonIso(:,:), &    !--- INDEX IN "isoTrac" AS f(tagging zone idx,  isotope idx)
+                                            iqTraPha(:,:)       !--- INDEX IN "qx"      AS f(isotopic tracer idx, phase idx)
+
+   !--- Aliases for older names + quantities to be removed soon
+   INTEGER,                 SAVE, POINTER ::  index_trac(:,:)   ! numero ixt en fn izone, indnum entre 1 et niso
+   INTEGER,                 SAVE, POINTER ::  iqiso(:,:)        ! donne indice iq en fn de (ixt,phase) 
+   INTEGER,                 SAVE, POINTER :: ntraciso, ntraceurs_zone
+   REAL,    SAVE :: qperemin, masseqmin, ratiomin
+   INTEGER, SAVE :: niso_possibles
+   LOGICAL, SAVE :: ok_isotopes, ok_iso_verif, ok_isotrac
    LOGICAL, SAVE, ALLOCATABLE ::       use_iso(:)
-   INTEGER, SAVE, ALLOCATABLE ::    iso_indnum(:)     !--- Gives 1<=idx<=niso_possibles as function(1<=iq <=nqtot)
-   INTEGER, SAVE, ALLOCATABLE :: indnum_fn_num(:)     !--- Gives 1<=idx<=niso           as function(1<=idx<=niso_possibles)
+   INTEGER, SAVE, ALLOCATABLE ::       iso_num(:)               !--- idx in [1,niso_possibles] = f(1<=iq <=nqtot)
+   INTEGER, SAVE, ALLOCATABLE ::    iso_indnum(:)               !--- idx in [1,niso]           = f(1<=iq <=nqtot)
+   INTEGER, SAVE, ALLOCATABLE :: indnum_fn_num(:)               !--- idx in [1,niso]           = f(1<=idx<=niso_possibles)
 
    !=== VARIABLES FOR ISOTOPES INITIALIZATION AND FOR INCA
-   REAL,               SAVE, ALLOCATABLE ::     tnat(:),  &     !--- Natural relative abundance of water isotope        (niso)
-                                         alpha_ideal(:)         !--- Ideal fractionning coefficient (for initial state) (niso)
-   INTEGER,            SAVE, ALLOCATABLE :: conv_flg(:),  &     !--- Convection     activation ; needed for INCA        (nbtr)
-                                             pbl_flg(:)         !--- Boundary layer activation ; needed for INCA        (nbtr)
-   CHARACTER(LEN=8),   SAVE, ALLOCATABLE ::   solsym(:)         !--- Names from INCA                                    (nbtr)
+   REAL,                SAVE, ALLOCATABLE ::     tnat(:), &     !--- Natural relative abundance of water isotope        (niso)
+                                          alpha_ideal(:)        !--- Ideal fractionning coefficient (for initial state) (niso)
+   INTEGER,             SAVE, ALLOCATABLE :: conv_flg(:), &     !--- Convection     activation ; needed for INCA        (nbtr)
+                                              pbl_flg(:)        !--- Boundary layer activation ; needed for INCA        (nbtr)
+   CHARACTER(LEN=8),    SAVE, ALLOCATABLE ::   solsym(:)        !--- Names from INCA                                    (nbtr)
    LOGICAL, PARAMETER :: lOldCode = .TRUE.
 
@@ -175,5 +170,5 @@
 !   05/94: F.Forget      Modif special traceur
 !   02/02: M-A Filiberti Lecture de traceur.def
-!   01/22: D. Cugnet     Nouveaux tracer.def et tracer_*.def + encapsulation (types tr et iso)
+!   01/22: D. Cugnet     Nouveaux tracer.def et tracer_*.def + encapsulation (types trac_type et isot_type)
 !
 !   Objet:
@@ -212,8 +207,7 @@
    TYPE(isot_type), POINTER             :: iso
 
-   CHARACTER(LEN=maxlen), ALLOCATABLE :: tnom_0(:), tnom_transp(:)        !--- Tracer short name + transporting fluid name
+   CHARACTER(LEN=maxlen), ALLOCATABLE :: tnom_0(:), tnom_transp(:)   !--- Tracer short name + transporting fluid name
    CHARACTER(LEN=maxlen)              :: tchaine
    INTEGER :: ierr
-   LOGICAL :: lINCA
 
    CHARACTER(LEN=*), PARAMETER :: modname="infotrac_init"
@@ -238,5 +232,5 @@
       !--- MESSAGE ABOUT THE CHOSEN CONFIGURATION
       msg1 = 'For type_trac = "'//TRIM(str(it))//'":'
-      SELECT CASE(type_trac)
+      SELECT CASE(str(it))
          CASE('inca'); CALL msg(TRIM(msg1)//' coupling with INCA chemistry model, config_inca='//config_inca, modname)
          CASE('inco'); CALL msg(TRIM(msg1)//' coupling jointly with INCA and CO2 cycle',  modname)
@@ -254,5 +248,5 @@
       !--- COHERENCE TEST BETWEEN "type_trac" AND PREPROCESSING KEYS
       SELECT CASE(str(it))
-         CASE('inca','inco')
+         CASE('inca', 'inco')
 #ifndef INCA
             CALL abort_gcm(modname, 'You must add cpp key INCA and compile with INCA code', 1)
@@ -283,5 +277,5 @@
 
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-   IF(lOldCode) THEN
+   IF(lOldCode) THEN         !--- "type_trac" is a single keyword => no need to loop on its parsed version "str(:)"
 !------------------------------------------------------------------------------------------------------------------------------
    !--- Determine nqtrue and (INCA only) nqo, nbtr
@@ -289,7 +283,6 @@
    IF(ierr /= 0) CALL abort_gcm(modname, 'file "traceur.def" not found !', 1)
    CALL msg('File "traceur.def" successfully opened.', modname)
-   lINCA = ANY(['inca','inco'] == type_trac)
-
-   IF(lINCA) THEN
+
+   IF(ANY(['inca','inco'] == type_trac)) THEN
 #ifdef INCA
       READ(90,*) nqo
@@ -299,18 +292,7 @@
       nqtrue = nbtr + nqo
       IF(ALL([2,3] /= nqo)) CALL abort_gcm(modname, 'Only 2 or 3 water phases allowed ; found nqo='//TRIM(int2str(nqo)), 1)
-      CALL msg('nqo    = '//TRIM(int2str(nqo)),    modname)
-      CALL msg('nbtr   = '//TRIM(int2str(nbtr)),   modname)
-      CALL msg('nqtrue = '//TRIM(int2str(nqtrue)), modname)
-      CALL msg('nqCO2  = '//TRIM(int2str(nqCO2)),  modname)
-      CALL msg('nqINCA = '//TRIM(int2str(nqINCA)), modname)
       ALLOCATE(hadv_inca(nqINCA), conv_flg_inca(nqINCA), solsym_inca(nqINCA))
       ALLOCATE(vadv_inca(nqINCA),  pbl_flg_inca(nqINCA))
       CALL init_transport(hadv_inca, vadv_inca, conv_flg_inca, pbl_flg_inca, solsym_inca)
-      ! DC passive CO2 tracer is at position 1: H2O was removed ; nqCO2/=0 in "inco" case only
-      ALLOCATE(conv_flg(nbtr),pbl_flg(nbtr),solsym(nbtr))
-      conv_flg = [(  1,        ic=1, nqCO2),conv_flg_inca]
-       pbl_flg = [(  1,        ic=1, nqCO2), pbl_flg_inca]
-       solsym  = [('CO2     ', ic=1, nqCO2), solsym_inca]
-      DEALLOCATE(conv_flg_inca, pbl_flg_inca)
 #endif
    ELSE
@@ -337,14 +319,10 @@
       END IF
 #endif
-      CALL msg('237: iq='//TRIM(int2str(iq)), modname)
       READ(90,'(I2,X,I2,X,A)',IOSTAT=ierr) hadv(iq),vadv(iq),tchaine
-      WRITE(msg1,'("hadv(",i0,"), vadv(",i0,") = ",i0,", ",i0)')iq, iq, hadv(iq), vadv(iq)
-      CALL msg(TRIM(msg1), modname)
-      CALL msg('tchaine = "'//TRIM(tchaine)//'"', modname)
-      CALL msg('infotrac 238: IOstatus='//TRIM(int2str(ierr)), modname)
       IF(ierr/=0) CALL abort_gcm('infotrac_init', 'Pb dans la lecture de traceur.def', 1)
       jq = INDEX(tchaine(1:LEN_TRIM(tchaine)),' ')
       CALL msg("Ancienne version de traceur.def: traceurs d'air uniquement", modname, iq==1 .AND. jq==0)
       CALL msg("Nouvelle version de traceur.def",                            modname, iq==1 .AND. jq/=0)
+      CALL msg('iq, hadv, vadv, tchaine ='//TRIM(strStack(int2str([iq, hadv(iq), vadv(iq)])))//', '//TRIM(tchaine), modname)
       IF(jq /= 0) THEN                                               !--- Space in the string chain => new format
          tnom_0     (iq) = tchaine(1:jq-1)
@@ -354,43 +332,34 @@
          tnom_transp(iq) = 'air'
       END IF
-      CALL msg(     'tnom_0(iq)=<'//TRIM(tnom_0(iq))     //'>', modname)
-      CALL msg('tnom_transp(iq)=<'//TRIM(tnom_transp(iq))//'>', modname)
-   END DO
-#ifdef INCA
-   DEALLOCATE(solsym_inca) 
-#endif
-
+   END DO
    CLOSE(90)
 
 #ifndef INCA
-   CALL msg('Valeur de traceur.def :', modname)
-   CALL msg('nombre total de traceurs '//TRIM(int2str(nqtrue)), modname)
-   DO iq = 1, nqtrue
-      CALL msg(strStack([int2str(hadv(iq)), int2str(vadv(iq)), tnom_0(iq), tnom_transp(iq)]), modname)
-   END DO
    IF(planet_type /= 'earth') nqo = 0                                !--- Same number of tracers in dynamics and physics
    IF(planet_type == 'earth') nqo = COUNT(delPhase(tnom_0) == 'H2O') !--- for all planets except for Earth
-   nbtr = nqtrue - nqo               
-   ALLOCATE(conv_flg(nbtr),pbl_flg(nbtr),solsym(nbtr))
-   conv_flg(1:nbtr) = 1                                     !--- Convection activated for all tracers
-   pbl_flg(1:nbtr) = 1                                     !--- Boundary layer activated for all tracers
-#endif
+   nbtr = nqtrue - nqo
+#endif
+
+   CALL msg('RAW CONTENT OF "traceur.def" FILE:', modname)
+   IF(dispTable('iiss', ['hadv  ', 'vadv  ', 'name  ', 'parent'], cat(tnom_0, tnom_transp), cat(hadv, vadv))) &
+      CALL abort_gcm(modname, "problem with the tracers table content", 1)
 
    !--- SET FIELDS %name, %parent, %phase, %component
-   tracers(:)%name      = tnom_0
-   tracers(:)%parent    = tnom_transp
-   tracers(:)%phase     = 'g'
+   tracers(:)%name      = old2newName(tnom_0)
+   tracers(:)%parent    = old2newName(tnom_transp)
+   tracers(:)%phase     = [( getPhase(tracers(iq)%name), iq=1, nqtrue )]
    tracers(:)%component = type_trac
-
-
    DO iq = 1, nqtrue
-      ip = strIdx([(addPhase('H2O',old_phases(ix:ix),''), ix=1, nphases)], strHead(tracers(iq)%name,'_',.TRUE.))
-      IF(ip == 0) CYCLE
-      tracers(iq)%phase = known_phases(ip:ip)
-      tracers(iq)%component = 'lmdz'
-   END DO
-   IF(lINCA) tracers(1+nqo:nqCO2+nqo)%component = 'co2i'
+      IF(addPhase('H2O',tracers(iq)%phase) == tracers(iq)%name) tracers(iq)%component = 'lmdz'
+   END DO
+   IF(ANY(['inca','inco'] == type_trac)) tracers(1+nqo:nqCO2+nqo)%component = 'co2i'
    CALL setGeneration(tracers)                                       !--- SET FIELDS %iGeneration, %gen0Name
-! manque "type"
+   WHERE(tracers(:)%iGeneration == 2) tracers(:)%type = 'tag'        !--- DEFAULT VALUE: "tracer"
+
+   !--- FINALIZE
+   DEALLOCATE(tnom_0, tnom_transp)
+#ifdef INCA
+   DEALLOCATE(hadv_inca, vadv_inca, conv_flag_inca, pbl_flag_inca, solsym_inca)
+#endif
 
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -400,5 +369,5 @@
    IF(fType == 0) CALL abort_gcm(modname, 'Missing tracers file: "traceur.def", "tracer.def" or "tracer_<keyword>.def file.',1)
    !---------------------------------------------------------------------------------------------------------------------------
-   IF(fType == 1) THEN                                               !=== FOUND AN OLD STYLE "traceur.def"
+   IF(fType == 1 .AND. ANY(['inca','inco'] == type_trac)) THEN       !=== FOUND OLD STYLE INCA "traceur.def" (single type_trac)
    !---------------------------------------------------------------------------------------------------------------------------
 #ifdef INCA
@@ -409,18 +378,11 @@
       nqtrue = nbtr + nqo                                            !--- Total number of "true" tracers
       IF(ALL([2,3] /= nqo)) CALL abort_gcm(modname, 'Only 2 or 3 water phases allowed ; found nqo='//TRIM(int2str(nqo)), 1)
-      CALL msg('nqo    = '//TRIM(int2str(nqo)),    modname)
-      CALL msg('nbtr   = '//TRIM(int2str(nbtr)),   modname)
-      CALL msg('nqtrue = '//TRIM(int2str(nqtrue)), modname)
-      CALL msg('nqCO2  = '//TRIM(int2str(nqCO2)),  modname)
-      CALL msg('nqINCA = '//TRIM(int2str(nqINCA)), modname)
-      ALLOCATE(hadv(nqtrue), hadv_inca(nqINCA), conv_flg_inca(nqINCA), solsym_inca(nqINCA))
-      ALLOCATE(vadv(nqtrue), vadv_inca(nqINCA),  pbl_flg_inca(nqINCA))
+      ALLOCATE(hadv(nqtrue), conv_flg(nbtr), hadv_inca(nqINCA), conv_flg_inca(nqINCA), solsym_inca(nqINCA))
+      ALLOCATE(vadv(nqtrue),  pbl_flg(nbtr), vadv_inca(nqINCA),  pbl_flg_inca(nqINCA), solsym(nbtr))
       CALL init_transport(hadv_inca, vadv_inca, conv_flg_inca, pbl_flg_inca, solsym_inca)
-      ! DC passive CO2 tracer is at position 1: H2O was removed ; nqCO2/=0 in "inco" case only
-      
-      conv_flg = [(  1       , k=1, nqCO2), conv_flg_inca]
-      pbl_flg  = [(  1       , k=1, nqCO2), pbl_flg_inca]
-      solsym   = [('CO2     ', k=1, nqCO2), solsym_inca]
-      DEALLOCATE(conv_flg_inca, pbl_flg_inca, solsym_inca)
+      !--- Passive CO2 tracer is at position 1 because: H2O has been removed ; nqCO2/=0 in "inco" case only
+      conv_flg(1:nbtr) = [(1,          k=1, nqCO2), conv_flg_inca]
+       pbl_flg(1:nbtr) = [(1,          k=1, nqCO2),  pbl_flg_inca]
+       solsym (1:nbtr) = [('CO2     ', k=1, nqCO2),   solsym_inca]
       ALLOCATE(ttr(nqtrue))
       ttr(1:nqo+nqCO2)                    = tracers
@@ -433,28 +395,21 @@
       lerr = getKey('hadv', had, ky=ttr(:)%keys); hadv(:) = [had, hadv_inca]
       lerr = getKey('vadv', vad, ky=ttr(:)%keys); vadv(:) = [vad, vadv_inca]
-      DEALLOCATE(had, hadv_inca, vad, vadv_inca)
       CALL MOVE_ALLOC(FROM=ttr, TO=tracers)
       CALL setGeneration(tracers)                                    !--- SET FIELDS %iGeneration, %gen0Name
-#else
-      nqo    = COUNT(delPhase(tracers(:)%name) == 'H2O')             !--- Number of water phases
+      DEALLOCATE(had, hadv_inca, vad, vadv_inca, conv_flg_inca, pbl_flg_inca, solsym_inca)
+#endif
+   !---------------------------------------------------------------------------------------------------------------------------
+   ELSE                                                              !=== FOUND NEW STYLE TRACERS CONFIGURATION FILE(S)
+   !---------------------------------------------------------------------------------------------------------------------------
+      nqo    =        COUNT(delPhase(tracers(:)%name)     == 'H2O' &
+                               .AND. tracers(:)%component == 'lmdz') !--- Number of water phases
       nqtrue = SIZE(tracers)                                         !--- Total number of "true" tracers
-      nbtr   = nqtrue - nqo                                          !--- Number of tracers passed to phytrac
+      nbtr   = nqtrue-COUNT(delPhase(tracers(:)%gen0Name) == 'H2O' &
+                               .AND. tracers(:)%component == 'lmdz') !--- Number of tracers passed to phytrac
       lerr = getKey('hadv', hadv, ky=tracers(:)%keys)
       lerr = getKey('vadv', vadv, ky=tracers(:)%keys)
-      ALLOCATE(solsym(nbtr))
-      conv_flg(1:nbtr)=1  !--- Convection activated for all tracers
-      pbl_flg(1:nbtr)=1   !--- Boundary layer activated for all tracers
-#endif
-   !---------------------------------------------------------------------------------------------------------------------------
-   ELSE                                                              !=== FOUND NEW STYLE TRACERS CONFIGURATION FILE(S)
-   !---------------------------------------------------------------------------------------------------------------------------
-      nqo    = COUNT(delPhase(tracers(:)%name) == 'H2O')             !--- Number of water phases
-      nqtrue = SIZE(tracers)                                         !--- Total number of "true" tracers
-      nbtr   = nqtrue - nqo                                          !--- Number of tracers passed to phytrac
-      lerr = getKey('hadv', hadv, ky=tracers(:)%keys)
-      lerr = getKey('vadv', vadv, ky=tracers(:)%keys)
-      ALLOCATE(solsym(nbtr))
-      conv_flg(1:nbtr)=1  !--- Convection activated for all tracers
-       pbl_flg(1:nbtr)=1  !--- Boundary layer activated for all tracers
+      ALLOCATE(conv_flg(nbtr), pbl_flg(nbtr), solsym(nbtr))
+      conv_flg(1:nbtr) = [(1, it=1, nbtr)]                           !--- Convection activated for all tracers
+       pbl_flg(1:nbtr) = [(1, it=1, nbtr)]                           !--- Boundary layer activated for all tracers
    !---------------------------------------------------------------------------------------------------------------------------
    END IF
@@ -488,8 +443,4 @@
       CALL msg('The total number of tracers needed is '//TRIM(int2str(nqtot)))
    END IF
-   CALL msg('nqo    = '//TRIM(int2str(nqo)),    modname)
-   CALL msg('nbtr   = '//TRIM(int2str(nbtr)),   modname)
-   CALL msg('nqtrue = '//TRIM(int2str(nqtrue)), modname)
-   CALL msg('nqtot  = '//TRIM(int2str(nqtot)),  modname)
 
 !==============================================================================================================================
@@ -527,12 +478,13 @@
       t1%iadv       = iad
       t1%isAdvected = iad >= 0
-      t1%isInPhysics= .not. (delPhase(t1%gen0Name) == 'H2O' .and. t1%component=='lmdz')  !=== TO BE COMPLETED WITH OTHER EXCEPTIONS: CO2i, SURSATURATED CLOUDS...
+      t1%isInPhysics= delPhase(t1%gen0Name) /= 'H2O' &
+                          .OR. t1%component /= 'lmdz' !=== OTHER EXCEPTIONS TO BE ADDED: CO2i, SURSATURATED WATER CLOUD...
       ttr(iq)       = t1
 
       !--- DEFINE THE HIGHER ORDER TRACERS, IF ANY
       nm = 0
-      IF(iad == 20) nm = 3                                             !--- 2nd order scheme
-      IF(iad == 30) nm = 9                                             !--- 3rd order scheme
-      IF(nm == 0) CYCLE                                                !--- No higher moments
+      IF(iad == 20) nm = 3                                           !--- 2nd order scheme
+      IF(iad == 30) nm = 9                                           !--- 3rd order scheme
+      IF(nm == 0) CYCLE                                              !--- No higher moments
       ttr(jq+1:jq+nm)             = t1
       ttr(jq+1:jq+nm)%name        = [(TRIM(t1%name)    //'-'//TRIM(suff(im)), im=1, nm) ]
@@ -549,7 +501,4 @@
    CALL indexUpdate(tracers)
 
-   CALL msg('Information stored in infotrac :', modname)
-   CALL msg('iadv  name  long_name :', modname)
-
    !=== TEST ADVECTION SCHEME
    DO iq=1,nqtot ; t1 => tracers(iq); iad = t1%iadv
@@ -568,5 +517,5 @@
 
       !--- ONLY VALID SCHEME NUMBER FOR WATER VAPOUR:            iadv = 14
-      ll = t1%name /= addPhase('H2O','g'); IF(lOldCode) ll = t1%name /= 'H2Ov'
+      ll = t1%name /= addPhase('H2O','g')
       IF(fmsg('WARNING ! iadv=14 is valid for water vapour only. Setting iadv=10 for "'//TRIM(t1%name)//'".', &
          modname, iad == 14 .AND. ll))                 t1%iadv = 10
@@ -578,71 +527,101 @@
    CALL infotrac_isoinit                    !--- SET FIELDS %type, %iso_iName, %iso_iZone, %iso_iPhase
    CALL getKey_init(tracers, isotopes)
-   IF(isoSelect('H2O')) RETURN                                    !--- Select water isotopes ; finished if no water isotopes
-   iH2O = ixIso                                                   !--- Keep track of water family index
+   IF(isoSelect('H2O')) RETURN              !--- Select water isotopes ; finished if no water isotopes
+   iH2O = ixIso                             !--- Keep track of water family index
 
    !--- Remove the isotopic tracers from the tracers list passed to phytrac
    nbtr    = nbtr -nqo*   ntiso             !--- ISOTOPIC TAGGING TRACERS ARE NOT PASSED TO THE PHYSICS
    nqtottr = nqtot-nqo*(1+ntiso)            !--- NO H2O-FAMILY    TRACER  IS      PASSED TO THE PHYSICS
-   CALL msg('702: nbtr, ntiso='//strStack(int2str([nbtr, ntiso])), modname)
-   CALL msg('704: nqtottr, nqtot, nqo = '//strStack(int2str([nqtottr, nqtot, nqo])), modname)
-   ! Rq: nqtottr n'est pas forcement egal a nbtr dans le cas ou nmom/=0
-   IF(COUNT(tracers%iso_iName == 0) - COUNT(delPhase(tracers(:)%name) == 'H2O' .AND. tracers(:)%component=='lmdz') /= nqtottr) &
-      CALL abort_gcm('infotrac_init', 'pb dans le calcul de nqtottr', 1)
-
-   !--- Finalize :
-   DEALLOCATE(tnom_0, tnom_transp)
+
+   ALLOCATE(conv_flg(nbtr), pbl_flg(nbtr), solsym(nbtr))
+#ifndef INCA
+   conv_flg(1:nbtr) = 1                                              !--- Convection activated for all tracers
+    pbl_flg(1:nbtr) = 1                                              !--- Boundary layer activated for all tracers
+#else
+   !--- Passive CO2 tracer is at position 1 because: H2O has been removed ; nqCO2/=0 in "inco" case only
+   conv_flg(1:nbtr) = [(  1,        ic=1, nqCO2),conv_flg_inca]
+    pbl_flg(1:nbtr) = [(  1,        ic=1, nqCO2), pbl_flg_inca]
+     solsym(1:nbtr) = [('CO2     ', ic=1, nqCO2),  solsym_inca]
+#endif
 
 ELSE
 
    CALL initIsotopes(tracers, isotopes)
-   nbIso = SIZE(isotopes); IF(nbIso==0) RETURN                    !--- No isotopes: finished.
-
-   !=== READ PHYSICAL PARAMETERS FROM "isotopes_params.def" FILE SPECIFIC TO WATER ISOTOPES
-   !    DONE HERE, AND NOT ONLY IN "infotrac_phy", BECAUSE SOME PHYSICAL PARAMS ARE NEEDED FOR RESTARTS (tnat, alpha_ideal)
-   CALL getKey_init(tracers, isotopes)
-   IF(isoSelect('H2O')) RETURN                                    !--- Select water isotopes ; finished if no water isotopes
-   iH2O = ixIso                                                   !--- Keep track of water family index
-   IF(getKey('tnat' , tnat,        isoName(1:niso))) CALL abort_gcm(modname, 'can''t read "tnat"', 1)
-   IF(getKey('alpha', alpha_ideal, isoName(1:niso))) CALL abort_gcm(modname, 'can''t read "alpha_ideal"', 1)
-
-   !=== ENSURE THE MEMBERS OF AN ISOTOPES FAMILY ARE PRESENT IN THE SAME PHASES
-   DO ix = 1, nbIso
-      iso => isotopes(ix)
-      !--- Check whether each isotope and tagging isotopic tracer is present in the same number of phases
-      DO it = 1, iso%ntiso
-         np = SUM([(COUNT(tracers(:)%name == addPhase(iso%trac(it), iso%phase(ip:ip))), ip=1, iso%nphas)])
-         IF(np == iso%nphas) CYCLE
-         WRITE(msg1,'("Found ",i0," phases for ",s," instead of ",i0)')np, iso%trac(it), iso%nphas
-         CALL abort_gcm(modname, msg1, 1)
+   nbIso = SIZE(isotopes)
+   nqtottr = nqtot - COUNT(tracers%gen0Name == 'H2O' .AND. tracers%component == 'lmdz')
+   IF(nbIso/=0) THEN                        !--- ISOTOPES FOUND
+
+      !=== READ PHYSICAL PARAMETERS FROM "isotopes_params.def" FILE SPECIFIC TO WATER ISOTOPES
+      !    DONE HERE, AND NOT ONLY IN "infotrac_phy", BECAUSE SOME PHYSICAL PARAMS ARE NEEDED FOR RESTARTS (tnat, alpha_ideal)
+      CALL getKey_init(tracers, isotopes)
+      IF(isoSelect('H2O')) RETURN           !--- Select water isotopes ; finished if no water isotopes
+      iH2O = ixIso                          !--- Keep track of water family index
+      IF(getKey('tnat' , tnat,        isoName(1:niso))) CALL abort_gcm(modname, 'can''t read "tnat"', 1)
+      IF(getKey('alpha', alpha_ideal, isoName(1:niso))) CALL abort_gcm(modname, 'can''t read "alpha_ideal"', 1)
+
+      !=== MAKE SURE THE MEMBERS OF AN ISOTOPES FAMILY ARE PRESENT IN THE SAME PHASES
+      DO ix = 1, nbIso
+         iso => isotopes(ix)
+         !--- Check whether each isotope and tagging isotopic tracer is present in the same number of phases
+         DO it = 1, iso%ntiso
+            np = SUM([(COUNT(tracers(:)%name == addPhase(iso%trac(it), iso%phase(ip:ip))), ip=1, iso%nphas)])
+            IF(np == iso%nphas) CYCLE
+            WRITE(msg1,'("Found ",i0," phases for ",a," instead of ",i0)')np, TRIM(iso%trac(it)), iso%nphas
+            CALL abort_gcm(modname, msg1, 1)
+         END DO
+         DO it = 1, iso%niso
+            nz = SUM([(COUNT(iso%trac == TRIM(iso%trac(it))//'_'//iso%zone(iz)), iz=1, iso%nzone)])
+            IF(nz == iso%nzone) CYCLE
+            WRITE(msg1,'("Found ",i0," tagging zones for ",a," instead of ",i0)')nz, TRIM(iso%trac(it)), iso%nzone
+            CALL abort_gcm(modname, msg1, 1)
+         END DO
       END DO
-      DO it = 1, iso%niso
-         nz = SUM([(COUNT(iso%trac == iso%trac(it)//'_'//iso%zone(iz)), iz=1, iso%nzone)])
-         IF(nz == iso%nzone) CYCLE
-         WRITE(msg1,'("Found ",i0," tagging zones for ",s," instead of ",i0)')nz, iso%trac(it), iso%nzone
-         CALL abort_gcm(modname, msg1, 1)
-      END DO
-   END DO
-   nqtottr = COUNT(tracers%iso_iName == 0)
+   END IF
 
 END IF
 
-   !=== DISPLAY THE RESULTING LIST
-   t => tracers
-   CALL msg('Information stored in infotrac :')
-   IF(dispTable('isssssssssiiiiiiiii', &
-      ['iq      ', 'name    ', 'longN.  ', 'gen0N.  ', 'parent  ', 'type    ', 'phase   ', 'compon. ', 'isAdv.  ', 'isPhy.  '&
-      ,'iadv    ', 'iGen.   ', 'iqPar.  ', 'nqDes.  ', 'nqChil. ', 'iso_iG. ', 'iso_iN. ', 'iso_iZ. ', 'iso_iP. '],          &
-      cat(t%name,  t%longName,  t%gen0Name,  t%parent,  t%type,  t%phase, &
-          t%component, bool2str(t%isAdvected), bool2str(t%isInPhysics)),  &
-      cat([(iq, iq=1, nqtot)],  t%iadv,  t%iGeneration, t%iqParent, t%nqDescen, &
-         t%nqChilds, t%iso_iGroup, t%iso_iName, t%iso_iZone, t%iso_iPhase))) &
-      CALL abort_gcm(modname, "problem with the tracers table content", 1)
+   !--- Note: nqtottr can differ from nbtr when nmom/=0
+!   IF(COUNT(tracers%iso_iName == 0) - COUNT(delPhase(tracers%name) == 'H2O' .AND. tracers%component == 'lmdz') /= nqtottr) &
+!      CALL abort_gcm('infotrac_init', 'pb dans le calcul de nqtottr', 1)
 
    !--- Some aliases to be removed later
-   ntraciso       => isotope%ntiso
-   ntraceurs_zone => isotope%nzone
+   ntraciso       => ntiso
+   ntraceurs_zone => nzone
    qperemin       =  min_qParent
    masseqmin      =  min_qMass
    ratiomin       =  min_ratio
+   iqiso          => iqTraPha
+   index_trac     => itZonIso
+
+   !=== DISPLAY THE RESULTS
+   CALL msg('nqo    = '//TRIM(int2str(nqo)),    modname)
+   CALL msg('nbtr   = '//TRIM(int2str(nbtr)),   modname)
+   CALL msg('nqtrue = '//TRIM(int2str(nqtrue)), modname)
+   CALL msg('nqtot  = '//TRIM(int2str(nqtot)),  modname)
+   CALL msg('niso   = '//TRIM(int2str(niso)),   modname)
+   CALL msg('ntiso  = '//TRIM(int2str(ntiso)),  modname)
+#ifdef INCA
+   CALL msg('nqCO2  = '//TRIM(int2str(nqCO2)),  modname)
+   CALL msg('nqINCA = '//TRIM(int2str(nqINCA)), modname)
+#endif
+   t => tracers
+   CALL msg('Information stored in infotrac :', modname)
+   IF(dispTable('isssssssssiiiiiiiii', &
+      ['iq    ', 'name  ', 'lName ', 'gen0N ', 'parent', 'type  ', 'phase ', 'compon', 'isAdv ', 'isPhy ', &
+       'iadv  ', 'iGen  ', 'iqPar ', 'nqDes ', 'nqChld', 'iGroup', 'iName ', 'iZone ', 'iPhase'],          &
+      cat(t%name, t%longName, t%gen0Name, t%parent, t%type, t%phase, t%component, bool2str(t%isAdvected),  &
+                                                                                  bool2str(t%isInPhysics)),&
+      cat([(iq, iq=1, nqtot)], t%iadv, t%iGeneration, t%iqParent, t%nqDescen, t%nqChilds, t%iso_iGroup,    &
+                                                    t%iso_iName, t%iso_iZone, t%iso_iPhase), sub=modname)) &
+      CALL abort_gcm(modname, "problem with the tracers table content", 1)
+   IF(niso > 0) THEN
+      CALL msg('Where, for isotopes family "'//TRIM(isotope%parent)//'":', modname)
+      CALL msg('  isoKeys = '//strStack(isoKeys%name), modname)
+      CALL msg('  isoName = '//strStack(isoName),      modname)
+      CALL msg('  isoZone = '//strStack(isoZone),      modname)
+      CALL msg('  isoPhas = '//TRIM(isoPhas),          modname)
+   ELSE
+      CALL msg('No isotopes identified.', modname)
+   END IF
    CALL msg('end', modname)
 
@@ -654,5 +633,5 @@
    !--- Purpose: Set fields %iqParent, %nqChilds, %iGeneration, %iqDescen, %nqDescen (old method)
    USE strings_mod, ONLY: strIdx
-   INTEGER               :: iq, ipere, ifils
+   INTEGER               :: iq, jq, ipere, ifils
    INTEGER, ALLOCATABLE  :: iqfils(:,:)
    CHARACTER(LEN=maxlen) :: msg1, modname='infotrac_init'
@@ -661,4 +640,6 @@
    !=== SET FIELDS %iqParent, %nqChilds
    ALLOCATE(iqfils(nqtot,nqtot)); iqfils(:,:) = 0
+   tracers(:)%nqChilds = 0
+   tracers(:)%iqParent = 0
 
    DO iq = 1, nqtot
@@ -684,7 +665,7 @@
    CALL msg('iqChilds = '//strStack(int2str(PACK(iqfils,MASK=.TRUE.))),modname)
 
-
    !=== SET FIELDS %iGeneration, %iqDescen, %nqDescen
    tracers(:)%iGeneration = 0
+   tracers(:)%nqDescen = 0
    DO iq = 1, nqtot
       ifils = iq
@@ -703,8 +684,13 @@
    END DO
 
+   !=== SET FIELD %gen0Name
+   DO iq = 1, nqtot
+      jq=iq; DO WHILE(tracers(jq)%iGeneration > 0); jq=tracers(jq)%iqParent; END DO
+      tracers(iq)%gen0Name = tracers(jq)%name
+   END DO
+
    CALL msg('nqDescen = '//TRIM(strStack(int2str(tracers(:)%nqDescen))), modname)
    CALL msg('nqDescen_tot = ' //TRIM(int2str(SUM(tracers(:)%nqDescen))), modname)
-   CALL msg('iqChilds = '//strStack(int2str(PACK(iqfils, MASK=.TRUE.))), modname)
-
+   CALL msg('iqDescen = '//strStack(int2str(PACK(iqfils, MASK=.TRUE.))), modname)
 
 END SUBROUTINE infotrac_setHeredity
@@ -719,158 +705,105 @@
    USE ioipsl_getincom
 #endif
+   USE readTracFiles_mod, ONLY: tnom_iso => newH2OIso
    IMPLICIT NONE
-   CHARACTER(LEN=3)      :: tnom_iso(niso_possibles)
-   INTEGER, ALLOCATABLE  :: nb_iso(:,:), nb_traciso(:,:), nb_isoind(:)
-   INTEGER               :: ii, ip, iq, it, iz, ixt, nzone_prec
+   INTEGER, ALLOCATABLE  :: nb_iso(:), nb_tiso(:), nb_zone(:), ix(:)
+   INTEGER               :: ii, ip, iq, it, iz, ixt
    TYPE(isot_type), POINTER :: i
-   TYPE(trac_type), POINTER :: t(:)
-   CHARACTER(LEN=maxlen)    :: tnom_trac
+   TYPE(trac_type), POINTER :: t(:), t1
+   CHARACTER(LEN=maxlen)    :: tnom_trac, modname, t0
    CHARACTER(LEN=maxlen), ALLOCATABLE :: str(:)
    LOGICAL, DIMENSION(:), ALLOCATABLE :: mask
+   REAL, ALLOCATABLE :: tnat0(:), alpha_ideal0(:)
    INCLUDE "iniprint.h"
 
-   tnom_iso = ['eau', 'HDO', 'O18', 'O17', 'HTO']
-   ALLOCATE(nb_iso       (niso_possibles,nqo))
-   ALLOCATE(nb_traciso   (niso_possibles,nqo))
-   ALLOCATE(use_iso      (niso_possibles))
-   ALLOCATE(indnum_fn_num(niso_possibles))
-   ALLOCATE(iso_indnum(nqtot))
-   ALLOCATE(nb_isoind(nqo))
-
-   iso_indnum   (:) = 0
-   use_iso      (:) = .FALSE.
-   indnum_fn_num(:) = 0
-   nb_iso     (:,:) = 0  
-   nb_traciso (:,:) = 0
-   nb_isoind    (:) = 0
-
-   DO iq=1, nqtot
-      IF(delPhase(tracers(iq)%name) == 'H2O' .OR. .NOT.tracers(iq)%isAdvected) CYCLE
-outer:DO ip = 1, nqo
-         DO ixt= 1,niso_possibles
-            tnom_trac = 'H2O'//old_phases(ip:ip)//'_'//TRIM(tnom_iso(ixt))
-            IF (tracers(iq)%name == tnom_trac) THEN
-               nb_iso(ixt,ip)         = nb_iso(ixt,ip)+1
-               nb_isoind (ip)         = nb_isoind (ip)+1
-               tracers(iq)%type       = 'tracer'
-               tracers(iq)%iso_iGroup = 1
-               tracers(iq)%iso_iName  = ixt
-               iso_indnum(iq)         = nb_isoind(ip)
-               indnum_fn_num(ixt)     = iso_indnum(iq)
-               tracers(iq)%iso_iPhase = ip
-               EXIT outer
-            ELSE IF(tracers(iq)%iqParent> 0) THEN
-               IF(tracers(tracers(iq)%iqParent)%name == tnom_trac) THEN
-                  nb_traciso(ixt,ip)     = nb_traciso(ixt,ip)+1
-                  iso_indnum(iq)         = indnum_fn_num(ixt)
-                  tracers(iq)%type       = 'tag'
-                  tracers(iq)%iso_iGroup = 1
-                  tracers(iq)%iso_iName  = ixt
-                  tracers(iq)%iso_iZone  = nb_traciso(ixt,ip)
-                  tracers(iq)%iso_iPhase = ip
-                  EXIT outer
-               END IF
-            END IF
-         END DO
-      END DO outer
-   END DO
-
-   niso = 0; nzone_prec = nb_traciso(1,1)
-   DO ixt = 1, niso_possibles
-      IF(nb_iso(ixt,1) == 0) CYCLE
-      IF(nb_iso(ixt,1) /= 1) CALL abort_gcm('infotrac_init', 'Isotopes are not well defined in traceur.def', 1)
-
-      ! on verifie que toutes les phases ont le meme nombre d'isotopes
-      IF(ANY(nb_iso(ixt,:) /= 1)) CALL abort_gcm('infotrac_init', 'Phases must have same number of isotopes', 1)
-
-      niso = niso+1
-      use_iso(ixt) = .TRUE.
-      nzone = nb_traciso(ixt,1)
-
-      ! on verifie que toutes les phases ont le meme nombre de traceurs d'isotopes
-      IF(ANY(nb_traciso(ixt,2:nqo) /= nzone)) CALL abort_gcm('infotrac_init','Phases must have same number of tracers',1)
-
-      ! on verifie que tous les isotopes ont le meme nombre de traceurs d'isotopes
-      IF(nzone /= nzone_prec) CALL abort_gcm('infotrac_init','Isotope tracers are not well defined in traceur.def',1)
-      nzone_prec = nzone
-   END DO
-
-   ! dimensions et flags isotopiques:
-   ntiso = niso*(nzone+1)
-   ok_isotopes = niso  > 0
-   ok_isotrac  = nzone > 0
- 
-   IF(ok_isotopes) THEN
-      ok_iso_verif = .FALSE.; CALL getin('ok_iso_verif', ok_iso_verif)
-      ok_init_iso  = .FALSE.; CALL getin('ok_init_iso',  ok_init_iso)
-   END IF
-      tnat        = [1.0, 155.76e-6, 2005.2e-6, 0.004/100., 0.0]
-      alpha_ideal = [1.0, 1.01,      1.006,     1.003,      1.0]
-!   END IF
-
-   ! remplissage du tableau iqiso(ntiso,phase)
-   ALLOCATE(iqiso(ntiso,nqo))   
-   iqiso(:,:)=0     
-   DO iq = 1, nqtot
-      IF(tracers(iq)%iso_iName <= 0) CYCLE
-      ixt = iso_indnum(iq) + tracers(iq)%iso_iZone*niso
-      iqiso(ixt, tracers(iq)%iso_iPhase) = iq
-   END DO
-
-   ! remplissage du tableau index_trac(nzone,niso)
-   ALLOCATE(index_trac(nzone, niso))  
-   IF(ok_isotrac) then
-      DO ii = 1, niso; index_trac(:, ii) = ii + niso*[(iz, iz=1, nzone)]; END DO
-   ELSE
-      index_trac(:,:)=0.0
-   END IF
-
+   modname = 'infotrac_isoinit'
+   tnat0       = [ 1.0 ,  155.76e-6, 2005.2e-6, 0.004/100., 0.0 ]    !--- Same length as tnom_iso
+   alpha_ideal0= [ 1.0 ,  1.01,      1.006,     1.003,      1.0 ]    !--- Same length as tnom_iso
    ALLOCATE(isotopes(1))                                             !--- Only water
-   nbIso = 1
+   nbIso = SIZE(isotopes)
    t => tracers
    i => isotopes(1)
    i%parent = 'H2O'
 
-   !--- Isotopes names list (embedded in the "keys" field)
-   i%niso  = niso
-   ALLOCATE(i%keys(i%niso))
-   mask = t%type=='tracer' .AND. delPhase(t%gen0Name)=='H2O' .AND. t%phase == 'g' .AND. t%iGeneration==1
-   str = strTail(PACK(delPhase(t%name), MASK=mask), '_')
-   CALL strReduce(str)
-   i%keys(:)%name = str
-
-   !--- Full isotopes list, with isotopes tagging tracers (if any) following the previous list
-   i%ntiso = ntiso
-   ALLOCATE(i%trac(i%ntiso))
-   mask = t%type=='tag'    .AND. delPhase(t%gen0Name)=='H2O' .AND. t%phase == 'g' .AND. t%iGeneration==2
+   !--- Effective isotopes names list (embedded in the "keys" field)
+   mask = t%type=='tracer' .AND. t%gen0Name==addPhase('H2O', 'g') .AND. t%iGeneration==1
    str = PACK(delPhase(t%name), MASK=mask)
    CALL strReduce(str)
+   i%niso = SIZE(str)
+   ALLOCATE(i%keys(i%niso))
+   i%keys(:)%name = str
+
+   !--- Check whether found isotopes are known
+   mask = [(ALL(tnom_iso /= str(ii)), ii=1, i%niso)]
+   IF(ANY(mask)) CALL abort_gcm(modname, 'The following isotopes are unknown: '//strStack(PACK(str, MASK=mask)), 1)
+
+   !--- Full isotopes list, with isotopes tagging tracers (if any) following the previous list
+   mask = t%type=='tag'    .AND. t%gen0Name==addPhase('H2O', 'g') .AND. t%iGeneration==2
+   str = PACK(delPhase(t%name), MASK=mask)
+   i%ntiso = i%niso + SIZE(str)
+   ALLOCATE(i%trac(i%ntiso))
    i%trac(:) = [i%keys(:)%name, str]
 
-   !--- Tagging zones names list
-   i%nzone = nzone
+   !--- Effective tagging zones names list
    i%zone = strTail(str, '_', .TRUE.)
+   CALL strReduce(i%zone)
+   i%nzone = SIZE(i%zone)
+   IF(i%ntiso /= i%niso*(i%nzone+1)) CALL abort_gcm(modname, 'Error in "ntiso" calculation', 1)
 
    !--- Effective phases list
-   i%nphas = nqo
    i%phase = ''
-   DO ip=1,nphases; IF(strIdx(t%name, addPhase('H2O',old_phases(ip:ip),''))/=0) i%phase=TRIM(i%phase)//known_phases(ip:ip); END DO
-
-   !--- Table: index in "qx" of an isotope, knowing its indices "it","ip" in "isotope%iName,%iPhase"
-   i%iTraPha = RESHAPE([((strIdx(t%name, TRIM(addPhase('H2O', new2oldPhase(i%phase(ip:ip)), ''))//'_'//TRIM(i%trac(it))), &
-                          it=1,i%ntiso), ip=1,i%nphas)], [i%ntiso,i%nphas])
+   DO ip=1,nphases; IF(strIdx(t%name, addPhase('H2O', ip))/=0) i%phase=TRIM(i%phase)//known_phases(ip:ip); END DO
+   i%nphas = LEN_TRIM(i%phase)
+
+   !--- Indexes related to isotopes
+   DO iq = 1, nqtot
+      t1 => tracers(iq)
+      t0 = t1%gen0Name
+      IF(t1%iGeneration==0 .OR. .NOT.t1%isAdvected .OR. delPhase(t0)/='H2O') CYCLE
+      t1%iso_iGroup = 1
+      t1%iso_iPhase = INDEX(i%phase, getPhase(t0))
+      t1%iso_iZone = strIdx(i%zone,  strTail(t1%name, '_'))
+      IF(t1%iso_iZone /= 0) t1%iso_iName = strIdx(i%keys(:)%name, delPhase(t1%parent))
+      IF(t1%iso_iZone == 0) t1%iso_iName = strIdx(i%keys(:)%name, delPhase(t1%name  ))
+   END DO
+
+   niso_possibles = SIZE(tnom_iso)
+!   ix = strIdx(tnom_iso, i%trac)
+!   tnat        = tnat0       (PACK(ix, MASK=ix/=0))
+!   alpha_ideal = alpha_ideal0(PACK(ix, MASK=ix/=0))
+   tnat        = tnat0
+   alpha_ideal = alpha_ideal0
+
+   !--- Tests
+   nb_iso  = [(COUNT(t%iso_iPhase == ip .AND. t%iGeneration == 1), ip=1, i%nphas)]
+   nb_tiso = [(COUNT(t%iso_iPhase == ip .AND. t%iGeneration == 2), ip=1, i%nphas)]
+   nb_zone = [(COUNT(t%iso_iZone  == iz), iz=1, i%nzone)]
+   IF(ANY(nb_iso (:) /= nb_iso (1))) CALL abort_gcm(modname, 'Phases must have same number of isotopes', 1)
+   IF(ANY(nb_tiso(:) /= nb_tiso(1))) CALL abort_gcm(modname, 'Phases must have same number of tagging tracers', 1)
+   IF(ANY(nb_zone(:) /= nb_zone(1))) CALL abort_gcm(modname, 'Isotopes must have the same number of tagging tracers', 1)
+
+   !--- Isotopic checking routines activation flag
+   i%check = .FALSE.; IF(i%niso > 0) CALL getin('ok_iso_verif', i%check)
+
+   !--- Table: index in "qx(:)" of an isotope, knowing its indices "it","ip" in "isotope%iName,%iPhase"
+   i%iqTraPha = RESHAPE([((strIdx(t%name, TRIM(addPhase(i%trac(it),ip,i%phase))),it=1,i%ntiso),ip=1,i%nphas)],[i%ntiso,i%nphas])
 
    !--- Table: index in "isotope%tracs(:)%name" of an isotopic tagging tracer, knowing its indices "iz","ip" in "isotope%iZone,%iName"
-   i%iZonIso = RESHAPE([((strIdx(i%trac,TRIM(i%trac(it))//'_'//TRIM(i%zone(iz))),iz=1,i%nzone),it=1,i%niso )],[i%nzone,i%niso ])
-   DO it=1,ntiso
-      WRITE(lunout,'(a,i0,a)')TRIM('infotrac_init')//': iqiso  (',it,',:) = '//strStack(int2str(iqiso(it,:)))
-      WRITE(lunout,'(a,i0,a)')TRIM('infotrac_init')//': iTraPha(',it,',:) = '//strStack(int2str(i%iTraPha(it,:)))
-   END DO
-   DO iz=1,nzone
-      WRITE(lunout,'(a,i0,a)')TRIM('infotrac_init')//': index_trac(',iz,',:) = '//strStack(int2str(index_trac(iz,:)))
-      WRITE(lunout,'(a,i0,a)')TRIM('infotrac_init')//': iZonIso   (',iz,',:) = '//strStack(int2str(i%iZonIso(iz,:)))
-   END DO
-
-   ! Finalize :
+   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])
+
+   DO it=1,i%ntiso; CALL msg('iqTraPha('//TRIM(int2str(it))//',:) = '//strStack(int2str(i%iqTraPha(it,:))), modname); END DO
+   DO iz=1,i%nzone; CALL msg('itZonIso('//TRIM(int2str(iz))//',:) = '//strStack(int2str(i%itZonIso(iz,:))), modname); END DO
+
+   !--- Isotopic quantities (to be removed soon)
+   ok_isotopes   = i%niso  > 0
+   ok_isotrac    = i%nzone > 0
+   ok_iso_verif  = i%check
+   niso_possibles = SIZE(tnom_iso)
+   iso_num       = [(strIdx(tnom_iso(:),    strHead(delPhase(tracers(iq)%name), '_')), iq=1, nqtot)]
+   iso_indnum    = [(strIdx(i%keys(:)%name, strHead(delPhase(tracers(iq)%name), '_')), iq=1, nqtot)]
+   indnum_fn_num = [(strIdx(i%keys(:)%name, tnom_iso(ixt)), ixt=1, niso_possibles)]
+   use_iso       = indnum_fn_num /= 0            !--- .TRUE. for the effectively used isotopes of the possible isotopes list
+
+   !--- Finalize :
    DEALLOCATE(nb_iso)
 
@@ -903,17 +836,17 @@
    lv = .FALSE.; IF(PRESENT(lVerbose)) lv = lVerbose
    lerr = .FALSE.
-   IF(iIso == ixIso) RETURN                                      !--- Nothing to do if the index is already OK
+   IF(iIso == ixIso) RETURN                                          !--- Nothing to do if the index is already OK
    lerr = iIso<=0 .OR. iIso>nbIso
    CALL msg('Inconsistent isotopes family index '//TRIM(int2str(iIso))//': should be > 0 and <= '//TRIM(int2str(nbIso))//'"',&
             ll=lerr .AND. lV)
    IF(lerr) RETURN
-   ixIso = iIso                                                  !--- Update currently selected family index
-   isotope => isotopes(ixIso)                                    !--- Select corresponding component
-   isoKeys => isotope%keys;    niso     = isotope%niso
-   isoName => isotope%trac;    ntiso    = isotope%ntiso
-   isoZone => isotope%zone;    nzone    = isotope%nzone
-   isoPhas => isotope%phase;   nphas    = isotope%nphas
-   iZonIso => isotope%iZonIso; isoCheck = isotope%check
-   iTraPha => isotope%iTraPha
+   ixIso = iIso                                                      !--- Update currently selected family index
+   isotope  => isotopes(ixIso)                                       !--- Select corresponding component
+   isoKeys  => isotope%keys;     niso     => isotope%niso
+   isoName  => isotope%trac;     ntiso    => isotope%ntiso
+   isoZone  => isotope%zone;     nzone    => isotope%nzone
+   isoPhas  => isotope%phase;    nphas    => isotope%nphas
+   itZonIso => isotope%itZonIso; isoCheck => isotope%check
+   iqTraPha => isotope%iqTraPha
 END FUNCTION isoSelectByIndex
 !==============================================================================================================================
Index: LMDZ6/trunk/libf/dyn3dmem/dynetat0_loc.F90
===================================================================
--- LMDZ6/trunk/libf/dyn3dmem/dynetat0_loc.F90	(revision 4119)
+++ LMDZ6/trunk/libf/dyn3dmem/dynetat0_loc.F90	(revision 4120)
@@ -7,19 +7,18 @@
 !-------------------------------------------------------------------------------
   USE parallel_lmdz
-  USE strings_mod, ONLY: maxlen
-  USE infotrac,    ONLY: nqtot, tracers, iqiso, iso_indnum, tnat, alpha_ideal, ok_isotopes
-  USE netcdf, ONLY: NF90_OPEN,  NF90_INQUIRE_DIMENSION, NF90_INQ_VARID,        &
-      NF90_NOWRITE, NF90_CLOSE, NF90_INQUIRE_VARIABLE,  NF90_GET_VAR, NF90_NoErr
+  USE infotrac,    ONLY: nqtot, tracers, niso, iqiso, iso_indnum, iso_num, tnat, alpha_ideal, ok_isotopes, iH2O
+  USE strings_mod, ONLY: maxlen, msg, strStack, real2str
+  USE netcdf,      ONLY: NF90_OPEN,  NF90_NOWRITE, NF90_INQUIRE_DIMENSION, NF90_INQ_VARID, &
+                         NF90_CLOSE, NF90_GET_VAR, NF90_INQUIRE_VARIABLE,  NF90_NoErr
+  USE readTracFiles_mod, ONLY: new2oldName
   USE control_mod, ONLY: planet_type
   USE assert_eq_m, ONLY: assert_eq
   USE comvert_mod, ONLY: pa,preff
-  USE comconst_mod, ONLY: cpp, daysec, dtvr, g, im, jm, kappa, lllm, &
-                          omeg, rad
+  USE comconst_mod, ONLY: cpp, daysec, dtvr, g, im, jm, kappa, lllm, omeg, rad
   USE logic_mod, ONLY: fxyhypb, ysinus
   USE serre_mod, ONLY: clon, clat, grossismx, grossismy
-  USE temps_mod, ONLY: annee_ref,day_ref,itau_dyn, &
-                       start_time,day_ini,hour_ini
+  USE temps_mod, ONLY: annee_ref, day_ini, day_ref, itau_dyn, start_time
   USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
-  
+
   IMPLICIT NONE
   include "dimensions.h"
@@ -40,5 +39,5 @@
 !===============================================================================
 ! Local variables:
-  CHARACTER(LEN=maxlen) :: msg, var, modname
+  CHARACTER(LEN=maxlen) :: mesg, var, modname, oldVar
   INTEGER, PARAMETER :: length=100
   INTEGER :: iq, fID, vID, idecal, ierr, iqParent, iName, iZone, iPhase
@@ -58,9 +57,9 @@
 !!!     .... while keeping everything OK for LMDZ EARTH
   IF(planet_type=="generic") THEN
-    WRITE(lunout,*)'NOTE NOTE NOTE : Planeto-like start files'
+    CALL msg('NOTE NOTE NOTE : Planeto-like start files', modname)
     idecal = 4
     annee_ref  = 2000
   ELSE
-    WRITE(lunout,*)'NOTE NOTE NOTE : Earth-like start files'
+    CALL msg('NOTE NOTE NOTE : Earth-like start files', modname)
     idecal = 5
     annee_ref  = tab_cntrl(5)
@@ -106,5 +105,5 @@
 
 !-------------------------------------------------------------------------------
-  WRITE(lunout,*)TRIM(modname)//': rad,omeg,g,cpp,kappa',rad,omeg,g,cpp,kappa
+  CALL msg('rad, omeg, g, cpp, kappa = '//TRIM(strStack(real2str([rad,omeg,g,cpp,kappa]))), modname)
   CALL check_dim(im,iim,'im','im')
   CALL check_dim(jm,jjm,'jm','jm')
@@ -120,6 +119,6 @@
   var="temps"
   IF(NF90_INQ_VARID(fID,var,vID)/=NF90_NoErr) THEN
-    WRITE(lunout,*)TRIM(modname)//": missing field <temps>"
-    WRITE(lunout,*)TRIM(modname)//": trying with <Time>"; var="Time"
+    CALL msg('missing field <temps> ; trying with <Time>', modname)
+    var="Time"
     CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
   END IF
@@ -153,33 +152,38 @@
   ALLOCATE(q_glo(ip1jmp1,llm))
   DO iq=1,nqtot
-    var=TRIM(tracers(iq)%name)
+    var = tracers(iq)%name
+    oldVar = new2oldName(var)
+    !--------------------------------------------------------------------------------------------------------------------------
+    IF(NF90_INQ_VARID(fID, var, vID) == NF90_NoErr) THEN                                 !=== REGULAR CASE
+      CALL get_var2(var,q_glo); q(ijb_u:ije_u,:,iq)=q_glo(ijb_u:ije_u,:)
+    !--------------------------------------------------------------------------------------------------------------------------
+    ELSE IF(NF90_INQ_VARID(fID, oldVar, vID) == NF90_NoErr) THEN                         !=== OLD NAME
+      CALL msg('Tracer <'//TRIM(var)//'> is missing => initialized to <'//TRIM(oldVar)//'>', modname)
+      CALL get_var2(oldVar, q_glo); q(ijb_u:ije_u,:,iq)=q_glo(ijb_u:ije_u,:)
+    !--------------------------------------------------------------------------------------------------------------------------
 #ifdef INCA
-    IF (var .eq. "O3" ) THEN 
-       IF(NF90_INQ_VARID(fID,var,vID) == NF90_NoErr) THEN
-          CALL get_var2(var,q_glo); q(ijb_u:ije_u,:,iq)=q_glo(ijb_u:ije_u,:); CYCLE
-       ELSE
-          WRITE(lunout,*) 'Tracer O3 is missing - it is initialized to OX'
-          IF(NF90_INQ_VARID(fID,"OX",vID) == NF90_NoErr) THEN
-             CALL get_var2("OX",q_glo); q(ijb_u:ije_u,:,iq)=q_glo(ijb_u:ije_u,:); CYCLE
-          ENDIF
-       ENDIF
-    ENDIF
+    ELSE IF(NF90_INQ_VARID(fID, 'OX',   vID) == NF90_NoErr .AND. var == 'O3') THEN       !=== INCA: OX INSTEAD OF O3
+      CALL msg('Tracer <O3> is missing => initialized to <OX>', modname)
+      CALL get_var2( 'OX' , q_glo); q(ijb_u:ije_u,:,iq)=q_glo(ijb_u:ije_u,:)
+    !--------------------------------------------------------------------------------------------------------------------------
 #endif
-    IF(NF90_INQ_VARID(fID,var,vID)==NF90_NoErr) THEN
-      CALL get_var2(var,q_glo); q(ijb_u:ije_u,:,iq)=q_glo(ijb_u:ije_u,:); CYCLE
+    ELSE IF(tracers(iq)%iso_iGroup == iH2O .AND. niso > 0) THEN                          !=== WATER ISOTOPES
+!     iName    = tracers(iq)%iso_iName  ! (next commit)
+      iName    = iso_num(iq)
+      iPhase   = tracers(iq)%iso_iPhase
+      iqParent = tracers(iq)%iqParent
+      IF(tracers(iq)%iso_iZone == 0) THEN
+         CALL msg('Tracer <'//TRIM(var)//'> is missing => initialized with a simplified Rayleigh distillation law.', modname)
+         q(ijb_u:ije_u,:,iq)= q(ijb_u:ije_u,:,iqParent)*tnat(iName)*(q(ijb_u:ije_u,:,iqParent)/30.e-3)**(alpha_ideal(iName)-1.)
+      ELSE
+         CALL msg('Tracer <'//TRIM(var)//'> is missing => initialized to its parent isotope concentration.', modname)
+         q(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqiso(iso_indnum(iq),iPhase))
+      END IF
+    !--------------------------------------------------------------------------------------------------------------------------
+    ELSE                                                                                 !=== MISSING: SET TO 0
+      CALL msg('Tracer <'//TRIM(var)//'> is missing => initialized to zero', modname)
+      q(ijb_u:ije_u,:,iq)=0.
+    !--------------------------------------------------------------------------------------------------------------------------
     END IF
-    WRITE(lunout,*)TRIM(modname)//": Tracer <"//TRIM(var)//"> is missing"
-    WRITE(lunout,*)"         It is hence initialized to zero"
-    q(ijb_u:ije_u,:,iq)=0.
-   !--- CRisi: for isotops, theoretical initialization using very simplified
-   !           Rayleigh distillation las.
-    iName = tracers(iq)%iso_iName
-    IF(.NOT.ok_isotopes .OR. iName <= 0) CYCLE
-    iZone = tracers(iq)%iso_iZone
-    iPhase= tracers(iq)%iso_iPhase
-    iqParent = tracers(iq)%iqParent
-    IF(iZone==0) q(:,:,iq) = q(:,:,iqParent)*tnat(iName)    &
-     &         *(q(:,:,iqParent)/30.e-3)**(alpha_ideal(iName)-1.)
-    IF(iZone==1) q(:,:,iq) = q(:,:,iqiso(iso_indnum(iq),iPhase))
   END DO
   DEALLOCATE(q_glo)
@@ -199,6 +203,6 @@
     s1='value of '//TRIM(str1)//' ='
     s2=' read in starting file differs from parametrized '//TRIM(str2)//' ='
-    WRITE(msg,'(10x,a,i4,2x,a,i4)'),s1,n1,s2,n2
-    CALL ABORT_gcm(TRIM(modname),TRIM(msg),1)
+    WRITE(mesg,'(10x,a,i4,2x,a,i4)')TRIM(ADJUSTL(s1)),n1,TRIM(ADJUSTL(s2)),n2
+    CALL ABORT_gcm(TRIM(modname),TRIM(mesg),1)
   END IF
 END SUBROUTINE check_dim
@@ -263,10 +267,10 @@
   IF(ierr==NF90_NoERR) RETURN
   SELECT CASE(typ)
-    CASE('inq');   msg="Field <"//TRIM(nam)//"> is missing"
-    CASE('get');   msg="Reading failed for <"//TRIM(nam)//">"
-    CASE('open');  msg="File opening failed for <"//TRIM(nam)//">"
-    CASE('close'); msg="File closing failed for <"//TRIM(nam)//">"
+    CASE('inq');   mesg="Field <"//TRIM(nam)//"> is missing"
+    CASE('get');   mesg="Reading failed for <"//TRIM(nam)//">"
+    CASE('open');  mesg="File opening failed for <"//TRIM(nam)//">"
+    CASE('close'); mesg="File closing failed for <"//TRIM(nam)//">"
   END SELECT
-  CALL ABORT_gcm(TRIM(modname),TRIM(msg),ierr)
+  CALL ABORT_gcm(TRIM(modname),TRIM(mesg),ierr)
 END SUBROUTINE err
 
Index: LMDZ6/trunk/libf/dyn3dmem/iniacademic_loc.F90
===================================================================
--- LMDZ6/trunk/libf/dyn3dmem/iniacademic_loc.F90	(revision 4119)
+++ LMDZ6/trunk/libf/dyn3dmem/iniacademic_loc.F90	(revision 4120)
@@ -6,5 +6,5 @@
   USE filtreg_mod, ONLY: inifilr
   USE infotrac,    ONLY: nqtot, niso_possibles, ok_isotopes, ok_iso_verif, tnat, alpha_ideal, &
-                         iqiso, tracers, iso_indnum
+                         iqiso, tracers, iso_indnum, iso_num
   USE control_mod, ONLY: day_step,planet_type
   use exner_hyb_m, only: exner_hyb
@@ -23,4 +23,5 @@
   USE temps_mod, ONLY: annee_ref, day_ini, day_ref
   USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
+  USE readTracFiles_mod, ONLY: addPhase
 
   !   Author:    Frederic Hourdin      original: 15/01/93
@@ -66,5 +67,5 @@
   real tetastrat ! potential temperature in the stratosphere, in K
   real tetajl(jjp1,llm)
-  INTEGER i,j,l,lsup,ij, iq, iName, iZone, iPhase, iqParent
+  INTEGER i,j,l,lsup,ij, iq, iName, iPhase, iqParent
 
   REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
@@ -280,19 +281,20 @@
            do iq=1,nqtot
               q(ijb_u:ije_u,:,iq)=0.
-!              IF(tracers(iq)%name == 'H2O'//phases_sep//'g') q(ijb_u:ije_u,:,iq)=1.e-10
-!              IF(tracers(iq)%name == 'H2O'//phases_sep//'l') q(ijb_u:ije_u,:,iq)=1.e-15
-              IF(tracers(iq)%name == 'H2Ov') q(ijb_u:ije_u,:,iq)=1.e-10
-              IF(tracers(iq)%name == 'H2Ol') q(ijb_u:ije_u,:,iq)=1.e-15
+              IF(tracers(iq)%name == addPhase('H2O', 'g')) q(ijb_u:ije_u,:,iq)=1.e-10
+              IF(tracers(iq)%name == addPhase('H2O', 'l')) q(ijb_u:ije_u,:,iq)=1.e-15
 
               ! CRisi: init des isotopes
               ! distill de Rayleigh très simplifiée
-              iName = tracers(iq)%iso_iName
+!             iName    = tracers(iq)%iso_iName  ! (next commit)
+              iName    = iso_num(iq)
               if (.NOT.ok_isotopes .OR. iName <= 0) CYCLE
-              iZone    = tracers(iq)%iso_iZone
               iPhase   = tracers(iq)%iso_iPhase
               iqParent = tracers(iq)%iqParent
-              if (iZone == 0) q(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqParent)*tnat(iName) &
-                                                  *(q(ijb_u:ije_u,:,iqParent)/30.e-3)**(alpha_ideal(iName)-1)
-              if (iZone == 1) q(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqiso(iso_indnum(iq),iPhase))
+              IF(tracers(iq)%iso_iZone == 0) THEN
+                 q(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqParent)*tnat(iName) &
+                                     *(q(ijb_u:ije_u,:,iqParent)/30.e-3)**(alpha_ideal(iName)-1.)
+              ELSE
+                 q(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqiso(iso_indnum(iq),iPhase))
+              END IF
            enddo
         else
Index: LMDZ6/trunk/libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90
===================================================================
--- LMDZ6/trunk/libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90	(revision 4119)
+++ LMDZ6/trunk/libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90	(revision 4120)
@@ -16,11 +16,5 @@
   USE mod_phys_lmdz_para, ONLY: klon_omp ! number of columns (on local omp grid)
   USE vertical_layers_mod, ONLY : init_vertical_layers
-  USE infotrac, ONLY: nqtot,nqo,nbtr,nqCO2,tracers,type_trac,&
-                      conv_flg,pbl_flg,solsym,&
-                      ok_isotopes,ok_iso_verif,ok_isotrac,&
-                      ok_init_iso,niso_possibles,tnat,&
-                      alpha_ideal,use_iso,iqiso,iso_indnum,&
-                      indnum_fn_num,index_trac,&
-                      niso,ntraceurs_zone,ntraciso,nqtottr
+  USE infotrac, ONLY: nbtr,nqCO2,tracers,isotopes,type_trac,conv_flg,pbl_flg,solsym,nqtottr
 #ifdef CPP_StratAer
   USE infotrac_phy, ONLY: nbtr_bin, nbtr_sulgas, id_OCS_strat, &
@@ -143,11 +137,5 @@
 
   ! Initialize tracer names, numbers, etc. for physics
-  CALL init_infotrac_phy(nqtot,nqo,nbtr,nqtottr,nqCO2,tracers,type_trac,&
-                         conv_flg,pbl_flg,solsym,&
-                         ok_isotopes,ok_iso_verif,ok_isotrac,&
-                         ok_init_iso,niso_possibles,tnat,&
-                         alpha_ideal,use_iso,iqiso,iso_indnum,&
-                         indnum_fn_num,index_trac,&
-                         niso,ntraceurs_zone,ntraciso)
+  CALL init_infotrac_phy(type_trac, tracers, isotopes, nqtottr, nqCO2, pbl_flg, conv_flg, solsym)
 
   ! Initializations for Reprobus
Index: LMDZ6/trunk/libf/misc/readTracFiles_mod.f90
===================================================================
--- LMDZ6/trunk/libf/misc/readTracFiles_mod.f90	(revision 4119)
+++ LMDZ6/trunk/libf/misc/readTracFiles_mod.f90	(revision 4120)
@@ -1,6 +1,6 @@
 MODULE readTracFiles_mod
 
-  USE strings_mod,    ONLY: msg, testFile,  strFind, strStack, strReduce,  strHead, strCount,   find, maxlen, fmsg, &
-             removeComment, cat, checkList, strIdx,  strParse, strReplace, strTail, reduceExpr, test, get_in, dispTable
+  USE strings_mod,    ONLY: msg, testFile,  strFind, strStack, strReduce,  strHead, strCount, find, fmsg, reduceExpr, &
+             removeComment, cat, checkList, str2int, strParse, strReplace, strTail, strIdx, maxlen, test, dispTable, get_in
   USE trac_types_mod, ONLY: trac_type, isot_type, keys_type
 
@@ -12,12 +12,13 @@
   PUBLIC :: readTracersFiles, indexUpdate, setGeneration             !--- TOOLS ASSOCIATED TO TRACERS  DESCRIPTORS
   PUBLIC :: readIsotopesFile                                         !--- TOOLS ASSOCIATED TO ISOTOPES DESCRIPTORS
-  PUBLIC :: getKey_init, getKey, setDirectKeys                       !--- GET/SET KEYS FROM/TO tracers & isotopes
-
-  PUBLIC :: known_phases, old_phases, nphases, phases_names, &       !--- VARIABLES RELATED TO THE PHASES
-            phases_sep, delPhase, addPhase, &                        !--- + ROUTINES TO ADD/REMOVE PHASE TO/FROM A NAME
-            old2newPhase, new2oldPhase
+  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
@@ -30,8 +31,11 @@
   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_1,    addPhase_m; END INTERFACE    addPhase
+  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
 !------------------------------------------------------------------------------------------------------------------------------
 
@@ -49,4 +53,10 @@
   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)
@@ -87,6 +97,6 @@
   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, oldH2O, newH2O
+  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(:)
@@ -142,8 +152,8 @@
         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 = TRIM(s(3))                                  !--- Set %name:   name of the tracer
-        tracs(it)%parent = tran0                                     !--- Set %parent: transporting fluid
-        IF(ns == 4) tracs(it)%parent = s(4)                          !---     default: 'air' or defined in the file
-        tracs(it)%phase = known_phases(1:1)                          !--- Set %phase:  tracer phase (default: "g"azeous)
+        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
@@ -151,15 +161,6 @@
       END DO
       CLOSE(90)
-      DO ip = 1, nphases                                             !--- Deal with old water names
-        oldH2O = 'H2O'//old_phases(ip:ip)
-        newH2O = 'H2O'//phases_sep//known_phases(ip:ip)
-        ix = strIdx(tracs(:)%name, oldH2O)
-        IF(ix == 0) CYCLE
-        tracs(ix)%name  = newH2O                                     !--- Set %name:   name of the tracer
-        WHERE(tracs(:)%parent == oldH2O) tracs(:)%parent = newH2O    !--- Set %parent: transporting fluid
-        tracs(ix)%phase = known_phases(ip:ip)                        !--- Set %phase:  tracer phase
-      END DO
       CALL setGeneration(tracs)                                      !--- Set %iGeneration and %gen0Name
-      WHERE(tracs%iGeneration == 3) tracs%type = 'tag'               !--- Set %type:        'tracer' or 'tag'
+      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
@@ -167,28 +168,29 @@
       tracs(:)%keys%name = tracs(:)%name                             !--- Copy tracers names in keys components
     !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-    CASE(2); IF(test(feedDBase(["tracer.def"],[type_trac]), lerr)) RETURN  !=== SINGLE FILE, COMA-SEPARATED SECTIONS LIST
+    CASE(2); IF(test(feedDBase(["tracer.def"], [type_trac], modname), lerr)) RETURN !=== SINGLE   FILE, MULTIPLE SECTIONS
     !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-    CASE(3); IF(test(feedDBase(  trac_files  , sections  ), lerr)) RETURN  !=== MULTIPLE FILES, ONE SECTION EACH FILE
+    CASE(3); IF(test(feedDBase(  trac_files  ,  sections,   modname), lerr)) RETURN !=== MULTIPLE FILES, SINGLE  SECTION
   !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   END SELECT
   !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
-  IF(ANY([2,3] == fType) .AND. nsec > 1) THEN
-    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
+  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) RESULT(lerr)
+LOGICAL FUNCTION feedDBase(fnames, snames, modname) RESULT(lerr)
 ! Purpose: Read the sections "snames(is)" (coma-separated list) from each "fnames(is)"
 !   file and create the corresponding tracers set descriptors in the database "dBase":
@@ -200,12 +202,12 @@
   CHARACTER(LEN=*), INTENT(IN)  :: fnames(:)                         !--- Files names
   CHARACTER(LEN=*), INTENT(IN)  :: snames(:)                         !--- Coma-deparated list of sections (one list each file)
-  INTEGER,  ALLOCATABLE :: ndb(:)                                    !--- Nuber of sections for each file
+  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, modname
+  CHARACTER(LEN=maxlen) :: fnm, snm
   INTEGER               :: idb, i
   LOGICAL :: ll
 !------------------------------------------------------------------------------------------------------------------------------
-  modname = 'feedDBase'
   !=== READ THE REQUIRED SECTIONS
   ll = strCount(snames, ',', ndb)                                    !--- Number of sections for each file
@@ -219,4 +221,5 @@
   !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     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
@@ -225,11 +228,8 @@
     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
   !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-  !=== DISPLAY BASIC INFORMATION
-  lerr = ANY([( dispTraSection('Expanded list for section "'//TRIM(dBase(idb)%name)//'"', dBase(idb)%name, modname), &
-                idb=1, SIZE(dBase) )])
 END FUNCTION feedDBase
 !------------------------------------------------------------------------------------------------------------------------------
@@ -406,14 +406,16 @@
   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)%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 = pa(ipr)
-        ttr(i)%keys = keys_type(ta(itr), tr(it)%keys%key, tr(it)%keys%val)
+        ttr(i)%name   = TRIM(ta(itr))
+        ttr(i)%parent = TRIM(pa(ipr))
+        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)%type      = tr(it)%type                     !--- Duplicating type
+    ttr(iq:iq+ntr*npr-1)%component = tr(it)%component                !--- Duplicating type
     iq = iq + ntr*npr
   !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -440,6 +442,6 @@
   tr(:)%iGeneration = -1                                             !--- error if -1
   nq = SIZE(tr, DIM=1)                                               !--- Number of tracers lines
-  lg = tr(:)%parent == tran0                                         !--- First generation tracers flag
-  WHERE(lg) tr(:)%iGeneration = 0                                    !--- First generation tracers
+  lg = tr(:)%parent == tran0                                         !--- Flag for generation 0 tracers
+  WHERE(lg) tr(:)%iGeneration = 0                                    !--- Generation 0 tracers
 
   !=== Determine generation for each tracer
@@ -511,6 +513,6 @@
     ll = tr(:)%name==TRIM(tnam)                                      !--- Mask for current tracer name
     IF(COUNT(ll)==1 ) CYCLE                                          !--- Tracer is not repeated
-    IF(tr(iq)%iGeneration>1) THEN
-      tdup(iq) = tnam                                                !--- gen>1: MUST be unique
+    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
@@ -531,5 +533,5 @@
 SUBROUTINE expandPhases(tr)
 !------------------------------------------------------------------------------------------------------------------------------
-! Purpose: Expand the phases in the tracers descriptor "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
@@ -538,4 +540,5 @@
   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
@@ -544,38 +547,39 @@
   nt = 0
   DO iq = 1, nq                                                      !--- GET THE NUMBER OF TRACERS
-    IF(tr(iq)%iGeneration /= 1) CYCLE
-    nc = COUNT(tr(:)%gen0Name==tr(iq)%name .AND. tr%iGeneration/=1)  !--- Number of childs of tr(iq)
+    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))
+  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                                                       !--- Need of a phase suffix
-    IF(lTg) lEx=lEx.AND.tr(iq)%iGeneration>1                         !--- No phase suffix for first generation tags
-    DO i=1,n                                                         !=== LOOP ON FIRST GENERATION ANCESTORS
-      jq=i0(i)                                                       !--- tr(jq): ith copy of 1st gen. ancestor of tr(iq)
-      IF(tr(iq)%iGeneration==1) jq=iq                                !--- Generation 1: current tracer phases only
+    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
-        trn=TRIM(tr(iq)%name); nam=trn                               !--- Tracer name (regular case)
+      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 = TRIM(nam)//phases_sep//pha(ip:ip)              !--- Phase extension needed
+        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)%name      = TRIM(nam)                                !--- Name with possibly phase suffix
         ttr(it)%keys%name = TRIM(nam)                                !--- Name inside the keys decriptor
-        ttr(it)%phase = pha(ip:ip)                                   !--- Single phase entry
-        IF(lEx.AND.tr(iq)%iGeneration>1) THEN
-          ttr(it)%parent   = TRIM(ttr(it)%parent)//phases_sep//pha(ip:ip)
-          ttr(it)%gen0Name = TRIM(ttr(it)%gen0Name)//phases_sep//pha(ip:ip)
+        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
+        it = it+1
       END DO
-      IF(tr(iq)%iGeneration==1) EXIT                                 !--- Break phase loop for gen 1
+      IF(tr(iq)%iGeneration==0) EXIT                                 !--- Break phase loop for gen 0
     END DO
   END DO
@@ -590,22 +594,23 @@
 !------------------------------------------------------------------------------------------------------------------------------
 ! Purpose: Sort tracers:
-!  * Put water at first places, in the "known_phases" order.
+!  * 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
+  TYPE(trac_type), ALLOCATABLE, INTENT(INOUT) :: tr(:)         !--- Tracer derived type vector
+!  TYPE(trac_type), ALLOCATABLE :: ttr(:)
+  INTEGER,         ALLOCATABLE :: iy(:), iz(:)
   INTEGER :: ig, ng, iq, jq, ip, nq, n, ix(SIZE(tr)), k
-  INTEGER, ALLOCATABLE :: iy(:), iz(:)
 !------------------------------------------------------------------------------------------------------------------------------
   nq = SIZE(tr)
-  iy = [(k, k=1, nq)]
   DO ip = nphases, 1, -1
-    iq = strIdx(tracers(:)%name, 'H2O'//phases_sep//known_phases(ip:ip))
-    IF(iq/=0) iy = [iq, iy(1:iq-1), iy(iq:nq)]
-  END DO
-  tr = tr(iy)                                                        !--- Water displaces at first positions
-  iq = 1
+    iq = strIdx(tr(:)%name, addPhase('H2O', ip))
+    IF(iq == 0) CYCLE
+    tr = tr([iq, 1:iq-1, iq+1:nq])
+!    tr(:)%name = nam
+  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
@@ -616,11 +621,12 @@
     END DO
   ELSE
-    DO jq = 1, nq                                                    !--- Loop on first generation tracers
-      IF(tr(jq)%iGeneration /= 1) CYCLE                              !--- Skip generations >= 1
-      ix(iq) = jq                                                    !--- First generation ancestor index first
-      iq = iq + 1
+    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)            !--- Generations number of the "tr(jq)" family
-      DO ig = 2, ng                                                  !--- Loop   on generations for the tr(jq) family
+      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(:)"
@@ -724,5 +730,5 @@
       tnam = TRIM(t1(iq)%name)                                       !--- Original name
       IF(COUNT(t1%name == tnam) == 1) CYCLE                          !--- Current tracer is not duplicated: finished
-      tnam_new = TRIM(tnam)//phases_sep//TRIM(sections(is)%name)     !--- Same with section extension
+      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
@@ -757,4 +763,5 @@
   INTEGER :: idb, iq, nq
   INTEGER, ALLOCATABLE :: hadv(:), vadv(:)
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: phas(:)
   TYPE(trac_type), POINTER :: tm(:)
   lerr = .FALSE.
@@ -762,9 +769,17 @@
   tm => dBase(idb)%trac
   nq = SIZE(tm)
-  IF(test(getKeyByName_im('hadv', hadv, tm(:)%name, tm(:)%keys),lerr)) RETURN
-  IF(test(getKeyByName_im('vadv', vadv, tm(:)%name, tm(:)%keys),lerr)) RETURN
+  !--- 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(test(dispTable('iiissis', ['iq        ','hadv      ','vadv      ','short name','parent    ','igen      ','phase     '], &
-    cat(tm(:)%name,  tm(:)%parent, tm(:)%phase), cat([(iq, iq=1, nq)],  hadv,  vadv, tm(:)%iGeneration)), lerr)) RETURN
+  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
 !==============================================================================================================================
@@ -825,15 +840,18 @@
 SUBROUTINE indexUpdate(tr)
   TYPE(trac_type), INTENT(INOUT) :: tr(:)
-  INTEGER :: iq, ig, ng, ngen
+  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)
-    ng = tr(iq)%iGeneration                                          !--- Generation of the current tracer
-    ix = idxAncestor(tr, igen = ng); ix = PACK(ix, ix/=0)            !--- Indexes of the tracers with ancestor tr(iq)
-    !--- Childs indexes in growing generation order
-    tr(iq)%iqDescen = [( PACK(ix, MASK = tr(ix)%iGeneration == ig), ig = ng+1, ngen)]
-    tr(iq)%nqDescen =     SUM(  [( COUNT(tr(ix)%iGeneration == ig), ig = ng+1, ngen)] )
-    tr(iq)%nqChilds =              COUNT(tr(ix)%iGeneration == ng+1)
+    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
@@ -847,5 +865,5 @@
 !=== NOTES:                                                                                                                ====
 !===  * Most of the "isot" components have been defined in the calling routine (initIsotopes):                             ====
-!===      parent,  nzone, zone(:),  niso, keys(:)%name,  ntiso, trac(:),  nphas, phas,  iTraPha(:,:),  iZonPhi(:,:)        ====
+!===      parent,  nzone, zone(:),  niso, keys(:)%name,  ntiso, trac(:),  nphas, phas,  iqTraPha(:,:),  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             ====
@@ -909,4 +927,8 @@
     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)
 
@@ -930,9 +952,9 @@
   LOGICAL, ALLOCATABLE :: ll(:)                                      !--- Mask
   TYPE(trac_type), POINTER   ::  t(:), t1
-  TYPE(isot_type), POINTER   ::  s
+  TYPE(isot_type), POINTER   ::  i
 
   t => trac
 
-  p = PACK(delPhase(t%parent), MASK = t%type=='tracer' .AND. t%iGeneration==2) !--- Parents of 2nd generation isotopes
+  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))
@@ -943,19 +965,19 @@
   isot(:)%parent = p
   DO ic = 1, SIZE(p)                                                 !--- Loop on isotopes classes
-    s => isot(ic)
-    iname = s%parent                                                 !--- Current isotopes class name (parent tracer name)
+    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"
-    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)
+    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 == 3
-    s%zone = PACK(strTail(t(:)%name,'_',lFirst=.TRUE.), MASK = ll)   !--- Tagging zones names  for isotopes category "iname"
-    CALL strReduce(s%zone)
-    s%nzone = SIZE(s%zone)                                           !--- Tagging zones number for isotopes category "iname"
+    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")
@@ -963,33 +985,32 @@
     str = PACK(delPhase(t(:)%name), MASK=ll)
     CALL strReduce(str)
-    s%ntiso = s%niso + SIZE(str)                                     !--- Number of isotopes + their geographic tracers [ntraciso]
-    ALLOCATE(s%trac(s%ntiso))
-    FORALL(it = 1:s%niso) s%trac(it) = s%keys(it)%name
-    FORALL(it = s%niso+1:s%ntiso) s%trac(it) = str(it-s%niso)
+    i%ntiso = i%niso + SIZE(str)                                     !--- Number of isotopes + their geographic tracers [ntraciso]
+    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"
-    s%phase = ''
-    DO ip = 1, nphases; ph = known_phases(ip:ip); IF(strIdx(t%name,addPhase(iname, ph)) /= 0) s%phase = TRIM(s%phase)//ph; END DO
-    s%nphas = LEN_TRIM(s%phase)                                       !--- Equal to "nqo" for water
+    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) CYCLE                       !--- Only deal with tracers descending on "iname"
+      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(s%trac, delPhase(strHead(t1%name,'_'))) !--- Current isotope       idx in effective isotopes list
-      t1%iso_iZone  = strIdx(s%zone,          strTail(t1%name,'_') ) !--- Current isotope zone  idx in effective zones    list
-      t1%iso_iPhase =  INDEX(s%phase,TRIM(t1%phase))                 !--- Current isotope phase idx in effective phases   list
-      IF(t1%iGeneration /= 3) t1%iso_iZone = 0                       !--- Skip possible generation 2 tagging tracers
+      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
-    s%iTraPha = RESHAPE( [( (strIdx(t(:)%name,  addPhase(s%trac(it),s%phase(ip:ip))),    it=1, s%ntiso), ip=1, s%nphas)], &
-                         [s%ntiso, s%nphas] )
-
+    i%iqTraPha = 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
-    s%iZonIso = RESHAPE( [( (strIdx(s%trac(:), TRIM(s%trac(it))//'_'//TRIM(s%zone(iz))), iz=1, s%nzone), it=1, s%niso )], &
-                         [s%nzone, s%niso] )
+    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
 
@@ -1023,6 +1044,6 @@
       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)')), &
-       lerr)) RETURN
+    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        
@@ -1078,6 +1099,6 @@
   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.)
+    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
@@ -1127,7 +1148,7 @@
 END SUBROUTINE getKey_init
 !==============================================================================================================================
-CHARACTER(LEN=maxlen) FUNCTION fgetKey(itr, keyn, ky, def_val) RESULT(out)
-!------------------------------------------------------------------------------------------------------------------------------
-! Purpose: Internal function ; get a key value in string format (this is the returned argument).
+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
@@ -1136,9 +1157,27 @@
   CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: def_val
 !------------------------------------------------------------------------------------------------------------------------------
-  INTEGER :: ik
-  ik = 0; IF(itr>0 .AND. itr<=SIZE(ky)) ik = strIdx(ky(itr)%key(:), keyn)
-  out = '';              IF(ik /= 0) out = ky(itr)%val(ik)           !--- Key was found
-  IF(PRESENT(def_val) .AND. ik == 0) out = def_val                   !--- Default value from arguments
-END FUNCTION fgetKey
+  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)
@@ -1151,99 +1190,108 @@
   CHARACTER(LEN=*),          INTENT(IN)  :: tname
   TYPE(keys_type), OPTIONAL, INTENT(IN)  :: ky(:)
-  INTEGER :: is
-  lerr = .FALSE.
+  CHARACTER(LEN=maxlen) :: tnam
+  INTEGER, ALLOCATABLE  :: is(:)
+  INTEGER :: i, itr
+  tnam = delPhase(strHead(tname,'_',.FALSE.))                        !--- Remove tag and phase
   IF(PRESENT(ky)) THEN
-    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
+    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 = getKeyByName_prv(keyn, tname, tracers(:)%keys); IF(val /= '') RETURN !--- "tracers" and "tnam"
+    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
-    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
+    !--- 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
-
-CONTAINS
-
-FUNCTION getKeyByName_prv(keyn, tname, ky) RESULT(val)
-  CHARACTER(LEN=maxlen)         :: val
-  CHARACTER(LEN=*), INTENT(IN)  :: keyn
-  CHARACTER(LEN=*), INTENT(IN)  :: tname
-  TYPE(keys_type),  INTENT(IN)  :: ky(:)
-  INTEGER :: itr, iky
-  val = ''; iky = 0
-  itr = strIdx(ky(:)%name, tname);                 IF(itr==0) RETURN           !--- Get the index of the wanted tracer
-  IF(itr /= 0) iky = strIdx(ky(itr)%key(:), keyn); IF(iky==0) RETURN           !--- Wanted key    index
-  val = ky(itr)%val(iky)
-END FUNCTION getKeyByName_prv
-
 END FUNCTION getKeyByName_s1
 !==============================================================================================================================
-LOGICAL FUNCTION getKeyByName_sm(keyn, val, tnam, ky) RESULT(lerr)
+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)  :: tnam(:)
-  TYPE(keys_type),  TARGET, OPTIONAL, INTENT(IN)  ::   ky(:)
-  CHARACTER(LEN=maxlen),    POINTER :: n(:)
-  INTEGER :: iq
-  n => tracers(:)%keys%name; IF(PRESENT(tnam)) n => tnam(:)
-  ALLOCATE(val(SIZE(n)))
-  IF(     PRESENT(ky)) lerr = ANY([(getKeyByName_s1(keyn, val(iq), n(iq), ky), iq=1, SIZE(n))])
-  IF(.NOT.PRESENT(ky)) lerr = ANY([(getKeyByName_s1(keyn, val(iq), n(iq)),     iq=1, SIZE(n))])
+  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, tnam, ky) RESULT(lerr)
+LOGICAL FUNCTION getKeyByName_i1(keyn, val, tname, ky) RESULT(lerr)
   CHARACTER(LEN=*),          INTENT(IN)  :: keyn
   INTEGER,                   INTENT(OUT) :: val
-  CHARACTER(LEN=*),          INTENT(IN)  :: tnam
+  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, tnam, ky)
-  IF(.NOT.PRESENT(ky)) lerr = getKeyByName_s1(keyn, sval, tnam)
-  IF(test(fmsg('key "'//TRIM(keyn)//'" or tracer "'//TRIM(tnam)//'" is missing',        modname, lerr), lerr)) RETURN
+  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(tnam)//'" is not an integer', modname, lerr), lerr)) RETURN
+  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, tnam, ky) RESULT(lerr)
+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)  :: tnam(:)
-  TYPE(keys_type),  TARGET, OPTIONAL, INTENT(IN)  ::   ky(:)
-  CHARACTER(LEN=maxlen), POINTER :: n(:)
-  INTEGER :: iq
-  n => tracers(:)%name; IF(PRESENT(tnam)) n => tnam(:)
-  ALLOCATE(val(SIZE(n)))
-  IF(     PRESENT(ky)) lerr = ANY([(getKeyByName_i1(keyn, val(iq), n(iq), ky), iq=1, SIZE(n))])
-  IF(.NOT.PRESENT(ky)) lerr = ANY([(getKeyByName_i1(keyn, val(iq), n(iq)),     iq=1, SIZE(n))])
+  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, tnam, ky) RESULT(lerr)
+LOGICAL FUNCTION getKeyByName_r1(keyn, val, tname, ky) RESULT(lerr)
   CHARACTER(LEN=*),          INTENT(IN)  :: keyn
   REAL,                      INTENT(OUT) :: val
-  CHARACTER(LEN=*),          INTENT(IN)  :: tnam
+  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, tnam, ky)
-  IF(.NOT.PRESENT(ky)) lerr = getKeyByName_s1(keyn, sval, tnam)
-  IF(test(fmsg('key "'//TRIM(keyn)//'" or tracer "'//TRIM(tnam)//'" is missing',    modname, lerr), lerr)) RETURN
+  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(tnam)//'" is not a real', modname, lerr), lerr)) RETURN
+  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, tnam, ky) RESULT(lerr)
+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)  :: tnam(:)
-  TYPE(keys_type),  TARGET, OPTIONAL, INTENT(IN)  ::   ky(:)
-  CHARACTER(LEN=maxlen), POINTER :: n(:)
-  INTEGER :: iq
-  n => tracers(:)%name;  IF(PRESENT(tnam)) n => tnam(:)
-  ALLOCATE(val(SIZE(n)))
-  IF(     PRESENT(ky)) lerr = ANY([(getKeyByName_r1(keyn, val(iq), n(iq), ky), iq=1, SIZE(n))])
-  IF(.NOT.PRESENT(ky)) lerr = ANY([(getKeyByName_r1(keyn, val(iq), n(iq)),     iq=1, SIZE(n))])
+  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
 !==============================================================================================================================
@@ -1276,40 +1324,135 @@
 END FUNCTION delPhase
 !------------------------------------------------------------------------------------------------------------------------------
-CHARACTER(LEN=maxlen) FUNCTION addPhase_1(s,pha,ph_sep) RESULT(out)
+CHARACTER(LEN=maxlen) FUNCTION addPhase_s1(s,pha) RESULT(out)
   CHARACTER(LEN=*),           INTENT(IN) :: s
   CHARACTER(LEN=1),           INTENT(IN) :: pha
-  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: ph_sep
-  CHARACTER(LEN=1) :: psep
   INTEGER :: l, i
   out = s
   IF(s == '') RETURN                                                           !--- Empty string: nothing to do
-  psep = phases_sep; IF(PRESENT(ph_sep)) psep = ph_sep
   i = INDEX(s, '_')                                                            !--- /=0 for <var>_<tag> tracers names
   l = LEN_TRIM(s)
-  IF(i == 0) out =  TRIM(s)//TRIM(psep)//pha                                   !--- <var>       => return <var><sep><pha>
-  IF(i /= 0) out = s(1:i-1)//TRIM(psep)//pha//'_'//s(i+1:l)                    !--- <var>_<tag> => return <var><sep><pha>_<tag>
-END FUNCTION addPhase_1
-!------------------------------------------------------------------------------------------------------------------------------
-FUNCTION addPhase_m(s,pha,ph_sep) RESULT(out)
+  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=*), OPTIONAL, INTENT(IN) :: ph_sep
   CHARACTER(LEN=maxlen),     ALLOCATABLE :: out(:)
-  CHARACTER(LEN=1) :: psep
   INTEGER :: k
-  psep = phases_sep; IF(PRESENT(ph_sep)) psep = ph_sep
-  out = [( addPhase_1(s(k), pha, psep), k=1, SIZE(s) )]
-END FUNCTION addPhase_m
-!------------------------------------------------------------------------------------------------------------------------------
-
-CHARACTER(LEN=1) FUNCTION old2newPhase(op) RESULT(np)
-  CHARACTER(LEN=1), INTENT(IN) :: op
-  np = known_phases(INDEX(old_phases,op):INDEX(old_phases,op))
-END FUNCTION old2newPhase
-
-CHARACTER(LEN=1) FUNCTION new2oldPhase(np) RESULT(op)
-  CHARACTER(LEN=1), INTENT(IN) :: np
-  op = old_phases(INDEX(known_phases,np):INDEX(known_phases,np))
-END FUNCTION new2oldPhase
+  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
+  newName = oldName
+  IF(PRESENT(iPhase)) iPhase = 1                                               !--- Default: gaseous phase
+  IF(oldName(1:MIN(3,LEN_TRIM(oldName))) /= 'H2O') RETURN                      !--- Not a water descendant
+  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
+!------------------------------------------------------------------------------------------------------------------------------
+
 
 !==============================================================================================================================
Index: LMDZ6/trunk/libf/misc/strings_mod.F90
===================================================================
--- LMDZ6/trunk/libf/misc/strings_mod.F90	(revision 4119)
+++ LMDZ6/trunk/libf/misc/strings_mod.F90	(revision 4120)
@@ -26,5 +26,5 @@
 !                 horzcat_d1,  horzcat_dm,
                                            horzcat_sm,  horzcat_im,  horzcat_rm; END INTERFACE cat
-  INTERFACE find;       MODULE PROCEDURE      strFind,    find_int,    find_boo; END INTERFACE find
+  INTERFACE find;         MODULE PROCEDURE    strFind,    find_int,    find_boo; END INTERFACE find
   INTERFACE dispOutliers; MODULE PROCEDURE dispOutliers_1, dispOutliers_2; END INTERFACE dispOutliers
   INTERFACE reduceExpr;   MODULE PROCEDURE   reduceExpr_1,   reduceExpr_m; END INTERFACE reduceExpr
@@ -105,12 +105,11 @@
   LOGICAL,          OPTIONAL, INTENT(IN) :: ll
   INTEGER,          OPTIONAL, INTENT(IN) :: unit
+  CHARACTER(LEN=maxlen) :: subn
   INTEGER :: unt
+  subn = '';    IF(PRESENT(modname)) subn = modname
   IF(PRESENT(ll)) THEN; IF(.NOT.ll) RETURN; END IF
   unt = lunout; IF(PRESENT(unit)) unt = unit
-  IF(PRESENT(modname)) THEN
-    WRITE(unt,'(a)') TRIM(modname)//': '//str              !--- Routine name provided
-  ELSE
-    WRITE(unt,'(a)') str                                   !--- Simple message
-  END IF
+  IF(subn == '') WRITE(unt,'(a)') str                                          !--- Simple message
+  IF(subn /= '') WRITE(unt,'(a)') TRIM(subn)//': '//str                        !--- Routine name provided
 END SUBROUTINE msg_1
 !==============================================================================================================================
@@ -123,15 +122,13 @@
   INTEGER,          OPTIONAL, INTENT(IN) :: nmax
   CHARACTER(LEN=maxlen), ALLOCATABLE :: s(:)
+  CHARACTER(LEN=maxlen) :: subn
   INTEGER :: unt, nmx, k
   LOGICAL :: l
+  subn = '';    IF(PRESENT(modname)) subn = modname
   l   = .TRUE.; IF(PRESENT(ll))     l = ll
   unt = lunout; IF(PRESENT(unit)) unt = unit
   nmx = 128;    IF(PRESENT(nmax)) nmx = nmax
   s = strStackm(str, ', ', nmx)
-  IF(PRESENT(modname)) THEN
-    DO k=1,SIZE(s); CALL msg_1(s(k), modname, l, unt); END DO
-  ELSE
-    DO k=1,SIZE(s); CALL msg_1(s(k), ll=l, unit=unt);  END DO
-  END IF
+  DO k=1,SIZE(s); CALL msg_1(s(k), subn,  l,   unt); END DO
 END SUBROUTINE msg_m
 !==============================================================================================================================
@@ -141,12 +138,10 @@
   LOGICAL,          OPTIONAL, INTENT(IN) :: ll
   INTEGER,          OPTIONAL, INTENT(IN) :: unit
+  CHARACTER(LEN=maxlen) :: subn
   INTEGER :: unt
+  subn = '';    IF(PRESENT(modname)) subn = modname
   l   = .TRUE.; IF(PRESENT(ll))     l = ll
   unt = lunout; IF(PRESENT(unit)) unt = unit
-  IF(PRESENT(modname)) THEN
-    CALL msg_1(str, modname, l, unt)
-  ELSE
-    CALL msg_1(str, ll=l, unit=unt)
-  END IF
+  CALL msg_1(str, subn, l, unt)
 END FUNCTION fmsg_1
 !==============================================================================================================================
@@ -157,13 +152,11 @@
   INTEGER,          OPTIONAL, INTENT(IN) :: unit
   INTEGER,          OPTIONAL, INTENT(IN)  :: nmax
+  CHARACTER(LEN=maxlen) :: subn
   INTEGER :: unt, nmx
+  subn = '';    IF(PRESENT(modname)) subn = modname
   l   = .TRUE.; IF(PRESENT(ll))     l = ll
   unt = lunout; IF(PRESENT(unit)) unt = unit
   nmx = 128;    IF(PRESENT(nmax)) nmx = nmax
-  IF(PRESENT(modname)) THEN
-    CALL msg_m(str, modname, l, unt, nmx)
-  ELSE
-    CALL msg_m(str, ll=l, unit=unt, nmax=nmx)
-  END IF
+  CALL msg_m(str, subn, l, unt, nmx)
 END FUNCTION fmsg_m
 !==============================================================================================================================
@@ -178,5 +171,5 @@
   out = str
   DO k=1,LEN_TRIM(str)
-    IF(str(k:k)>='A'.OR.str(k:k)<='Z') out(k:k)=ACHAR(IACHAR(str(k:k))+32)
+    IF(str(k:k)>='A' .AND. str(k:k)<='Z') out(k:k)=ACHAR(IACHAR(str(k:k))+32)
   END DO
 END FUNCTION strLower
@@ -187,5 +180,5 @@
   out = str
   DO k=1,LEN_TRIM(str)
-    IF(str(k:k)>='a'.OR.str(k:k)<='z') out(k:k)=ACHAR(IACHAR(str(k:k))-32)
+    IF(str(k:k)>='a' .AND. str(k:k)<='z') out(k:k)=ACHAR(IACHAR(str(k:k))-32)
   END DO
 END FUNCTION strUpper
@@ -222,7 +215,7 @@
   lf = .FALSE.; IF(PRESENT(lFirst)) lf=lFirst
   IF(PRESENT(sep)) THEN
-    out = [(strHead_1(str(k),sep,.NOT.lf),    k=1, SIZE(str))]
+    out = [(strHead_1(str(k), sep,   lf), k=1, SIZE(str))]
   ELSE
-    out = [(strHead_1(str(k),lFirst=.NOT.lf), k=1, SIZE(str))]
+    out = [(strHead_1(str(k), lFirst=lf), k=1, SIZE(str))]
   END IF
 END FUNCTION strHead_m
@@ -230,6 +223,6 @@
 !=== Extract the substring following the last (first if lFirst==TRUE) occurrence of separator "sep" in "str"   ================
 !=== Examples for str='a_b_c' and sep='_' and the bash command with the same effect:                           ================
-!===    * strHead(..,.FALSE.) = 'b_c'         ${str#*$sep}                                                     ================
-!===    * strHead(..,.TRUE.)  = 'c'           ${str##*$sep}                                                    ================
+!===    * strTail(..,.FALSE.) = 'c'           ${str#*$sep}                                                     ================
+!===    * strTail(..,.TRUE.)  = 'b_c'         ${str##*$sep}                                                    ================
 !==============================================================================================================================
 CHARACTER(LEN=maxlen) FUNCTION strTail_1(str,sep,lFirst) RESULT(out)
@@ -256,7 +249,7 @@
   lf = .FALSE.; IF(PRESENT(lFirst)) lf=lFirst
   IF(PRESENT(sep)) THEN
-    out = [(strTail_1(str(k),sep,.NOT.lf),    k=1, SIZE(str))]
+    out = [(strTail_1(str(k), sep,   lf), k=1, SIZE(str))]
   ELSE
-    out = [(strTail_1(str(k),lFirst=.NOT.lf), k=1, SIZE(str))]
+    out = [(strTail_1(str(k), lFirst=lf), k=1, SIZE(str))]
   END IF
 END FUNCTION strTail_m
@@ -861,5 +854,5 @@
 !=== The profile "p" describe in which order to pick up the columns from "s", "i" and "r" for display.
 !==============================================================================================================================
-LOGICAL FUNCTION dispTable(p, titles, s, i, r, rFmt, nmax, unit) RESULT(lerr)
+LOGICAL FUNCTION dispTable(p, titles, s, i, r, rFmt, nmax, unit, sub) RESULT(lerr)
   CHARACTER(LEN=*),           INTENT(IN)  :: p             !--- DISPLAY MAP: s/i/r
   CHARACTER(LEN=*),           INTENT(IN)  :: titles(:)     !--- TITLES (ONE EACH COLUMN)
@@ -870,7 +863,8 @@
   INTEGER,          OPTIONAL, INTENT(IN)  :: nmax          !--- Display less than "nrow" rows
   INTEGER,          OPTIONAL, INTENT(IN)  :: unit          !--- Output unit (default: screen)
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: sub           !--- Subroutine name
 
   CHARACTER(LEN=2048) :: row
-  CHARACTER(LEN=maxlen)  :: rFm, el
+  CHARACTER(LEN=maxlen)  :: rFm, el, subn
   CHARACTER(LEN=maxlen), ALLOCATABLE :: d(:,:)
   CHARACTER(LEN=1) :: s1, sp
@@ -881,5 +875,5 @@
   LOGICAL :: ls, li, lr
 
-!  modname = 'dispTable'
+  subn = '';    IF(PRESENT(sub)) subn = sub
   rFm = '*';    IF(PRESENT(rFmt)) rFm = rFmt               !--- Specified format for reals
   unt = lunout; IF(PRESENT(unit)) unt = unit               !--- Specified output unit
@@ -890,19 +884,19 @@
 
   !--- CHECK ARGUMENTS COHERENCE
-  lerr = np /= SIZE(titles); IF(fmsg('string "pattern" length and titles list mismatch', ll=lerr)) RETURN
-  IF(ls) THEN; ns = SIZE(s, DIM=1); ncol = ncol + SIZE(s, DIM=2)
-    lerr = COUNT([(p(ic:ic)=='s', ic=1, np)]) /= SIZE(s, DIM=2)
-  END IF
-  IF(li) THEN; ni = SIZE(i, DIM=1); ncol = ncol + SIZE(i, DIM=2)
-    lerr = COUNT([(p(ic:ic)=='i', ic=1, np)]) /= SIZE(i, DIM=2)
-  END IF
-  IF(lr) THEN; nr = SIZE(r, DIM=1); ncol = ncol + SIZE(r, DIM=2)
-    lerr = COUNT([(p(ic:ic)=='r', ic=1, np)]) /= SIZE(r, DIM=2)
-  END IF
-  IF(fmsg('string "pattern" length and arguments number mismatch', ll=lerr)) RETURN
-  lerr = ncol /= SIZE(titles); IF(fmsg('"titles" length and arguments number mismatch', ll=lerr)) RETURN
-  lerr = ls.AND.li.AND.ns/=ni; IF(fmsg('string and integer arguments lengths mismatch', ll=lerr)) RETURN
-  lerr = ls.AND.lr.AND.ns/=nr; IF(fmsg(   'string and real arguments lengths mismatch', ll=lerr)) RETURN
-  lerr = li.AND.lr.AND.ni/=nr; IF(fmsg(  'integer and real arguments lengths mismatch', ll=lerr)) RETURN
+  lerr = np /= SIZE(titles); IF(fmsg('string "pattern" length and titles list mismatch', subn, lerr)) RETURN
+  IF(ls) THEN
+    ns = SIZE(s, 1); ncol = ncol + SIZE(s, 2); lerr = COUNT([(p(ic:ic)=='s', ic=1, np)]) /= SIZE(s, 2)
+  END IF
+  IF(li) THEN
+    ni = SIZE(i, 1); ncol = ncol + SIZE(i, 2); lerr = COUNT([(p(ic:ic)=='i', ic=1, np)]) /= SIZE(i, 2)
+  END IF
+  IF(lr) THEN
+    nr = SIZE(r, 1); ncol = ncol + SIZE(r, 2); lerr = COUNT([(p(ic:ic)=='r', ic=1, np)]) /= SIZE(r, 2)
+  END IF
+  IF(fmsg('string "pattern" length and arguments number mismatch', subn, lerr)) RETURN
+  lerr = ncol /= SIZE(titles); IF(fmsg('"titles" length and arguments number mismatch', subn, lerr)) RETURN
+  lerr = ls.AND.li.AND.ns/=ni; IF(fmsg('string and integer arguments lengths mismatch', subn, lerr)) RETURN
+  lerr = ls.AND.lr.AND.ns/=nr; IF(fmsg(   'string and real arguments lengths mismatch', subn, lerr)) RETURN
+  lerr = li.AND.lr.AND.ni/=nr; IF(fmsg(  'integer and real arguments lengths mismatch', subn, lerr)) RETURN
   nrow = MAX(ns,ni,nr)+1
   nmx = nrow; IF(PRESENT(nmax)) nmx = MIN(nmx,nmax+1)
@@ -931,8 +925,8 @@
     END DO
     nr = LEN_TRIM(row)-1                                             !--- Final separator removed
-    CALL msg(row(1:nr), unit=unt)
+    CALL msg(row(1:nr), subn, unit=unt)
     IF(ir /= 1) CYCLE                                                !--- Titles are underlined
     row=''; DO ic=1,ncol; row=TRIM(row)//REPEAT('-',n(ic))//'+'; END DO
-    CALL msg(row(1:LEN_TRIM(row)-1), unit=unt)
+    CALL msg(row(1:LEN_TRIM(row)-1), subn, unit=unt)
   END DO
 
Index: LMDZ6/trunk/libf/misc/trac_types_mod.F90
===================================================================
--- LMDZ6/trunk/libf/misc/trac_types_mod.F90	(revision 4119)
+++ LMDZ6/trunk/libf/misc/trac_types_mod.F90	(revision 4120)
@@ -20,5 +20,5 @@
     CHARACTER(LEN=maxlen) :: type        = 'tracer'        !--- Type  (so far: 'tracer' / 'tag')
     CHARACTER(LEN=maxlen) :: phase       = 'g'             !--- Phase ('g'as / 'l'iquid / 's'olid)
-    CHARACTER(LEN=maxlen) :: component                     !--- Coma-separated list of components (Ex: lmdz,inca)
+    CHARACTER(LEN=maxlen) :: component   = ''              !--- Coma-separated list of components (Ex: lmdz,inca)
     INTEGER               :: iadv        = 10              !--- Advection scheme used
     INTEGER               :: iGeneration = -1              !--- Generation number (>=0)
@@ -47,8 +47,8 @@
     INTEGER                            :: ntiso = 0        !--- Number of isotopes, including tagging tracers
     INTEGER                            :: nphas = 0        !--- Number phases
-    INTEGER,               ALLOCATABLE :: iTraPha(:,:)     !--- Idx in "trac(1:niso)" = f(name(1:ntiso)),phas)
-                                                           !---        "iTraPha" former name: "iqiso"
-    INTEGER,               ALLOCATABLE :: iZonIso(:,:)     !--- Idx in "trac(1:ntiso)" = f(zone, name(1:niso))
-                                                           !---        "iZonIso" former name: "index_trac"
+    INTEGER,               ALLOCATABLE :: iqTraPha(:,:)    !--- Idx in "tracers(1:nqtot)" = f(name(1:ntiso)),phas)
+                                                           !---        "iqTraPha" former name: "iqiso"
+    INTEGER,               ALLOCATABLE :: itZonIso(:,:)    !--- Idx in "trac(1:ntiso)" = f(zone, name(1:niso))
+                                                           !---        "itZonIso" former name: "index_trac"
   END TYPE isot_type
 
Index: LMDZ6/trunk/libf/phylmd/infotrac_phy.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/infotrac_phy.F90	(revision 4119)
+++ LMDZ6/trunk/libf/phylmd/infotrac_phy.F90	(revision 4120)
@@ -4,29 +4,148 @@
 MODULE infotrac_phy
 
-! Infotrac for physics; for now contains the same information as infotrac for
-! the dynamics (could be further cleaned) and is initialized using values
-! provided by the dynamics
-
-  USE readTracFiles_mod, ONLY: trac_type, maxlen, delPhase
-
-! nqtot : total number of tracers and higher order of moment, water vapor and liquid included
-  INTEGER, SAVE :: nqtot
-!$OMP THREADPRIVATE(nqtot)
-
-!CR: on ajoute le nombre de traceurs de l eau
-  INTEGER, SAVE :: nqo
-!$OMP THREADPRIVATE(nqo)
-
-! nbtr : number of tracers not including higher order of moment or water vapor or liquid
-!        number of tracers used in the physics
-  INTEGER, SAVE :: nbtr
-!$OMP THREADPRIVATE(nbtr)
-
-  INTEGER, SAVE :: nqtottr
-!$OMP THREADPRIVATE(nqtottr)
-
-! ThL : number of CO2 tracers   ModThL
-  INTEGER, SAVE :: nqCO2
-!$OMP THREADPRIVATE(nqCO2)
+   USE       strings_mod, ONLY: msg, maxlen, strStack, strHead, strIdx, int2str
+   USE readTracFiles_mod, ONLY: trac_type, isot_type, keys_type, delPhase, getKey, tnom_iso => newH2OIso
+
+   IMPLICIT NONE
+
+   PRIVATE
+
+   !=== FOR TRACERS:
+   PUBLIC :: init_infotrac_phy                             !--- Initialization of the tracers
+   PUBLIC :: tracers, type_trac                            !--- Full tracers database, tracers type keyword
+   PUBLIC :: nqtot,   nbtr,   nqo,   nqCO2,   nqtottr      !--- Main dimensions
+   PUBLIC :: conv_flg, pbl_flg, solsym                     !--- Convection & boundary layer activation keys
+
+   !=== FOR ISOTOPES: General
+   PUBLIC :: isotopes,  nbIso                              !--- Derived type, full isotopes families database + nb of families
+   PUBLIC :: isoSelect, ixIso                              !--- Isotopes family selection tool + selected family index
+   !=== FOR ISOTOPES: Specific to water
+   PUBLIC :: iH2O                                          !--- H2O isotopes index
+   !=== FOR ISOTOPES: Depending on the selected isotopes family
+   PUBLIC :: isotope, isoKeys                              !--- Selected isotopes database + associated keys (cf. getKey)
+   PUBLIC :: isoName, isoZone, isoPhas                     !--- Isotopes and tagging zones names, phases
+   PUBLIC :: niso,    nzone,   nphas,   ntiso              !---  " " numbers + isotopes & tagging tracers number
+   PUBLIC :: itZonIso                                      !--- iq = function(tagging zone idx, isotope idx)
+   PUBLIC :: iqTraPha                                      !--- idx of tagging tracer in iName = function(isotope idx, phase idx)
+   PUBLIC :: isoCheck                                      !--- Run isotopes checking routines
+   !=== FOR BOTH TRACERS AND ISOTOPES
+   PUBLIC :: getKey                                        !--- Get a key from "tracers" or "isotope"
+
+   PUBLIC :: ntraciso, ntraceurs_zone, indnum_fn_num, use_iso, index_trac, iqiso
+   PUBLIC :: niso_possibles, ok_isotrac, ok_isotopes, ok_iso_verif
+
+   INTERFACE isoSelect; MODULE PROCEDURE isoSelectByIndex, isoSelectByName; END INTERFACE isoSelect
+
+!=== CONVENTIONS FOR TRACERS NUMBERS:
+!  |--------------------+-----------------------+-----------------+---------------+----------------------------|
+!  | water in different |    water tagging      |  water isotopes | other tracers | additional tracers moments |
+!  | phases: H2O_[gls]  |      isotopes         |                 |               |  for higher order schemes  |
+!  |--------------------+-----------------------+-----------------+---------------+----------------------------|
+!  |                    |                       |                 |               |                            |
+!  |<--     nqo      -->|<-- nqo*niso* nzone -->|<-- nqo*niso  -->|<--  nbtr   -->|<--        (nmom)        -->|         
+!  |                    |                                         |                                            |
+!  |                    |<-- nqo*niso*(nzone+1)  =   nqo*ntiso -->|<--    nqtottr = nbtr + nmom             -->|
+!  |                                                                              = nqtot - nqo*(ntiso+1)      |
+!  |                                                                                                           |
+!  |<--                        nqtrue  =  nbtr + nqo*(ntiso+1)                 -->|                            |
+!  |                                                                                                           |
+!  |<--                        nqtot   =  nqtrue + nmom                                                     -->|
+!  |                                                                                                           |
+!  |-----------------------------------------------------------------------------------------------------------|
+!  NOTES FOR THIS TABLE:
+!  * Used "niso", "nzone" and "ntiso" are components of "isotopes(ip)" for water (isotopes(ip)%parent == 'H2O'),
+!    since water is so far the sole tracers family, except passive CO2, removed from the main tracers table.
+!  * For water, "nqo" is equal to the more general field "isotopes(ip)%nphas".
+!  * "niso", "nzone", "ntiso", "nphas" are defined for other isotopic tracers families, if any.
+!
+!=== DERIVED TYPE EMBEDDING MOST OF THE TRACERS-RELATED QUANTITIES (LENGTH: nqtot)
+!    Each entry is accessible using "%" sign.
+!  |-------------+------------------------------------------------------+-------------+------------------------+
+!  |  entry      | Meaning                                              | Former name | Possible values        |
+!  |-------------+------------------------------------------------------+-------------+------------------------+
+!  | name        | Name (short)                                         | tname       |                        |
+!  | gen0Name    | Name of the 1st generation ancestor                  | /           |                        |
+!  | parent      | Name of the parent                                   | /           |                        |
+!  | longName    | Long name (with adv. scheme suffix) for outputs      | ttext       |                        |
+!  | type        | Type (so far: tracer or tag)                         | /           | tracer,tag             |
+!  | phase       | Phases list ("g"as / "l"iquid / "s"olid)             | /           | [g][l][s]              |
+!  | component   | Name(s) of the merged/cumulated section(s)           | /           | coma-separated names   |
+!  | iadv        | Advection scheme number                              | iadv        | 1-20,30 exc. 3-9,15,19 |
+!  | iGeneration | Generation (>=1)                                     | /           |                        |
+!  | isAdvected  | advected tracers flag (.TRUE. if iadv >= 0)          | /           | nqtrue  .TRUE. values  |
+!  | isInPhysics | tracers not extracted from the main table in physics | /           | nqtottr .TRUE. values  |
+!  | iqParent    | Index of the parent tracer                           | iqpere      | 1:nqtot                |
+!  | iqDescen    | Indexes of the childs       (all generations)        | iqfils      | 1:nqtot                |
+!  | nqDescen    | Number of the descendants   (all generations)        | nqdesc      | 1:nqtot                |
+!  | nqChilds    | Number of childs            (1st generation only)    | nqfils      | 1:nqtot                |
+!  | iso_iGroup  | Isotopes group index in isotopes(:)                  | /           | 1:nbIso                |
+!  | iso_iName   | Isotope  name  index in isotopes(iso_iGroup)%trac(:) | iso_indnum  | 1:niso                 |
+!  | iso_iZone   | Isotope  zone  index in isotopes(iso_iGroup)%zone(:) | zone_num    | 1:nzone                |
+!  | iso_iPhas   | Isotope  phase index in isotopes(iso_iGroup)%phas(:) | phase_num   | 1:nphas                |
+!  | keys        | key/val pairs accessible with "getKey" routine       | /           |                        |
+!  +-------------+------------------------------------------------------+-------------+------------------------+
+!
+!=== DERIVED TYPE EMBEDDING MOST OF THE ISOTOPES-RELATED QUANTITIES (LENGTH: nbIso, NUMBER OF ISOTOPES FAMILIES)
+!    Each entry is accessible using "%" sign.
+!  |-----------------+--------------------------------------------------+--------------------+-----------------+
+!  |  entry | length | Meaning                                          |    Former name     | Possible values |
+!  |-----------------+--------------------------------------------------+--------------------+-----------------+
+!  | parent          | Parent tracer (isotopes family name)             |                    |                 |
+!  | keys   | niso   | Isotopes keys/values pairs list + number         |                    |                 |
+!  | trac   | ntiso  | Isotopes + tagging tracers list + number         | / | ntraciso       |                 |
+!  | zone   | nzone  | Geographic tagging zones   list + number         | / | ntraceurs_zone |                 |
+!  | phase  | nphas  | Phases                     list + number         |                    | [g][l][s], 1:3  |
+!  | iqTraPha        | Index in "qx"           = f(name(1:ntiso)),phas) | iqiso              | 1:nqtot         |
+!  | itZonIso        | Index in "trac(1:ntiso)"= f(zone, name(1:niso))  | index_trac         | 1:ntiso         |
+!  +-----------------+--------------------------------------------------+--------------------+-----------------+
+
+   !=== DIMENSIONS OF THE TRACERS TABLES AND OTHER SCALAR VARIABLES
+   INTEGER,                 SAVE :: nqtot,  &                   !--- Tracers nb in dynamics (incl. higher moments + H2O)
+                                    nbtr,   &                   !--- Tracers nb in physics  (excl. higher moments + H2O)
+                                    nqo,    &                   !--- Number of water phases
+                                    nbIso,  &                   !--- Number of available isotopes family
+                                    nqtottr, &                  !--- Number of tracers passed to phytrac (TO BE DELETED ?)
+                                    nqCO2                       !--- Number of tracers of CO2  (ThL)
+   CHARACTER(LEN=maxlen),   SAVE :: type_trac                   !--- Keyword for tracers type
+!$OMP THREADPRIVATE(nqtot, nbtr, nqo, nbIso, nqtottr, nqCO2, type_trac)
+
+   !=== DERIVED TYPES EMBEDDING MOST INFORMATIONS ABOUT TRACERS AND ISOTOPES FAMILIES
+   TYPE(trac_type), TARGET, SAVE, ALLOCATABLE ::  tracers(:)    !=== TRACERS DESCRIPTORS VECTOR
+   TYPE(isot_type), TARGET, SAVE, ALLOCATABLE :: isotopes(:)    !=== ISOTOPES PARAMETERS VECTOR
+!$OMP THREADPRIVATE(tracers, isotopes)
+
+   !=== ALIASES FOR CURRENTLY SELECTED ISOTOPES FAMILY OF VARIABLES EMBEDDED IN THE VECTOR "isotopes"
+   TYPE(isot_type),         SAVE, POINTER :: isotope            !--- CURRENTLY SELECTED ISOTOPES FAMILY DESCRIPTOR
+   INTEGER,                 SAVE          :: ixIso, iH2O        !--- Index of the selected isotopes family and H2O family
+   LOGICAL,                 SAVE, POINTER :: isoCheck           !--- Flag to trigger the checking routines
+   TYPE(keys_type),         SAVE, POINTER :: isoKeys(:)         !--- ONE SET OF KEYS FOR EACH ISOTOPE (LISTED IN isoName)
+   CHARACTER(LEN=maxlen),   SAVE, POINTER :: isoName(:),   &    !--- ISOTOPES NAMES FOR THE CURRENTLY SELECTED FAMILY
+                                             isoZone(:),   &    !--- TAGGING ZONES  FOR THE CURRENTLY SELECTED FAMILY
+                                             isoPhas            !--- USED PHASES    FOR THE CURRENTLY SELECTED FAMILY
+   INTEGER,                 SAVE, POINTER ::  niso, nzone, &    !--- NUMBER OF ISOTOPES, TAGGING ZONES AND PHASES
+                                             nphas, ntiso, &    !--- NUMBER OF PHASES AND ISOTOPES + ISOTOPIC TAGGING TRACERS
+                                            itZonIso(:,:), &    !--- INDEX IN "isoTrac" AS f(tagging zone idx,  isotope idx)
+                                            iqTraPha(:,:)       !--- INDEX IN "qx"      AS f(isotopic tracer idx, phase idx)
+!$OMP THREADPRIVATE(isotope, ixIso,iH2O, isoCheck, isoKeys, isoName,isoZone,isoPhas, niso,nzone,nphas,ntiso, itZonIso,iqTraPha)
+
+   !=== VARIABLES FOR ISOTOPES INITIALIZATION AND FOR INCA
+   INTEGER,          SAVE,    ALLOCATABLE ::conv_flg(:),  &     !--- Convection     activation ; needed for INCA        (nbtr)
+                                             pbl_flg(:)         !--- Boundary layer activation ; needed for INCA        (nbtr)
+   CHARACTER(LEN=8), SAVE,    ALLOCATABLE ::  solsym(:)
+!$OMP THREADPRIVATE(conv_flg, pbl_flg, solsym)
+
+   !--- Aliases for older names + quantities to be removed             (will be replaced by:)
+   INTEGER, POINTER, SAVE :: ntraciso, ntraceurs_zone           !--- -> ntiso, nzone
+!$OMP THREADPRIVATE         (ntraciso, ntraceurs_zone)   
+   INTEGER, POINTER, SAVE :: index_trac(:,:), iqiso(:,:)        !--- -> itZonIso, iqTraPha
+!$OMP THREADPRIVATE         (index_trac,      iqiso)
+   INTEGER, SAVE :: niso_possibles                              !--- suppressed (use effective niso instead)
+!$OMP THREADPRIVATE(niso_possibles)
+   LOGICAL, SAVE :: ok_isotopes, ok_iso_verif, ok_isotrac       !--- -> niso>0, isoCheck, nzone>0
+!$OMP THREADPRIVATE(ok_isotopes, ok_iso_verif, ok_isotrac)
+   LOGICAL, SAVE, ALLOCATABLE :: use_iso(:)                     !--- suppressed
+!$OMP THREADPRIVATE             (use_iso)
+   INTEGER, SAVE, ALLOCATABLE :: indnum_fn_num(:)
+!$OMP THREADPRIVATE             (indnum_fn_num)
 
 #ifdef CPP_StratAer
@@ -38,94 +157,54 @@
 #endif
 
-! Tracers parameters
-  TYPE(trac_type), TARGET, ALLOCATABLE, SAVE :: tracers(:)
-!$OMP THREADPRIVATE(tracers)
-
-! conv_flg(it)=0 : convection desactivated for tracer number it 
-  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE  :: conv_flg
-!$OMP THREADPRIVATE(conv_flg)
-
-! pbl_flg(it)=0  : boundary layer diffusion desactivaded for tracer number it 
-  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE  :: pbl_flg
-!$OMP THREADPRIVATE(pbl_flg)
-
-  CHARACTER(len=4),SAVE :: type_trac
-!$OMP THREADPRIVATE(type_trac)
-  CHARACTER(len=8),DIMENSION(:),ALLOCATABLE, SAVE :: solsym
-!$OMP THREADPRIVATE(solsym)
-   
-    ! CRisi: cas particulier des isotopes
-    LOGICAL,SAVE :: ok_isotopes,ok_iso_verif,ok_isotrac,ok_init_iso
-!$OMP THREADPRIVATE(ok_isotopes,ok_iso_verif,ok_isotrac,ok_init_iso)
-    INTEGER :: niso_possibles   
-    PARAMETER ( niso_possibles=5)
-    real, DIMENSION (niso_possibles),SAVE :: tnat,alpha_ideal
-!$OMP THREADPRIVATE(tnat,alpha_ideal)
-    LOGICAL, DIMENSION(niso_possibles),SAVE ::  use_iso
-!$OMP THREADPRIVATE(use_iso)
-    INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE ::  iqiso ! donne indice iq en fn de (ixt,phase) 
-!$OMP THREADPRIVATE(iqiso)
-    INTEGER, ALLOCATABLE, DIMENSION(:), SAVE ::  iso_indnum ! donne numéro iso entre 1 et niso effectif en fn de nqtot
-!$OMP THREADPRIVATE(iso_indnum)
-    INTEGER, DIMENSION(niso_possibles), SAVE :: indnum_fn_num ! donne indice entre entre 1 et niso en fonction du numéro d isotope entre 1 et niso_possibles
-!$OMP THREADPRIVATE(indnum_fn_num)
-    INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE ::  index_trac ! numéro ixt en fn izone, indnum entre 1 et niso
-!$OMP THREADPRIVATE(index_trac)
-    INTEGER,SAVE :: niso,ntraceurs_zone,ntraciso
-!$OMP THREADPRIVATE(niso,ntraceurs_zone,ntraciso)
-
 CONTAINS
 
-  SUBROUTINE init_infotrac_phy(nqtot_,nqo_,nbtr_,nqtottr_,nqCO2_,tracers_,type_trac_,&
-                               conv_flg_,pbl_flg_,solsym_,&
-                               ok_isotopes_,ok_iso_verif_,ok_isotrac_,&
-                               ok_init_iso_,niso_possibles_,tnat_,&
-                               alpha_ideal_,use_iso_,iqiso_,iso_indnum_,&
-                               indnum_fn_num_,index_trac_,&
-                               niso_,ntraceurs_zone_,ntraciso_)
-
-    ! transfer information on tracers from dynamics to physics
-    USE print_control_mod, ONLY: prt_level, lunout
-    IMPLICIT NONE
-
-    INTEGER,INTENT(IN) :: nqtot_
-    INTEGER,INTENT(IN) :: nqo_
-    INTEGER,INTENT(IN) :: nbtr_
-    INTEGER,INTENT(IN) :: nqtottr_
-    INTEGER,INTENT(IN) :: nqCO2_
-    TYPE(trac_type), INTENT(IN) :: tracers_(nqtot_) ! tracers descriptors
-    CHARACTER(len=*),INTENT(IN) :: type_trac_
-    INTEGER,INTENT(IN) :: conv_flg_(nbtr_)
-    INTEGER,INTENT(IN) :: pbl_flg_(nbtr_)
-    CHARACTER(len=*),INTENT(IN) :: solsym_(nbtr_)
-    ! Isotopes:
-    LOGICAL,INTENT(IN) :: ok_isotopes_
-    LOGICAL,INTENT(IN) :: ok_iso_verif_
-    LOGICAL,INTENT(IN) :: ok_isotrac_
-    LOGICAL,INTENT(IN) :: ok_init_iso_
-    INTEGER,INTENT(IN) :: niso_possibles_
-    REAL,INTENT(IN) :: tnat_(niso_possibles_)
-    REAL,INTENT(IN) :: alpha_ideal_(niso_possibles_)
-    LOGICAL,INTENT(IN) :: use_iso_(niso_possibles_)
-    INTEGER,INTENT(IN) :: iqiso_(ntraciso_,nqo_)
-    INTEGER,INTENT(IN) :: iso_indnum_(nqtot_)
-    INTEGER,INTENT(IN) :: indnum_fn_num_(niso_possibles_)
-    INTEGER,INTENT(IN) :: index_trac_(ntraceurs_zone_,niso_)
-    INTEGER,INTENT(IN) :: niso_
-    INTEGER,INTENT(IN) :: ntraceurs_zone_
-    INTEGER,INTENT(IN) :: ntraciso_
-
-    INTEGER :: iq, itr
-    CHARACTER(LEN=maxlen), ALLOCATABLE :: tnames(:)
-    CHARACTER(LEN=maxlen) :: modname="init_infotrac_phy"
-
-    nqtot=nqtot_
-    nqo=nqo_
-    nbtr=nbtr_
-    nqCO2=nqCO2_
-    nqtottr=nqtottr_
-    ALLOCATE(tracers(nqtot)); tracers(:) = tracers_(:)
+SUBROUTINE init_infotrac_phy(type_trac_, tracers_, isotopes_, nqtottr_, nqCO2_, pbl_flg_, conv_flg_, solsym_)
+
+   USE print_control_mod, ONLY: prt_level, lunout
+
+   IMPLICIT NONE
+   CHARACTER(LEN=*),INTENT(IN) :: type_trac_
+   TYPE(trac_type), INTENT(IN) ::  tracers_(:)
+   TYPE(isot_type), INTENT(IN) :: isotopes_(:)
+   INTEGER,         INTENT(IN) :: nqtottr_
+   INTEGER,         INTENT(IN) :: nqCO2_
+   INTEGER,         INTENT(IN) :: conv_flg_(:)
+   INTEGER,         INTENT(IN) ::  pbl_flg_(:)
+   CHARACTER(LEN=*),INTENT(IN) ::   solsym_(:)
+
+   INTEGER :: iq, ixt
 #ifdef CPP_StratAer
-    IF (type_trac == 'coag') THEN
+   CHARACTER(LEN=maxlen), ALLOCATABLE :: tnames(:)
+#endif
+   CHARACTER(LEN=maxlen) :: modname="init_infotrac_phy"
+
+   type_trac = type_trac_
+   tracers   = tracers_
+   isotopes  = isotopes_
+   nqtottr   = nqtottr_
+   nqCO2     = nqCO2_
+   pbl_flg   =  pbl_flg_
+   conv_flg  = conv_flg_
+   solsym    = solsym_
+   nqtot     = SIZE(tracers_)
+   nqo       = COUNT(delPhase(tracers%name)=='H2O' .AND. tracers%iGeneration==0)
+   nbtr      = SIZE(conv_flg)
+   nbIso     = SIZE(isotopes_)
+
+   !=== Determine selected isotopes class related quantities:
+   !    ixIso, isotope, niso,isoKeys, ntiso,isoName, nzone,isoZone, nphas,isoPhas, itZonIso, iqTraPha, isoCheck
+   IF(.NOT.isoSelect('H2O')) iH2O = ixIso
+   IF(prt_level > 1) THEN
+      CALL msg('nqtot   = '//TRIM(int2str(nqtot)),   modname)
+      CALL msg('nbtr    = '//TRIM(int2str(nbtr )),   modname)
+      CALL msg('nqo     = '//TRIM(int2str(nqo  )),   modname)
+      CALL msg('niso    = '//TRIM(int2str(niso )),   modname)
+      CALL msg('ntiso   = '//TRIM(int2str(ntiso)),   modname)
+      CALL msg('nqtottr = '//TRIM(int2str(nqtottr)), modname)
+      CALL msg('nqCO2   = '//TRIM(int2str(nqCO2)),   modname)
+   END IF
+
+#ifdef CPP_StratAer
+   IF (type_trac == 'coag') THEN
       nbtr_bin    = COUNT([(tracers(iq)%name(1:3)=='BIN', iq=1, nqtot)])
       nbtr_sulgas = COUNT([(tracers(iq)%name(1:3)=='GAS', iq=1, nqtot)])
@@ -136,66 +215,73 @@
       id_H2SO4_strat = strIdx(tnames, 'GASH2SO4')
       id_TEST_strat  = strIdx(tnames, 'GASTEST' )
-      WRITE(lunout,*)'nbtr_bin       =', nbtr_bin
-      WRITE(lunout,*)'nbtr_sulgas    =', nbtr_sulgas
-      WRITE(lunout,*)'id_BIN01_strat =', id_BIN01_strat
-      WRITE(lunout,*)'id_OCS_strat   =',   id_OCS_strat
-      WRITE(lunout,*)'id_SO2_strat   =',   id_SO2_strat
-      WRITE(lunout,*)'id_H2SO4_strat =', id_H2SO4_strat
-      WRITE(lunout,*)'id_TEST_strat  =',  id_TEST_strat
-    END IF
-#endif
-    type_trac = type_trac_
-    ALLOCATE(conv_flg(nbtr))
-    conv_flg(:)=conv_flg_(:)
-    ALLOCATE(pbl_flg(nbtr))
-    pbl_flg(:)=pbl_flg_(:)
-    ALLOCATE(solsym(nbtr))
-    solsym(:)=solsym_(:)
-     
-    IF(prt_level.ge.1) THEN
-      write(lunout,*) TRIM(modname)//": nqtot,nqo,nbtr,nqCO2",nqtot,nqo,nbtr,nqCO2
-    ENDIF
-    
-    ! Isotopes:
-    
-    ! First check that the "niso_possibles" has the correct value
-    IF (niso_possibles.ne.niso_possibles_) THEN
-      CALL abort_physic(modname,&
-           "wrong value for parameter niso_possibles in infotrac_phy",1)
-    ENDIF
-    
-    ok_isotopes=ok_isotopes_
-    ok_iso_verif=ok_iso_verif_
-    ok_isotrac=ok_isotrac_
-    ok_init_iso=ok_init_iso_
-    
-    niso=niso_
-    ntraceurs_zone=ntraceurs_zone_
-    ntraciso=ntraciso_
-    
-    IF (ok_isotopes) THEN
-      tnat(:)=tnat_(:)
-      alpha_ideal(:)=alpha_ideal_(:)
-      use_iso(:)=use_iso_(:)
-    
-      ALLOCATE(iqiso(ntraciso,nqo))
-      iqiso(:,:)=iqiso_(:,:)
-      ALLOCATE(iso_indnum(nqtot))
-      iso_indnum(:)=iso_indnum_(:)
-      
-      indnum_fn_num(:)=indnum_fn_num_(:)
-      
-      ALLOCATE(index_trac(ntraceurs_zone,niso))
-      index_trac(:,:)=index_trac_(:,:)
-    ENDIF ! of IF(ok_isotopes)
-
-    WRITE(*,*) 'infotrac_phy 207: nqtottr=',nqtottr
-    WRITE(*,*) 'ntraciso,niso=',ntraciso,niso
+      CALL msg('nbtr_bin       ='//TRIM(int2str(nbtr_bin      )), modname)
+      CALL msg('nbtr_sulgas    ='//TRIM(int2str(nbtr_sulgas   )), modname)
+      CALL msg('id_BIN01_strat ='//TRIM(int2str(id_BIN01_strat)), modname)
+      CALL msg('id_OCS_strat   ='//TRIM(int2str(id_OCS_strat  )), modname)
+      CALL msg('id_SO2_strat   ='//TRIM(int2str(id_SO2_strat  )), modname)
+      CALL msg('id_H2SO4_strat ='//TRIM(int2str(id_H2SO4_strat)), modname)
+      CALL msg('id_TEST_strat  ='//TRIM(int2str(id_TEST_strat )), modname)
+   END IF
+#endif
+
+   !--- Isotopic quantities (to be removed soon)
+   ntraciso       => ntiso
+   ntraceurs_zone => nzone
+   iqiso          => iqTraPha
+   index_trac     => itZonIso
+   ok_isotopes    = niso  > 0
+   ok_isotrac     = nzone > 0
+   ok_iso_verif   = isoCheck
+   niso_possibles = SIZE(tnom_iso)
+   indnum_fn_num  = [(strIdx(isotope%keys(:)%name, tnom_iso(ixt)), ixt=1, niso_possibles)]
+   use_iso        = indnum_fn_num /= 0
 #ifdef ISOVERIF
-    ! DC: the "1" will be replaced by iH2O (H2O isotopes group index)
-    WRITE(*,*) 'iso_iName=',PACK(tracers(:)%iso_iName, MASK=tracers(:)%iso_iGroup==1)
-#endif
-
-  END SUBROUTINE init_infotrac_phy
+   CALL msg('iso_iName = '//strStack(int2str(PACK(tracers(:)%iso_iName, MASK=tracers(:)%iso_iGroup==iH2O))), modname)
+#endif
+
+END SUBROUTINE init_infotrac_phy
+
+
+!==============================================================================================================================
+!=== THE ROUTINE isoSelect IS USED TO SWITCH FROM AN ISOTOPE FAMILY TO ANOTHER: ISOTOPES DEPENDENT PARAMETERS ARE UPDATED
+!     Single generic "isoSelect" routine, using the predefined index of the parent (fast version) or its name (first call).
+!==============================================================================================================================
+LOGICAL FUNCTION isoSelectByName(iName, lVerbose) RESULT(lerr)
+   IMPLICIT NONE
+   CHARACTER(LEN=*),  INTENT(IN)  :: iName
+   LOGICAL, OPTIONAL, INTENT(IN) :: lVerbose
+   INTEGER :: iIso
+   LOGICAL :: lV
+   lV = .FALSE.; IF(PRESENT(lVerbose)) lV = lVerbose
+   iIso = strIdx(isotopes(:)%parent, iName)
+   lerr = iIso == 0
+   CALL msg('no isotope family named "'//TRIM(iName)//'"', ll=lerr .AND. lV)
+   IF(lerr) RETURN
+   lerr = isoSelectByIndex(iIso, lV)
+END FUNCTION isoSelectByName
+!==============================================================================================================================
+LOGICAL FUNCTION isoSelectByIndex(iIso, lVerbose) RESULT(lerr)
+   IMPLICIT NONE
+   INTEGER,           INTENT(IN) :: iIso
+   LOGICAL, OPTIONAL, INTENT(IN) :: lVerbose
+   LOGICAL :: lv
+   lv = .FALSE.; IF(PRESENT(lVerbose)) lv = lVerbose
+   lerr = .FALSE.
+   IF(iIso == ixIso) RETURN                                      !--- Nothing to do if the index is already OK
+   lerr = iIso<=0 .OR. iIso>nbIso
+   CALL msg('Inconsistent isotopes family index '//TRIM(int2str(iIso))//': should be > 0 and <= '//TRIM(int2str(nbIso))//'"',&
+            ll=lerr .AND. lV)
+   IF(lerr) RETURN
+   ixIso = iIso                                                  !--- Update currently selected family index
+   isotope  => isotopes(ixIso)                                   !--- Select corresponding component
+   isoKeys  => isotope%keys;     niso     => isotope%niso
+   isoName  => isotope%trac;     ntiso    => isotope%ntiso
+   isoZone  => isotope%zone;     nzone    => isotope%nzone
+   isoPhas  => isotope%phase;    nphas    => isotope%nphas
+   itZonIso => isotope%itZonIso; isoCheck => isotope%check
+   iqTraPha => isotope%iqTraPha
+END FUNCTION isoSelectByIndex
+!==============================================================================================================================
+
 
 END MODULE infotrac_phy
Index: LMDZ6/trunk/libf/phylmd/phys_output_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phys_output_mod.F90	(revision 4119)
+++ LMDZ6/trunk/libf/phylmd/phys_output_mod.F90	(revision 4120)
@@ -35,5 +35,6 @@
     USE iophy 
     USE dimphy
-    USE infotrac_phy, ONLY: nqtot, tracers, type_trac, niso, ntraciso, maxlen
+    USE infotrac_phy, ONLY: nqtot, tracers, type_trac, niso, ntraciso
+    USE strings_mod,  ONLY: maxlen
     USE ioipsl
     USE phys_cal_mod, only : hour, calend
Index: LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 4119)
+++ LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 4120)
@@ -25,5 +25,5 @@
 
     USE dimphy, ONLY: klon, klev, klevp1
-    USE infotrac_phy, ONLY: nbtr, nqtot, nqo, type_trac, tracers, niso, ntraciso, maxlen
+    USE infotrac_phy, ONLY: nbtr, nqtot, nqo, type_trac, tracers, niso, ntraciso
     USE strings_mod,  ONLY: maxlen
     USE mod_phys_lmdz_para, ONLY: is_north_pole_phy,is_south_pole_phy
Index: LMDZ6/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4119)
+++ LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4120)
@@ -1294,5 +1294,5 @@
        IF ((iflag_ice_thermo.gt.0).and.(nqo==2)) THEN
           WRITE (lunout, *) ' iflag_ice_thermo==1 requires 3 H2O tracers ', &
-               '(H2Ov, H2Ol, H2Oi) but nqo=', nqo, '. Might as well stop here.'
+               '(H2O_g, H2O_l, H2O_s) but nqo=', nqo, '. Might as well stop here.'
           abort_message='see above'
           CALL abort_physic(modname,abort_message,1)
@@ -1307,5 +1307,5 @@
        IF (ok_ice_sursat.AND.(nqo.NE.4)) THEN
           WRITE (lunout, *) ' ok_ice_sursat=y requires 4 H2O tracers ', &
-               '(H2Ov, H2Ol, H2Oi, H2Or) but nqo=', nqo, '. Might as well stop here.'
+               '(H2O_g, H2O_l, H2O_s, H2O_r) but nqo=', nqo, '. Might as well stop here.'
           abort_message='see above'
           CALL abort_physic(modname,abort_message,1)
@@ -2290,6 +2290,5 @@
     ELSE
 ! DC: make sure the final "1" index was meant for 1st H2O phase (vapor) !!!
-!       tr_seri(:,:,strIdx(tracers(:)%name,addPhase('H2O','g'))) = 0.0
-       tr_seri(:,:,strIdx(tracers(:)%name,addPhase('H2O','v',''))) = 0.0
+       tr_seri(:,:,strIdx(tracers(:)%name,addPhase('H2O','g'))) = 0.0
     ENDIF
 !
Index: LMDZ6/trunk/libf/phylmdiso/phys_output_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/phys_output_mod.F90	(revision 4119)
+++ LMDZ6/trunk/libf/phylmdiso/phys_output_mod.F90	(revision 4120)
@@ -35,5 +35,6 @@
     USE iophy 
     USE dimphy
-    USE infotrac_phy, ONLY: nqtot, tracers, type_trac, niso, ntraciso, maxlen
+    USE infotrac_phy, ONLY: nqtot, tracers, type_trac, niso, ntraciso
+    USE strings_mod,  ONLY: maxlen
     USE ioipsl
     USE phys_cal_mod, only : hour, calend
Index: LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 4119)
+++ LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 4120)
@@ -39,7 +39,8 @@
     USE ioipsl_getin_p_mod, ONLY : getin_p
     USE indice_sol_mod
-    USE infotrac_phy, ONLY: nqtot, nbtr, nqo, tracers, type_trac, nqCO2, ok_isotopes
-    USE readTracFiles_mod, ONLY: phases_sep
-    USE strings_mod,  ONLY: strIdx
+    USE infotrac, ONLY: iso_num, iso_indnum
+    USE infotrac_phy, ONLY: nqtot, nbtr, nqo, tracers, type_trac, nqCO2, ok_isotopes, indnum_fn_num
+    USE readTracFiles_mod, ONLY: addPhase
+    USE strings_mod,  ONLY: strIdx, strStack, int2str
     USE iophy
     USE limit_read_mod, ONLY : init_limit_read
@@ -126,5 +127,5 @@
 #ifdef ISO
     USE infotrac_phy, ONLY:  &
-        iqiso,iso_indnum,ok_isotrac,niso, ntraciso
+        iqiso,ok_isotrac,niso, ntraciso
      USE isotopes_mod, ONLY: iso_eau,iso_HDO,iso_O18,iso_O17,iso_HTO, &
         & bidouille_anti_divergence,ok_bidouille_wake, &
@@ -1392,5 +1393,5 @@
        IF ((iflag_ice_thermo.gt.0).and.(nqo==2)) THEN
           WRITE (lunout, *) ' iflag_ice_thermo==1 requires 3 H2O tracers ', &
-               '(H2Ov, H2Ol, H2Oi) but nqo=', nqo, '. Might as well stop here.'
+               '(H2O_g, H2O_l, H2O_s) but nqo=', nqo, '. Might as well stop here.'
           abort_message='see above'
           CALL abort_physic(modname,abort_message,1)
@@ -1405,5 +1406,5 @@
        IF (ok_ice_sursat.AND.(nqo.NE.4)) THEN
           WRITE (lunout, *) ' ok_ice_sursat=y requires 4 H2O tracers ', &
-               '(H2Ov, H2Ol, H2Oi, rnebi) but nqo=', nqo, '. Might as well stop here.'
+               '(H2O_g, H2O_l, H2O_s, H2O_r) but nqo=', nqo, '. Might as well stop here.'
           abort_message='see above'
           CALL abort_physic(modname,abort_message,1)
@@ -2433,7 +2434,5 @@
       endif !if (nqo.eq.3) then
 #endif
-      if (ixt.gt.niso) then
-      write(*,*) 'izone,iiso=',tracers(iqiso(ixt,ivap))%iso_iZone,iso_indnum(iqiso(ixt,ivap))  
-      endif
+      if (ixt.gt.niso) write(*,*) 'izone=',tracers(iqiso(ixt,ivap))%iso_iZone 
       DO k = 1, klev
        DO i = 1, klon
@@ -2494,6 +2493,5 @@
     ELSE
 ! DC: make sure the final "1" index was meant for 1st H2O phase (vapor) !!!
-!       tr_seri(:,:,strIdx(tracers(:)%name,'H2O'//phases_sep//'g')) = 0.0
-       tr_seri(:,:,strIdx(tracers(:)%name,'H2Ov')) = 0.0
+       tr_seri(:,:,strIdx(tracers(:)%name,addPhase('H2O','g'))) = 0.0
     ENDIF
 !