Index: /LMDZ6/trunk/libf/dyn3d/check_isotopes.F90
===================================================================
--- /LMDZ6/trunk/libf/dyn3d/check_isotopes.F90	(revision 5213)
+++ /LMDZ6/trunk/libf/dyn3d/check_isotopes.F90	(revision 5214)
@@ -2,5 +2,6 @@
    USE strings_mod, ONLY: maxlen, msg, strIdx, strStack, int2str, real2str
    USE infotrac,    ONLY: nqtot, niso, nphas, isotope, isoCheck, iqIsoPha, isoSelect, &
-                          ntiso, iH2O, nzone, tracers, isoName,  itZonIso, getKey
+                          ntiso, iH2O, nzone, tracers, isoName,  itZonIso
+   USE iso_params_mod,  ONLY: tnat_H216O, tnat_H217O, tnat_H218O, tnat_HDO, tnat_HTO
 #ifdef CPP_IOIPSL
    USE ioipsl,          ONLY: getin
@@ -15,7 +16,7 @@
    CHARACTER(LEN=maxlen) :: modname, msg1, nm(2)
    INTEGER :: ixt, ipha, k, i, iq, iiso, izon, ieau, iqeau, iqpar
-   INTEGER, ALLOCATABLE ::   ix(:)
+   INTEGER, ALLOCATABLE       ::   ix(:)
    REAL,    ALLOCATABLE, SAVE :: tnat(:)
-   REAL    :: xtractot, xiiso, deltaD, q1, q2
+   REAL :: xtractot, xiiso, deltaD, q1, q2
    REAL, PARAMETER :: borne     = 1e19,  &
                       errmax    = 1e-8,  &       !--- Max. absolute error
@@ -26,5 +27,6 @@
                       ridicule  = 1e-12
    INTEGER, SAVE :: iso_eau, iso_O17, iso_O18, iso_HDO, iso_HTO
-   LOGICAL, SAVE :: ltnat1, first=.TRUE.
+   LOGICAL       :: ltnat1
+   LOGICAL, SAVE :: first=.TRUE.
 
    modname='check_isotopes'
@@ -35,14 +37,10 @@
       ltnat1 = .TRUE.; CALL getin('tnateq1', ltnat1)
       ALLOCATE(tnat(niso))
-      iso_eau = strIdx(isoName,'H216O')
-      iso_O17 = strIdx(isoName,'H217O')
-      iso_O18 = strIdx(isoName,'H218O')
-      iso_HDO = strIdx(isoName,'HDO')
-      iso_HTO = strIdx(isoName,'HTO')
-      IF(ltnat1) THEN
-         tnat(:)=1.0
-      ELSE
-         IF(getKey('tnat', tnat)) CALL abort_gcm(modname, 'missing isotopic parameter', 1)
-      END IF
+      iso_eau = strIdx(isoName,'H216O'); IF(iso_eau /= 0) tnat(iso_eau) = tnat_H216O
+      iso_O17 = strIdx(isoName,'H217O'); IF(iso_O17 /= 0) tnat(iso_O17) = tnat_H217O
+      iso_O18 = strIdx(isoName,'H218O'); IF(iso_O18 /= 0) tnat(iso_O18) = tnat_H218O
+      iso_HDO = strIdx(isoName,'HDO');   IF(iso_HDO /= 0) tnat(iso_HDO) = tnat_HDO
+      iso_HTO = strIdx(isoName,'HTO');   IF(iso_HTO /= 0) tnat(iso_HTO) = tnat_HTO
+      IF(ltnat1) tnat(:) = 1.
       first = .FALSE.
    END IF
@@ -56,5 +54,5 @@
          DO k = 1, llm
             DO i = 1, ip1jmp1
-               IF(ABS(q(i,k,iq)) <= borne) CYCLE
+               IF(ABS(q(i,k,iq)) < borne) CYCLE
                WRITE(msg1,'(s,"(",i0,",",i0,",",i0,") = ",ES12.4)')TRIM(isoName(ixt)),i,k,iq,q(i,k,iq)
                CALL msg(msg1, modname)
Index: /LMDZ6/trunk/libf/dyn3d/dynetat0.F90
===================================================================
--- /LMDZ6/trunk/libf/dyn3d/dynetat0.F90	(revision 5213)
+++ /LMDZ6/trunk/libf/dyn3d/dynetat0.F90	(revision 5214)
@@ -7,5 +7,5 @@
 !-------------------------------------------------------------------------------
   USE infotrac,    ONLY: nqtot, tracers, niso, iqIsoPha, iH2O, isoName, &
-                         new2oldH2O, newHNO3, oldHNO3, getKey
+                         new2oldH2O, newHNO3, oldHNO3
   USE strings_mod, ONLY: maxlen, msg, strStack, real2str, int2str
   USE netcdf,      ONLY: NF90_OPEN,  NF90_NOWRITE, NF90_INQ_VARID, &
@@ -24,4 +24,5 @@
   USE ioipsl_getincom, ONLY: getin
 #endif
+  USE iso_params_mod   ! tnat_* and alpha_ideal_*
 
   IMPLICIT NONE
@@ -166,6 +167,12 @@
             CALL msg(' !!!  Beware: alpha_ideal put to 1  !!!', modname)
          ELSE
-          IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) &
-            CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1)
+            SELECT CASE(isoName(iName))
+              CASE('H216O'); tnat = tnat_H216O; alpha_ideal = alpha_ideal_H216O
+              CASE('H217O'); tnat = tnat_H217O; alpha_ideal = alpha_ideal_H217O
+              CASE('H218O'); tnat = tnat_H218O; alpha_ideal = alpha_ideal_H218O
+              CASE('HDO');   tnat = tnat_HDO;   alpha_ideal = alpha_ideal_HDO
+              CASE('HTO');   tnat = tnat_HTO;   alpha_ideal = alpha_ideal_HTO
+              CASE DEFAULT; CALL abort_gcm(TRIM(modname),'unknown isotope "'//TRIM(isoName(iName))//'" ; check tracer.def file',1)
+            END SELECT
          END IF
          CALL msg('Missing tracer <'//TRIM(var)//'> => initialized with a simplified Rayleigh distillation law.', modname)
Index: /LMDZ6/trunk/libf/dyn3d/iniacademic.F90
===================================================================
--- /LMDZ6/trunk/libf/dyn3d/iniacademic.F90	(revision 5213)
+++ /LMDZ6/trunk/libf/dyn3d/iniacademic.F90	(revision 5214)
@@ -5,5 +5,5 @@
 
   USE filtreg_mod, ONLY: inifilr
-  USE infotrac,    ONLY: nqtot, niso, iqIsoPha, tracers, getKey, isoName, addPhase
+  USE infotrac,    ONLY: nqtot, niso, iqIsoPha, tracers, isoName, addPhase
   USE control_mod, ONLY: day_step,planet_type
   use exner_hyb_m, only: exner_hyb
@@ -23,4 +23,5 @@
   use netcdf, only : NF90_NOWRITE,NF90_OPEN,NF90_NOERR,NF90_INQ_VARID
   use netcdf, only : NF90_CLOSE, NF90_GET_VAR
+  USE iso_params_mod   ! tnat_* and alpha_ideal_*
 
 
@@ -327,6 +328,13 @@
                     WRITE(lunout, *)'In '//TRIM(modname)//': !!!  Beware: alpha_ideal put to 1  !!!'
                  ELSE
-                    IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) &
-                    CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1)
+                    SELECT CASE(isoName(iName))
+                      CASE('H216O'); tnat = tnat_H216O; alpha_ideal = alpha_ideal_H216O
+                      CASE('H217O'); tnat = tnat_H217O; alpha_ideal = alpha_ideal_H217O
+                      CASE('H218O'); tnat = tnat_H218O; alpha_ideal = alpha_ideal_H218O
+                      CASE('HDO');   tnat = tnat_HDO;   alpha_ideal = alpha_ideal_HDO
+                      CASE('HTO');   tnat = tnat_HTO;   alpha_ideal = alpha_ideal_HTO
+                      CASE DEFAULT
+                         CALL abort_gcm(TRIM(modname),'unknown isotope "'//TRIM(isoName(iName))//'" ; check tracer.def file',1)
+                    END SELECT
                  END IF
                  q(:,:,iq) = q(:,:,iqParent)*tnat*(q(:,:,iqParent)/30.e-3)**(alpha_ideal-1.)
Index: DZ6/trunk/libf/dyn3d_common/iso_verif_dyn.F
===================================================================
--- /LMDZ6/trunk/libf/dyn3d_common/iso_verif_dyn.F	(revision 5213)
+++ 	(revision )
@@ -1,119 +1,0 @@
-        function iso_verif_noNaN_nostop(x,err_msg)
-        implicit none
-        ! si x est NaN, on affiche message
-        ! d'erreur et return 1 si erreur
-
-        ! input:
-        real x
-        character*(*) err_msg ! message d''erreur à afficher
-
-        ! output 
-        real borne
-        parameter (borne=1e19)
-        integer iso_verif_noNaN_nostop
-
-        if ((x.gt.-borne).and.(x.lt.borne)) then
-                iso_verif_noNAN_nostop=0
-        else
-            write(*,*) 'erreur detectee par iso_verif_nonNaN:'
-            write(*,*) err_msg
-            write(*,*) 'x=',x
-            iso_verif_noNaN_nostop=1
-        endif      
-
-        return
-        end
-
-        function iso_verif_egalite_nostop
-     :           (a,b,err_msg)
-        implicit none
-        ! compare a et b. Si pas egal, on affiche message
-        ! d'erreur et stoppe
-        ! pour egalite, on verifie erreur absolue et arreur relative
-
-        ! input:
-        real a, b
-        character*(*) err_msg ! message d''erreur à afficher
-
-        ! locals
-        real errmax ! erreur maximale en absolu.
-        real errmaxrel ! erreur maximale en relatif autorisée
-        parameter (errmax=1e-8)
-        parameter (errmaxrel=1e-3)
-
-        ! output
-        integer iso_verif_egalite_nostop
-
-        iso_verif_egalite_nostop=0
-
-        if (abs(a-b).gt.errmax) then
-          if (abs((a-b)/max(max(abs(b),abs(a)),1e-18))
-     :            .gt.errmaxrel) then
-            write(*,*) 'erreur detectee par iso_verif_egalite:'
-            write(*,*) err_msg
-            write(*,*) 'a=',a
-            write(*,*) 'b=',b
-            iso_verif_egalite_nostop=1
-          endif
-        endif     
-        
-        return
-        end        
-
-
-        function iso_verif_aberrant_nostop
-     :           (x,iso,q,err_msg)
-#ifdef CPP_IOIPSL
-        USE IOIPSL, ONLY: getin
-#else
-        USE ioipsl_getincom, ONLY: getin
-#endif
-        USE infotrac, ONLY: isoName, getKey
-        implicit none
-        
-        ! input:
-        real x,q
-        integer iso ! 2=HDO, 1=O18
-        character*(*) err_msg ! message d''erreur à afficher
-
-        ! locals
-        real qmin,deltaD
-        real deltaDmax,deltaDmin,tnat
-        parameter (qmin=1e-11)
-        parameter (deltaDmax=200.0,deltaDmin=-999.9)
-        LOGICAL, SAVE :: ltnat1
-        LOGICAL, SAVE :: lFirst=.TRUE.
-
-        ! output
-        integer iso_verif_aberrant_nostop
-
-        IF(lFirst) THEN
-           ltnat1 = .TRUE.; CALL getin('tnateq1', ltnat1)
-           lFirst = .FALSE.
-        END IF
-        iso_verif_aberrant_nostop=0
-
-        ! verifier que HDO est raisonable
-         if (q.gt.qmin) then
-             IF(ltnat1) THEN
-                tnat = 1.0
-             ELSE IF(getKey('tnat', tnat, isoName(iso))) THEN
-                  err_msg = 'Missing isotopic parameter "tnat"'
-                  iso_verif_aberrant_nostop=1
-                  RETURN
-             END IF
-             deltaD=(x/q/tnat-1)*1000
-             if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then
-                  write(*,*) 'erreur detectee par iso_verif_aberrant:'
-                  write(*,*) err_msg
-                  write(*,*) 'q=',q
-                  write(*,*) 'deltaD=',deltaD
-                  write(*,*) 'iso=',iso
-                  iso_verif_aberrant_nostop=1
-             endif !if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then
-          endif !if (q(i,k,iq).gt.qmin) then
-
-       
-        return
-        end        
-
Index: /LMDZ6/trunk/libf/dyn3d_common/iso_verif_dyn.F90
===================================================================
--- /LMDZ6/trunk/libf/dyn3d_common/iso_verif_dyn.F90	(revision 5214)
+++ /LMDZ6/trunk/libf/dyn3d_common/iso_verif_dyn.F90	(revision 5214)
@@ -0,0 +1,69 @@
+LOGICAL FUNCTION iso_verif_noNaN_nostop(x, err_msg) RESULT(lerr)
+   IMPLICIT NONE
+! If "x" is NaN, print an error message .
+   REAL,             INTENT(IN) :: x
+   CHARACTER(LEN=*), INTENT(IN) :: err_msg
+   REAL, PARAMETER :: borne = 1e19
+   lerr = x > -borne .AND. x < borne
+   IF(.NOT.lerr) RETURN
+   WRITE(*,*) 'erreur detectee par iso_verif_nonNaN:'
+   WRITE(*,*) err_msg
+   WRITE(*,*) 'x=',x
+END FUNCTION iso_verif_noNaN_nostop
+
+
+LOGICAL FUNCTION iso_verif_egalite_nostop(a, b, err_msg) RESULT(lerr)
+   IMPLICIT NONE
+!   Compare "a" and "b". If a/=b, print an error message.
+!   Both absolute and relative errors are checked for equality.
+   REAL,             INTENT(IN) :: a, b
+   CHARACTER(LEN=*), INTENT(IN) :: err_msg
+   REAL, PARAMETER :: errmax = 1e-8, &           ! max absolute error
+                   errmaxrel = 1e-3              ! max relative error
+   lerr = ABS(a-b) > errmax
+   IF(.NOT.lerr) RETURN
+   lerr = ABS( (a-b) / MAX(MAX(ABS(b), ABS(a)),1e-18) ) > errmaxrel
+   IF(.NOT.lerr) RETURN
+   WRITE(*,*) 'erreur detectee par iso_verif_egalite:'
+   WRITE(*,*) err_msg
+   WRITE(*,*) 'a=',a
+   WRITE(*,*) 'b=',b
+END FUNCTION iso_verif_egalite_nostop
+
+
+LOGICAL FUNCTION iso_verif_aberrant_nostop(x, iso, q, err_msg) RESULT(lerr)
+#ifdef CPP_IOIPSL
+   USE IOIPSL, ONLY: getin
+#else
+   USE ioipsl_getincom, ONLY: getin
+#endif
+   USE iso_params_mod, ONLY: tnat_HDO
+   IMPLICIT NONE
+   REAL,             INTENT(IN) :: x, q
+   INTEGER,          INTENT(IN) :: iso ! 2=HDO, 1=O18
+   CHARACTER(LEN=*), INTENT(IN) :: err_msg
+
+   REAL, PARAMETER :: qmin = 1e-11, &
+                 deltaDmax = 200.0, &
+                 deltaDmin =-999.9
+   LOGICAL       :: ltnat1
+   LOGICAL, SAVE :: lFirst=.TRUE.
+   REAL,    SAVE :: tnat
+   REAL          :: deltaD
+   IF(lFirst) THEN
+      ltnat1 = .TRUE.; CALL getin('tnateq1', ltnat1)
+      tnat = tnat_HDO; IF(ltnat1) tnat = 1.0
+      lFirst = .FALSE.
+   END IF
+   lerr = q > qmin
+   IF(.NOT.lerr) RETURN
+   deltaD = (x / q /tnat - 1.) * 1000.
+   lerr = deltaD > deltaDmax .OR. deltaD < deltaDmin
+   IF(.NOT.lerr) RETURN
+   WRITE(*,*) 'erreur detectee par iso_verif_aberrant:'
+   WRITE(*,*) err_msg
+   WRITE(*,*) 'q=',q
+   WRITE(*,*) 'deltaD=',deltaD
+   WRITE(*,*) 'iso=',iso
+END FUNCTION iso_verif_aberrant_nostop
+
Index: /LMDZ6/trunk/libf/dyn3dmem/check_isotopes_loc.F90
===================================================================
--- /LMDZ6/trunk/libf/dyn3dmem/check_isotopes_loc.F90	(revision 5213)
+++ /LMDZ6/trunk/libf/dyn3dmem/check_isotopes_loc.F90	(revision 5214)
@@ -3,5 +3,6 @@
    USE strings_mod, ONLY: maxlen, msg, strIdx, strStack, int2str, real2str
    USE infotrac,    ONLY: nqtot, niso, nphas, isotope, isoCheck, iqIsoPha, isoSelect, &
-                          ntiso, iH2O, nzone, tracers, isoName,  itZonIso, getKey
+                          ntiso, iH2O, nzone, tracers, isoName,  itZonIso
+   USE iso_params_mod, ONLY: tnat_H216O, tnat_H217O, tnat_H218O, tnat_HDO, tnat_HTO
 #ifdef CPP_IOIPSL
    USE ioipsl,          ONLY: getin
@@ -16,7 +17,7 @@
    CHARACTER(LEN=maxlen) :: modname, msg1, nm(2)
    INTEGER :: ixt, ipha, k, i, iq, iiso, izon, ieau, iqeau, iqpar
-   INTEGER, ALLOCATABLE ::   ix(:)
+   INTEGER, ALLOCATABLE       ::   ix(:)
    REAL,    ALLOCATABLE, SAVE :: tnat(:)         !--- OpenMP shared variable
-   REAL    :: xtractot, xiiso, deltaD, q1, q2
+   REAL :: xtractot, xiiso, deltaD, q1, q2
    REAL, PARAMETER :: borne     = 1e19,  &
                       errmax    = 1e-8,  &       !--- Max. absolute error
@@ -40,14 +41,10 @@
       ltnat1 = .TRUE.; CALL getin('tnateq1', ltnat1)
       ALLOCATE(tnat(niso))
-      iso_eau = strIdx(isoName,'H216O')
-      iso_O17 = strIdx(isoName,'H217O')
-      iso_O18 = strIdx(isoName,'H218O')
-      iso_HDO = strIdx(isoName,'HDO')
-      iso_HTO = strIdx(isoName,'HTO')
-      IF(ltnat1) THEN
-         tnat(:)=1.0
-      ELSE
-         IF(getKey('tnat', tnat)) CALL abort_gcm(modname, 'missing isotopic parameter', 1)
-      END IF
+      iso_eau = strIdx(isoName,'H216O'); IF(iso_eau /= 0) tnat(iso_eau) = tnat_H216O
+      iso_O17 = strIdx(isoName,'H217O'); IF(iso_O17 /= 0) tnat(iso_O17) = tnat_H217O
+      iso_O18 = strIdx(isoName,'H218O'); IF(iso_O18 /= 0) tnat(iso_O18) = tnat_H218O
+      iso_HDO = strIdx(isoName,'HDO');   IF(iso_HDO /= 0) tnat(iso_HDO) = tnat_HDO
+      iso_HTO = strIdx(isoName,'HTO');   IF(iso_HTO /= 0) tnat(iso_HTO) = tnat_HTO
+      IF(ltnat1) tnat(:) = 1.0
 !$OMP END MASTER
 !$OMP BARRIER
@@ -64,5 +61,5 @@
          DO k = 1, llm
             DO i = ijb, ije
-               IF(ABS(q(i,k,iq)) <= borne) CYCLE
+               IF(ABS(q(i,k,iq)) < borne) CYCLE
                WRITE(msg1,'(s,"(",i0,",",i0,",",i0,") = ",ES12.4)')TRIM(isoName(ixt)),i,k,iq,q(i,k,iq)
                CALL msg(msg1, modname)
Index: /LMDZ6/trunk/libf/dyn3dmem/dynetat0_loc.F90
===================================================================
--- /LMDZ6/trunk/libf/dyn3dmem/dynetat0_loc.F90	(revision 5213)
+++ /LMDZ6/trunk/libf/dyn3dmem/dynetat0_loc.F90	(revision 5214)
@@ -8,5 +8,5 @@
   USE parallel_lmdz
   USE infotrac,    ONLY: nqtot, tracers, niso, iqIsoPha, iH2O, isoName, &
-                         new2oldH2O, newHNO3, oldHNO3, getKey
+                         new2oldH2O, newHNO3, oldHNO3
   USE strings_mod, ONLY: maxlen, msg, strStack, real2str, int2str, strIdx
   USE netcdf,      ONLY: NF90_OPEN,  NF90_NOWRITE, NF90_INQUIRE_DIMENSION, NF90_INQ_VARID, &
@@ -25,4 +25,5 @@
   USE ioipsl_getincom, ONLY: getin
 #endif
+  USE iso_params_mod   ! tnat_* and alpha_ideal_*
 
   IMPLICIT NONE
@@ -191,6 +192,12 @@
             CALL msg(' !!!  Beware: alpha_ideal put to 1  !!!', modname)
          ELSE
-          IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) &
-            CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1)
+            SELECT CASE(isoName(iName))
+              CASE('H216O'); tnat = tnat_H216O; alpha_ideal = alpha_ideal_H216O
+              CASE('H217O'); tnat = tnat_H217O; alpha_ideal = alpha_ideal_H217O
+              CASE('H218O'); tnat = tnat_H218O; alpha_ideal = alpha_ideal_H218O
+              CASE('HDO');   tnat = tnat_HDO;   alpha_ideal = alpha_ideal_HDO
+              CASE('HTO');   tnat = tnat_HTO;   alpha_ideal = alpha_ideal_HTO
+              CASE DEFAULT; CALL abort_gcm(TRIM(modname),'unknown isotope "'//TRIM(isoName(iName))//'" ; check tracer.def file',1)
+            END SELECT
          END IF
          CALL msg('Missing tracer <'//TRIM(var)//'> => initialized with a simplified Rayleigh distillation law.', modname)
Index: /LMDZ6/trunk/libf/dyn3dmem/iniacademic_loc.F90
===================================================================
--- /LMDZ6/trunk/libf/dyn3dmem/iniacademic_loc.F90	(revision 5213)
+++ /LMDZ6/trunk/libf/dyn3dmem/iniacademic_loc.F90	(revision 5214)
@@ -5,5 +5,5 @@
 
   USE filtreg_mod, ONLY: inifilr
-  USE infotrac,    ONLY: nqtot, niso, iqIsoPha, tracers, getKey, isoName, addPhase
+  USE infotrac,    ONLY: nqtot, niso, iqIsoPha, tracers, isoName, addPhase
   USE control_mod, ONLY: day_step,planet_type
   use exner_hyb_m, only: exner_hyb
@@ -24,4 +24,5 @@
   use netcdf, only : NF90_NOWRITE,NF90_OPEN,NF90_NOERR,NF90_INQ_VARID
   use netcdf, only : NF90_CLOSE, NF90_GET_VAR
+  USE iso_params_mod   ! tnat_* and alpha_ideal_*
 
 
@@ -329,6 +330,13 @@
                     WRITE(lunout, *) 'In '//TRIM(modname)//': !!!  Beware: alpha_ideal put to 1  !!!'
                  ELSE
-                    IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) &
-                    CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1)
+                    SELECT CASE(isoName(iName))
+                      CASE('H216O'); tnat = tnat_H216O; alpha_ideal = alpha_ideal_H216O
+                      CASE('H217O'); tnat = tnat_H217O; alpha_ideal = alpha_ideal_H217O
+                      CASE('H218O'); tnat = tnat_H218O; alpha_ideal = alpha_ideal_H218O
+                      CASE('HDO');   tnat = tnat_HDO;   alpha_ideal = alpha_ideal_HDO
+                      CASE('HTO');   tnat = tnat_HTO;   alpha_ideal = alpha_ideal_HTO
+                      CASE DEFAULT
+                         CALL abort_gcm(TRIM(modname),'unknown isotope "'//TRIM(isoName(iName))//'" ; check tracer.def file',1)
+                    END SELECT
                  END IF
                  q(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqParent)*tnat*(q(ijb_u:ije_u,:,iqParent)/30.e-3)**(alpha_ideal-1.)
Index: /LMDZ6/trunk/libf/misc/iso_params_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/misc/iso_params_mod.F90	(revision 5214)
+++ /LMDZ6/trunk/libf/misc/iso_params_mod.F90	(revision 5214)
@@ -0,0 +1,10 @@
+MODULE iso_params_mod
+
+   REAL, PARAMETER :: tnat_H216O = 1.0000,    alpha_ideal_H216O = 1.000, &
+                      tnat_H217O = 40.000e-6, alpha_ideal_H217O = 1.003, &
+                      tnat_H218O = 2005.2e-6, alpha_ideal_H218O = 1.006, &
+                      tnat_HDO   = 155.76e-6, alpha_ideal_HDO   = 1.010, &
+                      tnat_HTO   = 0.0000,    alpha_ideal_HTO   = 1.000
+
+END MODULE iso_params_mod
+
Index: /LMDZ6/trunk/libf/phylmdiso/isotopes_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/isotopes_mod.F90	(revision 5213)
+++ /LMDZ6/trunk/libf/phylmdiso/isotopes_mod.F90	(revision 5214)
@@ -4,11 +4,42 @@
 MODULE isotopes_mod
    USE strings_mod,  ONLY: msg, real2str, int2str, bool2str, maxlen, strIdx, strStack
-   USE infotrac_phy, ONLY: isoName
+   USE infotrac_phy, ONLY: isoName, niso, ntiso
+   USE iso_params_mod
    IMPLICIT NONE
    INTERFACE get_in; MODULE PROCEDURE getinp_s, getinp_i, getinp_r, getinp_l;  END INTERFACE get_in
-   SAVE
 
   !--- Contains all isotopic variables + their initialization
   !--- Isotopes-specific routines are in isotopes_routines_mod to avoid circular dependencies with isotopes_verif_mod.
+
+   !--- Negligible lower thresholds: no need to check for absurd values under these lower limits
+   REAL, PARAMETER :: &
+      ridicule      = 1e-12,              & ! For mixing ratios
+      ridicule_rain = 1e-8,               & ! For rain fluxes (rain, zrfl...) in kg/s <-> 1e-3 mm/day
+      ridicule_evap = ridicule_rain*1e-2, & ! For evaporations                in kg/s <-> 1e-3 mm/day
+      ridicule_qsol = ridicule_rain,      & ! For qsol                        in kg <-> 1e-8 kg
+      ridicule_snow = ridicule_qsol         ! For snow                        in kg <-> 1e-8 kg
+   REAL, PARAMETER :: expb_max = 30.0
+
+   !--- Fractionation coefficients for H217O
+   REAL, PARAMETER :: fac_coeff_eq17_liq = 0.529,    &
+                      fac_coeff_eq17_ice = 0.529
+
+   !--- H218O reference
+   REAL, PARAMETER :: fac_enrichoce18 = 0.0005,  alpha_liq_sol_O18 = 1.00291,                    &
+                      talph1_O18 = 1137.,        talps1_O18 = 11.839,     tkcin0_O18 = 0.006,    &
+                      talph2_O18 = -0.4156,      talps2_O18 = -0.028244,  tkcin1_O18 = 0.000285, &
+                      talph3_O18 = -2.0667E-3,  tdifrel_O18 = 1./0.9723,  tkcin2_O18 = 0.00082
+
+   !--- Parameters that do not depend on the nature of water isotopes:
+   REAL, PARAMETER :: pxtmelt= 273.15,         & !--- temperature at which ice formation starts
+                      pxtice = 273.15 -  10.0, & !--- temperature at which all condensate is ice:
+                      pxtmin = 273.15 - 120.0, & !--- computation done only under -120°C
+                      pxtmax = 273.15 +  60.0, & !--- computation done only  over  +60°C
+                      tdifexp= 0.58,           & !--- a constant for alpha_eff for equilibrium below cloud base:
+                      tv0cin = 7.0,            & !--- wind threshold (m/s) for smooth/rough regime (Merlivat and Jouzel 1979)
+                      musi   = 1.0,            & !--- facteurs lambda et mu dans Si=musi-lambda*T
+                      Kd     = 2.5e-9,         & !--- diffusion in soil ; m2/s
+                      rh_cste_surf_cond = 0.6, & !--- cste_surf_cond case: rhs and/or Ts set to constants
+                      T_cste_surf_cond = 288.0
 
    !--- Isotopes indices (in [1,niso] ; non-existing => 0 index)
@@ -89,6 +120,6 @@
    !--- Vectors of length "niso"
    REAL, ALLOCATABLE, DIMENSION(:), SAVE :: &
-                    tnat, toce, tcorr, tdifrel
-!$OMP THREADPRIVATE(tnat, toce, tcorr, tdifrel)
+                    alpha, tnat, toce, tcorr, tdifrel
+!$OMP THREADPRIVATE(alpha, tnat, toce, tcorr, tdifrel)
    REAL, ALLOCATABLE, DIMENSION(:), SAVE :: &
                     talph1, talph2, talph3, talps1, talps2
@@ -100,43 +131,4 @@
                     alpha_liq_sol, Rdefault, Rmethox
 !$OMP THREADPRIVATE(alpha_liq_sol, Rdefault, Rmethox)
-!   REAL, SAVE ::    fac_coeff_eq17_liq, fac_coeff_eq17_ice
-!!$OMP THREADPRIVATE(fac_coeff_eq17_liq, fac_coeff_eq17_ice)
-
-   !--- H2[18]O reference
-   REAL, PARAMETER :: fac_enrichoce18=0.0005
-   REAL, PARAMETER :: alpha_liq_sol_O18=1.00291
-   REAL, PARAMETER :: talph1_O18=1137.
-   REAL, PARAMETER :: talph2_O18=-0.4156
-   REAL, PARAMETER :: talph3_O18=-2.0667E-3
-   REAL, PARAMETER :: talps1_O18=11.839
-   REAL, PARAMETER :: talps2_O18=-0.028244
-   REAL, PARAMETER :: tdifrel_O18=1./0.9723
-   REAL, PARAMETER :: tkcin0_O18=0.006
-   REAL, PARAMETER :: tkcin1_O18=0.000285
-   REAL, PARAMETER :: tkcin2_O18=0.00082
-   REAL, PARAMETER :: fac_coeff_eq17_liq=0.529
-   REAL, PARAMETER :: fac_coeff_eq17_ice=0.529
-
-   !---- Parameters that do not depend on the nature of water isotopes:
-   REAL, PARAMETER :: pxtmelt = 273.15 ! temperature at which ice formation starts
-   REAL, PARAMETER :: pxtice  = 273.15-10.0 ! -- temperature at which all condensate is ice:
-   REAL, PARAMETER :: pxtmin = 273.15 - 120.0   ! On ne calcule qu'au dessus de -120°C
-   REAL, PARAMETER :: pxtmax = 273.15 +  60.0   ! On ne calcule qu'au dessus de +60°C
-   REAL, PARAMETER :: tdifexp = 0.58 ! -- a constant for alpha_eff for equilibrium below cloud base:
-   REAL, PARAMETER :: tv0cin  = 7.0 ! wind threshold (m/s) for smooth/rough regime (Merlivat and Jouzel 1979)
-   REAL, PARAMETER :: musi=1.0  ! facteurs lambda et mu dans Si=musi-lambda*T
-   REAL, PARAMETER :: Kd=2.5e-9 ! m2/s ! diffusion dans le sol
-   REAL, PARAMETER :: rh_cste_surf_cond = 0.6 ! cas où cste_surf_cond: on met rhs ou/et Ts cste pour voir
-   REAL, PARAMETER :: T_cste_surf_cond = 288.0
-
-
-   !--- Negligible lower thresholds: no need to check for absurd values under these lower limits
-   REAL, PARAMETER :: &
-      ridicule      = 1e-12,              & ! For mixing ratios
-      ridicule_rain = 1e-8,               & ! For rain fluxes (rain, zrfl...) in kg/s <-> 1e-3 mm/day
-      ridicule_evap = ridicule_rain*1e-2, & ! For evaporations                in kg/s <-> 1e-3 mm/day
-      ridicule_qsol = ridicule_rain,      & ! For qsol                        in kg <-> 1e-8 kg
-      ridicule_snow = ridicule_qsol         ! For snow                        in kg <-> 1e-8 kg
-   REAL, PARAMETER :: expb_max = 30.0
 
    !--- Specific to HTO:
@@ -155,10 +147,8 @@
 
 SUBROUTINE iso_init()
-   USE infotrac_phy,       ONLY: ntiso, niso, getKey
-    USE strings_mod,       ONLY: maxlen
    IMPLICIT NONE
 
    !=== Local variables:
-   INTEGER :: ixt
+   INTEGER :: ixt, is
    LOGICAL :: ltnat1
    CHARACTER(LEN=maxlen) :: modname, sxt
@@ -185,230 +175,203 @@
    CALL msg('64: niso = '//TRIM(int2str(niso)), modname)
 
-   !--- Init de ntracisoOR: on ecrasera en cas de traceurs de tagging isotopiques
-   !                     (nzone>0) si complications avec ORCHIDEE
-   ntracisoOR = ntiso  
-
-   !--- Type of water isotopes:
-   iso_eau = strIdx(isoName, 'H216O'); CALL msg('iso_eau='//int2str(iso_eau), modname)
-   iso_HDO = strIdx(isoName, 'HDO');   CALL msg('iso_HDO='//int2str(iso_HDO), modname)
-   iso_O18 = strIdx(isoName, 'H218O'); CALL msg('iso_O18='//int2str(iso_O18), modname)
-   iso_O17 = strIdx(isoName, 'H217O'); CALL msg('iso_O17='//int2str(iso_O17), modname)
-   iso_HTO = strIdx(isoName, 'HTO');   CALL msg('iso_HTO='//int2str(iso_HTO), modname)
-
-   !--- Initialiaation: reading the isotopic parameters.
-   CALL get_in('lambda',     lambda_sursat, 0.004)
-   CALL get_in('thumxt1',    thumxt1,       0.75*1.2)
-   CALL get_in('ntot',       ntot,          20,  .FALSE.)
-   CALL get_in('h_land_ice', h_land_ice,    20., .FALSE.)
-   CALL get_in('P_veg',      P_veg,         1.0, .FALSE.)
-   CALL get_in('bidouille_anti_divergence', bidouille_anti_divergence, .FALSE.)
-   CALL get_in('essai_convergence',         essai_convergence,         .FALSE.)
-   CALL get_in('initialisation_iso',        initialisation_iso,        0)
-
-!  IF(nzone>0 .AND. initialisation_iso==0) &
-!      CALL get_in('initialisation_isotrac',initialisation_isotrac)
-   CALL get_in('modif_sst',      modif_sst,         0)
-   CALL get_in('deltaTtest',     deltaTtest,      0.0)     !--- For modif_sst>=1
-   CALL get_in('deltaTtestpoles',deltaTtestpoles, 0.0)     !--- For modif_sst>=2
-   CALL get_in( 'sstlatcrit',    sstlatcrit,     30.0)     !--- For modif_sst>=3
-   CALL get_in('dsstlatcrit',   dsstlatcrit,      0.0)     !--- For modif_sst>=3
+         !--- Init de ntracisoOR: on ecrasera en cas de traceurs de tagging isotopiques
+         !                     (nzone>0) si complications avec ORCHIDEE
+         ntracisoOR = ntiso
+
+         !--- Type of water isotopes:
+         iso_eau = strIdx(isoName, 'H216O'); CALL msg('iso_eau='//int2str(iso_eau), modname)
+         iso_HDO = strIdx(isoName, 'HDO');   CALL msg('iso_HDO='//int2str(iso_HDO), modname)
+         iso_O18 = strIdx(isoName, 'H218O'); CALL msg('iso_O18='//int2str(iso_O18), modname)
+         iso_O17 = strIdx(isoName, 'H217O'); CALL msg('iso_O17='//int2str(iso_O17), modname)
+         iso_HTO = strIdx(isoName, 'HTO');   CALL msg('iso_HTO='//int2str(iso_HTO), modname)
+
+         !--- Initialisation: reading the isotopic parameters.
+         CALL get_in('lambda',     lambda_sursat, 0.004); IF(ok_nocinsat) lambda_sursat = 0.
+         CALL get_in('thumxt1',    thumxt1,       0.75*1.2)
+         CALL get_in('ntot',       ntot,          20,  .FALSE.)
+         CALL get_in('h_land_ice', h_land_ice,    20., .FALSE.)
+         CALL get_in('P_veg',      P_veg,         1.0, .FALSE.)
+         CALL get_in('bidouille_anti_divergence', bidouille_anti_divergence, .FALSE.)
+         CALL get_in('essai_convergence',         essai_convergence,         .FALSE.)
+         CALL get_in('initialisation_iso',        initialisation_iso,        0)
+
+!        IF(nzone>0 .AND. initialisation_iso==0) &
+!           CALL get_in('initialisation_isotrac',initialisation_isotrac)
+         CALL get_in('modif_sst',      modif_sst,         0)
+         CALL get_in('deltaTtest',     deltaTtest,      0.0)     !--- For modif_sst>=1
+         CALL get_in('deltaTtestpoles',deltaTtestpoles, 0.0)     !--- For modif_sst>=2
+         CALL get_in( 'sstlatcrit',    sstlatcrit,     30.0)     !--- For modif_sst>=3
+         CALL get_in('dsstlatcrit',   dsstlatcrit,      0.0)     !--- For modif_sst>=3
 #ifdef ISOVERIF
-   CALL msg('iso_init 270:  sstlatcrit='//real2str( sstlatcrit), modname, sstlatcrit < 0.0) !--- For modif_sst>=2
-   CALL msg('iso_init 279: dsstlatcrit='//real2str(dsstlatcrit), modname, sstlatcrit < 0.0) !--- For modif_sst>=3
-   IF(modif_sst >= 2 .AND. sstlatcrit < 0.0) STOP
+         CALL msg('iso_init 270:  sstlatcrit='//real2str( sstlatcrit), modname, sstlatcrit < 0.0) !--- For modif_sst>=2
+         CALL msg('iso_init 279: dsstlatcrit='//real2str(dsstlatcrit), modname, sstlatcrit < 0.0) !--- For modif_sst>=3
+         IF(modif_sst >= 2 .AND. sstlatcrit < 0.0) STOP
 #endif              
 
-   CALL get_in('modif_sic', modif_sic,  0)
-   IF(modif_sic >= 1) &
-   CALL get_in('deltasic',  deltasic, 0.1)
-
-   CALL get_in('albedo_prescrit', albedo_prescrit, 0)
-   IF(albedo_prescrit == 1) THEN
-      CALL get_in('lon_min_albedo', lon_min_albedo, -200.)
-      CALL get_in('lon_max_albedo', lon_max_albedo,  200.)
-      CALL get_in('lat_min_albedo', lat_min_albedo, -100.)
-      CALL get_in('lat_max_albedo', lat_max_albedo,  100.)
-   END IF
-   CALL get_in('deltaP_BL',           deltaP_BL,     10.0)
-   CALL get_in('ruissellement_pluie', ruissellement_pluie, 0)
-   CALL get_in('alphak_stewart',      alphak_stewart,      1)
-   CALL get_in('tdifexp_sol',         tdifexp_sol,      0.67)
-   CALL get_in('calendrier_guide',    calendrier_guide,    0)
-   CALL get_in('cste_surf_cond',      cste_surf_cond,      0)
-   CALL get_in('mixlen',              mixlen,           35.0)
-   CALL get_in('evap_cont_cste',      evap_cont_cste,      0)
-   CALL get_in('deltaO18_evap_cont',  deltaO18_evap_cont,0.0)
-   CALL get_in('d_evap_cont',         d_evap_cont,       0.0)
-   CALL get_in('nudge_qsol',          nudge_qsol,          0)
-   CALL get_in('region_nudge_qsol',   region_nudge_qsol,   1)
-   nlevmaxO17 = 50
-   CALL msg('nlevmaxO17='//TRIM(int2str(nlevmaxO17)))
-   CALL get_in('no_pce',   no_pce,     0)
-   CALL get_in('A_satlim', A_satlim, 1.0)
-   CALL get_in('ok_restrict_A_satlim', ok_restrict_A_satlim, 0)
+         CALL get_in('modif_sic', modif_sic,  0)
+         IF(modif_sic >= 1) &
+         CALL get_in('deltasic',  deltasic, 0.1)
+
+         CALL get_in('albedo_prescrit', albedo_prescrit, 0)
+         IF(albedo_prescrit == 1) THEN
+            CALL get_in('lon_min_albedo', lon_min_albedo, -200.)
+            CALL get_in('lon_max_albedo', lon_max_albedo,  200.)
+            CALL get_in('lat_min_albedo', lat_min_albedo, -100.)
+            CALL get_in('lat_max_albedo', lat_max_albedo,  100.)
+         END IF
+         CALL get_in('deltaO18_oce',        deltaO18_oce,   0.0)
+         CALL get_in('deltaP_BL',           deltaP_BL,     10.0)
+         CALL get_in('ruissellement_pluie', ruissellement_pluie, 0)
+         CALL get_in('alphak_stewart',      alphak_stewart,      1)
+         CALL get_in('tdifexp_sol',         tdifexp_sol,      0.67)
+         CALL get_in('calendrier_guide',    calendrier_guide,    0)
+         CALL get_in('cste_surf_cond',      cste_surf_cond,      0)
+         CALL get_in('mixlen',              mixlen,           35.0)
+         CALL get_in('evap_cont_cste',      evap_cont_cste,      0)
+         CALL get_in('deltaO18_evap_cont',  deltaO18_evap_cont,0.0)
+         CALL get_in('d_evap_cont',         d_evap_cont,       0.0)
+         CALL get_in('nudge_qsol',          nudge_qsol,          0)
+         CALL get_in('region_nudge_qsol',   region_nudge_qsol,   1)
+         nlevmaxO17 = 50
+         CALL msg('nlevmaxO17='//TRIM(int2str(nlevmaxO17)))
+         CALL get_in('no_pce',   no_pce,     0)
+         CALL get_in('A_satlim', A_satlim, 1.0)
+         CALL get_in('ok_restrict_A_satlim', ok_restrict_A_satlim, 0)
 #ifdef ISOVERIF
-   CALL msg(' 315: A_satlim='//real2str(A_satlim), modname, A_satlim > 1.0)
-   IF(A_satlim > 1.0) STOP
+         CALL msg(' 315: A_satlim='//real2str(A_satlim), modname, A_satlim > 1.0)
+         IF(A_satlim > 1.0) STOP
 #endif
-!  CALL get_in('slope_limiterxy',   slope_limiterxy,  2.0) 
-!  CALL get_in('slope_limiterz',    slope_limiterz,   2.0) 
-   CALL get_in('modif_ratqs',       modif_ratqs,        0)
-   CALL get_in('Pcrit_ratqs',       Pcrit_ratqs,    500.0)
-   CALL get_in('ratqsbasnew',       ratqsbasnew,     0.05)
-   CALL get_in('fac_modif_evaoce',  fac_modif_evaoce, 1.0)
-   CALL get_in('ok_bidouille_wake', ok_bidouille_wake,  0)
-   ! si oui, la temperature de cond est celle de l'environnement, pour eviter
-   ! bugs quand temperature dans ascendances convs est mal calculee
-   CALL get_in('cond_temp_env',        cond_temp_env,        .FALSE.)
-   IF(ANY(isoName == 'HTO')) &
-   CALL get_in('ok_prod_nucl_tritium', ok_prod_nucl_tritium, .FALSE., .FALSE.)
-   CALL get_in('tnateq1', ltnat1, .TRUE.)
-
-   ! Ocean composition
-   CALL get_in('deltaO18_oce',  deltaO18_oce, 0.0)
-   
-   CALL msg('iso_O18, iso_HDO, iso_eau = '//TRIM(strStack(int2str([iso_O18, iso_HDO, iso_eau]))), modname)
-
-   !--------------------------------------------------------------
-   ! Isotope fractionation factors and a few isotopic constants
-   !--------------------------------------------------------------
-   ALLOCATE(tkcin0(niso))
-   ALLOCATE(tkcin1(niso))
-   ALLOCATE(tkcin2(niso))
-   ALLOCATE(tnat(niso))
-   ALLOCATE(tdifrel(niso))
-   ALLOCATE(toce(niso))
-   ALLOCATE(tcorr(niso))
-   ALLOCATE(talph1(niso))
-   ALLOCATE(talph2(niso))
-   ALLOCATE(talph3(niso))
-   ALLOCATE(talps1(niso))
-   ALLOCATE(talps2(niso))
-   ALLOCATE(alpha_liq_sol(niso))
-   ALLOCATE(Rdefault(niso))
-   ALLOCATE(Rmethox(niso))
-
-   do ixt=1,niso
-     if (ixt.eq.iso_HTO) then  ! Tritium
-       tkcin0(ixt) = 0.01056
-       tkcin1(ixt) = 0.0005016
-       tkcin2(ixt) = 0.0014432
-       tnat(ixt) = 0.0; IF(ltnat1) tnat(ixt)=1 
-       toce(ixt)=4.0E-19 ! rapport T/H = 0.2 TU Dreisigacker and Roether 1978
-       tcorr(ixt)=1.
-       tdifrel(ixt)=1./0.968
-       talph1(ixt)=46480.
-       talph2(ixt)=-103.87
-       talph3(ixt)=0.
-       talps1(ixt)=46480.
-       talps2(ixt)=-103.87
-       alpha_liq_sol(ixt)=1. 
-       Rmethox(ixt)=0.0
-     endif
-     if (ixt.eq.iso_O17) then  ! O17
-       pente_MWL=0.528 
-       tdifrel(ixt)=1./0.98555 ! valeur utilisée en 1D et dans modèle de LdG ! tdifrel(ixt)=1./0.985452 ! donné par Amaelle
-       fac_kcin= (tdifrel(ixt)-1.0)/(tdifrel_O18-1.0) ! fac_kcin=0.5145 ! donné par Amaelle
-       tkcin0(ixt) = tkcin0_O18*fac_kcin
-       tkcin1(ixt) = tkcin1_O18*fac_kcin
-       tkcin2(ixt) = tkcin2_O18*fac_kcin
-       tnat(ixt)=0.004/100. ! O17 représente 0.004% de l'oxygène
-       IF(ltnat1) tnat(ixt)=1 
-       toce(ixt)=tnat(ixt)*(1.0+deltaO18_oce/1000.0)**pente_MWL
-       tcorr(ixt)=1.0+fac_enrichoce18*pente_MWL ! donné par Amaelle            
-       talph1(ixt)=talph1_O18
-       talph2(ixt)=talph2_O18
-       talph3(ixt)=talph3_O18
-       talps1(ixt)=talps1_O18
-       talps2(ixt)=talps2_O18     
-       alpha_liq_sol(ixt)=(alpha_liq_sol_O18)**fac_coeff_eq17_liq 
-       Rdefault(ixt)=tnat(ixt)*(-3.15/1000.0+1.0)
-       Rmethox(ixt)=(230./1000.+1.)*tnat(ixt) !Zahn et al 2006
-     endif
-     if (ixt.eq.iso_O18) then  ! Oxygene18
-       tkcin0(ixt) = tkcin0_O18
-       tkcin1(ixt) = tkcin1_O18
-       tkcin2(ixt) = tkcin2_O18
-       tnat(ixt)=2005.2E-6; IF(ltnat1) tnat(ixt)=1 
-       toce(ixt)=tnat(ixt)*(1.0+deltaO18_oce/1000.0)
-       tcorr(ixt)=1.0+fac_enrichoce18
-       tdifrel(ixt)=tdifrel_O18
-       talph1(ixt)=talph1_O18
-       talph2(ixt)=talph2_O18
-       talph3(ixt)=talph3_O18
-       talps1(ixt)=talps1_O18
-       talps2(ixt)=talps2_O18
-       alpha_liq_sol(ixt)=alpha_liq_sol_O18   
-       Rdefault(ixt)=tnat(ixt)*(-6.0/1000.0+1.0)
-       Rmethox(ixt)=(130./1000.+1.)*tnat(ixt) !Zahn et al 2006  
-     endif
-     if (ixt.eq.iso_HDO) then ! Deuterium
-       pente_MWL=8.0
-       tdifrel(ixt)=1./0.9755 !          fac_kcin=0.88
-       fac_kcin= (tdifrel(ixt)-1)/(tdifrel_O18-1)
-       tkcin0(ixt) = tkcin0_O18*fac_kcin
-       tkcin1(ixt) = tkcin1_O18*fac_kcin
-       tkcin2(ixt) = tkcin2_O18*fac_kcin
-       tnat(ixt)=155.76E-6; IF(ltnat1) tnat(ixt)=1
-       toce(ixt)=tnat(ixt)*(1.0+pente_MWL*deltaO18_oce/1000.0)
-       tcorr(ixt)=1.0+fac_enrichoce18*pente_MWL          
-       talph1(ixt)=24844.
-       talph2(ixt)=-76.248
-       talph3(ixt)=52.612E-3
-       talps1(ixt)=16288.
-       talps2(ixt)=-0.0934
-       !ZXalpha_liq_sol=1.0192 ! Weston, Ralph, 1955
-       alpha_liq_sol(ixt)=1.0212 
-       ! valeur de Lehmann & Siegenthaler, 1991, Journal of
-       ! Glaciology, vol 37, p 23
-       Rdefault(ixt)=tnat(ixt)*((-6.0*pente_MWL+10.0)/1000.0+1.0)
-       Rmethox(ixt)=tnat(ixt)*(-25.0/1000.+1.) ! Zahn et al 2006
-     endif
-     if (ixt.eq.iso_eau) then ! Oxygene16
-       tkcin0(ixt) = 0.0
-       tkcin1(ixt) = 0.0
-       tkcin2(ixt) = 0.0
-       tnat(ixt)=1.
-       toce(ixt)=tnat(ixt)
-       tcorr(ixt)=1.0
-       tdifrel(ixt)=1.
-       talph1(ixt)=0.
-       talph2(ixt)=0.
-       talph3(ixt)=0.
-       talps1(ixt)=0.
-       talph3(ixt)=0.
-       alpha_liq_sol(ixt)=1.
-       Rdefault(ixt)=tnat(ixt)*1.0
-       Rmethox(ixt)=1.0
-     endif
-   enddo ! ixt=1,niso
-
-   IF(.NOT.Rdefault_smow) then
-        Rdefault(:) = 0.0
-        if (iso_eau.gt.0) Rdefault(iso_eau) = 1.0 ! correction Camille 30 mars 2023
-   ENDIF
-   write(*,*) 'Rdefault=',Rdefault
-   write(*,*) 'toce=',toce
-
-   !--- Sensitivity test: no kinetic effect in sfc evaporation
-   IF(ok_nocinsfc) THEN
-      tkcin0(1:niso) = 0.0
-      tkcin1(1:niso) = 0.0
-      tkcin2(1:niso) = 0.0
-   END IF
-
-   CALL msg('285: verif initialisation:', modname)
-   DO ixt=1,niso
-      sxt=int2str(ixt)
-      CALL msg(' * isoName('//TRIM(sxt)//') = <'//TRIM(isoName(ixt))//'>',  modname)
-      CALL msg(  '    tnat('//TRIM(sxt)//') = '//TRIM(real2str(tnat(ixt))), modname)
-!     CALL msg('    alpha_liq_sol('//TRIM(sxt)//') = '//TRIM(real2str(alpha_liq_sol(ixt))), modname)
-!     CALL msg(        '   tkcin0('//TRIM(sxt)//') = '//TRIM(real2str(tkcin0(ixt))),        modname)
-!     CALL msg(       '   tdifrel('//TRIM(sxt)//') = '//TRIM(real2str(tdifrel(ixt))),       modname)
-   END DO
-   CALL msg('69:     lambda = '//TRIM(real2str(lambda_sursat)), modname)
-   CALL msg('69:    thumxt1 = '//TRIM(real2str(thumxt1)),       modname)
-   CALL msg('69: h_land_ice = '//TRIM(real2str(h_land_ice)),    modname)
-   CALL msg('69:      P_veg = '//TRIM(real2str(P_veg)),         modname)
+!        CALL get_in('slope_limiterxy',   slope_limiterxy,  2.0) 
+!        CALL get_in('slope_limiterz',    slope_limiterz,   2.0) 
+         CALL get_in('modif_ratqs',       modif_ratqs,        0)
+         CALL get_in('Pcrit_ratqs',       Pcrit_ratqs,    500.0)
+         CALL get_in('ratqsbasnew',       ratqsbasnew,     0.05)
+         CALL get_in('fac_modif_evaoce',  fac_modif_evaoce, 1.0)
+         CALL get_in('ok_bidouille_wake', ok_bidouille_wake,  0)
+         ! si oui, la temperature de cond est celle de l'environnement, pour eviter
+         ! bugs quand temperature dans ascendances convs est mal calculee
+         CALL get_in('cond_temp_env',        cond_temp_env,        .FALSE.)
+         IF(ANY(isoName == 'HTO')) &
+         CALL get_in('ok_prod_nucl_tritium', ok_prod_nucl_tritium, .FALSE., .FALSE.)
+         CALL get_in('tnateq1', ltnat1, .TRUE.)
+
+         CALL msg('iso_O18, iso_HDO, iso_eau = '//TRIM(strStack(int2str([iso_O18, iso_HDO, iso_eau]))), modname)
+
+         !--------------------------------------------------------------
+         ! Parameters that depend on the nature of water isotopes:
+         !--------------------------------------------------------------
+         ALLOCATE(tnat (niso), talph1(niso),  talps1(niso), tkcin0(niso), tdifrel (niso), alpha        (niso))
+         ALLOCATE(toce (niso), talph2(niso),  talps2(niso), tkcin1(niso), Rdefault(niso), alpha_liq_sol(niso))
+         ALLOCATE(tcorr(niso), talph3(niso),                tkcin2(niso), Rmethox (niso))
+
+         !=== H216O
+         is = iso_eau
+         IF(is /= 0) THEN
+            tdifrel (is) = 1.0
+            alpha   (is) = alpha_ideal_H216O
+            tnat    (is) = tnat_H216O; IF(ltnat1) tnat(is) = 1.0
+            toce    (is) = tnat(is)
+            tcorr   (is) = 1.0
+            talph1  (is) = 0.0;  talps1(is) = 0.0;  tkcin0(is)  = 0.0
+            talph2  (is) = 0.0;  talps2(is) = 0.0;  tkcin1(is)  = 0.0
+            talph3  (is) = 0.0;                     tkcin2(is)  = 0.0
+            Rdefault(is) = tnat(is)*1.0
+            Rmethox (is) = 1.0
+            alpha_liq_sol(is) = 1.0
+         END IF
+
+         !=== H217O
+         is = iso_O17
+         IF(is /= 0) THEN; pente_MWL = 0.528
+            tdifrel (is) = 1./0.98555  ! used in 1D and in LdG model ; tdifrel=1./0.985452: from Amaelle
+            alpha   (is) = alpha_ideal_H217O
+            tnat    (is) = tnat_H217O; IF(ltnat1) tnat(is) = 1.0
+            toce    (is) = tnat(is)*(1.0+deltaO18_oce/1000.0)**pente_MWL
+            tcorr   (is) = 1.0+fac_enrichoce18*pente_MWL
+            fac_kcin = (tdifrel(is)-1.0)/(tdifrel_O18-1.0)           ! fac_kcin=0.5145:     from Amaelle
+            talph1  (is) = talph1_O18;  talps1(is) = talps1_O18;  tkcin0(is) = tkcin0_O18*fac_kcin
+            talph2  (is) = talph2_O18;  talps2(is) = talps2_O18;  tkcin1(is) = tkcin1_O18*fac_kcin
+            talph3  (is) = talph3_O18;                            tkcin2(is) = tkcin2_O18*fac_kcin
+            Rdefault(is) = tnat(is)*(1.0-3.15/1000.)
+            Rmethox (is) = tnat(is)*(1.0+230./1000.)
+            alpha_liq_sol(is) = alpha_liq_sol_O18**fac_coeff_eq17_liq
+         END IF
+
+         !=== H218O
+         is = iso_O18
+         IF(is /= 0) THEN
+            tdifrel (is) = tdifrel_O18
+            alpha   (is) = alpha_ideal_H218O
+            tnat    (is) = tnat_H218O; IF(ltnat1) tnat(is) = 1.0
+            toce    (is) = tnat(is)*(1.0+deltaO18_oce/1000.0)
+            tcorr   (is) = 1.0+fac_enrichoce18
+            talph1  (is) = talph1_O18;  talps1(is) = talps1_O18;  tkcin0(is) = tkcin0_O18
+            talph2  (is) = talph2_O18;  talps2(is) = talps2_O18;  tkcin1(is) = tkcin1_O18
+            talph3  (is) = talph3_O18;                            tkcin2(is) = tkcin2_O18
+            Rdefault(is) = tnat(is)*(1.0-6.00/1000.)
+            Rmethox (is) = tnat(is)*(1.0+130./1000.)  ! Zahn & al. 2006
+            alpha_liq_sol(is) = alpha_liq_sol_O18
+         END IF
+
+         !=== HDO
+         is = iso_HDO
+         IF(is /= 0) THEN; pente_MWL = 8.0
+            tdifrel (is) = 1./0.9755                  ! fac_kcin=0.88
+            alpha   (is) = alpha_ideal_HDO
+            tnat    (is) = tnat_HDO; IF(ltnat1) tnat(is) = 1.0
+            toce    (is) = tnat(is)*(1.0+deltaO18_oce/1000.0*pente_MWL)
+            tcorr   (is) = 1.0+fac_enrichoce18*pente_MWL
+            fac_kcin = (tdifrel(is)-1.0)/(tdifrel_O18-1.0)
+            talph1  (is) = 24844.;      talps1(is) = 16288.;      tkcin0(is) = tkcin0_O18*fac_kcin
+            talph2  (is) = -76.248;     talps2(is) = -0.0934;     tkcin1(is) = tkcin1_O18*fac_kcin
+            talph3  (is) = 52.612E-3;                             tkcin2(is) = tkcin2_O18*fac_kcin
+            Rdefault(is) = tnat(is)*(1.0+(10.0-6.0*pente_MWL)/1000.)
+            Rmethox (is) = tnat(is)*(1.0-25.0/1000.)
+            alpha_liq_sol(is) = 1.0212      ! Lehmann & Siegenthaler, 1991, Jo. of Glaciology, vol 37, p 23
+                                            ! alpha_liq_sol=1.0192: Weston, Ralph, 1955
+         END IF
+
+         !=== HTO
+         is = iso_HTO
+         IF(is /= 0) THEN
+            tdifrel (is) = 1./0.968
+            alpha   (is) = alpha_ideal_HTO
+            tnat    (is) = tnat_HTO; IF(ltnat1) tnat(is) = 1.0
+            toce    (is) = 4.0E-19          ! ratio T/H = 0.2 TU Dreisigacker & Roether 1978
+            tcorr   (is) = 1.0
+            talph1  (is) = 46480.;      talps1(is) = 46480.;      tkcin0(is) = 0.01056
+            talph2  (is) = -103.87;     talps2(is) = -103.87;     tkcin1(is) = 0.0005016
+            talph3  (is) = 0.0;                                   tkcin2(is) = 0.0014432
+            Rdefault(is) = 0.0
+            Rmethox (is) = 0.0
+            alpha_liq_sol(is) = 1.0
+         END IF
+
+         IF(.NOT. Rdefault_smow) THEN
+            Rdefault(:) = 0.0; IF(iso_eau > 0) Rdefault(iso_eau) = 1.0
+         END IF
+         WRITE(*,*) 'Rdefault = ',Rdefault
+         WRITE(*,*) 'toce = ', toce
+
+         !--- Sensitivity test: no kinetic effect in sfc evaporation
+         IF(ok_nocinsfc) THEN
+            tkcin0(1:niso) = 0.0
+            tkcin1(1:niso) = 0.0
+            tkcin2(1:niso) = 0.0
+         END IF
+
+         CALL msg('285: verif initialisation:', modname)
+         DO ixt=1,niso
+            sxt=int2str(ixt)
+            CALL msg(' * isoName('//TRIM(sxt)//') = <'//TRIM(isoName(ixt))//'>',  modname)
+            CALL msg(  '    tnat('//TRIM(sxt)//') = '//TRIM(real2str(tnat(ixt))), modname)
+!           CALL msg('    alpha_liq_sol('//TRIM(sxt)//') = '//TRIM(real2str(alpha_liq_sol(ixt))), modname)
+!           CALL msg(        '   tkcin0('//TRIM(sxt)//') = '//TRIM(real2str(tkcin0(ixt))),        modname)
+!           CALL msg(       '   tdifrel('//TRIM(sxt)//') = '//TRIM(real2str(tdifrel(ixt))),       modname)
+         END DO
+         CALL msg('69:     lambda = '//TRIM(real2str(lambda_sursat)), modname)
+         CALL msg('69:    thumxt1 = '//TRIM(real2str(thumxt1)),       modname)
+         CALL msg('69: h_land_ice = '//TRIM(real2str(h_land_ice)),    modname)
+         CALL msg('69:      P_veg = '//TRIM(real2str(P_veg)),         modname)
 
 END SUBROUTINE iso_init
