Index: LMDZ6/branches/Portage_acc/libf/phylmdiso/isotopes_mod.F90
===================================================================
--- LMDZ6/branches/Portage_acc/libf/phylmdiso/isotopes_mod.F90	(revision 4128)
+++ LMDZ6/branches/Portage_acc/libf/phylmdiso/isotopes_mod.F90	(revision 4446)
@@ -3,686 +3,404 @@
 
 MODULE isotopes_mod
-USE infotrac_phy, ONLY: ntraciso,niso,indnum_fn_num,use_iso, &
-&       niso_possibles
-IMPLICIT NONE
-SAVE
-
-! contient toutes les variables isotopiques et leur initialisation
-! les routines specifiquement isotopiques sont dans
-! isotopes_routines_mod pour éviter dépendance circulaire avec
-! isotopes_verif_mod.
-
-
-! indices des isotopes
-integer, save :: iso_eau,iso_HDO,iso_O18,iso_O17,iso_HTO ! indices de 1 à niso: les isos n'existant pas sont mis à 0
-!$OMP THREADPRIVATE(iso_eau,iso_HDO,iso_O18,iso_O17,iso_HTO)
-
-integer :: iso_eau_possible,iso_HDO_possible,iso_O18_possible,iso_O17_possible,iso_HTO_possible ! indices de 1 à niso_possibles: ils correspondent aux tableaux définis dans infotrac:
-! tnom_iso=(/'eau','HDO','O18','O17','HTO'/)
-! ce sont ces indices qui doivent être utilisés avec use_iso, puisque use_iso est défini comme DIMENSION(niso_possibles)
-parameter (iso_eau_possible=1)
-parameter (iso_HDO_possible=2)
-parameter (iso_O18_possible=3)
-parameter (iso_O17_possible=4)
-parameter (iso_HTO_possible=5)
-
-integer, save :: ntracisoOR
+   USE strings_mod,  ONLY: msg, real2str, int2str, bool2str, maxlen, strIdx, strStack
+   USE infotrac_phy, ONLY: isoName
+   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.
+
+   !--- Isotopes indices (in [1,niso] ; non-existing => 0 index)
+   INTEGER, SAVE :: iso_eau, iso_HDO, iso_O18, iso_O17, iso_HTO
+!$OMP THREADPRIVATE(iso_eau, iso_HDO, iso_O18, iso_O17, iso_HTO)
+
+   INTEGER, SAVE :: ntracisoOR
 !$OMP THREADPRIVATE(ntracisoOR)
 
-! variables indépendantes des isotopes
-
-real, save :: pxtmelt,pxtice,pxtmin,pxtmax
-!$OMP THREADPRIVATE(pxtmelt,pxtice,pxtmin,pxtmax)
-real, save ::  tdifexp, tv0cin, thumxt1
+   !--- Variables not depending on isotopes
+   REAL,    SAVE :: pxtmelt, pxtice, pxtmin, pxtmax
+!$OMP THREADPRIVATE(pxtmelt, pxtice, pxtmin, pxtmax)
+   REAL,    SAVE :: tdifexp, tv0cin, thumxt1
 !$OMP THREADPRIVATE(tdifexp, tv0cin, thumxt1)
-integer, save ::  ntot 
+   INTEGER, SAVE :: ntot 
 !$OMP THREADPRIVATE(ntot)
-real, save ::  h_land_ice
+   REAL,    SAVE :: h_land_ice
 !$OMP THREADPRIVATE(h_land_ice)
-real, save ::  P_veg
+   REAL,    SAVE :: P_veg
 !$OMP THREADPRIVATE(P_veg)
-real, save ::  musi,lambda_sursat
-!$OMP THREADPRIVATE(lambda_sursat)
-real, save ::  Kd
+   REAL,    SAVE :: musi, lambda_sursat
+!$OMP THREADPRIVATE(musi, lambda_sursat)
+   REAL,    SAVE :: Kd
 !$OMP THREADPRIVATE(Kd)
-real, save ::  rh_cste_surf_cond,T_cste_surf_cond
-!$OMP THREADPRIVATE(rh_cste_surf_cond,T_cste_surf_cond)
-
-logical, save ::   bidouille_anti_divergence
-                ! si true, rappel régulier de xteau vers q, pour éviter dérives lentes
+   REAL,    SAVE :: rh_cste_surf_cond, T_cste_surf_cond
+!$OMP THREADPRIVATE(rh_cste_surf_cond, T_cste_surf_cond)
+   LOGICAL, SAVE :: bidouille_anti_divergence    ! T: regularly, xteau <- q to avoid slow drifts
 !$OMP THREADPRIVATE(bidouille_anti_divergence)
-logical, save ::   essai_convergence 
-                ! si false, on fait rigoureusement comme dans LMDZ sans isotopes, 
-                ! meme si c'est génant pour les isotopes
+   LOGICAL, SAVE :: essai_convergence            ! F: as in LMDZ without isotopes (bad for isotopes)
 !$OMP THREADPRIVATE(essai_convergence)
-integer, save ::   initialisation_iso
-                ! 0: dans fichier
-                ! 1: R=0
-                ! 2: R selon distill rayleigh
-                ! 3: R=Rsmow
+   INTEGER, SAVE :: initialisation_iso           ! 0: file ; 1: R=0 ; 2: R=distill. Rayleigh ; 3: R=Rsmow
 !$OMP THREADPRIVATE(initialisation_iso)
-integer, save ::  modif_SST ! 0 par defaut, 1 si on veut modifier la sst
-                ! 2 et 3: profils de SST
+   INTEGER, SAVE :: modif_SST                    ! 0: default ; 1: modified SST ; 2, 3: SST profiles
 !$OMP THREADPRIVATE(modif_SST)
-real, save ::  deltaTtest ! modif de la SST, uniforme.
+   REAL,    SAVE :: deltaTtest                   ! Uniform modification of the SST
 !$OMP THREADPRIVATE(deltaTtest)
-integer, save ::  modif_sic ! on met des trous dans glace de mer
+   INTEGER, SAVE :: modif_sic                    ! Holes in the Sea Ice
 !$OMP THREADPRIVATE(modif_sic)
-real, save ::  deltasic ! fraction de trous minimale
+   REAL,    SAVE :: deltasic                     ! Minimal holes fraction
 !$OMP THREADPRIVATE(deltasic)
-real, save ::  deltaTtestpoles
+   REAL,    SAVE :: deltaTtestpoles
 !$OMP THREADPRIVATE(deltaTtestpoles)
-real, save ::  sstlatcrit
-!$OMP THREADPRIVATE(sstlatcrit)
-real, save ::  dsstlatcrit
-!$OMP THREADPRIVATE(dsstlatcrit)
-real, save ::  deltaO18_oce
+   REAL,    SAVE :: sstlatcrit, dsstlatcrit
+!$OMP THREADPRIVATE(sstlatcrit, dsstlatcrit)
+   REAL,    SAVE :: deltaO18_oce
 !$OMP THREADPRIVATE(deltaO18_oce)
-integer, save ::  albedo_prescrit ! 0 par defaut
-                        ! 1 si on veut garder albedo constant
+   INTEGER, SAVE :: albedo_prescrit              ! 0: default ; 1: constant albedo
 !$OMP THREADPRIVATE(albedo_prescrit)
-real, save ::  lon_min_albedo,lon_max_albedo
-!$OMP THREADPRIVATE(lon_min_albedo,lon_max_albedo)
-real, save :: lat_min_albedo,lat_max_albedo
-!$OMP THREADPRIVATE(lat_min_albedo,lat_max_albedo)
-real, save ::  deltaP_BL,tdifexp_sol
+   REAL,    SAVE :: lon_min_albedo, lon_max_albedo, lat_min_albedo, lat_max_albedo
+!$OMP THREADPRIVATE(lon_min_albedo, lon_max_albedo, lat_min_albedo, lat_max_albedo)
+   REAL,    SAVE :: deltaP_BL,tdifexp_sol
 !$OMP THREADPRIVATE(deltaP_BL,tdifexp_sol)
-integer, save ::  ruissellement_pluie,alphak_stewart
-!$OMP THREADPRIVATE(ruissellement_pluie,alphak_stewart)
-integer, save ::  calendrier_guide
+   INTEGER, SAVE :: ruissellement_pluie, alphak_stewart
+!$OMP THREADPRIVATE(ruissellement_pluie, alphak_stewart)
+   INTEGER, SAVE :: calendrier_guide
 !$OMP THREADPRIVATE(calendrier_guide)
-integer, save ::  cste_surf_cond
+   INTEGER, SAVE :: cste_surf_cond
 !$OMP THREADPRIVATE(cste_surf_cond)
-real, save ::  mixlen
+   REAL,    SAVE :: mixlen
 !$OMP THREADPRIVATE(mixlen)
-integer, save ::  evap_cont_cste
+   INTEGER, SAVE :: evap_cont_cste
 !$OMP THREADPRIVATE(evap_cont_cste)
-real, save ::  deltaO18_evap_cont,d_evap_cont
-!$OMP THREADPRIVATE(deltaO18_evap_cont,d_evap_cont)
-integer, save ::  nudge_qsol,region_nudge_qsol
-!$OMP THREADPRIVATE(nudge_qsol,region_nudge_qsol)
-integer, save :: nlevmaxO17
+   REAL,    SAVE :: deltaO18_evap_cont, d_evap_cont
+!$OMP THREADPRIVATE(deltaO18_evap_cont, d_evap_cont)
+   INTEGER, SAVE :: nudge_qsol, region_nudge_qsol
+!$OMP THREADPRIVATE(nudge_qsol, region_nudge_qsol)
+   INTEGER, SAVE :: nlevmaxO17
 !$OMP THREADPRIVATE(nlevmaxO17)
-integer, save ::  no_pce
-!	real, save :: slope_limiterxy,slope_limiterz
+   INTEGER, SAVE :: no_pce
 !$OMP THREADPRIVATE(no_pce)
-real, save ::  A_satlim
+   REAL,    SAVE :: A_satlim
 !$OMP THREADPRIVATE(A_satlim)
-integer, save ::  ok_restrict_A_satlim,modif_ratqs
-!$OMP THREADPRIVATE(ok_restrict_A_satlim,modif_ratqs)
-real, save ::  Pcrit_ratqs,ratqsbasnew
-!$OMP THREADPRIVATE(Pcrit_ratqs,ratqsbasnew)
-real, save ::  fac_modif_evaoce
+   INTEGER, SAVE :: ok_restrict_A_satlim, modif_ratqs
+!$OMP THREADPRIVATE(ok_restrict_A_satlim, modif_ratqs)
+   REAL,    SAVE :: Pcrit_ratqs, ratqsbasnew
+!$OMP THREADPRIVATE(Pcrit_ratqs, ratqsbasnew)
+   REAL,    SAVE :: fac_modif_evaoce
 !$OMP THREADPRIVATE(fac_modif_evaoce)
-integer, save ::  ok_bidouille_wake
+   INTEGER, SAVE :: ok_bidouille_wake
 !$OMP THREADPRIVATE(ok_bidouille_wake)
-logical ::  cond_temp_env
+   LOGICAL, SAVE :: cond_temp_env
 !$OMP THREADPRIVATE(cond_temp_env)
 
-
-! variables tableaux fn de niso
-real, ALLOCATABLE, DIMENSION(:), save :: tnat, toce, tcorr
-!$OMP THREADPRIVATE(tnat, toce, tcorr)
-real, ALLOCATABLE, DIMENSION(:), save :: tdifrel
-!$OMP THREADPRIVATE(tdifrel)
-real, ALLOCATABLE, DIMENSION(:), save :: talph1, talph2, talph3
-!$OMP THREADPRIVATE(talph1, talph2, talph3)
-real, ALLOCATABLE, DIMENSION(:), save :: talps1, talps2
-!$OMP THREADPRIVATE(talps1, talps2)
-real, ALLOCATABLE, DIMENSION(:), save :: tkcin0, tkcin1, tkcin2
+   !--- Vectors of length "niso"
+   REAL, ALLOCATABLE, DIMENSION(:), SAVE :: &
+                    tnat, toce, tcorr, tdifrel
+!$OMP THREADPRIVATE(tnat, toce, tcorr, tdifrel)
+   REAL, ALLOCATABLE, DIMENSION(:), SAVE :: &
+                    talph1, talph2, talph3, talps1, talps2
+!$OMP THREADPRIVATE(talph1, talph2, talph3, talps1, talps2)
+   REAL, ALLOCATABLE, DIMENSION(:), SAVE :: &
+                    tkcin0, tkcin1, tkcin2
 !$OMP THREADPRIVATE(tkcin0, tkcin1, tkcin2)
-real, ALLOCATABLE, DIMENSION(:), save :: alpha_liq_sol
-!$OMP THREADPRIVATE(alpha_liq_sol)
-real, ALLOCATABLE, DIMENSION(:), save :: Rdefault, Rmethox
-!$OMP THREADPRIVATE(Rdefault, Rmethox)
-character*3, ALLOCATABLE, DIMENSION(:), save :: striso
-!$OMP THREADPRIVATE(striso)
-real, save :: fac_coeff_eq17_liq, fac_coeff_eq17_ice
+   REAL, ALLOCATABLE, DIMENSION(:), SAVE :: &
+                    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)
 
-      real ridicule ! valeur maximale pour qu'une variable de type
-                    ! rapoport de mélange puisse être considérée comme négligeable. Si
-                    ! négligeable, alors on ne verifie pas si sa compo iso esta bérrante.
-      parameter (ridicule=1e-12)      
-!      parameter (ridicule=1)
-!
-      real ridicule_rain ! valeur limite de ridicule pour les flux de pluies (rain, zrfl...)
-      parameter (ridicule_rain=1e-8) ! en kg/s <-> 1e-3mm/day
-
-      real ridicule_evap ! valeur limite de ridicule pour les evap
-      parameter (ridicule_evap=ridicule_rain*1e-2) ! en kg/s <-> 1e-3mm/day
-
-      real ridicule_qsol ! valeur limite de ridicule pour les qsol
-      parameter (ridicule_qsol=ridicule_rain) ! en kg <-> 1e-8kg
-
-      real ridicule_snow ! valeur limite de ridicule pour les snow
-      parameter (ridicule_snow=ridicule_qsol) ! en kg/s <-> 1e-8kg
-      
-        real expb_max
-        parameter (expb_max=30.0)
-
-        ! spécifique au tritium:
-        
-
-logical, save :: ok_prod_nucl_tritium ! si oui, production de tritium par essais nucleaires
+   !--- 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:
+   LOGICAL, SAVE :: ok_prod_nucl_tritium    !--- TRUE => HTO production by nuclear tests
 !$OMP THREADPRIVATE(ok_prod_nucl_tritium)
-        integer nessai
-        parameter (nessai=486)
-        integer, save :: day_nucl(nessai)
-!$OMP THREADPRIVATE(day_nucl)
-        integer, save :: month_nucl(nessai)
-!$OMP THREADPRIVATE(month_nucl)
-        integer, save :: year_nucl(nessai)
-!$OMP THREADPRIVATE(year_nucl)
-        real, save :: lat_nucl(nessai)
-!$OMP THREADPRIVATE(lat_nucl)
-        real, save :: lon_nucl(nessai)
-!$OMP THREADPRIVATE(lon_nucl)
-        real, save :: zmin_nucl(nessai)
-!$OMP THREADPRIVATE(zmin_nucl)
-        real, save :: zmax_nucl(nessai)
-!$OMP THREADPRIVATE(zmax_nucl)
-        real, save :: HTO_nucl(nessai)
-!$OMP THREADPRIVATE(HTO_nucl)
-
+   INTEGER, PARAMETER :: nessai = 486
+   INTEGER, DIMENSION(nessai), SAVE :: &
+                    day_nucl, month_nucl, year_nucl
+!$OMP THREADPRIVATE(day_nucl, month_nucl, year_nucl)
+   REAL,    DIMENSION(nessai), SAVE :: &
+                    lat_nucl, lon_nucl, zmin_nucl, zmax_nucl, HTO_nucl
+!$OMP THREADPRIVATE(lat_nucl, lon_nucl, zmin_nucl, zmax_nucl, HTO_nucl)
+ 
  
 CONTAINS
 
-  SUBROUTINE iso_init()
-      use ioipsl_getin_p_mod, ONLY : getin_p
-      implicit none
-
-! -- local variables:
-
-      integer ixt
-      ! référence O18
-      real fac_enrichoce18
-      real alpha_liq_sol_O18, &
-     &     talph1_O18,talph2_O18,talph3_O18, &
-     &     talps1_O18,talps2_O18, &
-     &     tkcin0_O18,tkcin1_O18,tkcin2_O18, &
-     &     tdifrel_O18  
+SUBROUTINE iso_init()
+   USE infotrac_phy,       ONLY: ntiso, niso, getKey
+    USE strings_mod,       ONLY: maxlen
+   IMPLICIT NONE
+
+   !=== Local variables:
+   INTEGER :: ixt
+
+   !--- H2[18]O reference
+   REAL :: fac_enrichoce18, alpha_liq_sol_O18, &
+           talph1_O18, talph2_O18, talph3_O18, talps1_O18, talps2_O18, &
+           tkcin0_O18, tkcin1_O18, tkcin2_O18, tdifrel_O18  
+
+   !--- For H2[17]O
+   REAL    :: fac_kcin, pente_MWL
       
-      ! cas de l'O17
-      real fac_kcin
-      real pente_MWL
-      integer ierr
-      
-      logical ok_nocinsat, ok_nocinsfc !sensi test
-      parameter (ok_nocinsfc=.FALSE.)  ! if T: no kinetic effect in sfc evap
-      parameter (ok_nocinsat=.FALSE.)  ! if T: no sursaturation effect for ice
-      logical Rdefault_smow
-      parameter (Rdefault_smow=.FALSE.) ! si T: Rdefault=smow; si F: nul
-      ! pour le tritium
-      integer iessai
-
-    write(*,*) 'iso_init 219: entree'
-
-! allocations mémoire
-allocate (tnat(niso))
-allocate (toce(niso))
-allocate (tcorr(niso))
-allocate (tdifrel(niso))
-allocate (talph1(niso))
-allocate (talph2(niso))
-allocate (talph3(niso))
-allocate (talps1(niso))
-allocate (talps2(niso))
-allocate (tkcin0(niso))
-allocate (tkcin1(niso))
-allocate (tkcin2(niso))
-allocate (alpha_liq_sol(niso))
-allocate (Rdefault(niso))
-allocate (Rmethox(niso))
-allocate (striso(niso))
-
-
-!--------------------------------------------------------------
-! General:
-!--------------------------------------------------------------
-
-! -- verif du nombre d'isotopes:
-      write(*,*) 'iso_init 64: niso=',niso
-
-! init de ntracisoOR: on écrasera en cas de nzone>0 si complications avec
-! ORCHIDEE
-      ntracisoOR=ntraciso  
-              
-! -- Type of water isotopes:
-
-        iso_eau=indnum_fn_num(1)
-        iso_HDO=indnum_fn_num(2)
-        iso_O18=indnum_fn_num(3)
-        iso_O17=indnum_fn_num(4)
-        iso_HTO=indnum_fn_num(5)
-        write(*,*) 'iso_init 59: iso_eau=',iso_eau
-        write(*,*) 'iso_HDO=',iso_HDO
-        write(*,*) 'iso_O18=',iso_O18
-        write(*,*) 'iso_O17=',iso_O17
-        write(*,*) 'iso_HTO=',iso_HTO
-        write(*,*) 'iso_init 251: use_iso=',use_iso
-
-      ! initialisation
-        lambda_sursat=0.004
-        thumxt1=0.75*1.2
-        ntot=20
-        h_land_ice=20. ! à comparer aux 3000mm de snow_max
-        P_veg=1.0
-        bidouille_anti_divergence=.false.
-        essai_convergence=.false.
-        initialisation_iso=0
-        modif_sst=0
-        modif_sic=0
-        deltaTtest=0.0
-        deltasic=0.1
-        deltaTtestpoles=0.0
-        sstlatcrit=30.0
-        deltaO18_oce=0.0
-        albedo_prescrit=0
-        lon_min_albedo=-200
-        lon_max_albedo=200
-        lat_min_albedo=-100
-        lat_max_albedo=100
-        deltaP_BL=10.0
-        ruissellement_pluie=0
-        alphak_stewart=1
-        tdifexp_sol=0.67
-        calendrier_guide=0
-        cste_surf_cond=0
-mixlen=35.0        
-evap_cont_cste=0.0
-deltaO18_evap_cont=0.0
-d_evap_cont=0.0
-nudge_qsol=0
-region_nudge_qsol=1
-nlevmaxO17=50
-no_pce=0
-A_satlim=1.0
-ok_restrict_A_satlim=0
-!        slope_limiterxy=2.0
-!        slope_limiterz=2.0
-modif_ratqs=0
-Pcrit_ratqs=500.0
-ratqsbasnew=0.05
-
-fac_modif_evaoce=1.0
-ok_bidouille_wake=0
-cond_temp_env=.false.
-! si oui, la temperature de cond est celle de l'environnement,
-! pour eviter bugs quand temperature dans ascendances convs est
-! mal calculee
-ok_prod_nucl_tritium=.false.
-
-! lecture des paramètres isotopiques:
-! pour que ça marche en openMP, il faut utiliser getin_p. Car le getin ne peut
-! être appelé que par un thread à la fois, et ça pose tout un tas de problème,
-! d'où tout un tas de magouilles comme dans conf_phys_m. A terme, tout le monde
-! lira par getin_p.
-call getin_p('lambda',lambda_sursat)
-call getin_p('thumxt1',thumxt1)
-call getin_p('ntot',ntot)
-call getin_p('h_land_ice',h_land_ice)
-call getin_p('P_veg',P_veg)
-call getin_p('bidouille_anti_divergence',bidouille_anti_divergence)
-call getin_p('essai_convergence',essai_convergence)
-call getin_p('initialisation_iso',initialisation_iso)
-!if (nzone>0) then      
-!if (initialisation_iso.eq.0) then
-!  call getin_p('initialisation_isotrac',initialisation_isotrac)
-!endif !if (initialisation_iso.eq.0) then
-!endif !if (nzone>0)
-call getin_p('modif_sst',modif_sst)
-if (modif_sst.ge.1) then
-call getin_p('deltaTtest',deltaTtest)
-if (modif_sst.ge.2) then
-  call getin_p('deltaTtestpoles',deltaTtestpoles)
-  call getin_p('sstlatcrit',sstlatcrit)
+   !--- Sensitivity tests
+   LOGICAL, PARAMETER ::   ok_nocinsfc = .FALSE. ! if T: no kinetic effect in sfc evap
+   LOGICAL, PARAMETER ::   ok_nocinsat = .FALSE. ! if T: no sursaturation effect for ice
+   LOGICAL, PARAMETER :: Rdefault_smow = .FALSE. ! if T: Rdefault=smow; if F: nul
+
+   !--- For [3]H
+   INTEGER :: iessai
+
+   CHARACTER(LEN=maxlen) :: modname, sxt
+   REAL, ALLOCATABLE :: tmp(:)
+
+   modname = 'iso_init'
+   CALL msg('219: entree', modname)
+
+   !--------------------------------------------------------------
+   ! General:
+   !--------------------------------------------------------------
+
+   !--- Check number of isotopes
+   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
 #ifdef ISOVERIF
-      !call iso_verif_positif(sstlatcrit,'iso_init 107')
-      if (sstlatcrit.lt.0.0) then
-        write(*,*) 'iso_init 270: sstlatcrit=',sstlatcrit
-        stop
-      endif
+   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              
-  if (modif_sst.ge.3) then  
-      call getin_p('dsstlatcrit',dsstlatcrit)
+
+   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
+   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 iso_verif_positif(dsstlatcrit,'iso_init 110')
-      if (sstlatcrit.lt.0.0) then
-        write(*,*) 'iso_init 279: dsstlatcrit=',dsstlatcrit
-        stop
-      endif
-#endif              
-  endif !if (modif_sst.ge.3) then
-endif !if (modif_sst.ge.2) then
-endif !  if (modif_sst.ge.1) then   
-call getin_p('modif_sic',modif_sic)
-if (modif_sic.ge.1) then
-call getin_p('deltasic',deltasic)
-endif !if (modif_sic.ge.1) then
-
-call getin_p('albedo_prescrit',albedo_prescrit)
-call getin_p('lon_min_albedo',lon_min_albedo)
-call getin_p('lon_max_albedo',lon_max_albedo)
-call getin_p('lat_min_albedo',lat_min_albedo)
-call getin_p('lat_max_albedo',lat_max_albedo)
-call getin_p('deltaO18_oce',deltaO18_oce)
-call getin_p('deltaP_BL',deltaP_BL)
-call getin_p('ruissellement_pluie',ruissellement_pluie)
-call getin_p('alphak_stewart',alphak_stewart)
-call getin_p('tdifexp_sol',tdifexp_sol)
-call getin_p('calendrier_guide',calendrier_guide)
-call getin_p('cste_surf_cond',cste_surf_cond)
-call getin_p('mixlen',mixlen)
-call getin_p('evap_cont_cste',evap_cont_cste)
-call getin_p('deltaO18_evap_cont',deltaO18_evap_cont)
-call getin_p('d_evap_cont',d_evap_cont)  
-call getin_p('nudge_qsol',nudge_qsol)
-call getin_p('region_nudge_qsol',region_nudge_qsol)
-call getin_p('no_pce',no_pce)
-call getin_p('A_satlim',A_satlim)
-call getin_p('ok_restrict_A_satlim',ok_restrict_A_satlim)
-#ifdef ISOVERIF      
-!call iso_verif_positif(1.0-A_satlim,'iso_init 158')
-      if (A_satlim.gt.1.0) then
-        write(*,*) 'iso_init 315: A_satlim=',A_satlim
-        stop
-      endif
-#endif          
-!      call getin_p('slope_limiterxy',slope_limiterxy) 
-!      call getin_p('slope_limiterz',slope_limiterz) 
-call getin_p('modif_ratqs',modif_ratqs)
-call getin_p('Pcrit_ratqs',Pcrit_ratqs)
-call getin_p('ratqsbasnew',ratqsbasnew)
-call getin_p('fac_modif_evaoce',fac_modif_evaoce)
-call getin_p('ok_bidouille_wake',ok_bidouille_wake)
-call getin_p('cond_temp_env',cond_temp_env)
-if (use_iso(iso_HTO_possible)) then
-  ok_prod_nucl_tritium=.true.
-  call getin_p('ok_prod_nucl_tritium',ok_prod_nucl_tritium)
-endif
-
-write(*,*) 'lambda,thumxt1=',lambda_sursat,thumxt1
-write(*,*) 'bidouille_anti_divergence=',bidouille_anti_divergence
-write(*,*) 'essai_convergence=',essai_convergence
-write(*,*) 'initialisation_iso=',initialisation_iso
-write(*,*) 'modif_sst=',modif_sst
-if (modif_sst.ge.1) then
-write(*,*) 'deltaTtest=',deltaTtest
-if (modif_sst.ge.2) then  
-write(*,*) 'deltaTtestpoles,sstlatcrit=', &
-&           deltaTtestpoles,sstlatcrit
-if (modif_sst.ge.3) then   
- write(*,*) 'dsstlatcrit=',dsstlatcrit
-endif !if (modif_sst.ge.3) then
-endif !if (modif_sst.ge.2) then
-endif !if (modif_sst.ge.1) then
-write(*,*) 'modif_sic=',modif_sic
-if (modif_sic.ge.1) then  
-write(*,*) 'deltasic=',deltasic
-endif !if (modif_sic.ge.1) then
-write(*,*) 'deltaO18_oce=',deltaO18_oce
-write(*,*) 'albedo_prescrit=',albedo_prescrit
-if (albedo_prescrit.eq.1) then
- write(*,*) 'lon_min_albedo,lon_max_albedo=', &
-&           lon_min_albedo,lon_max_albedo
- write(*,*) 'lat_min_albedo,lat_max_albedo=', &
-&           lat_min_albedo,lat_max_albedo
-endif !if (albedo_prescrit.eq.1) then
-write(*,*) 'deltaP_BL,ruissellement_pluie,alphak_stewart=', &
-&       deltaP_BL,ruissellement_pluie,alphak_stewart
-write(*,*) 'cste_surf_cond=',cste_surf_cond
-write(*,*) 'mixlen=',mixlen
-write(*,*) 'tdifexp_sol=',tdifexp_sol
-write(*,*) 'calendrier_guide=',calendrier_guide
-write(*,*) 'evap_cont_cste=',evap_cont_cste
-write(*,*) 'deltaO18_evap_cont,d_evap_cont=', &
-&           deltaO18_evap_cont,d_evap_cont
-write(*,*) 'nudge_qsol,region_nudge_qsol=', &
-&  nudge_qsol,region_nudge_qsol  
-write(*,*) 'nlevmaxO17=',nlevmaxO17
-write(*,*) 'no_pce=',no_pce
-write(*,*) 'A_satlim=',A_satlim
-write(*,*) 'ok_restrict_A_satlim=',ok_restrict_A_satlim
-!      write(*,*) 'slope_limiterxy=',slope_limiterxy
-!      write(*,*) 'slope_limiterz=',slope_limiterz
-write(*,*) 'modif_ratqs=',modif_ratqs
-write(*,*) 'Pcrit_ratqs=',Pcrit_ratqs
-write(*,*) 'ratqsbasnew=',ratqsbasnew
-write(*,*) 'fac_modif_evaoce=',fac_modif_evaoce
-write(*,*) 'ok_bidouille_wake=',ok_bidouille_wake
-write(*,*) 'cond_temp_env=',cond_temp_env
-write(*,*) 'ok_prod_nucl_tritium=',ok_prod_nucl_tritium
-          
-
-!--------------------------------------------------------------
-! Parameters that do not depend on the nature of water isotopes:
-!--------------------------------------------------------------
-
-! -- temperature at which ice condensate starts to form (valeur ECHAM?):
-pxtmelt=273.15
-!      pxtmelt=273.15-10.0 ! test PHASE
-
-! -- temperature at which all condensate is ice:
-pxtice=273.15-10.0
-!      pxtice=273.15-30.0 ! test PHASE
-
-! -- minimum temperature to calculate fractionation coeff
-pxtmin=273.15-120.0 ! On ne calcule qu'au dessus de -120°C
-pxtmax=273.15+60.0 ! On ne calcule qu'au dessus de +60°C
-! remarque: les coeffs ont été mesurés seulement jusq'à -40!
-
-! -- a constant for alpha_eff for equilibrium below cloud base:
-tdifexp=0.58
-tv0cin=7.0
-
-! facteurs lambda et mu dans Si=musi-lambda*T
-musi=1.0
-if (ok_nocinsat) then
-lambda_sursat = 0.0 ! no sursaturation effect
-endif            
-
-
-! diffusion dans le sol
-Kd=2.5e-9 ! m2/s   
-
-! cas où cste_surf_cond: on met rhs ou/et Ts cste pour voir
-rh_cste_surf_cond=0.6
-T_cste_surf_cond=288.0
-
-!--------------------------------------------------------------
-! Parameters that depend on the nature of water isotopes:
-!--------------------------------------------------------------
-! ** constantes locales
-fac_enrichoce18=0.0005
-! on a alors tcor018=1+fac_enrichoce18
-! tcorD=1+fac_enrichoce18*8
-! tcorO17=1+fac_enrichoce18*0.528
-alpha_liq_sol_O18=1.00291 ! valeur de Lehmann & Siegenthaler, 1991, 
-  ! Journal of Glaciology, vol 37, p 23
-talph1_O18=1137.
-talph2_O18=-0.4156
-talph3_O18=-2.0667E-3
-talps1_O18=11.839
-talps2_O18=-0.028244
-tkcin0_O18 = 0.006
-tkcin1_O18 = 0.000285
-tkcin2_O18 = 0.00082
-tdifrel_O18= 1./0.9723
-
-! rapport des ln(alphaeq) entre O18 et O17
-fac_coeff_eq17_liq=0.529 ! donné par Amaelle
-!      fac_coeff_eq17_ice=0.528 ! slope MWL 
-fac_coeff_eq17_ice=0.529
-
-
-write(*,*) 'iso_O18,iso_HDO,iso_eau=',iso_O18,iso_HDO,iso_eau
-do 999 ixt = 1, niso
-write(*,*) 'iso_init 80: ixt=',ixt
-
-
-! -- kinetic factor for surface evaporation:
-! (cf: kcin = tkcin0                  if |V|<tv0cin
-!      kcin = tkcin1*|Vsurf| + tkcin2 if |V|>tv0cin )
-! (Rq: formula discontinuous for |V|=tv0cin... )        
-
-! -- main:
-if (ixt.eq.iso_HTO) then ! Tritium
-  tkcin0(ixt) = 0.01056
-  tkcin1(ixt) = 0.0005016
-  tkcin2(ixt) = 0.0014432
-  tnat(ixt)=0.
-  !toce(ixt)=2.2222E-8 ! corrigé par Alex Cauquoin
-  !toce(ixt)=1.0E-18 ! rapport 3H/1H ocean
-  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. 
-  Rdefault(ixt)=0.0
-  Rmethox(ixt)=0.0
-  striso(ixt)='HTO'
-endif     
-if (ixt.eq.iso_O17) then ! Deuterium
-  pente_MWL=0.528
-!          tdifrel(ixt)=1./0.985452 ! donné par Amaelle 
-  tdifrel(ixt)=1./0.98555 ! valeur utilisée en 1D et dans modèle de LdG
-!          fac_kcin=0.5145 ! donné par Amaelle
-  fac_kcin= (tdifrel(ixt)-1.0)/(tdifrel_O18-1.0)
-  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
-  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 
-  if (Rdefault_smow) then    
-        Rdefault(ixt)=tnat(ixt)*(-3.15/1000.0+1.0)
-  else
-        Rdefault(ixt)=0.0
-  endif
-  Rmethox(ixt)=(230./1000.+1.)*tnat(ixt) !Zahn et al 2006
-  striso(ixt)='O17'
-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
-  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   
-  if (Rdefault_smow) then    
-        Rdefault(ixt)=tnat(ixt)*(-6.0/1000.0+1.0)
-  else
-        Rdefault(ixt)=0.0
-  endif
-  Rmethox(ixt)=(130./1000.+1.)*tnat(ixt) !Zahn et al 2006   
-!       write(*,*) 'iso_init 163: ZXalpha_liq_sol=',ZXalpha_liq_sol
-  striso(ixt)='O18'
-  write(*,*) 'isotopes_mod 519: ixt,striso(ixt)=',ixt,striso(ixt)
-endif
-
-if (ixt.eq.iso_HDO) then ! Deuterium
-  pente_MWL=8.0
-!          fac_kcin=0.88
-  tdifrel(ixt)=1./0.9755
-  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
-  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
-  if (Rdefault_smow) then    
-    Rdefault(ixt)=tnat(ixt)*((-6.0*pente_MWL+10.0)/1000.0+1.0)
-  else
-    Rdefault(ixt)=0.0
-  endif
-  Rmethox(ixt)=tnat(ixt)*(-25.0/1000.+1.) ! Zahn et al 2006
-  striso(ixt)='HDO'
-  write(*,*) 'isotopes_mod 548: ixt,striso(ixt)=',ixt,striso(ixt)
-endif
-
-!       write(*,*) 'iso_init 163: ZXalpha_liq_sol=',ZXalpha_liq_sol
-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.
-  if (Rdefault_smow) then
-        Rdefault(ixt)=tnat(ixt)*1.0
-  else
-        Rdefault(ixt)=1.0
-  endif
-  Rmethox(ixt)=1.0
-  striso(ixt)='eau'
-endif
-
-999   continue
-
-! test de sensibilité: 
-if (ok_nocinsfc) then ! no kinetic effect in sfc evaporation
- do ixt=1,niso
-  tkcin0(ixt) = 0.0 
-  tkcin1(ixt) = 0.0 
-  tkcin2(ixt) = 0.0
- enddo 
-endif
-
-! nom des isotopes
-
-! verif
-write(*,*) 'iso_init 285: verif initialisation:'
-
-do ixt=1,niso
-  write(*,*) '* striso(',ixt,')=<'//striso(ixt)//'>'
-  write(*,*) 'tnat(',ixt,')=',tnat(ixt)
-!          write(*,*) 'alpha_liq_sol(',ixt,')=',alpha_liq_sol(ixt)
-!          write(*,*) 'tkcin0(',ixt,')=',tkcin0(ixt)
-!          write(*,*) 'tdifrel(',ixt,')=',tdifrel(ixt)
-enddo
-write(*,*) 'iso_init 69: lambda=',lambda_sursat
-write(*,*) 'iso_init 69: thumxt1=',thumxt1
-write(*,*) 'iso_init 69: h_land_ice=',h_land_ice
-write(*,*) 'iso_init 69: P_veg=',P_veg   
-
-    return
+   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.)
+
+   !--------------------------------------------------------------
+   ! Parameters that do not depend on the nature of water isotopes:
+   !--------------------------------------------------------------
+   ! -- temperature at which ice condensate starts to form (valeur ECHAM?):
+   pxtmelt = 273.15
+
+   ! -- temperature at which all condensate is ice:
+   pxtice  = 273.15-10.0
+
+   !- -- test PHASE
+!   pxtmelt = 273.15 - 10.0
+!   pxtice  = 273.15 - 30.0
+
+   ! -- minimum temperature to calculate fractionation coeff
+   pxtmin = 273.15 - 120.0   ! On ne calcule qu'au dessus de -120°C
+   pxtmax = 273.15 +  60.0   ! On ne calcule qu'au dessus de +60°C
+   !    Remarque: les coeffs ont ete mesures seulement jusq'à -40!
+
+   ! -- a constant for alpha_eff for equilibrium below cloud base:
+   tdifexp = 0.58
+   tv0cin  = 7.0
+
+   ! facteurs lambda et mu dans Si=musi-lambda*T
+   musi=1.0
+   if (ok_nocinsat) lambda_sursat = 0.0          ! no sursaturation effect
+
+   ! diffusion dans le sol
+   Kd=2.5e-9 ! m2/s   
+
+   ! cas où cste_surf_cond: on met rhs ou/et Ts cste pour voir
+   rh_cste_surf_cond = 0.6
+    T_cste_surf_cond = 288.0
+   
+   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:
+   !--------------------------------------------------------------
+   IF(getKey('tnat',    tnat,    isoName)) CALL abort_physic(modname, 'can''t get tnat',    1)
+   IF(getKey('toce',    toce,    isoName)) CALL abort_physic(modname, 'can''t get toce',    1)
+   IF(getKey('tcorr',   tcorr,   isoName)) CALL abort_physic(modname, 'can''t get tcorr',   1)
+   IF(getKey('talph1',  talph1,  isoName)) CALL abort_physic(modname, 'can''t get talph1',  1)
+   IF(getKey('talph2',  talph2,  isoName)) CALL abort_physic(modname, 'can''t get talph2',  1)
+   IF(getKey('talph3',  talph3,  isoName)) CALL abort_physic(modname, 'can''t get talph3',  1)
+   IF(getKey('talps1',  talps1,  isoName)) CALL abort_physic(modname, 'can''t get talps1',  1)
+   IF(getKey('talps2',  talps2,  isoName)) CALL abort_physic(modname, 'can''t get talps2',  1)
+   IF(getKey('tkcin0',  tkcin0,  isoName)) CALL abort_physic(modname, 'can''t get tkcin0',  1)
+   IF(getKey('tkcin1',  tkcin1,  isoName)) CALL abort_physic(modname, 'can''t get tkcin1',  1)
+   IF(getKey('tkcin2',  tkcin2,  isoName)) CALL abort_physic(modname, 'can''t get tkcin2',  1)
+   IF(getKey('tdifrel', tdifrel, isoName)) CALL abort_physic(modname, 'can''t get tdifrel', 1)
+   IF(getKey('alpha_liq_sol', alpha_liq_sol, isoName)) CALL abort_physic(modname, 'can''t get alpha_liq_sol',  1)
+   IF(getKey('Rdefault',Rdefault,isoName)) CALL abort_physic(modname, 'can''t get Rdefault',1)
+   IF(getKey('Rmethox', Rmethox, isoName)) CALL abort_physic(modname, 'can''t get Rmethox', 1)
+   IF(.NOT.Rdefault_smow) Rdefault(:) = 0.0
+
+   !--- 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
 
+
+SUBROUTINE getinp_s(nam, val, def, lDisp)
+   USE ioipsl_getincom, ONLY: getin
+   USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
+   USE mod_phys_lmdz_omp_data, ONLY :  is_omp_root
+   USE mod_phys_lmdz_transfert_para, ONLY : bcast
+   CHARACTER(LEN=*),           INTENT(IN)    :: nam
+   CHARACTER(LEN=*),           INTENT(INOUT) :: val
+   CHARACTER(LEN=*), OPTIONAL, INTENT(IN)    :: def
+   LOGICAL,          OPTIONAL, INTENT(IN)    :: lDisp
+   LOGICAL :: lD
+!$OMP BARRIER
+   IF(is_mpi_root.AND.is_omp_root) THEN
+      IF(PRESENT(def)) val=def; CALL getin(nam,val)
+      lD=.TRUE.; IF(PRESENT(lDisp)) lD=lDisp
+      IF(lD) CALL msg(TRIM(nam)//' = '//TRIM(val))
+  END IF
+  CALL bcast(val)
+END SUBROUTINE getinp_s
+
+SUBROUTINE getinp_i(nam, val, def, lDisp)
+   USE ioipsl_getincom, ONLY: getin
+   USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
+   USE mod_phys_lmdz_omp_data, ONLY :  is_omp_root
+   USE mod_phys_lmdz_transfert_para, ONLY : bcast
+   CHARACTER(LEN=*),  INTENT(IN)    :: nam
+   INTEGER,           INTENT(INOUT) :: val
+   INTEGER, OPTIONAL, INTENT(IN)    :: def
+   LOGICAL, OPTIONAL, INTENT(IN)    :: lDisp
+   LOGICAL :: lD
+!$OMP BARRIER
+   IF(is_mpi_root.AND.is_omp_root) THEN
+      IF(PRESENT(def)) val=def; CALL getin(nam,val)
+      lD=.TRUE.; IF(PRESENT(lDisp)) lD=lDisp
+      IF(lD) CALL msg(TRIM(nam)//' = '//TRIM(int2str(val)))
+  END IF
+  CALL bcast(val)
+END SUBROUTINE getinp_i
+
+SUBROUTINE getinp_r(nam, val, def, lDisp)
+   USE ioipsl_getincom, ONLY: getin
+   USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
+   USE mod_phys_lmdz_omp_data, ONLY :  is_omp_root
+   USE mod_phys_lmdz_transfert_para, ONLY : bcast
+   CHARACTER(LEN=*),  INTENT(IN)    :: nam
+   REAL,              INTENT(INOUT) :: val
+   REAL,    OPTIONAL, INTENT(IN)    :: def
+   LOGICAL, OPTIONAL, INTENT(IN)    :: lDisp
+   LOGICAL :: lD
+!$OMP BARRIER
+   IF(is_mpi_root.AND.is_omp_root) THEN
+      IF(PRESENT(def)) val=def; CALL getin(nam,val)
+      lD=.TRUE.; IF(PRESENT(lDisp)) lD=lDisp
+      IF(lD) CALL msg(TRIM(nam)//' = '//TRIM(real2str(val)))
+  END IF
+  CALL bcast(val)
+END SUBROUTINE getinp_r
+
+SUBROUTINE getinp_l(nam, val, def, lDisp)
+   USE ioipsl_getincom, ONLY: getin
+   USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
+   USE mod_phys_lmdz_omp_data, ONLY :  is_omp_root
+   USE mod_phys_lmdz_transfert_para, ONLY : bcast
+   CHARACTER(LEN=*),  INTENT(IN)    :: nam
+   LOGICAL,           INTENT(INOUT) :: val
+   LOGICAL, OPTIONAL, INTENT(IN)    :: def
+   LOGICAL, OPTIONAL, INTENT(IN)    :: lDisp
+   LOGICAL :: lD
+!$OMP BARRIER
+   IF(is_mpi_root.AND.is_omp_root) THEN
+      IF(PRESENT(def)) val=def; CALL getin(nam,val)
+      lD=.TRUE.; IF(PRESENT(lDisp)) lD=lDisp
+      IF(lD) CALL msg(TRIM(nam)//' = '//TRIM(bool2str(val)))
+  END IF
+  CALL bcast(val)
+END SUBROUTINE getinp_l
 
 END MODULE isotopes_mod
