Index: LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/atke_exchange_coeff_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/atke_exchange_coeff_mod.F90	(revision 4482)
+++ LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/atke_exchange_coeff_mod.F90	(revision 4482)
@@ -0,0 +1,1 @@
+link ../phylmd/atke_exchange_coeff_mod.F90
Index: LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/atke_turbulence_ini_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/atke_turbulence_ini_mod.F90	(revision 4482)
+++ LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/atke_turbulence_ini_mod.F90	(revision 4482)
@@ -0,0 +1,1 @@
+link ../phylmd/atke_turbulence_ini_mod.F90
Index: LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/fisrtilp.F90
===================================================================
--- LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/fisrtilp.F90	(revision 4175)
+++ LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/fisrtilp.F90	(revision 4482)
@@ -143,7 +143,7 @@
   ! Coeffients de fraction lessivee : pour OFF-LINE
   !
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfrac_nucl
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfrac_1nucl
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfrac_impa
+  REAL, DIMENSION(klon,klev),      INTENT(INOUT)  :: pfrac_nucl
+  REAL, DIMENSION(klon,klev),      INTENT(INOUT)  :: pfrac_1nucl
+  REAL, DIMENSION(klon,klev),      INTENT(INOUT)  :: pfrac_impa
   !
   ! Fraction d'aerosols lessivee par impaction et par nucleation
Index: LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/isotopes_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/isotopes_mod.F90	(revision 4175)
+++ LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/isotopes_mod.F90	(revision 4482)
@@ -133,5 +133,4 @@
 
 SUBROUTINE iso_init()
-   USE ioipsl_getin_p_mod, ONLY: getin_p
    USE infotrac_phy,       ONLY: ntiso, niso, getKey
     USE strings_mod,       ONLY: maxlen
@@ -175,16 +174,11 @@
 
    !--- Type of water isotopes:
-   iso_eau = strIdx(isoName, 'H2[16]O'); CALL msg('iso_eau='//int2str(iso_eau), modname)
-   iso_HDO = strIdx(isoName, 'H[2]HO');  CALL msg('iso_HDO='//int2str(iso_HDO), modname)
-   iso_O18 = strIdx(isoName, 'H2[18]O'); CALL msg('iso_O18='//int2str(iso_O18), modname)
-   iso_O17 = strIdx(isoName, 'H2[17]O'); CALL msg('iso_O17='//int2str(iso_O17), modname)
-   iso_HTO = strIdx(isoName, 'H[3]HO');  CALL msg('iso_HTO='//int2str(iso_HTO), modname)
-
-   ! initialisation
-   ! lecture des parametres isotopiques:
-   ! pour que ca marche en openMP, il faut utiliser getin_p. Car le getin ne peut
-   ! etre appele que par un thread a la fois, et ca pose tout un tas de problemes,
-   ! d'ou tout un tas de magouilles comme dans conf_phys_m. A terme, tout le monde
-   ! lira par getin_p.
+   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)
@@ -252,5 +246,5 @@
    ! bugs quand temperature dans ascendances convs est mal calculee
    CALL get_in('cond_temp_env',        cond_temp_env,        .FALSE.)
-   IF(ANY(isoName == 'H[3]HO')) &
+   IF(ANY(isoName == 'HTO')) &
    CALL get_in('ok_prod_nucl_tritium', ok_prod_nucl_tritium, .FALSE., .FALSE.)
 
@@ -339,14 +333,16 @@
    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=*),  INTENT(IN)    :: def
-   LOGICAL, OPTIONAL, INTENT(IN)    :: lDisp
+   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(.NOT.(is_mpi_root.AND.is_omp_root)) RETURN
-   val=def; CALL getin(nam,val); CALL bcast(val)
-   lD=.TRUE.; IF(PRESENT(lDisp)) lD=lDisp
-   IF(lD) CALL msg(TRIM(nam)//' = '//TRIM(val))
+   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
 
@@ -358,12 +354,14 @@
    CHARACTER(LEN=*),  INTENT(IN)    :: nam
    INTEGER,           INTENT(INOUT) :: val
-   INTEGER,           INTENT(IN)    :: def
+   INTEGER, OPTIONAL, INTENT(IN)    :: def
    LOGICAL, OPTIONAL, INTENT(IN)    :: lDisp
    LOGICAL :: lD
 !$OMP BARRIER
-   IF(.NOT.(is_mpi_root.AND.is_omp_root)) RETURN
-   val=def; CALL getin(nam,val); CALL bcast(val)
-   lD=.TRUE.; IF(PRESENT(lDisp)) lD=lDisp
-   IF(lD) CALL msg(TRIM(nam)//' = '//TRIM(int2str(val)))
+   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
 
@@ -375,12 +373,14 @@
    CHARACTER(LEN=*),  INTENT(IN)    :: nam
    REAL,              INTENT(INOUT) :: val
-   REAL,              INTENT(IN)    :: def
+   REAL,    OPTIONAL, INTENT(IN)    :: def
    LOGICAL, OPTIONAL, INTENT(IN)    :: lDisp
    LOGICAL :: lD
 !$OMP BARRIER
-   IF(.NOT.(is_mpi_root.AND.is_omp_root)) RETURN
-   val=def; CALL getin(nam,val); CALL bcast(val)
-   lD=.TRUE.; IF(PRESENT(lDisp)) lD=lDisp
-   IF(lD) CALL msg(TRIM(nam)//' = '//TRIM(real2str(val)))
+   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
 
@@ -392,12 +392,14 @@
    CHARACTER(LEN=*),  INTENT(IN)    :: nam
    LOGICAL,           INTENT(INOUT) :: val
-   LOGICAL,           INTENT(IN)    :: def
+   LOGICAL, OPTIONAL, INTENT(IN)    :: def
    LOGICAL, OPTIONAL, INTENT(IN)    :: lDisp
    LOGICAL :: lD
 !$OMP BARRIER
-   IF(.NOT.(is_mpi_root.AND.is_omp_root)) RETURN
-   val=def; CALL getin(nam,val); CALL bcast(val)
-   lD=.TRUE.; IF(PRESENT(lDisp)) lD=lDisp
-   IF(lD) CALL msg(TRIM(nam)//' = '//TRIM(bool2str(val)))
+   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
 
Index: LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/isotopes_routines_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/isotopes_routines_mod.F90	(revision 4175)
+++ LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/isotopes_routines_mod.F90	(revision 4482)
@@ -1731,7 +1731,7 @@
        else
            write(*,*) 'iso_fractcal 1734: non prévu: ixt=',ixt
-#ifdef ISOVERIF
-           stop
-#endif
+!#ifdef ISOVERIF
+           CALL abort_physic('isotopes_routines_mod', 'iso_fractcal 1734', 1)
+!#endif
        endif
 
@@ -1949,7 +1949,7 @@
        else
            write(*,*) 'iso_fractcal 1954: non prévu: ixt=',ixt
-#ifdef ISOVERIF
-           stop
-#endif
+!#ifdef ISOVERIF
+          CALL abort_physic('isotopes_routines_mod', 'iso_fractcal 1954', 1)
+!#endif
        endif
 
@@ -2185,5 +2185,5 @@
         write(*,*) 'qsol(i)=',qsol(i),' mais evap(i)=',evap(i)
         write(*,*) 'q1lay(i)=',q1lay(i)
-        stop
+        CALL abort_physic('isotopes_routines_mod', 'iso_surf 2187', 1)
  endif !if (q1lay(i).gt.ridicule) then
 endif !if (qsol(i).gt.ridicule_rain) then
@@ -2414,5 +2414,5 @@
     write(*,*) 'iso_surf>iso_rosee_givre 3189: evap=',evap(i)
     write(*,*) 'q1lay(i)=',q1lay(i)
-    stop
+    CALL abort_physic('isotopes_routines_mod', 'iso_surf 2416', 1)
 endif !if (q1lay(i).gt.0) then
 
@@ -2552,5 +2552,5 @@
   write(*,*) 'il,Pqisup(il),Eqi(il),Pqiinf(il)=', &
 &          il,Pqisup(il),Eqi(il),Pqiinf(il) 
-  stop
+  CALL abort_physic('isotopes_routines_mod', 'stewart 2554', 1)
 endif !if (iso_verif_egalite
 enddo !do il=1,ncas 
@@ -16508,246 +16508,162 @@
       end subroutine phyiso_etat0_dur
 
-      subroutine phyiso_etat0_fichier( &
-     &           snow,run_off_lic_0, &
-     &           xtsnow,xtrun_off_lic_0, &
-     &           Rland_ice)
-      USE dimphy, only: klon,klev
-      !USE mod_grid_phy_lmdz
-      !USE mod_phys_lmdz_para
-      USE iophy
-      USE phys_state_var_mod, ONLY: q_ancien,xt_ancien,wake_deltaq,wake_deltaxt, &
-#ifdef ISOVERIF
-        rain_fall,snow_fall,fevap,qsol, &
-#endif
-        xtrain_fall,xtsnow_fall,ql_ancien,xtl_ancien,qs_ancien,xts_ancien, &
-        fxtevap,xtsol
-      !USE iostart
-      !USE write_field_phy
-      USE indice_sol_mod, only: nbsrf  
-  USE isotopes_mod, ONLY: isoName,iso_HDO,iso_eau
-#ifdef ISOVERIF
-  USE isotopes_verif_mod
+SUBROUTINE phyiso_etat0_fichier(snow, run_off_lic_0, xtsnow, xtrun_off_lic_0, Rland_ice)
+   USE dimphy,             ONLY: klon,klev
+   USE iophy
+   USE phys_state_var_mod, ONLY: q_ancien, xt_ancien, wake_deltaq, wake_deltaxt, &
+#ifdef ISOVERIF
+     rain_fall, snow_fall, fevap,qsol, &
+#endif
+     xtrain_fall, xtsnow_fall, ql_ancien, xtl_ancien, qs_ancien, xts_ancien, fxtevap, xtsol
+   USE indice_sol_mod,    ONLY: nbsrf
+   USE isotopes_mod,      ONLY: isoName,iso_HDO,iso_eau
+   USE phyetat0_get_mod,  ONLY: phyetat0_get, phyetat0_srf
+   USE readTracFiles_mod, ONLY: new2oldH2O
+   USE strings_mod,       ONLY: strIdx, strTail, maxlen, msg, int2str
+#ifdef ISOVERIF
+   USE isotopes_verif_mod
 #endif
 #ifdef ISOTRAC
- USE isotrac_mod, ONLY: strtrac,initialisation_isotrac,index_iso, &
-&       index_zone,izone_init
-#endif
-        implicit none
+   USE isotrac_mod, ONLY: strtrac, initialisation_isotrac, index_iso, index_zone, izone_init
+#endif
+   IMPLICIT NONE
 
 #include "netcdf.inc"
 #include "dimsoil.h"
 #include "clesphys.h"
-! #include "thermcell.h"
 #include "compbl.h"    
 
-        ! inputs 
-        !REAL qsol(klon)
-        REAL snow(klon,nbsrf)
-        !REAL evap(klon,nbsrf)
-        REAL run_off_lic_0(klon)
-        ! outputs   
-        !REAL xtsol(niso,klon)
-        REAL xtsnow(niso,klon,nbsrf) 
-        !REAL xtevap(ntraciso,klon,nbsrf)     
-        REAL xtrun_off_lic_0(niso,klon)
-        REAL Rland_ice(niso,klon)
-
-        ! locals
-        real iso_tmp(klon)
-        real iso_tmp_lonlev(klon,klev)
-        real iso_tmp_lonsrf(klon,nbsrf)
-        INTEGER ierr
-        integer i,ixt,k,nsrf
-        INTEGER nid, nvarid
-        CHARACTER*2 str2
-        CHARACTER*5 str5
-        real xmin,xmax    
-        CHARACTER*50 outiso  
-        integer lnblnk
-        LOGICAL :: found,phyetat0_get,phyetat0_srf
-
-!#ifdef ISOVERIF
-!      integer iso_verif_egalite_nostop
-!#endif
-!#ifdef ISOVERIF
-!        real deltaD
-!        integer iso_verif_noNaN_nostop
-!#endif
+   REAL, INTENT(IN) ::             snow     (klon,nbsrf)
+   REAL, INTENT(IN) ::    run_off_lic_0     (klon)
+   REAL, INTENT(OUT) ::          xtsnow(niso,klon,nbsrf) 
+   REAL, INTENT(OUT) :: xtrun_off_lic_0(niso,klon)
+   REAL, INTENT(OUT) ::       Rland_ice(niso,klon)
+
+   INTEGER :: ierr, i, ixt, k, nsrf, nid, nvarid, lnblnk
+   CHARACTER(LEN=2) :: str2
+   CHARACTER(LEN=5) :: str5
+   CHARACTER(LEN=maxlen) :: outiso, oldIso, modname, nam(2)
+   REAL :: xmin, xmax
+   LOGICAL :: found
 #ifdef ISOTRAC
-        integer iiso,izone
-#endif
-
-
-   write(*,*) 'phyiso_etat0_fichier 3'
-   write(*,*) 'niso=',niso
-   write(*,*) 'isoName(1)='//TRIM(isoName(1))
-
-   do ixt=1,ntraciso
-
-      outiso=TRIM(isoName(ixt))
-      i = INDEX(outiso, '_', .TRUE.)
-      outiso = outiso(1:i-1)//outiso(i+1:LEN_TRIM(outiso))
-      write(*,*) 'phyiso_etat0_fichier 16621: ixt,outiso=',ixt,TRIM(outiso)
-
-           
-      ! on lit seulement si ixt<=niso ou si on initialise les traceurs d'après
-      ! fichier:
+   INTEGER :: iiso, izone
+#endif
+
+   modname = 'phyiso_etat0_fichier'
+   CALL msg('3', modname)
+   CALL msg('niso = '//TRIM(int2str(niso)), modname)
+   CALL msg('isoName(1) = '//TRIM(isoName(1)), modname)
+
+   DO ixt = 1, ntraciso
+
+      outiso = isoName(ixt)
+      oldIso = strTail(new2oldH2O(outiso), '_')            !--- Remove "H2O_" from "H2O_<iso>[_<tag>]"
+      ! on lit seulement si ixt<=niso ou si on initialise les traceurs d'après fichier:
 #ifdef ISOTRAC
-      if ((ixt.le.niso).or.(initialisation_isotrac.eq.0)) then
-#endif
-
-      found=phyetat0_srf(1,iso_tmp_lonsrf,"XTSNOW"//TRIM(outiso),"Surface snow",0.)
-      if (.NOT.found) CALL abort_physic('isotopes_routines_mod', &
-                            'phyiso_etat0_fichier 16581: variable isotopique not found',1)
-      xtsnow(ixt,:,:)=iso_tmp_lonsrf(:,:)
-      
-      found=phyetat0_srf(1,iso_tmp_lonsrf,"XTEVAP"//TRIM(outiso),"evaporation",0.)
-      fxtevap(ixt,:,:)=iso_tmp_lonsrf(:,:)
-
-      found=phyetat0_get(1,iso_tmp,"xtrain_f"//TRIM(outiso),"xrain fall",0.)
-      xtrain_fall(ixt,:)=iso_tmp(:)
-
-      found=phyetat0_get(1,iso_tmp,"xtsnow_f"//TRIM(outiso),"snow fall",0.)
-      xtsnow_fall(ixt,:)=iso_tmp(:)
-
-      found=phyetat0_get(klev,iso_tmp_lonlev,"XTANCIEN"//TRIM(outiso),"QANCIEN",0.)
-      xt_ancien(ixt,:,:)=iso_tmp_lonlev(:,:)
-
-      found=phyetat0_get(klev,iso_tmp_lonlev,"XTLANCIEN"//TRIM(outiso),"QLANCIEN",0.)
-      xtl_ancien(ixt,:,:)=iso_tmp_lonlev(:,:)
-
-      found=phyetat0_get(klev,iso_tmp_lonlev,"XTSANCIEN"//TRIM(outiso),"QSANCIEN",0.)
-      xts_ancien(ixt,:,:)=iso_tmp_lonlev(:,:)
-
-      found=phyetat0_get(1,iso_tmp,"XTRUNOFFLIC0"//TRIM(outiso),"RUNOFFLIC0",0.)  
-      xtrun_off_lic_0(ixt,:)=iso_tmp(:)
-
-      found=phyetat0_get(klev,iso_tmp_lonlev,"WAKE_DELTAXT"//TRIM(outiso),"Delta hum. wake/env",0.)  
-      wake_deltaxt(ixt,:,:)=iso_tmp_lonlev(:,:)
-
-#ifdef ISOVERIF           
-      if ((ixt.eq.iso_eau).and.(iso_eau.gt.0)) then
-        do i=1,klon
-         call iso_verif_egalite(xtrain_fall(iso_eau,i),rain_fall(i), &
-     &           'phyisoetat0_fichier 231a')
-         call iso_verif_egalite(xtsnow_fall(iso_eau,i),snow_fall(i), &
-     &           'phyisoetat0_fichier 231b')
-         DO nsrf = 1, nbsrf
-         call iso_verif_egalite(fxtevap(iso_eau,i,nsrf),fevap(i,nsrf), &
-     &           'phyisoetat0_fichier 231c')
-         call iso_verif_egalite(xtsnow(iso_eau,i,nsrf),snow(i,nsrf), &
-     &           'phyisoetat0_fichier 231d')
-         enddo !DO nsrf = 1, nbsrf
-        enddo !do i=1,klon        
-      endif !if (iso_eau.gt.0) then  
-        if ((iso_HDO.gt.0).and.(ixt.eq.iso_HDO)) then
-              do k=1,klev
-               do i=1,klon
-                if (q_ancien(i,k).gt.2e-3) then
-                call iso_verif_aberrant(xt_ancien(iso_hdo,i,k) &
-     &           /q_ancien(i,k),'phyisoetat0_fichier 312')
-                endif !if (q_ancien(i,k).gt.2e-3) then
-               enddo !do i=1,klon
-              enddo !do k=1,klev
-      endif !if ((iso_HDO.gt.0).and.(ixt.eq.iso_HDO)) then 
-      if (iso_eau.gt.0) then   
-        do i=1,klon
-          if (iso_verif_egalite_nostop(run_off_lic_0(i), &
-     &           xtrun_off_lic_0(iso_eau,i), &
-     &          'phyiso_etat0_fichier 326').eq.1) then
-            write(*,*) 'i=',i
-            stop
-          endif !if (iso_verif_egalite_nostop(run_off_lic_0(i),
-        enddo !do i=1,klon
-      endif !if (iso_eau.gt.0) then 
-#endif
-
-       ! ces variables n'ont pas de traceurs:
-       if (ixt.le.niso) then
-        found=phyetat0_get(1,iso_tmp,"XTSOL"//TRIM(outiso),"Surface hmidity / bucket",0.)  
-        xtsol(ixt,:)=iso_tmp(:)
-
-        found=phyetat0_get(1,iso_tmp,"Rland_ice"//TRIM(outiso),"R land ice",0.)
-        Rland_ice(ixt,:)=iso_tmp(:)
-
-#ifdef ISOVERIF
-      do i=1,klon
-          if (iso_verif_noNaN_nostop(xtsol(ixt,i), &
-     &          'phyiso_etat0_fichier 95').eq.1) then
-            write(*,*) 'ixt,i=',ixt,i
-            stop
-          endif       
-      enddo !do i=1,klon
-#endif
-
-       endif
-
+      IF(ixt <= niso .OR. initialisation_isotrac == 0) THEN
+#endif
+      found = phyetat0iso_srf3(xtsnow,      "XTSNOW", "Surface snow", 0.)
+      if (.NOT.found) CALL abort_physic('isotopes_routines_mod', 'phyiso_etat0_fichier 16581: unfound isotopic variable',1)
+      found = phyetat0iso_srf3(fxtevap,     "XTEVAP", "evaporation",  0.)
+      found = phyetat0iso_get2(xtrain_fall, "xtrain_f", "xrain fall", 0.)
+      found = phyetat0iso_get2(xtsnow_fall, "xtsnow_f", "xsnow fall", 0.)
+      found = phyetat0iso_get3(xt_ancien,   "XTANCIEN",  "QANCIEN",   0.)
+      found = phyetat0iso_get3(xtl_ancien,  "XTLANCIEN", "QLANCIEN",  0.)
+      found = phyetat0iso_get3(xts_ancien,  "XTSANCIEN", "QSANCIEN",  0.)
+      found = phyetat0iso_get2(xtrun_off_lic_0, "XTRUNOFFLIC0", "RUNOFFLIC0", 0.)
+      found = phyetat0iso_get3(wake_deltaxt,  "WAKE_DELTAXT", "Delta hum. wake/env",  0.)
+#ifdef ISOVERIF
+      IF(ixt == iso_eau .AND. iso_eau > 0) THEN
+         DO i=1,klon
+            CALL iso_verif_egalite(xtrain_fall(iso_eau,i),rain_fall(i),TRIM(modname)//' 231a')
+            CALL iso_verif_egalite(xtsnow_fall(iso_eau,i),snow_fall(i),TRIM(modname)//' 231b')
+            DO nsrf = 1, nbsrf
+               CALL iso_verif_egalite(fxtevap(iso_eau,i,nsrf),fevap(i,nsrf),TRIM(modname)//' 231c')
+               CALL iso_verif_egalite( xtsnow(iso_eau,i,nsrf), snow(i,nsrf),TRIM(modname)//' 231d')
+            END DO
+         END DO
+      END IF
+      IF(ixt == iso_HDO .AND. iso_HDO > 0) THEN
+         DO k=1,klev
+            DO i=1,klon
+               IF(q_ancien(i,k) > 2e-3) &
+                  CALL iso_verif_aberrant(xt_ancien(iso_hdo,i,k)/q_ancien(i,k),TRIM(modname)//' 312')
+            END DO
+         END DO
+      END IF
+      IF(iso_eau > 0 .AND. ixt == iso_eau) THEN
+         DO i=1,klon
+            IF(iso_verif_egalite_nostop(run_off_lic_0(i),xtrun_off_lic_0(iso_eau,i),TRIM(modname)//' 326') == 1) THEN
+               WRITE(*,*) 'i=',i
+               STOP
+            END IF
+         END DO
+      END IF
+#endif
+      ! ces variables n'ont pas de traceurs:
+      IF(ixt <= niso) THEN
+         found = phyetat0iso_get2(xtsol, "XTSOL", "Surface humidity / bucket", 0.)
+         found = phyetat0iso_get2(Rland_ice, "Rland_ice", "SR land ice", 0.)
+#ifdef ISOVERIF
+
+         DO i=1,klon
+            IF(iso_verif_noNaN_nostop(xtsol(ixt,i),TRIM(modname)//' 95') == 1) THEN
+               WRITE(*,*) 'ixt,i=',ixt,i
+               STOP
+            END IF
+         END DO
+#endif
+      END IF
 #ifdef ISOTRAC
-     endif !if ((ixt.le.niso).or.(initialisation_isotrac.eq.0)) then
-#endif
-
-  enddo !do ixt=1,ntraciso
+      END IF ! IF(ixt > niso .OR. initialisation_isotrac == 0))
+#endif
+
+   END DO
 
 #ifdef ISOTRAC
-        if (initialisation_isotrac.ne.0) then
-        ! on n'initialise pas d'après le fichier
-        ! l'eau normale est mise dans la zone izone_init
-
-        do ixt=niso+1,ntraciso
-
-             iiso=index_iso(ixt)
-
-             if (index_zone(ixt).eq.izone_init) then
-                do i=1,klon
-                 do nsrf = 1, nbsrf
-                  fxtevap(ixt,i,nsrf)=fxtevap(iiso,i,nsrf)
-                 enddo !do nsrf = 1, nbsrf
-                 xtsnow_fall(ixt,i)=xtsnow_fall(iiso,i)
-                 xtrain_fall(ixt,i)=xtrain_fall(iiso,i)
-                 do k=1,klev
-                   xt_ancien(ixt,i,k)=xt_ancien(iiso,i,k)
-                   xtl_ancien(ixt,i,k)=xtl_ancien(iiso,i,k)
-                   xts_ancien(ixt,i,k)=xts_ancien(iiso,i,k)
-                   wake_deltaxt(ixt,i,k)= wake_deltaxt(iiso,i,k)    
-                 enddo
-                enddo !do i=1,klon
-             else !if (index_zone(ixt).eq.izone_init) then
-                do i=1,klon
-                 do nsrf = 1, nbsrf
-                  fxtevap(ixt,i,nsrf)=0.0
-                 enddo !do nsrf = 1, nbsrf
-                 xtsnow_fall(ixt,i)=0.0
-                 xtrain_fall(ixt,i)=0.0
-                 do k=1,klev
-                   xt_ancien(ixt,i,k)=0.0
-                   xtl_ancien(ixt,i,k)=0.0
-                   xts_ancien(ixt,i,k)=0.0
-                 enddo
-                enddo !do i=1,klon
-             endif !if (index_zone(ixt).eq.izone_init) then
-
-         enddo  !do ixt=1,niso
-      endif !if (initialisation_isotrac.eq.0) then
-
-
-#ifdef ISOVERIF
-        DO nsrf = 1, nbsrf
-         do i=1,klon
-               call iso_verif_traceur(fxtevap(1,i,nsrf), &
-     &                   'phyiso_etat0_fichier 426')
-         enddo !do i=1,klon 
-        enddo !DO nsrf = 1, nbsrf
-        do i=1,klon
-           call iso_verif_traceur(xtrain_fall(1,i), &
-     &                   'phyiso_etat0_fichier 466')
-           call iso_verif_traceur(xtsnow_fall(1,i), &
-     &                   'phyiso_etat0_fichier 468')
-        enddo !do i=1,klon 
-        do k=1,klev
-          do i=1,klon
-               call iso_verif_traceur(xt_ancien(1,i,k), &
-     &                   'phyiso_etat0_fichier 591')
-          enddo !do i=1,klon
-        enddo !do k=1,klev             
+   IF(initialisation_isotrac /= 0) THEN
+      ! On n'initialise pas d'apres le fichier. L'eau normale est mise dans la zone izone_init
+      DO ixt=niso+1,ntraciso
+         iiso=index_iso(ixt)
+         IF(index_zone(ixt) == izone_init) THEN
+            DO i = 1, klon
+               fxtevap(ixt,i,1:nbsrf) = fxtevap(iiso,i,1:nbsrf)
+               xtsnow_fall(ixt,i) = xtsnow_fall(iiso,i)
+               xtrain_fall(ixt,i) = xtrain_fall(iiso,i)
+               DO k = 1, klev
+                  xt_ancien   (ixt,i,k) = xt_ancien   (iiso,i,k)
+                  xtl_ancien  (ixt,i,k) = xtl_ancien  (iiso,i,k)
+                  xts_ancien  (ixt,i,k) = xts_ancien  (iiso,i,k)
+                  wake_deltaxt(ixt,i,k) = wake_deltaxt(iiso,i,k)    
+               END DO
+            END DO
+         ELSE
+            DO i = 1, klon
+               fxtevap(ixt,i,1:nbsrf)=0.0
+               xtsnow_fall(ixt,i)=0.0
+               xtrain_fall(ixt,i)=0.0
+               xt_ancien (ixt,i,1:klev) = 0.0
+               xtl_ancien(ixt,i,1:klev) = 0.0
+               xts_ancien(ixt,i,1:klev) = 0.0
+            END DO
+         END IF
+      END DO
+   END IF
+
+#ifdef ISOVERIF
+   DO nsrf = 1, nbsrf
+      DO i = 1, klon
+         CALL iso_verif_traceur(fxtevap(1,i,nsrf), 'phyiso_etat0_fichier 426')
+      END DO
+   END DO
+   DO i=1,klon
+      CALL iso_verif_traceur(xtrain_fall(1,i), 'phyiso_etat0_fichier 466')
+      CALL iso_verif_traceur(xtsnow_fall(1,i), 'phyiso_etat0_fichier 468')
+   END DO
+   DO k = 1, klev
+      DO i = 1, klon
+         CALL iso_verif_traceur(xt_ancien(1,i,k), 'phyiso_etat0_fichier 591')
+      END DO
+   END DO
 #endif 
         ! endif ISOVERIF       
@@ -16755,9 +16671,43 @@
         ! endif ISOTRAC      
 
-! on ferme le fichier
-!      CALL close_startphy
-! déjà fermé dans phyetat0
+CONTAINS
+
+LOGICAL FUNCTION phyetat0iso_get2(field, pref, descr, default) RESULT(lFound)
+  REAL,             INTENT(INOUT) :: field(:,:)
+  CHARACTER(LEN=*), INTENT(IN)    :: pref, descr
+  REAL,             INTENT(IN)    :: default
+  REAL :: iso_tmp(klon)
+  nam(1) = TRIM(pref)//TRIM(outiso)
+  nam(2) = TRIM(pref)//TRIM(oldIso)
+  lFound = phyetat0_get(iso_tmp, nam, descr, default)
+  field(ixt,:) = iso_tmp
+END FUNCTION phyetat0iso_get2
+
+
+LOGICAL FUNCTION phyetat0iso_get3(field, pref, descr, default) RESULT(lFound)
+  REAL,             INTENT(INOUT) :: field(:,:,:)
+  CHARACTER(LEN=*), INTENT(IN)    :: pref, descr
+  REAL,             INTENT(IN)    :: default
+  REAL :: iso_tmp_lonlev(klon,klev)
+  nam(1) = TRIM(pref)//TRIM(outiso)
+  nam(2) = TRIM(pref)//TRIM(oldIso)
+  lFound = phyetat0_get(iso_tmp_lonlev, nam, descr, default)
+  field(ixt,:,:) = iso_tmp_lonlev(:,:)
+END FUNCTION phyetat0iso_get3
+
+LOGICAL FUNCTION phyetat0iso_srf3(field, pref, descr, default) RESULT(lFound)
+  REAL,             INTENT(INOUT) :: field(:,:,:)
+  CHARACTER(LEN=*), INTENT(IN)    :: pref, descr
+  REAL,             INTENT(IN)    :: default
+  REAL :: iso_tmp_lonsrf(klon,nbsrf)
+  nam(1) = TRIM(pref)//TRIM(outiso)
+  nam(2) = TRIM(pref)//TRIM(oldIso)
+  lFound = phyetat0_srf(iso_tmp_lonsrf, nam, descr, default)
+  field(ixt,:,:) = iso_tmp_lonsrf
+END FUNCTION phyetat0iso_srf3
 
         end subroutine phyiso_etat0_fichier
+
+
 
 
Index: LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/isotopes_verif_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/isotopes_verif_mod.F90	(revision 4175)
+++ LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/isotopes_verif_mod.F90	(revision 4482)
@@ -5,8 +5,8 @@
 MODULE isotopes_verif_mod
 !use isotopes_mod, ONLY: 
-#ifdef ISOTRAC
-   USE isotrac_mod, ONLY: nzone
-#endif
-USE infotrac_phy, ONLY: ntraciso=>ntiso, niso, itZonIso
+!#ifdef ISOTRAC
+!   USE isotrac_mod, ONLY: nzone
+!#endif
+USE infotrac_phy, ONLY: ntraciso=>ntiso, niso, itZonIso, nzone
 implicit none
 save
Index: LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/isotrac_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/isotrac_mod.F90	(revision 4175)
+++ LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/isotrac_mod.F90	(revision 4482)
@@ -1,104 +1,27 @@
 #ifdef ISO
 #ifdef ISOTRAC
-! $Id: $
 
 MODULE isotrac_mod
-use infotrac_phy, ONLY: niso,ntiso,ntraceurs_zone=>nzone
-use isotopes_mod, only: ridicule
-
-IMPLICIT NONE
-SAVE
-
-! contient toutes les variables traceurs isotopiques + les routines specifiquement
-! traceurs isotopiques
-
-      real ridicule_trac
-      parameter (ridicule_trac=ridicule*1e4)
-
-integer, save ::  option_traceurs 
-integer, save ::  ntraceurs_zone_opt ! ntraceurs_zone propre à l'option
-! on vérifie que ça correspond bien à ntraceurs_zone d'infotrac
-integer, save ::  ntraceurs_zoneOR
-!$OMP THREADPRIVATE(option_traceurs,ntraceurs_zone_opt,ntraceurs_zoneOR)
-integer, save ::  initialisation_isotrac
-                ! 1 pour idéalisé
-                ! 0 pour lecture dans fichier
-!$OMP THREADPRIVATE(initialisation_isotrac)
-
-! variables spécifiques aux différentes options, mais necessaires au
-! calcul du nombre de zones de traceurs
-! si option=3
-integer, save :: use_bassin_atlantic
-!$OMP THREADPRIVATE(use_bassin_atlantic)
-integer, save :: use_bassin_medit
-!$OMP THREADPRIVATE(use_bassin_medit)
-integer, save :: use_bassin_indian
-!$OMP THREADPRIVATE(use_bassin_indian)
-integer, save :: use_bassin_austral
-!$OMP THREADPRIVATE(use_bassin_austral)
-integer, save :: use_bassin_pacific
-!$OMP THREADPRIVATE(use_bassin_pacific)
-integer, save :: use_bassin_merarabie
-!$OMP THREADPRIVATE(use_bassin_merarabie)
-integer, save :: use_bassin_golfebengale
-!$OMP THREADPRIVATE(use_bassin_golfebengale)
-integer, save :: use_bassin_indiansud
-!$OMP THREADPRIVATE(use_bassin_indiansud)
-integer, save :: use_bassin_tropics
-!$OMP THREADPRIVATE(use_bassin_tropics)
-integer, save :: use_bassin_midlats
-!$OMP THREADPRIVATE(use_bassin_midlats)
-integer, save :: use_bassin_hauteslats
-!$OMP THREADPRIVATE(use_bassin_hauteslats)
-integer, save :: bassin_atlantic
-!$OMP THREADPRIVATE(bassin_atlantic)
-integer, save :: bassin_medit
-!$OMP THREADPRIVATE(bassin_medit)
-integer, save :: bassin_indian
-!$OMP THREADPRIVATE(bassin_indian)
-integer, save :: bassin_austral
-!$OMP THREADPRIVATE(bassin_austral)
-integer, save :: bassin_pacific
-!$OMP THREADPRIVATE(bassin_pacific)
-integer, save :: bassin_merarabie
-!$OMP THREADPRIVATE(bassin_merarabie)
-integer, save :: bassin_golfebengale
-!$OMP THREADPRIVATE(bassin_golfebengale)
-integer, save :: bassin_indiansud
-!$OMP THREADPRIVATE(bassin_indiansud)
-integer, save :: bassin_tropics
-!$OMP THREADPRIVATE(bassin_tropics)
-integer, save :: bassin_midlats
-!$OMP THREADPRIVATE(bassin_midlats)
-integer, save :: bassin_hauteslats
-!$OMP THREADPRIVATE(bassin_hauteslats)
-! si option=4
-integer nzone_temp
-parameter (nzone_temp=1)
-real, save :: zone_temp1,zone_tempf,zone_tempa  
-!$OMP THREADPRIVATE(zone_temp1,zone_tempf,zone_tempa)
-! si option 14
-integer nzone_lat
-parameter (nzone_lat=4)
-integer nzone_pres
-parameter (nzone_pres=3)
-real, save :: zone_pres1,zone_presf,zone_presa 
-!$OMP THREADPRIVATE(zone_pres1,zone_presf,zone_presa) 
-real, save :: dlattag,lattag_min
-!$OMP THREADPRIVATE(dlattag,lattag_min)
-
-
-! option 1: on trace evap ocean et continent séparement  
+  USE infotrac_phy,      ONLY: niso, ntiso, nzone
+  USE readTracFiles_mod, ONLY: delPhase
+  USE isotopes_mod,      ONLY: ridicule, get_in
+
+  IMPLICIT NONE
+  SAVE
+
+!=== CONTENT: ALL THE ISOTOPIC TRACERS RELATED VARIABLES ===
+!
+! option 1: on trace evap ocean et continent separement  
 ! option 2: on trace evap ocean, continent et evap precip
-! option 3: on trace evap différents bassins océaniques 
-!       + continents + résidu
-!       attention, choisir dans ce cas les bassins océaniques
+! option 3: on trace evap differents bassins oceaniques 
+!       + continents + residu
+!       attention, choisir dans ce cas les bassins oceaniques
 !       dans iso_traceurs_opt3F90.h 
-! option 4: tracage par température minimale
-!       dans ce cas, on définit des bins dans iso_traceurs_opt4.h
-! option 5: pour AMMA: on taggue résidu/AEJ/flux mousson/Harmattan
+! option 4: tracage par temperature minimale
+!       dans ce cas, on definit des bins dans iso_traceurs_opt4.h
+! option 5: pour AMMA: on taggue residu/AEJ/flux mousson/Harmattan
 ! option 6: taggage des ddfts
-! option 7: pour Sandrine: taggage de la vapeur à 700hPa pour omega500<-20 TODO
-! option 8: pour Sandrine: taggage de la vapeur entre 950 et 800hPa, omega de 0 à 25 hPa et de l'évaoration en omega<-20. TODO
+! option 7: pour Sandrine: taggage de la vapeur a 700hPa pour omega500<-20 TODO
+! option 8: pour Sandrine: taggage de la vapeur entre 950 et 800hPa, omega de 0 a 25 hPa et de l'evaoration en omega<-20. TODO
 ! option 9: taggage du condensat et de la revap precip
 ! option 10: taggage evap oce, transpiration et evaporation
@@ -107,910 +30,645 @@
 ! option 12: taggage evap oce, sol nu, canop et reste evap cont. 
 ! A utiliser quand on couple avec ORCHIDEE
-! option 13: taggage température minimale + revap precip
-! option 14: taggage lat et altitude de dernière saturation (niveaux de pression) + evap surf
+! option 13: taggage temperature minimale + revap precip
+! option 14: taggage lat et altitude de derniere saturation (niveaux de pression) + evap surf
 ! otion 15: taggage irrigation
 ! option 16: taggage precip selon saisons et fonte neige: seulement pour ORCHIDEE
-! option 17: taggage température minimum de condensation directement dans la convection et la cond LS, + evap sfc, condensat et precipitation
+! option 17: taggage temperature minimum de condensation directement dans la convection et la cond LS, + evap sfc, condensat et precipitation
 ! option 18: idem 17 mais on tague qsmin au lieu de Tmin
 ! option 19: on tag vap residuelle, vap residuelle dans ddfts, sfc, cond, rev
 ! option 20: on taggue vapeur tropicale vs vapeur extratropicale
 ! option 21: taggage de 2 boites 3D: extratropiques (>35°) et UT tropicale (15-15°, > 500hPa)
-! option 22: tagage de la vapeur proccessée dans les zones très convectives
+! option 22: tagage de la vapeur proccessee dans les zones tres convectives
                 
-        ! ces variables sont initialisées dans traceurs_init
+   !--- nzone_opt (value of nzone for the selected option) must be equal to nzone as defined in onfotrac
+   REAL, PARAMETER :: ridicule_trac = ridicule * 1e4
+   INTEGER, SAVE :: option_traceurs, nzone_opt, nzoneOR
+!$OMP THREADPRIVATE(option_traceurs,nzone_opt,nzoneOR)
+   INTEGER, SAVE :: initialisation_isotrac
+!$OMP THREADPRIVATE(initialisation_isotrac)      
+                ! 1 pour idealise
+                ! 0 pour lecture dans fichier
+
+   !=== VARIABLES SPECIFIC TO THE SELECTED OPTION, BUT NEEDED FOR THE COMPUTATION OF THE NUMBER OF ZONES ; TO BE INITIALIZED IN traceurs_init
+
+   !--- option 3
+   LOGICAL, SAVE :: use_bassin_Austral, use_bassin_Atlantic, use_bassin_MidLats, use_bassin_SouthIndian, use_bassin_MerArabie
+!$OMP THREADPRIVATE(use_bassin_Austral, use_bassin_Atlantic, use_bassin_MidLats, use_bassin_SouthIndian, use_bassin_MerArabie)
+   INTEGER, SAVE ::     bassin_Austral,     bassin_Atlantic,     bassin_MidLats,     bassin_SouthIndian,     bassin_MerArabie
+!$OMP THREADPRIVATE(    bassin_Austral,     bassin_Atlantic,     bassin_MidLats,     bassin_SouthIndian,     bassin_MerArabie)
+   LOGICAL, SAVE :: use_bassin_Pacific, use_bassin_Indian,   use_bassin_Tropics, use_bassin_BengalGolf,  use_bassin_HighLats, use_bassin_Medit
+!$OMP THREADPRIVATE(use_bassin_Pacific, use_bassin_Indian,   use_bassin_Tropics, use_bassin_BengalGolf,  use_bassin_HighLats, use_bassin_Medit)
+   INTEGER, SAVE ::     bassin_Pacific,     bassin_Indian,       bassin_Tropics,     bassin_BengalGolf,      bassin_HighLats,     bassin_Medit
+!$OMP THREADPRIVATE(    bassin_Pacific,     bassin_Indian,       bassin_Tropics,     bassin_BengalGolf,      bassin_HighLats,     bassin_Medit)
+
+   !--- option 4
+   INTEGER, PARAMETER :: nzone_temp = 1
+   REAL,   SAVE ::  zone_temp1, zone_tempf, zone_tempa  
+!$OMP THREADPRIVATE(zone_temp1, zone_tempf, zone_tempa)
+   REAL,   SAVE ::  zone_temp(nzone_temp-1)
+!$OMP THREADPRIVATE(zone_temp)
+
+   !--- option 5
+   INTEGER, SAVE :: izone_aej, izone_harmattan, izone_mousson
+!$OMP THREADPRIVATE(izone_aej, izone_harmattan, izone_mousson)
+
+   !--- option 6
+   INTEGER, SAVE :: izone_ddft
+!$OMP THREADPRIVATE(izone_ddft)
+
+   !--- option 10
+   INTEGER, SAVE :: izone_contfrac
+!$OMP THREADPRIVATE(izone_contfrac)
+
+   !--- option 12	
+   INTEGER, SAVE :: izone_contcanop
+!$OMP THREADPRIVATE(izone_contcanop)
+
+   !--- option 13
+   INTEGER, PARAMETER :: nzone_pres = 3
+   REAL, SAVE ::  zone_pres(nzone_pres-1)
+!$OMP THREADPRIVATE(zone_pres)
+
+   !--- option 14
+   INTEGER, PARAMETER :: nzone_lat = 4
+   REAL,    SAVE :: zone_pres1, zone_presf, zone_presa 
+!$OMP THREADPRIVATE(zone_pres1, zone_presf, zone_presa) 
+   REAL,    SAVE :: dlattag, lattag_min, zone_lat(nzone_lat-1)
+!$OMP THREADPRIVATE(dlattag, lattag_min, zone_lat)
+
+   !--- option 15
+   INTEGER, SAVE :: izone_irrig
+!$OMP THREADPRIVATE(izone_irrig)
+
+   !--- option 17
+   REAL,    SAVE :: seuil_tag_tmin, seuil_tag_tmin_ls
+!$OMP THREADPRIVATE(seuil_tag_tmin, seuil_tag_tmin_ls)
+  INTEGER,  SAVE :: option_seuil_tag_tmin
+!$OMP THREADPRIVATE(option_seuil_tag_tmin)
+
+   !--- option 20
+   INTEGER, SAVE :: izone_trop, izone_extra
+!$OMP THREADPRIVATE(izone_trop, izone_extra)
+   REAL,    SAVE :: lim_tag20 
+!$OMP THREADPRIVATE(lim_tag20)
+
+   !--- option 21: on garde izone_trop, izone_extra  
+
+   !--- option 22
+   INTEGER, SAVE :: izone_conv_BT, izone_conv_UT
+!$OMP THREADPRIVATE(izone_conv_BT, izone_conv_UT)
+   REAL,    SAVE :: lim_precip_tag22 
+!$OMP THREADPRIVATE(lim_precip_tag22)
+
         
-integer, ALLOCATABLE, DIMENSION(:), save :: index_iso
-!$OMP THREADPRIVATE(index_iso)
-integer, ALLOCATABLE, DIMENSION(:), save ::  index_zone
-!$OMP THREADPRIVATE(index_zone)
-integer, ALLOCATABLE, DIMENSION(:,:), save ::  itZonIso_loc ! il y a déjà un itZonIso dans infotrac: vérifier que c'est le même
-!$OMP THREADPRIVATE(itZonIso_loc)
-character*3, ALLOCATABLE, DIMENSION(:), save :: strtrac
-!$OMP THREADPRIVATE(strtrac)
-! -> tout ça passe maintenant par infotrac
-
-integer, ALLOCATABLE, DIMENSION(:), save :: bassin_map
-integer, ALLOCATABLE, DIMENSION(:,:), save :: boite_map
-!$OMP THREADPRIVATE(bassin_map,boite_map)
-
-
-        ! traitement recyclage et evap
-integer, save :: izone_cont ! pour le recyclage continental
-!$OMP THREADPRIVATE(izone_cont)
-integer, save :: izone_oce ! pour l'océan
-!$OMP THREADPRIVATE(izone_oce)
-integer, save :: izone_poubelle ! pour les petits résidus numériques
+  INTEGER, ALLOCATABLE, SAVE :: index_iso(:), index_zone(:), itZonIso_loc(:,:)
+!$OMP THREADPRIVATE(            index_iso,    index_zone,    itZonIso_loc)
+  CHARACTER(LEN=3), ALLOCATABLE :: strtrac(:)
+!$OMP THREADPRIVATE(               strtrac)
+  INTEGER, ALLOCATABLE, SAVE :: bassin_map(:), boite_map(:,:)
+!$OMP THREADPRIVATE(            bassin_map,    boite_map)
+
+   !=== RECYCLING AND EVAPORATION TREATMENT
+   INTEGER, SAVE :: izone_cont, izone_oce        !--- For land and ocean recycling
+!$OMP THREADPRIVATE(izone_cont, izone_oce)
+   INTEGER, SAVE :: izone_poubelle               !--- For small numerical residues
 !$OMP THREADPRIVATE(izone_poubelle)
-integer, save :: izone_init ! pour l'initialisation par défaut
+   INTEGER, SAVE :: izone_init                   !--- For default initialization
 !$OMP THREADPRIVATE(izone_init)
-integer, save :: izone_revap ! pour l'évap des gouttes
+   INTEGER, SAVE :: izone_revap                  !--- For droplets evaporation
 !$OMP THREADPRIVATE(izone_revap)
-integer, save :: option_revap
-!$OMP THREADPRIVATE(option_revap)
-integer, save :: option_tmin
-!$OMP THREADPRIVATE(option_tmin)
-integer, save :: option_cond
-!$OMP THREADPRIVATE(option_cond)
-integer, save :: izone_cond
-!$OMP THREADPRIVATE(izone_cond)
-real evap_franche
-parameter (evap_franche=1e-6) ! en kg/m2/s
-
-! specifique à option 4:
-real, save ::  zone_temp(nzone_temp-1)
-!$OMP THREADPRIVATE(zone_temp)
-! si option 5
-integer, save :: izone_aej
-!$OMP THREADPRIVATE(izone_aej)
-integer, save :: izone_harmattan
-!$OMP THREADPRIVATE(izone_harmattan)
-integer, save :: izone_mousson
-!$OMP THREADPRIVATE(izone_mousson)
-! si option 6
-integer, save :: izone_ddft
-!$OMP THREADPRIVATE(izone_ddft)
-! si option 10
-integer, save :: izone_contfrac
-!$OMP THREADPRIVATE(izone_contfrac)
-! si option 12	
-integer, save :: izone_contcanop
-!$OMP THREADPRIVATE(izone_contcanop)
-! specifique à option 13:
-real, save ::  zone_pres(nzone_pres-1)
-!$OMP THREADPRIVATE(zone_pres)
-! si option 14
-real, save ::  zone_lat(nzone_lat-1)
-!$OMP THREADPRIVATE(zone_lat)
-! si option 15
-integer, save :: izone_irrig
-!$OMP THREADPRIVATE(izone_irrig)
-! si option 17
-real, save ::  seuil_tag_tmin
-!$OMP THREADPRIVATE(seuil_tag_tmin)
-real, save ::  seuil_tag_tmin_ls
-!$OMP THREADPRIVATE(seuil_tag_tmin_ls)
-integer, save :: option_seuil_tag_tmin
-!$OMP THREADPRIVATE(option_seuil_tag_tmin)
-! si option 20
-integer, save :: izone_trop,izone_extra
-real, save ::  lim_tag20 
-!$OMP THREADPRIVATE(izone_trop,izone_extra,lim_tag20)
-! si option 21: on garde izone_trop,izone_extra  
-! si opt 22
-integer, save :: izone_conv_BT,izone_conv_UT
-real, save ::  lim_precip_tag22 
-!$OMP THREADPRIVATE(izone_conv_BT,izone_conv_UT,lim_precip_tag22)
-
+   INTEGER, SAVE :: option_revap, option_tmin, option_cond, izone_cond
+!$OMP THREADPRIVATE(option_revap, option_tmin, option_cond, izone_cond)
+   REAL, PARAMETER :: evap_franche = 1e-6        !--- In kg/m2/s
 
 CONTAINS
 
-      subroutine iso_traceurs_init()
-
-      use IOIPSL ! getin
-      USE infotrac_phy, ONLY: itZonIso
-      USE isotopes_mod, ONLY: iso_eau,ntracisoOR,initialisation_iso
-      USE dimphy, only: klon,klev
-
-        implicit none
-
-
-        ! définition de quelles zones et quelles isotopes représentent
-        ! les traceurs
-
-        ! inputs, outputs
-        ! ! c'est les variables dans traceurs.h qui sont modifiées
-
-        ! locals
-        integer itrac,izone,ixt,k
-        integer izone_pres,izone_lat
-        character*2 strz,strz_preslat
-        character*1 strz_pres,strz_lat
-        integer ntraceurs_zone_opt
-
-        ! vérifier que on a bien l'eau comme traceurs
-        if (iso_eau.eq.0) then
-            write(*,*) 'traceurs_init 18: isotrac ne marche que si ', &
-     &            'on met l''eau comme isotope'
-            stop
-        endif
-
-        ! initialiser
-        option_traceurs=0
-        initialisation_isotrac=0
-
-        ! allouer
-        allocate (index_iso(ntiso))
-        allocate (index_zone(ntiso))
-        allocate (itZonIso_loc(ntraceurs_zone,niso))
-        allocate (strtrac(ntraceurs_zone))
-        allocate (bassin_map(klon))
-        allocate (boite_map(klon,klev))
-
-        if (initialisation_iso.eq.0) then
-          call getin('initialisation_isotrac',initialisation_isotrac)
-          write(*,*) 'initialisation_isotrac=',initialisation_isotrac
-        endif !if (initialisation_iso.eq.0) then
-
-        ! lire l'option de traçage
-        call getin('option_traceurs',option_traceurs)
-        write(*,*) 'option_traceurs=',option_traceurs
-
-        ! cas général: pas de traceurs dans ORCHIDEE
-        ntracisoOR=niso
-
-        ! partie à éditer ! pour définir les différentes zones
-        if (option_traceurs.eq.1) then
-          ! on trace continents/ocean  
-
-          ntraceurs_zone_opt=2
-          izone_cont=1
-          izone_oce=2          
-          izone_poubelle=2 ! zone où on met les flux non physiques, de
-                ! réajustement
-          izone_init=2 ! zone d'initialisation par défaut          
-          option_revap=0
-          option_tmin=0
-          izone_revap=0
-          option_cond=0
-
-          strtrac(izone_cont)='con'
-          strtrac(izone_oce)='oce'
-
-        elseif (option_traceurs.eq.2) then
-          ! on trace continent/ ocean/reevap des gouttes
-
-          ntraceurs_zone_opt=3
-          izone_cont=1
-          izone_oce=2
-          izone_poubelle=2 ! zone où on met les flux non physiques, de
-                ! réajustement
-          izone_init=2 ! zone d'initialisation par défaut
-          option_revap=1
-          option_tmin=0
-          izone_revap=3
-          option_cond=0
-
-          strtrac(izone_cont)='con'
-          strtrac(izone_oce)='oce'
-          strtrac(izone_revap)='rev'
-          
-
-        else if (option_traceurs.eq.3) then
-            ! on trace des bassins océaniques + un résidu. On ne trace
-            ! pas l'évap des gouttes à part
-            ! le résidu est la dernère dimension
-            
-          ! lire les use_bassin
-          call getin('use_bassin_atlantic',use_bassin_atlantic)     
-          call getin('use_bassin_medit',use_bassin_medit)      
-          call getin('use_bassin_indian',use_bassin_indian)     
-          call getin('use_bassin_austral',use_bassin_austral)     
-          call getin('use_bassin_pacific',use_bassin_pacific)     
-          call getin('use_bassin_merarabie',use_bassin_merarabie)     
-          call getin('use_bassin_golfebengale',use_bassin_golfebengale)     
-          call getin('use_bassin_indiansud',use_bassin_indiansud)     
-          call getin('use_bassin_tropics',use_bassin_tropics)     
-          call getin('use_bassin_midlats',use_bassin_midlats)     
-          call getin('use_bassin_hauteslats',use_bassin_hauteslats)
-
-          write(*,*) 'use_bassin_atlantic=' ,use_bassin_atlantic  
-          write(*,*) 'use_bassin_medit=' ,use_bassin_medit
-          write(*,*) 'use_bassin_indian=' ,use_bassin_indian
-          write(*,*) 'use_bassin_austral=' ,use_bassin_austral
-          write(*,*) 'use_bassin_merarabie=' ,use_bassin_merarabie
-          write(*,*) 'use_bassin_golfebengale=' ,use_bassin_golfebengale
-          write(*,*) 'use_bassin_indiansud=' ,use_bassin_indiansud
-          write(*,*) 'use_bassin_tropics=' ,use_bassin_tropics
-          write(*,*) 'use_bassin_midlats=' ,use_bassin_midlats
-          write(*,*) 'use_bassin_hauteslats=' ,use_bassin_hauteslats
-
-        
-          ntraceurs_zone_opt=2 &
-     &                   +use_bassin_atlantic &
-     &                   +use_bassin_medit &
-     &                   +use_bassin_indian &
-     &                   +use_bassin_austral &
-     &                   +use_bassin_pacific &
-     &                   +use_bassin_merarabie &
-     &                   +use_bassin_golfebengale &
-     &                   +use_bassin_indiansud &
-     &                   +use_bassin_tropics &
-     &                   +use_bassin_midlats &
-     &                   +use_bassin_hauteslats
-
-          izone_cont=ntraceurs_zone
-          izone_oce=0 ! pas de sens car séparée en bassins          
-          izone_poubelle=ntraceurs_zone-1 ! zone où on met les flux non physiques, de
-                ! réajustement
-          izone_init=ntraceurs_zone-1 ! zone d'initialisation par défaut
-          option_revap=0 ! on ne trace pas les gouttes 
-          option_tmin=0
-          izone_revap=0 ! pas de sens car on taggue pas les gouttes séparemment  
-          option_cond=0
-
-          ! si on a use_bassin_indian, on n'a pas le découpage détaillé
-          ! de l'indian:
+   SUBROUTINE iso_traceurs_init()
+
+   USE infotrac_phy, ONLY: itZonIso, isoName, isoZone
+   USE isotopes_mod, ONLY: iso_eau, ntracisoOR, initialisation_iso
+   USE dimphy,       ONLY: klon, klev
+   USE  strings_mod, ONLY: int2str, strStack, strTail, strHead, strIdx, fmsg
+
+   IMPLICIT NONE
+   ! Define which zones and isotopes correspond to isotopic tagging tracers
+   ! Modify traceurs.h variables
+   INTEGER :: izone, ixt, k
+   INTEGER :: izone_pres, izone_lat
+   INTEGER :: nzone_opt
+
+   IF(fmsg("traceurs_init 18: isotrac ne marche que si on met l'eau comme isotope", 'iso_traceurs_init', iso_eau==0)) STOP
+
+   !--- Initialize
+   option_traceurs = 0
+   initialisation_isotrac = 0
+
+   !--- Allocate
+   ALLOCATE(index_iso (ntiso))
+   ALLOCATE(index_zone(ntiso))
+   ALLOCATE(itZonIso_loc(nzone,niso))
+   ALLOCATE(strtrac(nzone))
+   ALLOCATE(bassin_map(klon))
+   ALLOCATE( boite_map(klon,klev))
+
+   IF(initialisation_iso == 0) CALL get_in('initialisation_isotrac', initialisation_isotrac)
+
+   !--- Read tracing option
+   CALL get_in('option_traceurs', option_traceurs)
+
+   !--- Genral case: no traceurs in ORCHIDEE
+   ntracisoOR=niso
+
+   ! partie a editer ! pour definir les differentes zones
+   SELECT CASE(option_traceurs)
+      !========================================================================================================================
+      CASE(1)      !=== TRACING LAND/OCEAN
+      !========================================================================================================================
+         nzone_opt=2
+         izone_cont=1
+         izone_oce=2
+         izone_poubelle=2    ! zone ou on met les flux non physiques, de reajustement
+         izone_init=2        ! zone d'initialisation par defaut          
+         option_revap=0
+         option_tmin=0
+         izone_revap=0
+         option_cond=0
+         strtrac(izone_cont) = 'con'
+         strtrac(izone_oce)  = 'oce'
+      !========================================================================================================================
+      CASE(2)      !=== TRACING LAND/OCEAN/DROPLETS REEVAPORATION
+      !========================================================================================================================
+         nzone_opt=3
+         izone_cont=1
+         izone_oce=2
+         izone_poubelle=2    ! zone ou on met les flux non physiques, de reajustement
+         izone_init=2        ! zone d'initialisation par defaut          
+         option_revap=1
+         option_tmin=0
+         izone_revap=3
+         option_cond=0
+         strtrac(izone_cont) = 'con'
+         strtrac(izone_oce)  = 'oce'
+         strtrac(izone_revap)= 'rev'
+      !========================================================================================================================
+      CASE(3)      !=== TRACING OCEANS BASINS + RESIDUE (LAST DIMENSION). NO DROPLETS EVAPORATION TRACING.
+      !========================================================================================================================
+         ! lire les use_bassin
+         CALL get_in('use_bassin_Atlantic',   use_bassin_Atlantic)
+         CALL get_in('use_bassin_Medit',      use_bassin_Medit)
+         CALL get_in('use_bassin_Indian',     use_bassin_Indian)
+         CALL get_in('use_bassin_Austral',    use_bassin_Austral)
+         CALL get_in('use_bassin_Pacific',    use_bassin_Pacific)
+         CALL get_in('use_bassin_MerArabie',  use_bassin_MerArabie)
+         CALL get_in('use_bassin_BengalGolf', use_bassin_BengalGolf)
+         CALL get_in('use_bassin_SouthIndian',use_bassin_SouthIndian)
+         CALL get_in('use_bassin_Tropics',    use_bassin_Tropics)
+         CALL get_in('use_bassin_Midlats',    use_bassin_Midlats)
+         CALL get_in('use_bassin_HighLats',   use_bassin_HighLats)
+         nzone_opt  =  2  +  COUNT([use_bassin_Atlantic, use_bassin_Medit,     use_bassin_Indian,     &
+            use_bassin_Austral,     use_bassin_Pacific,  use_bassin_MerArabie, use_bassin_BengalGolf, &
+            use_bassin_SouthIndian, use_bassin_Tropics,  use_bassin_Midlats,   use_bassin_HighLats])
+         izone_cont=nzone
+         izone_oce=0             ! pas de sens car separee en bassins          
+         izone_poubelle=nzone-1  ! zone ou on met les flux non physiques, de reajustement
+         izone_init=nzone-1      ! zone d'initialisation par defaut
+         option_revap=0          ! on ne trace pas les gouttes 
+         option_tmin=0
+         izone_revap=0           ! pas de sens car on taggue pas les gouttes separemment  
+         option_cond=0
 #ifdef ISOVERIF
-          if (use_bassin_indian.eq.1) then
-!              call iso_verif_egalite(float(use_bassin_merarabie), &
-!     &            0.0,'iso_traceurs_init 73: revoir def des bassins')
-               if ((use_bassin_merarabie.ne.0).or. &
-      &            (use_bassin_indiansud.ne.0).or. &
-      &            (use_bassin_golfebengale.ne.0)) then
-                write(*,*) 'traceurs_init 73'
-                stop
-               endif 
-!              call iso_verif_egalite(float(use_bassin_golfebengale), &
-!     &            0.0,'iso_traceurs_init 73: revoir def des bassins')
-!              call iso_verif_egalite(float(use_bassin_indiansud), &
-!     &            0.0,'iso_traceurs_init 73: revoir def des bassins')
-          endif
+         IF(use_bassin_Indian) THEN   !=== NON COMPATIBLE WITH A DETAILED INDIAN CUTTING
+            IF(use_bassin_MerArabie .OR. use_bassin_SouthIndian .OR. use_bassin_BengalGolf) THEN
+               WRITE(*,*) 'traceurs_init 73'; STOP
+            END IF
+!           CALL iso_verif_egalite(float(use_bassin_MerArabie),   0.0, 'iso_traceurs_init 73: revoir def des bassins')
+!           CALL iso_verif_egalite(float(use_bassin_BengalGolf),  0.0, 'iso_traceurs_init 73: revoir def des bassins')
+!           CALL iso_verif_egalite(float(use_bassin_SouthIndian), 0.0, 'iso_traceurs_init 73: revoir def des bassins')
+         END IF
 #endif    
-          
-          bassin_atlantic= max(use_bassin_atlantic,1)
-          bassin_medit=max(use_bassin_atlantic &
-     &           +use_bassin_medit,1)
-          bassin_indian=max(use_bassin_atlantic &
-     &           +use_bassin_medit &
-     &           +use_bassin_indian,1)
-          bassin_austral=max(use_bassin_atlantic &
-     &           +use_bassin_medit &
-     &           +use_bassin_indian &
-     &           +use_bassin_austral,1)
-          bassin_pacific=max(use_bassin_atlantic &
-     &           +use_bassin_medit &
-     &           +use_bassin_indian &
-     &           +use_bassin_austral &
-     &           +use_bassin_pacific,1)
-          bassin_merarabie=max(use_bassin_atlantic &
-     &           +use_bassin_medit &
-     &           +use_bassin_indian &
-     &           +use_bassin_austral &
-     &           +use_bassin_pacific &
-     &           +use_bassin_merarabie,1)
-          bassin_golfebengale=max(use_bassin_atlantic& 
-     &           +use_bassin_medit &
-     &           +use_bassin_indian & 
-     &           +use_bassin_austral &
-     &           +use_bassin_pacific &
-     &           +use_bassin_merarabie &
-     &           +use_bassin_golfebengale,1)
-          bassin_indiansud=max(use_bassin_atlantic &
-     &           +use_bassin_medit &
-     &           +use_bassin_indian &
-     &           +use_bassin_austral &
-     &           +use_bassin_pacific &
-     &           +use_bassin_merarabie &
-     &           +use_bassin_golfebengale &
-     &           +use_bassin_indiansud,1) 
-          bassin_tropics=max(use_bassin_atlantic &
-     &                       +use_bassin_medit &
-     &                       +use_bassin_indian &
-     &                       +use_bassin_austral &
-     &                       +use_bassin_pacific &
-     &                       +use_bassin_merarabie &
-     &                       +use_bassin_golfebengale &
-     &                       +use_bassin_indiansud, &
-     &                       +use_bassin_tropics,1) 
-          bassin_midlats=max(use_bassin_atlantic &
-     &                       +use_bassin_medit &
-     &                       +use_bassin_indian &
-     &                       +use_bassin_austral &
-     &                       +use_bassin_pacific &
-     &                       +use_bassin_merarabie &
-     &                       +use_bassin_golfebengale &
-     &                       +use_bassin_indiansud &
-     &                       +use_bassin_tropics &
-     &                       +use_bassin_midlats,1) 
-          bassin_hauteslats=max(use_bassin_atlantic & 
-     &                       +use_bassin_medit &
-     &                       +use_bassin_indian &
-     &                       +use_bassin_austral &
-     &                       +use_bassin_pacific &
-     &                       +use_bassin_merarabie &
-     &                       +use_bassin_golfebengale &
-     &                       +use_bassin_indiansud &
-     &                       +use_bassin_tropics &
-     &                       +use_bassin_midlats &
-     &                       +use_bassin_hauteslats,1) 
-
-          write(*,*) 'bassin_atlantic=' ,bassin_atlantic  
-          write(*,*) 'bassin_medit=' ,bassin_medit
-          write(*,*) 'bassin_indian=' ,bassin_indian
-          write(*,*) 'bassin_austral=' ,bassin_austral
-          write(*,*) 'bassin_merarabie=' ,bassin_merarabie
-          write(*,*) 'bassin_golfebengale=' ,bassin_golfebengale
-          write(*,*) 'bassin_indiansud=' ,bassin_indiansud
-          write(*,*) 'bassin_tropics=' ,bassin_tropics
-          write(*,*) 'bassin_midlats=' ,bassin_midlats
-          write(*,*) 'bassin_hauteslats=' ,bassin_hauteslats
-
-          if (use_bassin_atlantic.eq.1) then
-            strtrac(bassin_atlantic)='atl'
-          endif
-          if (use_bassin_medit.eq.1) then
-            strtrac(bassin_medit)='med'
-          endif
-          if (use_bassin_indian.eq.1) then
-            strtrac(bassin_indian)='ind'
-          endif
-          if (use_bassin_austral.eq.1) then
-            strtrac(bassin_austral)='aus'
-          endif
-          if (use_bassin_pacific.eq.1) then
-            strtrac(bassin_pacific)='pac'
-          endif
-          if (use_bassin_merarabie.eq.1) then
-            strtrac(bassin_merarabie)='ara'
-          endif
-          if (use_bassin_golfebengale.eq.1) then
-            strtrac(bassin_golfebengale)='ben'
-          endif
-          if (use_bassin_indiansud.eq.1) then
-            strtrac(bassin_indiansud)='ins'
-          endif
-          if (use_bassin_tropics.eq.1) then
-            strtrac(bassin_tropics)='tro'
-          endif
-          if (use_bassin_midlats.eq.1) then
-            strtrac(bassin_midlats)='mid'
-          endif
-          if (use_bassin_hauteslats.eq.1) then
-            strtrac(bassin_hauteslats)='hau'
-          endif
-          strtrac(ntraceurs_zone-1)='res'
-          strtrac(ntraceurs_zone)='con'
-
-        else if (option_traceurs.eq.4) then
-          ! on trace les température minimales vécues
-          ! comme dans article sur LdG sauf pas de revap
-            
-          zone_temp1=293.0 ! en K
-!          zone_tempf=223.0 ! en K
-          zone_tempf=243.0 ! en K
- ! courbure de la relation entre l'indice et la température: 0 pour linéaire, <0 pour plus de détal en bas
-
+         bassin_Atlantic   = 1
+         bassin_Medit      = bassin_Atlantic    + COUNT([use_bassin_Medit]);       WRITE(*,*) 'bassin_Atlantic    =' ,bassin_Atlantic
+         bassin_Indian     = bassin_Medit       + COUNT([use_bassin_Indian]);      WRITE(*,*) 'bassin_Medit       =' ,bassin_Medit
+         bassin_Austral    = bassin_Indian      + COUNT([use_bassin_Austral]);     WRITE(*,*) 'bassin_Indian      =' ,bassin_Indian
+         bassin_Pacific    = bassin_Austral     + COUNT([use_bassin_Pacific]);     WRITE(*,*) 'bassin_Austral     =' ,bassin_Austral
+         bassin_MerArabie  = bassin_Pacific     + COUNT([use_bassin_MerArabie]);   WRITE(*,*) 'bassin_MerArabie   =' ,bassin_MerArabie
+         bassin_BengalGolf = bassin_MerArabie   + COUNT([use_bassin_BengalGolf]);  WRITE(*,*) 'bassin_BengalGolf  =' ,bassin_BengalGolf
+         bassin_SouthIndian= bassin_BengalGolf  + COUNT([use_bassin_SouthIndian]); WRITE(*,*) 'bassin_SouthIndian =' ,bassin_SouthIndian
+         bassin_Tropics    = bassin_SouthIndian + COUNT([use_bassin_Tropics]);     WRITE(*,*) 'bassin_Tropics     =' ,bassin_Tropics
+         bassin_MidLats    = bassin_Tropics     + COUNT([use_bassin_MidLats]);     WRITE(*,*) 'bassin_MidLats     =' ,bassin_MidLats
+         bassin_HighLats   = bassin_MidLats     + COUNT([use_bassin_HighLats]);    WRITE(*,*) 'bassin_HighLats    =' ,bassin_HighLats
+         IF(use_bassin_atlantic   ) strtrac(bassin_atlantic)   = 'atl'
+         IF(use_bassin_medit      ) strtrac(bassin_medit)      = 'med'
+         IF(use_bassin_indian     ) strtrac(bassin_indian)     = 'ind'
+         IF(use_bassin_austral    ) strtrac(bassin_austral)    = 'aus'
+         IF(use_bassin_pacific    ) strtrac(bassin_pacific)    = 'pac'
+         IF(use_bassin_merarabie  ) strtrac(bassin_merarabie)  = 'ara'
+         IF(use_bassin_BengalGolf ) strtrac(bassin_BengalGolf) = 'ben'
+         IF(use_bassin_SouthIndian) strtrac(bassin_SouthIndian)= 'ins'
+         IF(use_bassin_tropics    ) strtrac(bassin_tropics)    = 'tro'
+         IF(use_bassin_midlats    ) strtrac(bassin_midlats)    = 'mid'
+         IF(use_bassin_HighLats   ) strtrac(bassin_HighLats)   = 'hau'
+         strtrac(nzone-1)='res'
+         strtrac(nzone)='con'
+      !========================================================================================================================
+      CASE(4)      !=== TRACING MINIMAL EXPERIENCED TEMPERATURE AS IN THE ARTICLE ON LfG, EXCEPT NO REVAPORATION
+      !========================================================================================================================
+         zone_temp1 = 293.0  ! en K
+!        zone_tempf = 223.0  ! en K
+         zone_tempf = 243.0  ! en K
+        ! courbure de la relation entre l'indice et la temperature: 0 pour lineaire, <0 pour plus de detal en bas
         ! zone 1: >= zone_temp1
-        ! zone 2 à 4: intermédiaire, 
+        ! zone 2 a 4: intermediaire, 
         ! zone 5: <zone_tempf
-        
-          ntraceurs_zone_opt=nzone_temp+1
-
-          zone_tempa=-4.0 ! en K
-          izone_cont=ntraceurs_zone
-          izone_oce=ntraceurs_zone  
-          izone_poubelle=ntraceurs_zone
-          izone_init=ntraceurs_zone ! zone d'initialisation par défaut
-          option_revap=0 
-          option_tmin=0  
-          izone_revap=0
-          option_cond=0
-          do izone=1,nzone_temp
-            write(strz,'(i2.2)') izone
-            strtrac(izone)='t'//strz
-            write(*,*) 'izone,strz,strtrac=',izone,strz,strtrac(izone)
-          enddo
-          strtrac(izone_poubelle)='pou'
-
-          ! initialisation des zones de tempéarture
-          do izone=1,nzone_temp-1
-            zone_temp(izone)=zone_temp1+float(izone-1) &
-     &                      *(zone_tempa*float(izone-nzone_temp+1) &
-     &                      +(zone_tempf-zone_temp1)/float(nzone_temp-2))
-          enddo
-          write(*,*) 'iso_trac_init 183: zone_temp=',zone_temp          
-
-        elseif (option_traceurs.eq.5) then
-          ! on trace AEJ/flux de mousson/Harmattan
-!          write(*,*) 'iso_traceurs_init 129'
-
-          ntraceurs_zone_opt=4
-          izone_cont=1
-          izone_oce=1
-          izone_poubelle=1 ! zone où on met les flux non physiques, de
-                ! réajustement
-          izone_init=1 ! zone d'initialisation par défaut
-          option_revap=0
-          option_tmin=0
-          izone_revap=0
-          izone_aej=2
-          izone_mousson=3
-          izone_harmattan=4 
-          option_cond=0
-
-          strtrac(izone_poubelle)='res'
-          strtrac(izone_aej)='aej'
-          strtrac(izone_mousson)='mou'
-          strtrac(izone_harmattan)='sah'
-
-        elseif (option_traceurs.eq.6) then
-          ! on trace les ddfts
-
-          ntraceurs_zone_opt=2
-          izone_cont=1
-          izone_oce=1
-          izone_poubelle=1 ! zone où on met les flux non physiques, de
-                ! réajustement
-          izone_init=1 ! zone d'initialisation par défaut
-          option_revap=0
-          option_tmin=0
-          izone_revap=0
-          izone_ddft=2
-          option_cond=0
-
-          strtrac(izone_poubelle)='res'
-          strtrac(izone_ddft)='dft'
-
-        elseif (option_traceurs.eq.9) then
-          ! on trace le condensat
-
-          ntraceurs_zone_opt=3
-          izone_cont=1
-          izone_oce=1
-          izone_poubelle=1 ! zone où on met les flux non physiques, de
-                ! réajustement
-          izone_init=1 ! zone d'initialisation par défaut
-          option_revap=1
-          option_tmin=0
-          izone_revap=2
-          izone_cond=3
-          option_cond=1
-
-          ! 1 par défaut pour colorier à la fois condensat LS et
-          ! condensat convectif. Mais on peut mettre 2 si on ne veut que
-          ! collorier que le condensat convectif.
-          call getin('option_cond',option_cond)
-          write(*,*) 'option_cond=',option_cond
-
-          strtrac(izone_poubelle)='res'
-          strtrac(izone_cond)='con'
-          strtrac(izone_revap)='rev'
-
-        elseif (option_traceurs.eq.10) then
-          ! on trace l'évap venant de ocean/continent no frac/continent frac
-          !  utilse seulement si couplé avec ORCHIDEE
-#ifdef CPP_VEGET
-#else
-          write(*,*) 'iso_traceurs_init 219: option_traceurs=10 ', &
-     &                      'inutile si on ne couple pas avec ORCHIDEE'
-          stop
+         nzone_opt=nzone_temp+1
+         zone_tempa=-4.0     ! en K
+         izone_cont=nzone
+         izone_oce=nzone  
+         izone_poubelle=nzone
+         izone_init=nzone    ! zone d'initialisation par defaut
+         option_revap=0 
+         option_tmin=0  
+         izone_revap=0
+         option_cond=0
+         DO izone=1,nzone_temp
+            strtrac(izone) = 't'//TRIM(int2str(izone))
+            WRITE(*,*) 'izone, strtrac=', izone, strtrac(izone)
+         END DO
+         strtrac(izone_poubelle)='pou'
+         ! Initialization of temperatures zones
+         DO izone=1,nzone_temp-1
+            zone_temp(izone) = zone_temp1+float(izone-1)            &
+                            * (zone_tempa*float(izone-nzone_temp+1) &
+                            + (zone_tempf-zone_temp1)/float(nzone_temp-2))
+         END DO
+         WRITE(*,*) 'iso_trac_init 183: zone_temp=', zone_temp
+      !========================================================================================================================
+      CASE(5)      !=== TRACING AEJ/MOONSOON FLUX/Harmattan
+      !========================================================================================================================
+!        WRITE*,*) 'iso_traceurs_init 129'
+         nzone_opt=4
+         izone_cont=1
+         izone_oce=1
+         izone_poubelle=1    ! zone ou on met les flux non physiques, de reajustement
+         izone_init=1        ! zone d'initialisation par defaut          
+         option_revap=0
+         option_tmin=0
+         izone_revap=0
+         izone_aej=2
+         izone_mousson=3
+         izone_harmattan=4 
+         option_cond=0
+         strtrac(izone_poubelle) = 'res'
+         strtrac(izone_aej)      = 'aej'
+         strtrac(izone_mousson)  = 'mou'
+         strtrac(izone_harmattan)= 'sah'
+      !========================================================================================================================
+      CASE(6)      !=== TRACING DDFTS
+      !========================================================================================================================
+         nzone_opt=2
+         izone_cont=1
+         izone_oce=1
+         izone_poubelle=1    ! zone ou on met les flux non physiques, de reajustement
+         izone_init=1        ! zone d'initialisation par defaut          
+         option_revap=0
+         option_tmin=0
+         izone_revap=0
+         izone_ddft=2
+         option_cond=0
+         strtrac(izone_poubelle)='res'
+         strtrac(izone_ddft)='dft'
+      !========================================================================================================================
+      CASE(9)      !=== TRACING CONDENSATION
+      !========================================================================================================================
+         nzone_opt=3
+         izone_cont=1
+         izone_oce=1
+         izone_poubelle=1    ! zone ou on met les flux non physiques, de reajustement
+         izone_init=1        ! zone d'initialisation par defaut          
+         option_revap=1
+         option_tmin=0
+         izone_revap=2
+         izone_cond=3
+         option_cond=1
+         ! 1 par defaut pour colorier a la fois condensat LS et condensat convectif.
+         ! Mais on peut mettre 2 si on ne veut que colorier que le condensat convectif.
+         CALL get_in('option_cond',option_cond)
+         strtrac(izone_poubelle)='res'
+         strtrac(izone_cond)='con'
+         strtrac(izone_revap)='rev'
+      !========================================================================================================================
+      CASE(10)     !=== TRACING EVAPORATION FROM OCEAN/LAND, NON FRAC/LAND FRAC ; ONLY WHEN COUPLED WITH ORCHIDEE
+      !========================================================================================================================
+#ifndef CPP_VEGET
+         WRITE(*,*) 'iso_traceurs_init 219: option_traceurs=10 inutile si on ne couple pas avec ORCHIDEE'; STOP
 #endif          
-
-          ntraceurs_zone_opt=3
-          izone_cont=1 ! sous-entendu non fractionnant
-          izone_oce=2
-          izone_poubelle=2 ! zone où on met les flux non physiques, de
-                ! réajustement
-          izone_init=2 ! zone d'initialisation par défaut
-          option_revap=0
-          option_tmin=0
-          izone_revap=0
-          izone_contfrac=3
-          izone_contcanop=3
-          izone_irrig=0
-          option_cond=0
-
-          strtrac(izone_oce)='oce'
-          strtrac(izone_cont)='con'  
-          strtrac(izone_contfrac)='enu'  ! evap sol nu
-
-        elseif (option_traceurs.eq.11) then
-          ! on trace reevap des gouttes et le reste
-
-          ntraceurs_zone_opt=2
-          izone_cont=1
-          izone_oce=1
-          izone_poubelle=1 ! zone où on met les flux non physiques, de
-                ! réajustement
-          izone_init=1 ! zone d'initialisation par défaut
-          option_revap=1
-          option_tmin=0
-          izone_revap=2
-          izone_irrig=0
-          option_cond=0
-
-          strtrac(izone_poubelle)='res'
-          strtrac(izone_revap)='rev'
-
-        elseif (option_traceurs.eq.12) then
-          ! on trace evap du sol nu, evap de la canopée, reste de l'evap cont et
-          ! evap oce
-#ifdef CPP_VEGET
-#else
-          write(*,*) 'iso_traceurs_init 257: option_traceurs=10 ', &
-     &                      'inutile si on ne couple pas avec ORCHIDEE'
-          stop
+         nzone_opt=3
+         izone_cont=1        ! sous-entendu non fractionnant
+         izone_oce=2
+         izone_poubelle=2    ! zone ou on met les flux non physiques, de reajustement
+         izone_init=2        ! zone d'initialisation par defaut
+         option_revap=0
+         option_tmin=0
+         izone_revap=0
+         izone_contfrac=3
+         izone_contcanop=3
+         izone_irrig=0
+         option_cond=0
+         strtrac(izone_oce)='oce'
+         strtrac(izone_cont)='con'  
+         strtrac(izone_contfrac)='enu'  ! evap sol nu
+      !========================================================================================================================
+      CASE(11)     !=== TRACING DROPLETS REEVAPORATION + REST
+      !========================================================================================================================
+         nzone_opt=2
+         izone_cont=1
+         izone_oce=1
+         izone_poubelle=1    ! zone ou on met les flux non physiques, de reajustement
+         izone_init=1        ! zone d'initialisation par defaut
+         option_revap=1
+         option_tmin=0
+         izone_revap=2
+         izone_irrig=0
+         option_cond=0
+         strtrac(izone_poubelle)='res'
+         strtrac(izone_revap)='rev'
+      !========================================================================================================================
+      CASE(12)     !=== TRACING NAKED GROUND EVAPORATION, CANOPY EVAPORATION, REST OF LAND EVAPORATION AND OCEAN EVAPORATION
+      !========================================================================================================================
+#ifndef CPP_VEGET
+         WRITE(*,*) 'iso_traceurs_init 257: option_traceurs=10 inutile si on ne couple pas avec ORCHIDEE'; STOP
 #endif           
-
-          ntraceurs_zone_opt=2
-          izone_cont=1
-          izone_oce=2
-          izone_poubelle=2 ! zone où on met les flux non physiques, de
-                ! réajustement
-          izone_init=2 ! zone d'initialisation par défaut
-          option_revap=0
-          option_tmin=0
-          izone_revap=0
-          izone_contfrac=3
-          izone_contcanop=4
-          izone_irrig=0   
-          option_cond=0
-
-          strtrac(izone_oce)='oce'
-          strtrac(izone_cont)='con'
-          strtrac(izone_contfrac)='enu'  ! evap sol nu
-          strtrac(izone_contcanop)='eca'  ! evap canop
-
-       else if (option_traceurs.eq.13) then
-          ! on trace les température minimales vécues + la revap
-          ! comme dans article sur LdG
-            
-        zone_temp1=293.0         ! en K        
-!        parameter (zone_tempf=223.0) ! en K
-        zone_tempf=243.0 ! en K
-        zone_tempa=-4.0 ! courbure de la relation entre l'indice et la température: 0 pour linéaire, <0 pour plus de détal en bas
-
-        ! zone 1: >= zone_temp1
-        ! zone 2 à 4: intermédiaire, 
-        ! zone 5: <zone_tempf
-        
-          ntraceurs_zone_opt=nzone_temp+1
-          
-          izone_cont=1
-          izone_oce=1  
-          izone_poubelle=1
-          izone_init=1 ! zone d'initialisation par défaut
-          option_revap=1   
-          option_tmin=0
-          izone_revap=ntraceurs_zone
-          izone_irrig=0
-          option_cond=0
-          do izone=1,nzone_temp
-            write(strz,'(i2.2)') izone
-            strtrac(izone)='t'//strz
-            write(*,*) 'izone,strz,strtrac=',izone,strz,strtrac(izone)
-          enddo
-          strtrac(izone_revap)='rev'
-
-          ! initialisation des zones de tempéarture
-          do izone=1,nzone_temp-1
-            zone_temp(izone)=zone_temp1+float(izone-1) &
-     &                      *(zone_tempa*float(izone-nzone_temp+1) &
-     &                      +(zone_tempf-zone_temp1)/float(nzone_temp-2))
-          enddo
-          write(*,*) 'zone_temp=',zone_temp
-
-       else if (option_traceurs.eq.14) then
-          ! on trace les pres et lat de dernière saturation définies
-          ! comme rh>90%
-            
-        zone_pres1=600.0*100.0 ! en Pa       
-        zone_presf=300.0*100.0 ! en Pa
-        zone_presa=0.0 ! courbure de la relation entre l'indice et la température: 0 pour linéaire, <0 pour plus de détal en bas
-
-        lattag_min=10.0 ! en degrès
-        dlattag=15.0
-
-        ! zone 1: >= zone_pres1
-        ! zone 2 à 4: intermédiaire, 
-        ! zone 5: <zone_presf
-        
-         ntraceurs_zone_opt=nzone_pres*nzone_lat+1          
-          izone_cont=ntraceurs_zone
-          izone_oce=ntraceurs_zone
-          izone_poubelle=ntraceurs_zone
-          izone_init=ntraceurs_zone ! zone d'initialisation par défaut
-          option_revap=0  
-          option_tmin=0 
-          izone_revap=0
-          izone_irrig=0
-          option_cond=0
-          do izone_pres=1,nzone_pres
-           do izone_lat=1,nzone_lat
-            write(strz_pres,'(i1.1)') izone_pres
-            write(strz_lat,'(i1.1)') izone_lat
-            strz_preslat=strz_pres//strz_lat
-            izone=izone_lat+(izone_pres-1)*nzone_lat
-            strtrac(izone)='t'//strz_preslat
-            write(*,*) 'izone_pres,izone_lat,strtrac=', &
-     &                        izone_pres,izone_lat,izone,strtrac(izone)
-           enddo !do izone_lat=1,nzone_lat
-          enddo !do izone_pres=1,nzone_pres
-          strtrac(ntraceurs_zone)='sfc'
-
-          ! initialisation des zones de tempéarture
-          do izone=1,nzone_pres-1
-            zone_pres(izone)=zone_pres1+float(izone-1) &
-     &                      *(zone_presa*float(izone-nzone_pres+1) &
-     &                      +(zone_presf-zone_pres1)/float(nzone_pres-2))
-          enddo !do izone=1,nzone_pres-1
-          write(*,*) 'traceurs_init 332: zone_pres=',zone_pres
-!          stop
-!
-       elseif (option_traceurs.eq.15) then
-          ! on trace l'irrigation dans ORCHIDEE
-#ifdef CPP_VEGET
-#else
-          write(*,*) 'iso_traceurs_init 257: option_traceurs=15 ', &
-     &                      'inutile si on ne couple pas avec ORCHIDEE'
-          stop
+         nzone_opt=2
+         izone_cont=1
+         izone_oce=2
+         izone_poubelle=2    ! zone ou on met les flux non physiques, de reajustement
+         izone_init=2        ! zone d'initialisation par defaut
+         option_revap=0
+         option_tmin=0
+         izone_revap=0
+         izone_contfrac=3
+         izone_contcanop=4
+         izone_irrig=0   
+         option_cond=0
+         strtrac(izone_oce)='oce'
+         strtrac(izone_cont)='con'
+         strtrac(izone_contfrac)='enu' ! evap sol nu
+         strtrac(izone_contcanop)='eca'! evap canop
+      !========================================================================================================================
+      CASE(13)     !=== TRACING MINIMUM EXPERIENCED TEMPERATIRES + REEVAPORATION AS IN THE ARTICLE ON LdG
+      !========================================================================================================================
+         zone_temp1=293.0    ! en K        
+!        zone_tempf=223.0    ! en K
+         zone_tempf=243.0    ! en K
+         zone_tempa=-4.0     ! courbure de la relation entre l'indice et la temperature: 0 pour lineaire, <0 pour plus de detal en bas
+         ! zone 1: >= zone_temp1
+         ! zone 2 a 4: intermediaire, 
+         ! zone 5: <zone_tempf
+         nzone_opt=nzone_temp+1
+         izone_cont=1
+         izone_oce=1  
+         izone_poubelle=1
+         izone_init=1        ! zone d'initialisation par defaut
+         option_revap=1   
+         option_tmin=0
+         izone_revap=nzone
+         izone_irrig=0
+         option_cond=0
+         DO izone=1,nzone_temp
+            strtrac(izone) = 't'//TRIM(int2str(izone))
+            WRITE(*,*) 'izone, strtrac = ', izone, strtrac(izone)
+         END DO
+         strtrac(izone_revap)='rev'
+         ! initialisation des zones de tempearture
+         DO izone=1,nzone_temp-1
+            zone_temp(izone) = zone_temp1+float(izone-1) &
+                             *(zone_tempa*float(izone-nzone_temp+1) &
+                             +(zone_tempf-zone_temp1)/float(nzone_temp-2))
+         END DO
+         WRITE(*,*) 'zone_temp=',zone_temp
+      !========================================================================================================================
+      CASE(14)     !=== TRACING PRES AND LAT OF LAST SATURATION DEFINED AS rh>90%
+      !========================================================================================================================
+         zone_pres1=600.0*100.0   ! en Pa       
+         zone_presf=300.0*100.0   ! en Pa
+         zone_presa=0.0           ! courbure de la relation entre l'indice et la temperature: 0 pour lineaire
+         lattag_min=10.0          ! en degres
+         dlattag=15.0
+         ! zone 1: >= zone_pres1
+         ! zone 2 a 4: intermediaire, 
+         ! zone 5: <zone_presf
+         nzone_opt=nzone_pres*nzone_lat+1          
+         izone_cont=nzone
+         izone_oce=nzone
+         izone_poubelle=nzone
+         izone_init=nzone         ! zone d'initialisation par defaut
+         option_revap=0  
+         option_tmin=0 
+         izone_revap=0
+         izone_irrig=0
+         option_cond=0
+         DO izone_pres=1,nzone_pres
+            DO izone_lat=1,nzone_lat
+               izone=izone_lat+(izone_pres-1)*nzone_lat
+               strtrac(izone) = 't'//TRIM(int2str(izone_pres))//TRIM(int2str(izone_lat))
+               write(*,*) 'izone_pres, izone_lat, izone, strtrac = ',izone_pres, izone_lat, izone, strtrac(izone)
+            END DO
+         END DO
+         strtrac(nzone)='sfc'
+         ! initialisation des zones de temperature
+         DO izone=1,nzone_pres-1
+            zone_pres(izone) = zone_pres1+float(izone-1) &
+                             *(zone_presa*float(izone-nzone_pres+1) &
+                             +(zone_presf-zone_pres1)/float(nzone_pres-2))
+         END DO
+         WRITE(*,*) 'traceurs_init 332: zone_pres=',zone_pres
+      !========================================================================================================================
+      CASE(15)     !=== TRACING IRRIGATION IN ORCHIDEE
+      !========================================================================================================================
+#ifndef CPP_VEGET
+         WRITE(*,*) 'iso_traceurs_init 257: option_traceurs=15 inutile si on ne couple pas avec ORCHIDEE'; STOP
 #endif
-
-          ntraceurs_zone_opt=1
-          izone_cont=1
-          izone_oce=1
-          izone_poubelle=1 ! zone où on met les flux non physiques, de
-                ! réajustement
-          izone_init=1 ! zone d'initialisation par défaut
-          option_revap=0
-          option_tmin=0
-          izone_revap=0
-          izone_contfrac=0
-          izone_contcanop=0 
-          izone_irrig=2 
-          option_cond=0
-          
-          strtrac(izone_poubelle)='res'
-          strtrac(izone_irrig)='irrig'
-
-          ! dans ce cas particulier, il y a des traceurs dans ORCHIDEE
-          ntracisoOR=ntiso
-
-        else if ((option_traceurs.eq.17).or. &
-     &           (option_traceurs.eq.18)) then
-          ! on trace les température minimales vécues
-          ! comme dans article sur LdG sauf pas de revap
-            
-        zone_temp1=12.0e-3 ! en kg/kg        
-        zone_tempf=0.2e-3 ! en kg/kg
-        zone_tempa=1.2e-3 ! courbure de la relation entre l'indice et la température: 0 pour linéaire, <0 pour plus de détail en bas
-
-!	parameter (zone_temp1=14.0e-3) ! en kg/kg        
-!   	parameter (zone_tempf=0.2e-3) ! en kg/kg
-!	parameter (zone_tempa=0.5e-3)       
-
-!        parameter (zone_temp1=10.0e-3) ! en kg/kg
-!   	parameter (zone_tempf=0.5e-3) ! en kg/kg
-!	parameter (zone_tempa=0.5e-3)
-
-        ! zone 1: >= zone_temp1
-        ! zone 2 à 4: intermédiaire, 
-        ! zone 5: <zone_tempf
-        
-        ntraceurs_zone_opt=nzone_temp+3
-        
-          izone_cont=nzone_temp+1
-          izone_oce=nzone_temp+1 
-          izone_poubelle=nzone_temp+1
-          izone_init=nzone_temp+1 ! zone d'initialisation par défaut
-          option_revap=1  
-          option_tmin=1 
-          option_cond=1
-
-          izone_revap=nzone_temp+3
-          izone_cond=nzone_temp+2
-          do izone=1,nzone_temp
-            write(strz,'(i2.2)') izone
-            strtrac(izone)='t'//strz
-            write(*,*) 'izone,strz,strtrac=',izone,strz,strtrac(izone)
-          enddo !do izone=1,nzone_temp
-          strtrac(izone_poubelle)='sfc'
-          strtrac(izone_cond)='con'
-          strtrac(izone_revap)='rev'
-
-          ! initialisation des zones de tempéarture
-          do izone=1,nzone_temp-1
-            zone_temp(izone)=zone_temp1+float(izone-1) &
-     &                      *(zone_tempa*float(izone-nzone_temp+1) &
-     &             +(zone_tempf-zone_temp1)/float(nzone_temp-2))
-          enddo
-         write(*,*) 'zone_temp1,zone_tempf,zone_tempa=', &
-     &              zone_temp1,zone_tempf,zone_tempa
-          write(*,*) 'zone_temp=',zone_temp
-!          stop          
-
-        else if (option_traceurs.eq.19) then
-
-        zone_temp1=12.0e-3 ! en kg/kg        
-        zone_tempf=0.2e-3 ! en kg/kg
-        zone_tempa=1.2e-3 ! courbure de la relation entre l'indice et la température: 0 pour linéaire, <0 pour plus de détail en bas
-
-!	parameter (zone_temp1=14.0e-3) ! en kg/kg        
-!   	parameter (zone_tempf=0.2e-3) ! en kg/kg
-!	parameter (zone_tempa=0.5e-3)       
-
-!        parameter (zone_temp1=10.0e-3) ! en kg/kg
-!   	parameter (zone_tempf=0.5e-3) ! en kg/kg
-!	parameter (zone_tempa=0.5e-3)
-
-        ! zone 1: >= zone_temp1
-        ! zone 2 à 4: intermédiaire, 
-        ! zone 5: <zone_tempf
-        
-        ntraceurs_zone_opt=nzone_temp+4
-        
-          izone_cont=nzone_temp+1
-          izone_oce=nzone_temp+1 
-          izone_poubelle=nzone_temp+1
-          if (option_seuil_tag_tmin.eq.1) then
-            izone_init=nzone_temp+1 ! zone d'initialisation par défaut
-          else
+         nzone_opt=1
+         izone_cont=1
+         izone_oce=1
+         izone_poubelle=1    ! zone ou on met les flux non physiques, de reajustement
+         izone_init=1        ! zone d'initialisation par defaut
+         option_revap=0
+         option_tmin=0
+         izone_revap=0
+         izone_contfrac=0
+         izone_contcanop=0 
+         izone_irrig=2 
+         option_cond=0
+         strtrac(izone_poubelle)='res'
+         strtrac(izone_irrig)='irrig'
+         ! dans ce cas particulier, il y a des traceurs dans ORCHIDEE
+         ntracisoOR=ntiso
+      !========================================================================================================================
+      CASE(17,18)  !=== TRACING MINIMAL EXPERIENCES TEMPERATURES AS IN THE ARTICLE ABOUT LdG, BUT NO EVAPORATION
+      !========================================================================================================================
+         zone_temp1=12.0e-3  ! en kg/kg        
+         zone_tempf=0.2e-3   ! en kg/kg
+         zone_tempa=1.2e-3   ! courbure de la relation entre l'indice et la temperature: 0 pour lineaire
+!	 zone_temp1=14.0e-3  ! en kg/kg        
+!   	 zone_tempf=0.2e-3   ! en kg/kg
+!	 zone_tempa=0.5e-3       
+!        zone_temp1=10.0e-3  ! en kg/kg
+!   	 zone_tempf=0.5e-3   ! en kg/kg
+!	 zone_tempa=0.5e-3
+         ! zone 1: >= zone_temp1
+         ! zone 2 a 4: intermediaire, 
+         ! zone 5: <zone_tempf
+         nzone_opt=nzone_temp+3
+         izone_cont=nzone_temp+1
+         izone_oce=nzone_temp+1 
+         izone_poubelle=nzone_temp+1
+         izone_init=nzone_temp+1 ! zone d'initialisation par defaut
+         option_revap=1  
+         option_tmin=1 
+         option_cond=1
+         izone_revap=nzone_temp+3
+         izone_cond=nzone_temp+2
+         DO izone=1,nzone_temp
+            strtrac(izone) = 't'//TRIM(int2str(izone))
+            WRITE(*,*) 'izone, strtrac = ', izone, strtrac(izone)
+         END DO !do izone=1,nzone_temp
+         strtrac(izone_poubelle)='sfc'
+         strtrac(izone_cond)='con'
+         strtrac(izone_revap)='rev'
+         ! initialisation des zones de tempearture
+         DO izone=1,nzone_temp-1
+            zone_temp(izone) = zone_temp1+float(izone-1) &
+                             *(zone_tempa*float(izone-nzone_temp+1) &
+                             +(zone_tempf-zone_temp1)/float(nzone_temp-2))
+         END DO
+         WRITE(*,*) 'zone_temp1,zone_tempf,zone_tempa=',zone_temp1,zone_tempf,zone_tempa
+         WRITE(*,*) 'zone_temp=',zone_temp
+!        STOP          
+      !========================================================================================================================
+      CASE(19)     !=== TRACING TROPICAL AND EXTRATROPICAL VAPOUR
+      !========================================================================================================================
+         zone_temp1=12.0e-3  ! en kg/kg        
+         zone_tempf=0.2e-3   ! en kg/kg
+         zone_tempa=1.2e-3   ! courbure de la relation entre l'indice et la temperature: 0 pour lineaire, <0 pour plus de detail en bas
+!        zone_temp1=14.0e-3  ! en kg/kg        
+!        zone_tempf=0.2e-3   ! en kg/kg
+!        zone_tempa=0.5e-3
+!        zone_temp1=10.0e-3  ! en kg/kg        
+!        zone_tempf=0.5e-3   ! en kg/kg
+!        zone_tempa=0.5e-3
+         ! zone 1: >= zone_temp1
+         ! zone 2 a 4: intermediaire, 
+         ! zone 5: <zone_tempf
+         nzone_opt=nzone_temp+4
+         izone_cont=nzone_temp+1
+         izone_oce=nzone_temp+1 
+         izone_poubelle=nzone_temp+1
+         IF(option_seuil_tag_tmin == 1) THEN
+            izone_init=nzone_temp+1 ! zone d'initialisation par defaut
+         ELSE
             izone_init=nzone_temp
-          endif
-          option_revap=1   
-          izone_revap=nzone_temp+3
-          izone_cond=nzone_temp+2
-          izone_ddft=nzone_temp+4
-          option_tmin=1          
-          option_cond=1
-          do izone=1,nzone_temp
-            write(strz,'(i2.2)') izone
-            strtrac(izone)='t'//strz
-            write(*,*) 'izone,strz,strtrac=',izone,strz,strtrac(izone)
-          enddo !do izone=1,nzone_temp
-          strtrac(izone_poubelle)='sfc'
-          strtrac(izone_cond)='con'
-          strtrac(izone_revap)='rev'
-          strtrac(izone_ddft)='dft'
-
-        elseif (option_traceurs.eq.20) then
-          ! on vapeur tropical/extractropicale/recyclage extractropical
-          ! pour comprendre controles humidité et isotopes subtropicaux.        
-          
-          lim_tag20=35.0
-          call getin('lim_tag20',lim_tag20)
-          write(*,*) 'lim_tag20=',lim_tag20
-
-          ntraceurs_zone_opt=3
-          izone_cont=1
-          izone_oce=1
-          izone_poubelle=2 ! zone où on met les flux non physiques, de
-                ! réajustement
-          izone_init=2 ! zone d'initialisation par défaut
-          option_revap=0
-          option_tmin=0
-          izone_revap=0
-          izone_trop=2
-          izone_extra=3
-
-          strtrac(izone_trop)='tro' ! vapeur tropicale
-          strtrac(izone_extra)='ext' ! vapeur extractropicale evaporée 
-                ! dans les tropiques
-          strtrac(izone_cont)='rec' ! recyclage
-
-        elseif (option_traceurs.eq.21) then
-          ! on trace 2 boites 3D: UT tropicale et extratropiques 
-          ! fonctionnement similaire à option 5 pour taggage des zones
-          ! AMMA
-!          write(*,*) 'iso_traceurs_init 129'
-
-          ntraceurs_zone_opt=3
-          izone_cont=1
-          izone_oce=1
-          izone_poubelle=1 ! zone où on met les flux non physiques, de
-                ! réajustement
-          izone_init=1 ! zone d'initialisation par défaut
-          option_revap=0
-          option_tmin=0
-          izone_revap=0
-          izone_trop=2
-          izone_extra=3
-          option_cond=0
-
-          strtrac(izone_poubelle)='res'
-          strtrac(izone_trop)='tro'
-          strtrac(izone_extra)='ext'
-
-        elseif (option_traceurs.eq.22) then
-          ! on trace la vapeur qui a été processée dans zones de
-          ! convections à 3 niveaux: BT, MT et UT
-
-          lim_precip_tag22=20.0
-          call getin('lim_precip_tag22',lim_precip_tag22)
-          write(*,*) 'lim_precip_tag22=',lim_precip_tag22
-
-          ntraceurs_zone_opt=3
-          izone_cont=1
-          izone_oce=1
-          izone_poubelle=1 ! zone où on met les flux non physiques, de
-                ! réajustement
-          izone_init=1 ! zone d'initialisation par défaut
-          option_revap=0
-          option_tmin=0
-          izone_revap=0
-          izone_conv_BT=2
-          izone_conv_UT=3
-          option_cond=0
-
-          strtrac(izone_poubelle)='res'
-          strtrac(izone_conv_BT)='cbt'
-          strtrac(izone_conv_UT)='cut'
-
-        else
-            write(*,*) 'traceurs_init 36: option pas encore prévue'
-            stop
-        endif
-
-        
-          if (ntraceurs_zone_opt.ne.ntraceurs_zone) then
-                write(*,*) 'ntraceurs_zone_opt,ntraceurs_zone=', &
-                        & ntraceurs_zone_opt,ntraceurs_zone
-                call abort_physic ('isotrac_mod','ntraceurs_zone incoherent',1)
-          endif
-
-       
-        ! seuil sur le taux de condensation
-        if (option_tmin.eq.1) then
-          seuil_tag_tmin=0.01
-          call getin('seuil_tag_tmin',seuil_tag_tmin)
-          write(*,*) 'seuil_tag_tmin=',seuil_tag_tmin
-
-          seuil_tag_tmin_ls=seuil_tag_tmin
-          call getin('seuil_tag_tmin_ls',seuil_tag_tmin_ls)
-          write(*,*) 'seuil_tag_tmin_ls=',seuil_tag_tmin_ls
-
-          option_seuil_tag_tmin=1
-          call getin('option_seuil_tag_tmin',option_seuil_tag_tmin)
-          write(*,*) 'option_seuil_tag_tmin=',option_seuil_tag_tmin
-        endif
-
-
-        do ixt=1,niso
-           index_zone(ixt)=0
-           index_iso(ixt)=ixt
-        enddo
-        itrac=niso        
-        do izone=1,ntraceurs_zone
-          do ixt=1,niso
-            itrac=itrac+1
-            index_zone(itrac)=izone
-            index_iso(itrac)=ixt
-            itZonIso_loc(izone,ixt)=itrac
-            if (itZonIso(izone,ixt).ne.itZonIso_loc(izone,ixt)) then
-                write(*,*) 'isotrac 989: izone,ixt,itrac=',izone,ixt,itrac
-                CALL abort_physic ('isotrac','isotrac 989',1)
-            endif
-          enddo
-        enddo
+         END IF
+         option_revap=1   
+         izone_revap=nzone_temp+3
+         izone_cond=nzone_temp+2
+         izone_ddft=nzone_temp+4
+         option_tmin=1          
+         option_cond=1
+         DO izone=1,nzone_temp
+            strtrac(izone) = 't'//TRIM(int2str(izone))
+            WRITE(*,*) 'izone, strtrac = ', izone, strtrac(izone)
+         END DO
+         strtrac(izone_poubelle)='sfc'
+         strtrac(izone_cond)='con'
+         strtrac(izone_revap)='rev'
+         strtrac(izone_ddft)='dft'
+      !========================================================================================================================
+      CASE(20)     !=== TRACING TROPICAL/EXTRATROPICAL/EXTRATROPICAL RECYCLING TO STUDY HUMIDITY AND SUBTROPICAL ISOTOPES CONTROL
+      !========================================================================================================================
+         CALL get_in('lim_tag20', lim_tag20, 35.0)
+         nzone_opt=3
+         izone_cont=1
+         izone_oce=1
+         izone_poubelle=2    ! zone ou on met les flux non physiques, de reajustement
+         izone_init=2        ! zone d'initialisation par defaut
+         option_revap=0
+         option_tmin=0
+         izone_revap=0
+         izone_trop=2
+         izone_extra=3
+         strtrac(izone_trop)='tro'     ! tropical vapour
+         strtrac(izone_extra)='ext'    ! extratropical vapour evaporated in the tropics
+         strtrac(izone_cont)='rec'     ! recycling
+      !========================================================================================================================
+      CASE(21)     !=== TRACING TWO 3D BOXES: TROPICAL UT AND EXTRATROPICS ; SIMILAR TO 5 FOR AMMA ZONES TAGGING
+      !========================================================================================================================
+!        WRITE(*,*) 'iso_traceurs_init 129'
+         nzone_opt=3
+         izone_cont=1
+         izone_oce=1
+         izone_poubelle=1    ! zone ou on met les flux non physiques, de reajustement
+         izone_init=1        ! zone d'initialisation par defaut
+         option_revap=0
+         option_tmin=0
+         izone_revap=0
+         izone_trop=2
+         izone_extra=3
+         option_cond=0
+         strtrac(izone_poubelle)='res'
+         strtrac(izone_trop)='tro'
+         strtrac(izone_extra)='ext'
+      !========================================================================================================================
+      CASE(22)     !=== TRACING WATER VAPOUR PROCESSED IN THE 3-LEVELS SCONVECTION ZONES BT, MT AND UT
+      !========================================================================================================================
+         CALL get_in('lim_precip_tag22', lim_precip_tag22, 20.0)
+         nzone_opt=3
+         izone_cont=1
+         izone_oce=1
+         izone_poubelle=1    ! zone ou on met les flux non physiques, de reajustement
+         izone_init=1        ! zone d'initialisation par defaut
+         option_revap=0
+         option_tmin=0
+         izone_revap=0
+         izone_conv_BT=2
+         izone_conv_UT=3
+         option_cond=0
+         strtrac(izone_poubelle)='res'
+         strtrac(izone_conv_BT)='cbt'
+         strtrac(izone_conv_UT)='cut'
+      CASE DEFAULT
+         WRITE(*,*) 'traceurs_init 36: option pas encore prevue' ; STOP
+   END SELECT
+
+   IF(nzone_opt /= nzone) THEN
+      WRITE(*,*) 'nzone_opt, nzone=', nzone_opt, nzone
+      CALL abort_physic ('isotrac_mod','nzone incoherent',1)
+   END IF
+
+   !--- Condensation rate threshold
+   IF(option_tmin == 1) THEN
+      seuil_tag_tmin = 0.01
+      CALL get_in('seuil_tag_tmin',        seuil_tag_tmin,        0.01)
+      CALL get_in('seuil_tag_tmin_ls',     seuil_tag_tmin_ls,     seuil_tag_tmin)
+      CALL get_in('option_seuil_tag_tmin', option_seuil_tag_tmin, 1)
+   END IF
+
+   index_zone = [(strIdx(isoZone, strTail(isoName(ixt) ,'_',.TRUE.)), ixt=1, ntiso)]
+   index_iso  = [(strIdx(isoName, strHead(isoName(ixt) ,'_',.TRUE.)), ixt=1, ntiso)]
+   itZonIso_loc = itZonIso(:,:)
 #ifdef ISOVERIF
-!        call iso_verif_egalite(float(itrac),float(ntiso), &
-!     &           'traceurs_init 50')
-        if (itrac.ne.ntiso) then
-          write(*,*) 'traceurs_init 50'
-          stop
-        endif
-     
-        write(*,*) 'traceurs_init 65: bilan de l''init:'
-        write(*,*) 'index_zone=',index_zone(1:ntiso)
-        write(*,*) 'index_iso=',index_iso(1:ntiso)
-        write(*,*) 'itZonIso=',itZonIso(1:ntraceurs_zone,1:niso)
-        do izone=1,ntraceurs_zone
-          write(*,*) 'strtrac(',izone,')=',strtrac(izone)
-        enddo !do izone=1,ntraceurs_zone
-        write(*,*) 'ntracisoOR=',ntracisoOR
+   WRITE(*,*) 'traceurs_init 65: bilan de l''init:'
+   WRITE(*,*) 'index_zone = '//TRIM(strStack(int2str(index_zone(1:ntiso))))
+   WRITE(*,*) 'index_iso  = '//TRIM(strStack(int2str(index_iso (1:ntiso))))
+   DO izone=1,nzone
+      WRITE(*,*)'itZonIso('//TRIM(int2str(izone))//',:) = '//strStack(int2str(itZonIso(izone,:)))
+   END DO
+   DO izone=1,nzone
+      WRITE(*,*)'strtrac('//TRIM(int2str(izone))//',:) = '//TRIM(strtrac(izone))
+   END DO
+   WRITE(*,*) 'ntracisoOR=',ntracisoOR
 #endif  
 
-        end subroutine iso_traceurs_init
-
+END SUBROUTINE iso_traceurs_init
 
 END MODULE isotrac_mod
Index: LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/isotrac_routines_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/isotrac_routines_mod.F90	(revision 4175)
+++ LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/isotrac_routines_mod.F90	(revision 4482)
@@ -1113,10 +1113,10 @@
       USE isotrac_mod, only: use_bassin_atlantic,use_bassin_medit, &
 &       use_bassin_indian,use_bassin_austral,use_bassin_pacific, &
-&       use_bassin_merarabie,use_bassin_golfebengale,use_bassin_indiansud, &
-&       use_bassin_tropics,use_bassin_midlats,use_bassin_hauteslats, &
+&       use_bassin_MerArabie,use_bassin_BengalGolf,use_bassin_SouthIndian, &
+&       use_bassin_tropics,use_bassin_midlats,use_bassin_HighLats, &
 &       bassin_atlantic,bassin_medit, &
 &       bassin_indian,bassin_austral,bassin_pacific, &
-&       bassin_merarabie,bassin_golfebengale,bassin_indiansud, &
-&       bassin_tropics,bassin_midlats,bassin_hauteslats
+&       bassin_MerArabie,bassin_BengalGolf,bassin_SouthIndian, &
+&       bassin_tropics,bassin_midlats,bassin_HighLats
       implicit none
       ! répond true si lat,lon se trouve dans le bassin numéroté bassin
@@ -1137,6 +1137,5 @@
       write(*,*) 'is_in_basin 84: entree,bassin=',bassin
 #endif
-      if ((use_bassin_atlantic.eq.1).and. &
-     &           (bassin.eq.bassin_atlantic)) then
+      if (use_bassin_atlantic .and. bassin==bassin_atlantic) then
 #ifdef ISOVERIF           
           write(*,*) 'bassin Atlantique?'
@@ -1169,6 +1168,5 @@
           endif
 
-      else if ((use_bassin_medit.eq.1).and. &
-     &           (bassin.eq.bassin_medit)) then
+      else if (use_bassin_medit .and. bassin==bassin_medit) then
 #ifdef ISOVERIF           
           write(*,*) 'bassin Medit?'
@@ -1183,6 +1181,5 @@
           endif
 
-      else if ((use_bassin_indian.eq.1).and. &
-     &           (bassin.eq.bassin_indian)) then
+      else if (use_bassin_indian .and. bassin==bassin_indian) then
 #ifdef ISOVERIF           
           write(*,*) 'bassin indian?'
@@ -1199,6 +1196,5 @@
           endif   
 
-      else if ((use_bassin_indiansud.eq.1).and. &
-     &           (bassin.eq.bassin_indiansud)) then
+      else if (use_bassin_SouthIndian .and. bassin==bassin_SouthIndian) then
 #ifdef ISOVERIF           
           write(*,*) 'bassin indian hemisphere Sud?'
@@ -1209,6 +1205,5 @@
           endif
           
-      else if ((use_bassin_merarabie.eq.1).and. &
-     &           (bassin.eq.bassin_merarabie)) then
+      else if (use_bassin_MerArabie .and. bassin==bassin_MerArabie) then
 #ifdef ISOVERIF           
           write(*,*) 'bassin Mer d''Arabie?'
@@ -1219,6 +1214,5 @@
           endif 
 
-      else if ((use_bassin_golfebengale.eq.1).and. &
-     &           (bassin.eq.bassin_golfebengale)) then
+      else if (use_bassin_BengalGolf .and. bassin==bassin_BengalGolf) then
 #ifdef ISOVERIF           
           write(*,*) 'bassin Golfe du Bengale?'
@@ -1229,6 +1223,5 @@
           endif         
 
-      else if ((use_bassin_pacific.eq.1).and. &
-     &           (bassin.eq.bassin_pacific)) then
+      else if (use_bassin_pacific .and. bassin==bassin_pacific) then
 #ifdef ISOVERIF           
           write(*,*) 'bassin Pacific?'
@@ -1267,6 +1260,5 @@
           endif
 
-      else if ((use_bassin_austral.eq.1).and. &
-     &           (bassin.eq.bassin_austral)) then  
+      else if (use_bassin_austral .and. bassin==bassin_austral) then  
 #ifdef ISOVERIF           
           write(*,*) 'bassin austral?'
@@ -1277,6 +1269,5 @@
           endif  
 
-      else if ((use_bassin_hauteslats.eq.1).and. &
-     &           (bassin.eq.bassin_hauteslats)) then  
+      else if (use_bassin_HighLats .and. bassin==bassin_HighLats) then  
 #ifdef ISOVERIF           
           write(*,*) 'bassin hautes lats?'
@@ -1287,6 +1278,5 @@
           endif
 
-      else if ((use_bassin_tropics.eq.1).and. &
-     &           (bassin.eq.bassin_tropics)) then  
+      else if (use_bassin_tropics .and. bassin==bassin_tropics) then  
 #ifdef ISOVERIF           
           write(*,*) 'bassin tropics?'
@@ -1297,6 +1287,5 @@
           endif
 
-       else if ((use_bassin_midlats.eq.1).and. &
-     &           (bassin.eq.bassin_midlats)) then  
+       else if (use_bassin_midlats .and. bassin==bassin_midlats) then  
 #ifdef ISOVERIF           
           write(*,*) 'bassin mid lats?'
@@ -1314,14 +1303,14 @@
           write(*,*) 'bassin_indian=' ,bassin_indian
           write(*,*) 'bassin_austral=' ,bassin_austral
-          write(*,*) 'bassin_merarabie=' ,bassin_merarabie
-          write(*,*) 'bassin_golfebengale=' ,bassin_golfebengale
-          write(*,*) 'bassin_indiansud=' ,bassin_indiansud
+          write(*,*) 'bassin_MerArabie=' ,bassin_MerArabie
+          write(*,*) 'bassin_BengalGolf=' ,bassin_BengalGolf
+          write(*,*) 'bassin_SouthIndian=' ,bassin_SouthIndian
           write(*,*) 'use_bassin_atlantic=' ,use_bassin_atlantic  
           write(*,*) 'use_bassin_medit=' ,use_bassin_medit
           write(*,*) 'use_bassin_indian=' ,use_bassin_indian
           write(*,*) 'use_bassin_austral=' ,use_bassin_austral
-          write(*,*) 'use_bassin_merarabie=' ,use_bassin_merarabie
-          write(*,*) 'use_bassin_golfebengale=' ,use_bassin_golfebengale
-          write(*,*) 'use_bassin_indiansud=' ,use_bassin_indiansud
+          write(*,*) 'use_bassin_MerArabie=' ,use_bassin_MerArabie
+          write(*,*) 'use_bassin_BengalGolf=' ,use_bassin_BengalGolf
+          write(*,*) 'use_bassin_SouthIndian=' ,use_bassin_SouthIndian
           stop
       endif
@@ -2333,5 +2322,6 @@
         subroutine iso_verif_traceur_jbid_vect(x,n,m)
         USE isotopes_mod, ONLY: bidouille_anti_divergence,iso_eau,ridicule
-        use isotrac_mod, only: ntraceurs_zone=>nzone
+        !use isotrac_mod, only: ntraceurs_zone=>nzone
+        USE infotrac_phy, ONLY: ntraceurs_zone=>nzone
         implicit none
         
Index: LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/lscp_ini_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/lscp_ini_mod.F90	(revision 4482)
+++ LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/lscp_ini_mod.F90	(revision 4482)
@@ -0,0 +1,1 @@
+link ../phylmd/lscp_ini_mod.F90
Index: LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/pbl_surface_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/pbl_surface_mod.F90	(revision 4175)
+++ LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/pbl_surface_mod.F90	(revision 4482)
@@ -14,5 +14,5 @@
   USE mod_grid_phy_lmdz,   ONLY : klon_glo
   USE ioipsl
-  USE surface_data,        ONLY : type_ocean, ok_veget
+  USE surface_data,        ONLY : type_ocean, ok_veget, landice_opt
   USE surf_land_mod,       ONLY : surf_land
   USE surf_landice_mod,    ONLY : surf_landice
@@ -23,4 +23,5 @@
   USE climb_wind_mod,      ONLY : climb_wind_down, climb_wind_up
   USE coef_diff_turb_mod,  ONLY : coef_diff_turb
+  USE atke_exchange_coeff_mod, ONLY :  atke_compute_km_kh
   USE ioipsl_getin_p_mod,  ONLY : getin_p
   USE cdrag_mod
@@ -406,6 +407,7 @@
 #endif
     USE ioipsl_getin_p_mod, ONLY : getin_p
-    use phys_state_var_mod, only: ds_ns, dt_ns, delta_sst, delta_sal, zsig, zmea
-    use phys_output_var_mod, only: dter, dser, tkt, tks, taur, sss
+    use phys_state_var_mod, only: ds_ns, dt_ns, delta_sst, delta_sal, dter, &
+         dser, dt_ds, zsig, zmea
+    use phys_output_var_mod, only: tkt, tks, taur, sss
 #ifdef CPP_XIOS
     USE wxios, ONLY: missing_val
@@ -1028,9 +1030,9 @@
     REAL, DIMENSION(klon)              :: yrmu0
     ! Martin
-
-    REAL, DIMENSION(klon):: ydelta_sst, ydelta_sal, yds_ns, ydt_ns, ydter, ydser, &
-         ytkt, ytks, ytaur, ysss
-    ! compression of delta_sst, delta_sal, ds_ns, dt_ns, dter, dser, tkt, tks,
-    ! taur, sss on ocean points
+    REAL, DIMENSION(klon)              :: yri0
+    REAL, DIMENSION(klon):: ydelta_sst, ydelta_sal, yds_ns, ydt_ns, ydter, &
+         ydser, ydt_ds, ytkt, ytks, ytaur, ysss
+    ! compression of delta_sst, delta_sal, ds_ns, dt_ns, dter, dser,
+    ! dt_ds, tkt, tks, taur, sss on ocean points
 
 #ifdef ISO
@@ -1277,4 +1279,5 @@
 
     ytke=0.
+    yri0(:)=0.
 !FC 
     y_treedrg=0.
@@ -1813,4 +1816,7 @@
              ydelta_sal(:knon) = delta_sal(ni(:knon))
              ydelta_sst(:knon) = delta_sst(ni(:knon))
+             ydter(:knon) = dter(ni(:knon))
+             ydser(:knon) = dser(ni(:knon))
+             ydt_ds(:knon) = dt_ds(ni(:knon))
           end if
           
@@ -1844,5 +1850,5 @@
         CALL cdrag(knon, nsrf, &
             speed, yt(:,1), yq(:,1), zgeo1, ypaprs(:,1),&
-            yts, yqsurf, yz0m, yz0h, &
+            yts, yqsurf, yz0m, yz0h, yri0, 0, &
             ycdragm, ycdragh, zri1, pref )
 
@@ -1878,5 +1884,5 @@
             CALL cdrag(knon, nsrf, &
             speed_x, yt_x(:,1), yq_x(:,1), zgeo1_x, ypaprs(:,1),&
-            yts_x, yqsurf_x, yz0m, yz0h, &
+            yts_x, yqsurf_x, yz0m, yz0h, yri0, 0, &
             ycdragm_x, ycdragh_x, zri1_x, pref_x )
 
@@ -1905,5 +1911,5 @@
         CALL cdrag(knon, nsrf, &
             speed_w, yt_w(:,1), yq_w(:,1), zgeo1_w, ypaprs(:,1),&
-            yts_w, yqsurf_w, yz0m, yz0h, &
+            yts_w, yqsurf_w, yz0m, yz0h, yri0, 0, &
             ycdragm_w, ycdragh_w, zri1_w, pref_w )
 !
@@ -1947,4 +1953,13 @@
 
        ENDIF
+
+        IF (iflag_pbl>=50) THEN
+
+        CALL atke_compute_km_kh(knon,klev,yu,yv,yt, &
+             ypplay,ypaprs,ytke,ycoefm, ycoefh)
+
+        ELSE
+
+
         CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
             ypaprs, ypplay, yu, yv, yq, yt, yts, yqsurf, ycdragm, &
@@ -1962,4 +1977,7 @@
        ENDDO
        ENDIF
+
+       ENDIF ! iflag_pbl >= 50
+
         IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh ',ycoefh
 !
@@ -1977,4 +1995,13 @@
 
        ENDIF
+
+        IF (iflag_pbl>=50) THEN
+
+        CALL atke_compute_km_kh(knon,klev,yu_x,yv_x,yt_x, &
+             ypplay,ypaprs,ytke_x,ycoefm_x, ycoefh_x)
+
+        ELSE
+
+
         CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
             ypaprs, ypplay, yu_x, yv_x, yq_x, yt_x, yts_x, yqsurf_x, ycdragm_x, &
@@ -1992,4 +2019,7 @@
        ENDDO
        ENDIF
+
+       ENDIF ! iflag_pbl >= 50
+
         IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh_x ',ycoefh_x
 !
@@ -2005,4 +2035,13 @@
       print *,' args coef_diff_turb: ytke_w ', ytke_w
        ENDIF
+
+        IF (iflag_pbl>=50) THEN
+
+        CALL atke_compute_km_kh(knon,klev,yu_w,yv_w,yt_w, &
+             ypplay,ypaprs,ytke_w,ycoefm_w, ycoefh_w)
+
+        ELSE
+
+
         CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
             ypaprs, ypplay, yu_w, yv_w, yq_w, yt_w, yts_w, yqsurf_w, ycdragm_w, &
@@ -2019,4 +2058,7 @@
        ENDDO
        ENDIF
+
+       ENDIF ! iflag_pbl >= 50
+
         IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh_w ',ycoefh_w
 !
@@ -2382,67 +2424,70 @@
        CASE(is_lic)
           ! Martin
-          CALL surf_landice(itap, dtime, knon, ni, &
-               rlon, rlat, debut, lafin, &
-               yrmu0, ylwdown, yalb, zgeo1, &
-               ysolsw, ysollw, yts, ypplay(:,1), &
-!!jyg               ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
-               ycdragh, ycdragm, yrain_f, ysnow_f, yt1, yq1,&
-               AcoefH, AcoefQ, BcoefH, BcoefQ, &
-               AcoefU, AcoefV, BcoefU, BcoefV, &
-               ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
-               ysnow, yqsurf, yqsol, yagesno, &
-               ytsoil, yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap,yfluxsens,yfluxlat, &
-               ytsurf_new, y_dflux_t, y_dflux_q, &
-               yzmea, yzsig, ycldt, &
-               ysnowhgt, yqsnow, ytoice, ysissnow, &
-               yalb3_new, yrunoff, &
-               y_flux_u1, y_flux_v1 &
-#ifdef ISO
-           &    ,yxtrain_f, yxtsnow_f,yxt1,yRland_ice &
-           &    ,yxtsnow,yxtsol,yxtevap &
+          IF (landice_opt .LT. 2) THEN
+             ! Land ice is treated by LMDZ and not by ORCHIDEE
+             
+             CALL surf_landice(itap, dtime, knon, ni, &
+                  rlon, rlat, debut, lafin, &
+                  yrmu0, ylwdown, yalb, zgeo1, &
+                  ysolsw, ysollw, yts, ypplay(:,1), &
+                  !!jyg               ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
+                  ycdragh, ycdragm, yrain_f, ysnow_f, yt1, yq1,&
+                  AcoefH, AcoefQ, BcoefH, BcoefQ, &
+                  AcoefU, AcoefV, BcoefU, BcoefV, &
+                  ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
+                  ysnow, yqsurf, yqsol, yagesno, &
+                  ytsoil, yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap,yfluxsens,yfluxlat, &
+                  ytsurf_new, y_dflux_t, y_dflux_q, &
+                  yzmea, yzsig, ycldt, &
+                  ysnowhgt, yqsnow, ytoice, ysissnow, &
+                  yalb3_new, yrunoff, &
+                  y_flux_u1, y_flux_v1 &
+#ifdef ISO
+                  &    ,yxtrain_f, yxtsnow_f,yxt1,yRland_ice &
+                  &    ,yxtsnow,yxtsol,yxtevap &
 #endif              
-           &    )
-
-!jyg<
-!!          alb3_lic(:)=0.
-!>jyg
-          DO j = 1, knon
-             i = ni(j)
-             alb3_lic(i) = yalb3_new(j)
-             snowhgt(i)   = ysnowhgt(j)
-             qsnow(i)     = yqsnow(j)
-             to_ice(i)    = ytoice(j)
-             sissnow(i)   = ysissnow(j)
-             runoff(i)    = yrunoff(j)
-          ENDDO
-          ! Martin
-! Special DICE MPL 05082013 puis BOMEX MPL 20150410
-       IF (ok_prescr_ust) THEN
-          DO j=1,knon
-          y_flux_u1(j)=ycdragm(j)*(1.+sqrt(yu(j,1)*yu(j,1)+yv(j,1)*yv(j,1)))*yu(j,1)*ypplay(j,1)/RD/yt(j,1)
-          y_flux_v1(j)=ycdragm(j)*(1.+sqrt(yu(j,1)*yu(j,1)+yv(j,1)*yv(j,1)))*yv(j,1)*ypplay(j,1)/RD/yt(j,1)
-          ENDDO
-      ENDIF
-
+                  &    )
+             
+             !jyg<
+             !!          alb3_lic(:)=0.
+             !>jyg
+             DO j = 1, knon
+                i = ni(j)
+                alb3_lic(i) = yalb3_new(j)
+                snowhgt(i)   = ysnowhgt(j)
+                qsnow(i)     = yqsnow(j)
+                to_ice(i)    = ytoice(j)
+                sissnow(i)   = ysissnow(j)
+                runoff(i)    = yrunoff(j)
+             ENDDO
+             ! Martin
+             ! Special DICE MPL 05082013 puis BOMEX MPL 20150410
+             IF (ok_prescr_ust) THEN
+                DO j=1,knon
+                   y_flux_u1(j)=ycdragm(j)*(1.+sqrt(yu(j,1)*yu(j,1)+yv(j,1)*yv(j,1)))*yu(j,1)*ypplay(j,1)/RD/yt(j,1)
+                   y_flux_v1(j)=ycdragm(j)*(1.+sqrt(yu(j,1)*yu(j,1)+yv(j,1)*yv(j,1)))*yv(j,1)*ypplay(j,1)/RD/yt(j,1)
+                ENDDO
+             ENDIF
+             
 #ifdef ISOVERIF
-        do j=1,knon
-          do ixt=1,ntraciso
-            call iso_verif_noNaN(yxtevap(ixt,j), &
-         &      'pbl_surface 1095a: apres surf_landice')
-            call iso_verif_noNaN(yxtsol(ixt,j), &
-         &      'pbl_surface 1095b: apres surf_landice')
-          enddo
-        enddo
+             do j=1,knon
+                do ixt=1,ntraciso
+                   call iso_verif_noNaN(yxtevap(ixt,j), &
+                        &      'pbl_surface 1095a: apres surf_landice')
+                   call iso_verif_noNaN(yxtsol(ixt,j), &
+                        &      'pbl_surface 1095b: apres surf_landice')
+                enddo
+             enddo
 #endif
 #ifdef ISOVERIF
-        !write(*,*) 'pbl_surface_mod 1060: sortie surf_landice'
-        do j=1,knon
-          if (iso_eau.gt.0) then     
-                 call iso_verif_egalite(yxtsnow(iso_eau,j), &
-     &                  ysnow(j),'pbl_surf_mod 1064')
-           endif !if (iso_eau.gt.0) then
-        enddo !do i=1,klon
-#endif
-          
+             !write(*,*) 'pbl_surface_mod 1060: sortie surf_landice'
+             do j=1,knon
+                if (iso_eau.gt.0) then     
+                   call iso_verif_egalite(yxtsnow(iso_eau,j), &
+                        &                  ysnow(j),'pbl_surf_mod 1064')
+                endif !if (iso_eau.gt.0) then
+             enddo !do i=1,klon
+#endif
+          END IF
        CASE(is_oce)
            CALL surf_ocean(rlon, rlat, ysolsw, ysollw, yalb_vis, &
@@ -2459,5 +2504,5 @@
                y_flux_u1, y_flux_v1, ydelta_sst(:knon), ydelta_sal(:knon), &
                yds_ns(:knon), ydt_ns(:knon), ydter(:knon), ydser(:knon), &
-               ytkt(:knon), ytks(:knon), ytaur(:knon), ysss &
+               ydt_ds(:knon), ytkt(:knon), ytks(:knon), ytaur(:knon), ysss &
 #ifdef ISO
          &      ,yxtrain_f, yxtsnow_f,yxt1,Roce, &
@@ -3391,4 +3436,9 @@
           taur(ni(:knon)) = ytaur(:knon)
           sss(ni(:knon)) = ysss(:knon)
+
+          if (activate_ocean_skin == 2 .and. type_ocean == "couple") then
+             dt_ds = missing_val
+             dt_ds(ni(:knon)) = ydt_ds(:knon)
+          end if
        end if
 
@@ -4128,5 +4178,6 @@
 
     USE indice_sol_mod
-    use phys_state_var_mod, only: delta_sal, ds_ns, dt_ns, delta_sst
+    use phys_state_var_mod, only: delta_sal, ds_ns, dt_ns, delta_sst, dter, &
+         dser, dt_ds
     use config_ocean_skin_m, only: activate_ocean_skin
 
@@ -4234,4 +4285,7 @@
                          delta_sal(i) = 0.
                          delta_sst(i) = 0.
+                         dter(i) = 0.
+                         dser(i) = 0.
+                         dt_ds(i) = 0.
                       end if
                       
Index: LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/phyetat0.F90
===================================================================
--- LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/phyetat0.F90	(revision 4175)
+++ 	(revision )
@@ -1,706 +1,0 @@
-! $Id: phyetat0.F90 3890 2021-05-05 15:15:06Z jyg $
-
-SUBROUTINE phyetat0 (fichnom, clesphy0, tabcntr0)
-
-  USE dimphy, only: klon, zmasq, klev
-  USE iophy, ONLY : init_iophy_new
-  USE ocean_cpl_mod,    ONLY : ocean_cpl_init
-  USE fonte_neige_mod,  ONLY : fonte_neige_init
-  USE pbl_surface_mod,  ONLY : pbl_surface_init
-#ifdef ISO
-  USE fonte_neige_mod,  ONLY : fonte_neige_init_iso
-  USE pbl_surface_mod,  ONLY : pbl_surface_init_iso
-#endif
-  USE surface_data,     ONLY : type_ocean, version_ocean
-  USE phys_state_var_mod, ONLY : ancien_ok, clwcon, detr_therm, phys_tstep, &
-       qsol, fevap, z0m, z0h, agesno, &
-       du_gwd_rando, du_gwd_front, entr_therm, f0, fm_therm, &
-       falb_dir, falb_dif, prw_ancien, prlw_ancien, prsw_ancien, &
-       ftsol, pbl_tke, pctsrf, q_ancien, ql_ancien, qs_ancien, radpas, radsol, rain_fall, ratqs, &
-       rnebcon, rugoro, sig1, snow_fall, solaire_etat0, sollw, sollwdown, &
-       solsw, solswfdiff, t_ancien, u_ancien, v_ancien, w01, wake_cstar, wake_deltaq, &
-       wake_deltat, wake_delta_pbl_TKE, delta_tsurf, beta_aridity, wake_fip, wake_pe, &
-       wake_s, wake_dens, awake_dens, cv_gen, zgam, zmax0, zmea, zpic, zsig, &
-#ifdef ISO
-       fxtevap, xtsol, xt_ancien, xtl_ancien, xts_ancien, wake_deltaxt, &
-       xtrain_fall,xtsnow_fall, &
-#endif
-       zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl, u10m, v10m, treedrg, &
-       ale_wake, ale_bl_stat, ds_ns, dt_ns, delta_sst, delta_sal, ratqs_inter
-!FC
-  USE geometry_mod,     ONLY: longitude_deg, latitude_deg
-  USE iostart,          ONLY: close_startphy, get_field, get_var, open_startphy
-  USE infotrac_phy,     ONLY: nqtot, nbtr, types_trac, tracers
-  USE traclmdz_mod,     ONLY: traclmdz_from_restart
-  USE carbon_cycle_mod, ONLY: carbon_cycle_tr, carbon_cycle_cpl, co2_send
-  USE indice_sol_mod,   ONLY: nbsrf, is_ter, epsfra, is_lic, is_oce, is_sic
-  USE ocean_slab_mod,   ONLY: nslay, tslab, seaice, tice, ocean_slab_init
-  USE time_phylmdz_mod, ONLY: init_iteration, pdtphys, itau_phy
-#ifdef CPP_XIOS
-  USE wxios, ONLY: missing_val
-#else
-  use netcdf, only: missing_val => nf90_fill_real
-#endif
-  use config_ocean_skin_m, only: activate_ocean_skin
-#ifdef ISO
-  USE infotrac_phy, ONLY: niso
-  USE isotopes_routines_mod, ONLY: phyisoetat0
-  USE isotopes_mod, ONLY: iso_eau
-#ifdef ISOVERIF
-  USE isotopes_verif_mod, ONLY: iso_verif_egalite_vect2D,iso_verif_egalite
-#endif
-#endif
-
-  IMPLICIT none
-  !======================================================================
-  ! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
-  ! Objet: Lecture de l'etat initial pour la physique
-  !======================================================================
-  include "dimsoil.h"
-  include "clesphys.h"
-  include "alpale.h"
-  include "compbl.h"
-  include "YOMCST.h"
-  !======================================================================
-  CHARACTER*(*) fichnom
-
-  ! les variables globales lues dans le fichier restart
-
-  REAL tsoil(klon, nsoilmx, nbsrf)
-  REAL qsurf(klon, nbsrf)
-  REAL snow(klon, nbsrf)
-  real fder(klon)
-  REAL run_off_lic_0(klon)
-  REAL fractint(klon)
-  REAL trs(klon, nbtr)
-  REAL zts(klon)
-  ! pour drag arbres FC
-  REAL drg_ter(klon,klev)
-
-  CHARACTER*6 ocean_in
-  LOGICAL ok_veget_in
-
-  INTEGER        longcles
-  PARAMETER    ( longcles = 20 )
-  REAL clesphy0( longcles )
-
-  REAL xmin, xmax
-
-  INTEGER nid, nvarid
-  INTEGER ierr, i, nsrf, isoil , k
-  INTEGER length
-  PARAMETER (length=100)
-  INTEGER it, iq, isw
-  REAL tab_cntrl(length), tabcntr0(length)
-  CHARACTER*7 str7
-  CHARACTER*2 str2
-  LOGICAL :: found,phyetat0_get,phyetat0_srf
-  REAL :: lon_startphy(klon), lat_startphy(klon)
-
-#ifdef ISO
-  REAL xtsnow(niso,klon, nbsrf)
-  REAL xtrun_off_lic_0(niso,klon)
-  REAL Rland_ice(niso,klon)
-#endif
-  ! FH1D
-  !     real iolat(jjm+1)
-  !real iolat(jjm+1-1/(iim*jjm))
-
-  ! Ouvrir le fichier contenant l'etat initial:
-
-  CALL open_startphy(fichnom)
-
-  ! Lecture des parametres de controle:
-
-  CALL get_var("controle", tab_cntrl)
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-  ! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
-  ! Les constantes de la physiques sont lues dans la physique seulement.
-  ! Les egalites du type
-  !             tab_cntrl( 5 )=clesphy0(1)
-  ! sont remplacees par
-  !             clesphy0(1)=tab_cntrl( 5 )
-  ! On inverse aussi la logique.
-  ! On remplit les tab_cntrl avec les parametres lus dans les .def
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-  DO i = 1, length
-     tabcntr0( i ) = tab_cntrl( i )
-  ENDDO
-
-  tab_cntrl(1)=pdtphys
-  tab_cntrl(2)=radpas
-
-  ! co2_ppm : value from the previous time step
-  IF (carbon_cycle_tr .OR. carbon_cycle_cpl) THEN
-     co2_ppm = tab_cntrl(3)
-     RCO2    = co2_ppm * 1.0e-06 * RMCO2 / RMD
-     ! ELSE : keep value from .def
-  ENDIF
-
-! co2_ppm0 : initial value of atmospheric CO2 (from create_etat0_limit.e .def)
-! co2_ppm0   = tab_cntrl(16)
-! initial value for interactive CO2 run when there is no tracer field for CO2 in restart
-  co2_ppm0=284.32
-
-  solaire_etat0      = tab_cntrl(4)
-  tab_cntrl(5)=iflag_con
-  tab_cntrl(6)=nbapp_rad
-
-  IF (iflag_cycle_diurne.GE.1) tab_cntrl( 7) = iflag_cycle_diurne
-  IF (soil_model) tab_cntrl( 8) =1.
-  IF (new_oliq) tab_cntrl( 9) =1.
-  IF (ok_orodr) tab_cntrl(10) =1.
-  IF (ok_orolf) tab_cntrl(11) =1.
-  IF (ok_limitvrai) tab_cntrl(12) =1.
-
-  itau_phy = tab_cntrl(15)
-
-  clesphy0(1)=tab_cntrl( 5 )
-  clesphy0(2)=tab_cntrl( 6 )
-  clesphy0(3)=tab_cntrl( 7 )
-  clesphy0(4)=tab_cntrl( 8 )
-  clesphy0(5)=tab_cntrl( 9 )
-  clesphy0(6)=tab_cntrl( 10 )
-  clesphy0(7)=tab_cntrl( 11 )
-  clesphy0(8)=tab_cntrl( 12 )
-
-  ! set time iteration
-   CALL init_iteration(itau_phy)
-
-  ! read latitudes and make a sanity check (because already known from dyn)
-  CALL get_field("latitude",lat_startphy)
-  DO i=1,klon
-    IF (ABS(lat_startphy(i)-latitude_deg(i))>=1) THEN
-      WRITE(*,*) "phyetat0: Error! Latitude discrepancy wrt startphy file:",&
-                 " i=",i," lat_startphy(i)=",lat_startphy(i),&
-                 " latitude_deg(i)=",latitude_deg(i)
-      ! This is presumably serious enough to abort run
-      CALL abort_physic("phyetat0","discrepancy in latitudes!",1)
-    ENDIF
-    IF (ABS(lat_startphy(i)-latitude_deg(i))>=0.0001) THEN
-      WRITE(*,*) "phyetat0: Warning! Latitude discrepancy wrt startphy file:",&
-                 " i=",i," lat_startphy(i)=",lat_startphy(i),&
-                 " latitude_deg(i)=",latitude_deg(i)
-    ENDIF
-  ENDDO
-
-  ! read longitudes and make a sanity check (because already known from dyn)
-  CALL get_field("longitude",lon_startphy)
-  DO i=1,klon
-    IF (ABS(lon_startphy(i)-longitude_deg(i))>=1) THEN
-      IF (ABS(360-ABS(lon_startphy(i)-longitude_deg(i)))>=1) THEN
-        WRITE(*,*) "phyetat0: Error! Longitude discrepancy wrt startphy file:",&
-                   " i=",i," lon_startphy(i)=",lon_startphy(i),&
-                   " longitude_deg(i)=",longitude_deg(i)
-        ! This is presumably serious enough to abort run
-        CALL abort_physic("phyetat0","discrepancy in longitudes!",1)
-      ENDIF
-    ENDIF
-    IF (ABS(lon_startphy(i)-longitude_deg(i))>=1) THEN
-      IF (ABS(360-ABS(lon_startphy(i)-longitude_deg(i))) > 0.0001) THEN
-        WRITE(*,*) "phyetat0: Warning! Longitude discrepancy wrt startphy file:",&
-                   " i=",i," lon_startphy(i)=",lon_startphy(i),&
-                   " longitude_deg(i)=",longitude_deg(i)
-      ENDIF 
-    ENDIF
-  ENDDO
-
-  ! Lecture du masque terre mer
-
-  CALL get_field("masque", zmasq, found)
-  IF (.NOT. found) THEN
-     PRINT*, 'phyetat0: Le champ <masque> est absent'
-     PRINT *, 'fichier startphy non compatible avec phyetat0'
-  ENDIF
-
-  ! Lecture des fractions pour chaque sous-surface
-
-  ! initialisation des sous-surfaces
-
-  pctsrf = 0.
-
-  ! fraction de terre
-
-  CALL get_field("FTER", pctsrf(:, is_ter), found)
-  IF (.NOT. found) PRINT*, 'phyetat0: Le champ <FTER> est absent'
-
-  ! fraction de glace de terre
-
-  CALL get_field("FLIC", pctsrf(:, is_lic), found)
-  IF (.NOT. found) PRINT*, 'phyetat0: Le champ <FLIC> est absent'
-
-  ! fraction d'ocean
-
-  CALL get_field("FOCE", pctsrf(:, is_oce), found)
-  IF (.NOT. found) PRINT*, 'phyetat0: Le champ <FOCE> est absent'
-
-  ! fraction glace de mer
-
-  CALL get_field("FSIC", pctsrf(:, is_sic), found)
-  IF (.NOT. found) PRINT*, 'phyetat0: Le champ <FSIC> est absent'
-
-  !  Verification de l'adequation entre le masque et les sous-surfaces
-
-  fractint( 1 : klon) = pctsrf(1 : klon, is_ter)  &
-       + pctsrf(1 : klon, is_lic)
-  DO i = 1 , klon
-     IF ( abs(fractint(i) - zmasq(i) ) .GT. EPSFRA ) THEN
-        WRITE(*, *) 'phyetat0: attention fraction terre pas ',  &
-             'coherente ', i, zmasq(i), pctsrf(i, is_ter) &
-             , pctsrf(i, is_lic)
-        WRITE(*, *) 'Je force la coherence zmasq=fractint'
-        zmasq(i) = fractint(i)
-     ENDIF
-  ENDDO
-  fractint (1 : klon) =  pctsrf(1 : klon, is_oce)  &
-       + pctsrf(1 : klon, is_sic)
-  DO i = 1 , klon
-     IF ( abs( fractint(i) - (1. - zmasq(i))) .GT. EPSFRA ) THEN
-        WRITE(*, *) 'phyetat0 attention fraction ocean pas ',  &
-             'coherente ', i, zmasq(i) , pctsrf(i, is_oce) &
-             , pctsrf(i, is_sic)
-        WRITE(*, *) 'Je force la coherence zmasq=1.-fractint'
-        zmasq(i) = 1. - fractint(i)
-     ENDIF
-  ENDDO
-
-!===================================================================
-! Lecture des temperatures du sol:
-!===================================================================
-
-  found=phyetat0_get(1,ftsol(:,1),"TS","Surface temperature",283.)
-  IF (found) THEN
-     DO nsrf=2,nbsrf
-        ftsol(:,nsrf)=ftsol(:,1)
-     ENDDO
-  ELSE
-     found=phyetat0_srf(1,ftsol,"TS","Surface temperature",283.)
-  ENDIF
-
-!===================================================================
-  ! Lecture des albedo difus et direct
-!===================================================================
-
-  DO nsrf = 1, nbsrf
-     DO isw=1, nsw
-        IF (isw.GT.99) THEN
-           PRINT*, "Trop de bandes SW"
-           call abort_physic("phyetat0", "", 1)
-        ENDIF
-        WRITE(str2, '(i2.2)') isw
-        found=phyetat0_srf(1,falb_dir(:, isw,:),"A_dir_SW"//str2//"srf","Direct Albedo",0.2)
-        found=phyetat0_srf(1,falb_dif(:, isw,:),"A_dif_SW"//str2//"srf","Direct Albedo",0.2)
-     ENDDO
-  ENDDO
-
-  found=phyetat0_srf(1,u10m,"U10M","u a 10m",0.)
-  found=phyetat0_srf(1,v10m,"V10M","v a 10m",0.)
-
-!===================================================================
-  ! Lecture des temperatures du sol profond:
-!===================================================================
-
-   DO isoil=1, nsoilmx
-        IF (isoil.GT.99) THEN
-           PRINT*, "Trop de couches "
-           call abort_physic("phyetat0", "", 1)
-        ENDIF
-        WRITE(str2,'(i2.2)') isoil
-        found=phyetat0_srf(1,tsoil(:, isoil,:),"Tsoil"//str2//"srf","Temp soil",0.)
-        IF (.NOT. found) THEN
-           PRINT*, "phyetat0: Le champ <Tsoil"//str7//"> est absent"
-           PRINT*, "          Il prend donc la valeur de surface"
-           tsoil(:, isoil, :)=ftsol(:, :)
-        ENDIF
-   ENDDO
-
-!=======================================================================
-! Lecture precipitation/evaporation
-!=======================================================================
-
-  found=phyetat0_srf(1,qsurf,"QS","Near surface hmidity",0.)
-  found=phyetat0_get(1,qsol,"QSOL","Surface hmidity / bucket",0.)
-  found=phyetat0_srf(1,snow,"SNOW","Surface snow",0.)
-  found=phyetat0_srf(1,fevap,"EVAP","evaporation",0.)
-  found=phyetat0_get(1,snow_fall,"snow_f","snow fall",0.)
-  found=phyetat0_get(1,rain_fall,"rain_f","rain fall",0.)
-
-!=======================================================================
-! Radiation
-!=======================================================================
-
-  found=phyetat0_get(1,solsw,"solsw","net SW radiation surf",0.)
-  found=phyetat0_get(1,solswfdiff,"solswfdiff","fraction of SW radiation surf that is diffuse",1.)
-  found=phyetat0_get(1,sollw,"sollw","net LW radiation surf",0.)
-  found=phyetat0_get(1,sollwdown,"sollwdown","down LW radiation surf",0.)
-  IF (.NOT. found) THEN
-     sollwdown(:) = 0. ;  zts(:)=0.
-     DO nsrf=1,nbsrf
-        zts(:)=zts(:)+ftsol(:,nsrf)*pctsrf(:,nsrf)
-     ENDDO
-     sollwdown(:)=sollw(:)+RSIGMA*zts(:)**4
-  ENDIF
-
-  found=phyetat0_get(1,radsol,"RADS","Solar radiation",0.)
-  found=phyetat0_get(1,fder,"fder","Flux derivative",0.) 
-
-
-  ! Lecture de la longueur de rugosite 
-  found=phyetat0_srf(1,z0m,"RUG","Z0m ancien",0.001)
-  IF (found) THEN
-     z0h(:,1:nbsrf)=z0m(:,1:nbsrf)
-  ELSE
-     found=phyetat0_srf(1,z0m,"Z0m","Roughness length, momentum ",0.001)
-     found=phyetat0_srf(1,z0h,"Z0h","Roughness length, enthalpy ",0.001)
-  ENDIF
-!FC
-  IF (ifl_pbltree>0) THEN
-!CALL get_field("FTER", pctsrf(:, is_ter), found)
-    treedrg(:,1:klev,1:nbsrf)= 0.0
-    CALL get_field("treedrg_ter", drg_ter(:,:), found)
-!  found=phyetat0_srf(1,treedrg,"treedrg","drag from vegetation" , 0.)
-    !lecture du profile de freinage des arbres
-    IF (.not. found ) THEN
-      treedrg(:,1:klev,1:nbsrf)= 0.0
-    ELSE
-      treedrg(:,1:klev,is_ter)= drg_ter(:,:)
-!     found=phyetat0_srf(klev,treedrg,"treedrg","freinage arbres",0.)
-    ENDIF
-  ELSE
-    ! initialize treedrg(), because it will be written in restartphy.nc
-    treedrg(:,:,:) = 0.0 
-  ENDIF
-
-  ! Lecture de l'age de la neige:
-  found=phyetat0_srf(1,agesno,"AGESNO","SNOW AGE",0.001)
-
-  ancien_ok=.true.
-  ancien_ok=ancien_ok.AND.phyetat0_get(klev,t_ancien,"TANCIEN","TANCIEN",0.)
-  ancien_ok=ancien_ok.AND.phyetat0_get(klev,q_ancien,"QANCIEN","QANCIEN",0.)
-  ancien_ok=ancien_ok.AND.phyetat0_get(klev,ql_ancien,"QLANCIEN","QLANCIEN",0.)
-  ancien_ok=ancien_ok.AND.phyetat0_get(klev,qs_ancien,"QSANCIEN","QSANCIEN",0.)
-  ancien_ok=ancien_ok.AND.phyetat0_get(klev,u_ancien,"UANCIEN","UANCIEN",0.)
-  ancien_ok=ancien_ok.AND.phyetat0_get(klev,v_ancien,"VANCIEN","VANCIEN",0.)
-  ancien_ok=ancien_ok.AND.phyetat0_get(1,prw_ancien,"PRWANCIEN","PRWANCIEN",0.)
-  ancien_ok=ancien_ok.AND.phyetat0_get(1,prlw_ancien,"PRLWANCIEN","PRLWANCIEN",0.)
-  ancien_ok=ancien_ok.AND.phyetat0_get(1,prsw_ancien,"PRSWANCIEN","PRSWANCIEN",0.)
-
-  ! Ehouarn: addtional tests to check if t_ancien, q_ancien contain
-  !          dummy values (as is the case when generated by ce0l,
-  !          or by iniaqua)
-  IF ( (maxval(q_ancien).EQ.minval(q_ancien))       .OR. &
-       (maxval(ql_ancien).EQ.minval(ql_ancien))     .OR. &
-       (maxval(qs_ancien).EQ.minval(qs_ancien))     .OR. &
-       (maxval(prw_ancien).EQ.minval(prw_ancien))   .OR. &
-       (maxval(prlw_ancien).EQ.minval(prlw_ancien)) .OR. &
-       (maxval(prsw_ancien).EQ.minval(prsw_ancien)) .OR. &
-       (maxval(t_ancien).EQ.minval(t_ancien)) ) THEN
-    ancien_ok=.false.
-  ENDIF
-
-  found=phyetat0_get(klev,clwcon,"CLWCON","CLWCON",0.)
-  found=phyetat0_get(klev,rnebcon,"RNEBCON","RNEBCON",0.)
-  found=phyetat0_get(klev,ratqs,"RATQS","RATQS",0.)
-
-  found=phyetat0_get(1,run_off_lic_0,"RUNOFFLIC0","RUNOFFLIC0",0.)
-
-!==================================
-!  TKE
-!==================================
-!
-  IF (iflag_pbl>1) then
-     found=phyetat0_srf(klev+1,pbl_tke,"TKE","Turb. Kinetic. Energ. ",1.e-8)
-  ENDIF
-
-  IF (iflag_pbl>1 .AND. iflag_wake>=1  .AND. iflag_pbl_split >=1 ) then
-    found=phyetat0_srf(klev+1,wake_delta_pbl_tke,"DELTATKE","Del TKE wk/env",0.)
-!!    found=phyetat0_srf(1,delta_tsurf,"DELTA_TSURF","Delta Ts wk/env ",0.)
-    found=phyetat0_srf(1,delta_tsurf,"DELTATS","Delta Ts wk/env ",0.)
-!!    found=phyetat0_srf(1,beta_aridity,"BETA_S","Aridity factor ",1.)
-    found=phyetat0_srf(1,beta_aridity,"BETAS","Aridity factor ",1.)
-  ENDIF   !(iflag_pbl>1 .AND. iflag_wake>=1 .AND. iflag_pbl_split >=1 )
-
-!==================================
-!  thermiques, poches, convection
-!==================================
-
-! Emanuel
-  found=phyetat0_get(klev,sig1,"sig1","sig1",0.)
-  found=phyetat0_get(klev,w01,"w01","w01",0.)
-
-! Wake
-  found=phyetat0_get(klev,wake_deltat,"WAKE_DELTAT","Delta T wake/env",0.)
-  found=phyetat0_get(klev,wake_deltaq,"WAKE_DELTAQ","Delta hum. wake/env",0.)
-  found=phyetat0_get(1,wake_s,"WAKE_S","Wake frac. area",0.)
-!jyg<
-!  Set wake_dens to -1000. when there is no restart so that the actual
-!  initialization is made in calwake.
-!!  found=phyetat0_get(1,wake_dens,"WAKE_DENS","Wake num. /unit area",0.)
-  found=phyetat0_get(1,wake_dens,"WAKE_DENS","Wake num. /unit area",-1000.)
-  found=phyetat0_get(1,awake_dens,"AWAKE_DENS","Active Wake num. /unit area",0.)
-  found=phyetat0_get(1,cv_gen,"CV_GEN","CB birth rate",0.)
-!>jyg
-  found=phyetat0_get(1,wake_cstar,"WAKE_CSTAR","WAKE_CSTAR",0.)
-  found=phyetat0_get(1,wake_pe,"WAKE_PE","WAKE_PE",0.)
-  found=phyetat0_get(1,wake_fip,"WAKE_FIP","WAKE_FIP",0.)
-
-! Thermiques
-  found=phyetat0_get(1,zmax0,"ZMAX0","ZMAX0",40.)
-  found=phyetat0_get(1,f0,"F0","F0",1.e-5)
-  found=phyetat0_get(klev+1,fm_therm,"FM_THERM","Thermals mass flux",0.)
-  found=phyetat0_get(klev,entr_therm,"ENTR_THERM","Thermals Entrain.",0.)
-  found=phyetat0_get(klev,detr_therm,"DETR_THERM","Thermals Detrain.",0.)
-
-! ALE/ALP
-  found=phyetat0_get(1,ale_bl,"ALE_BL","ALE BL",0.)
-  found=phyetat0_get(1,ale_bl_trig,"ALE_BL_TRIG","ALE BL_TRIG",0.)
-  found=phyetat0_get(1,alp_bl,"ALP_BL","ALP BL",0.)
-  found=phyetat0_get(1,ale_wake,"ALE_WAKE","ALE_WAKE",0.)
-  found=phyetat0_get(1,ale_bl_stat,"ALE_BL_STAT","ALE_BL_STAT",0.)
-
-! fisrtilp/Clouds 0.002 could be ratqsbas. But can stay like this as well
-  found=phyetat0_get(klev,ratqs_inter,"RATQS_INTER","Relative width of the lsc sugrid scale water",0.002)
-
-!===========================================
-  ! Read and send field trs to traclmdz
-!===========================================
-
-!--OB now this is for co2i - ThL: and therefore also for inco
-  IF (ANY(types_trac == 'co2i') .OR. ANY(types_trac == 'inco')) THEN
-     IF (carbon_cycle_cpl) THEN
-        ALLOCATE(co2_send(klon), stat=ierr)
-        IF (ierr /= 0) CALL abort_physic('phyetat0', 'pb allocation co2_send', 1)
-        found=phyetat0_get(1,co2_send,"co2_send","co2 send",co2_ppm)
-     ENDIF
-  ELSE IF (ANY(types_trac == 'lmdz')) THEN
-     it = 0
-     DO iq = 1, nqtot
-        IF(.NOT.(tracers(iq)%isAdvected .AND. tracers(iq)%isInPhysics)) CYCLE
-        it = it+1
-        found=phyetat0_get(1,trs(:,it),"trs_"//TRIM(tracers(iq)%name), &
-                                  "Surf trac"//TRIM(tracers(iq)%name),0.)
-     END DO
-     CALL traclmdz_from_restart(trs)
-  ENDIF
-
-#ifdef ISO
-        ! initialise les isotopes       
-        write(*,*) 'phyetat0 1069' 
-         CALL phyisoetat0 (snow,run_off_lic_0, &
-     &           xtsnow,xtrun_off_lic_0, &
-     &           Rland_ice)
-#ifdef ISOVERIF 
-      write(*,*) 'phyetat0 1074'
-      if (iso_eau.gt.0) then
-      call iso_verif_egalite_vect2D(  &
-     &           xtsnow,snow, &
-     &           'phyetat0 1101a',niso,klon,nbsrf)
-        do i=1,klon  
-              call iso_verif_egalite(Rland_ice(iso_eau,i),1.0, &
-     &         'phyetat0 1101b')
-         enddo
-      endif
-      write(*,*) 'phyetat0 1102'
-#endif
-#endif
-
-!===========================================
-!  ondes de gravite / relief
-!===========================================
-
-!  ondes de gravite non orographiques
-  IF (ok_gwd_rando) found = &
-       phyetat0_get(klev,du_gwd_rando,"du_gwd_rando","du_gwd_rando",0.)
-  IF (.NOT. ok_hines .AND. ok_gwd_rando) found &
-       = phyetat0_get(klev,du_gwd_front,"du_gwd_front","du_gwd_front",0.)
-
-!  prise en compte du relief sous-maille
-  found=phyetat0_get(1,zmea,"ZMEA","sub grid orography",0.)
-  found=phyetat0_get(1,zstd,"ZSTD","sub grid orography",0.)
-  found=phyetat0_get(1,zsig,"ZSIG","sub grid orography",0.)
-  found=phyetat0_get(1,zgam,"ZGAM","sub grid orography",0.)
-  found=phyetat0_get(1,zthe,"ZTHE","sub grid orography",0.)
-  found=phyetat0_get(1,zpic,"ZPIC","sub grid orography",0.)
-  found=phyetat0_get(1,zval,"ZVAL","sub grid orography",0.)
-  found=phyetat0_get(1,zmea,"ZMEA","sub grid orography",0.)
-  found=phyetat0_get(1,rugoro,"RUGSREL","sub grid orography",0.)
-
-!===========================================
-! Initialize ocean
-!===========================================
-
-  IF ( type_ocean == 'slab' ) THEN
-      CALL ocean_slab_init(phys_tstep, pctsrf)
-      IF (nslay.EQ.1) THEN
-        found=phyetat0_get(1,tslab,"tslab01","tslab",0.)
-        IF (.NOT. found) THEN
-            found=phyetat0_get(1,tslab,"tslab","tslab",0.)
-        ENDIF
-      ELSE
-          DO i=1,nslay
-            WRITE(str2,'(i2.2)') i
-            found=phyetat0_get(1,tslab(:,i),"tslab"//str2,"tslab",0.)  
-          ENDDO
-      ENDIF
-      IF (.NOT. found) THEN 
-          PRINT*, "phyetat0: Le champ <tslab> est absent"
-          PRINT*, "Initialisation a tsol_oce"
-          DO i=1,nslay
-              tslab(:,i)=MAX(ftsol(:,is_oce),271.35)
-          ENDDO
-      ENDIF 
-
-      ! Sea ice variables
-      IF (version_ocean == 'sicINT') THEN
-          found=phyetat0_get(1,tice,"slab_tice","slab_tice",0.)
-          IF (.NOT. found) THEN 
-              PRINT*, "phyetat0: Le champ <tice> est absent"
-              PRINT*, "Initialisation a tsol_sic"
-                  tice(:)=ftsol(:,is_sic)
-          ENDIF 
-          found=phyetat0_get(1,seaice,"seaice","seaice",0.)
-          IF (.NOT. found) THEN
-              PRINT*, "phyetat0: Le champ <seaice> est absent"
-              PRINT*, "Initialisation a 0/1m suivant fraction glace"
-              seaice(:)=0.
-              WHERE (pctsrf(:,is_sic).GT.EPSFRA)
-                  seaice=917.
-              ENDWHERE
-          ENDIF
-      ENDIF !sea ice INT
-  ENDIF ! Slab        
-
-  if (activate_ocean_skin >= 1) then
-     if (activate_ocean_skin == 2 .and. type_ocean == 'couple') then
-        found = phyetat0_get(1, delta_sal, "delta_sal", &
-             "ocean-air interface salinity minus bulk salinity", 0.)
-        found = phyetat0_get(1, delta_sst, "delta_SST", &
-             "ocean-air interface temperature minus bulk SST", 0.)
-     end if
-     
-     found = phyetat0_get(1, ds_ns, "dS_ns", "delta salinity near surface", 0.)
-     found = phyetat0_get(1, dt_ns, "dT_ns", "delta temperature near surface", &
-          0.)
-
-     where (pctsrf(:, is_oce) == 0.)
-        ds_ns = missing_val
-        dt_ns = missing_val
-        delta_sst = missing_val
-        delta_sal = missing_val
-     end where
-  end if
-
-  ! on ferme le fichier
-  CALL close_startphy
-
-  ! Initialize module pbl_surface_mod 
-
-#ifdef ISOVERIF
-        write(*,*) 'phyetat0 572: snow(994,:)=',snow(994,2)
-        write(*,*) 'xtsnow(:,994,2)=',xtsnow(:,994,2)
-#endif
-
-  CALL pbl_surface_init(fder, snow, qsurf, tsoil)
-#ifdef ISO
-  CALL pbl_surface_init_iso(xtsnow,Rland_ice)
-#endif
-
-  ! Initialize module ocean_cpl_mod for the case of coupled ocean
-  IF ( type_ocean == 'couple' ) THEN
-     CALL ocean_cpl_init(phys_tstep, longitude_deg, latitude_deg)
-  ENDIF
-
-!  CALL init_iophy_new(latitude_deg, longitude_deg)
-
-  ! Initilialize module fonte_neige_mod      
-  CALL fonte_neige_init(run_off_lic_0)
-#ifdef ISO
-   CALL fonte_neige_init_iso(xtrun_off_lic_0)
-#endif
-
-END SUBROUTINE phyetat0
-
-!===================================================================
-FUNCTION phyetat0_get(nlev,field,name,descr,default)
-!===================================================================
-! Lecture d'un champ avec contrôle
-! Function logique dont le resultat indique si la lecture
-! s'est bien passée
-! On donne une valeur par defaut dans le cas contraire
-!===================================================================
-
-USE iostart, ONLY : get_field
-USE dimphy, only: klon 
-USE print_control_mod, ONLY: lunout
-
-IMPLICIT NONE
-
-LOGICAL phyetat0_get
-
-! arguments
-INTEGER,INTENT(IN) :: nlev
-CHARACTER*(*),INTENT(IN) :: name,descr
-REAL,INTENT(IN) :: default
-REAL,DIMENSION(klon,nlev),INTENT(INOUT) :: field
-
-! Local variables
-LOGICAL found
-
-   CALL get_field(name, field, found)
-   IF (.NOT. found) THEN
-     WRITE(lunout,*) "phyetat0: Le champ <",TRIM(name),"> est absent"
-     WRITE(lunout,*) "Depart legerement fausse. Mais je continue"
-     field(:,:)=default
-   ENDIF
-   WRITE(lunout,*) name, descr, MINval(field),MAXval(field)
-   phyetat0_get=found
-
-RETURN
-END FUNCTION phyetat0_get
-
-!================================================================
-FUNCTION phyetat0_srf(nlev,field,name,descr,default)
-!===================================================================
-! Lecture d'un champ par sous-surface avec contrôle
-! Function logique dont le resultat indique si la lecture
-! s'est bien passée
-! On donne une valeur par defaut dans le cas contraire
-!===================================================================
-
-USE iostart, ONLY : get_field
-USE dimphy, only: klon 
-USE indice_sol_mod, only: nbsrf
-USE print_control_mod, ONLY: lunout
-
-IMPLICIT NONE
-
-LOGICAL phyetat0_srf
-! arguments
-INTEGER,INTENT(IN) :: nlev
-CHARACTER*(*),INTENT(IN) :: name,descr
-REAL,INTENT(IN) :: default
-REAL,DIMENSION(klon,nlev,nbsrf),INTENT(INOUT) :: field
-
-! Local variables
-LOGICAL found,phyetat0_get
-INTEGER nsrf
-CHARACTER*2 str2
-  
-     IF (nbsrf.GT.99) THEN
-        WRITE(lunout,*) "Trop de sous-mailles"
-        call abort_physic("phyetat0", "", 1)
-     ENDIF
-
-     DO nsrf = 1, nbsrf
-        WRITE(str2, '(i2.2)') nsrf
-        found= phyetat0_get(nlev,field(:,:, nsrf), &
-        name//str2,descr//" srf:"//str2,default)
-     ENDDO
-
-     phyetat0_srf=found
-
-RETURN
-END FUNCTION phyetat0_srf
-
Index: LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/phyetat0_get_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/phyetat0_get_mod.F90	(revision 4482)
+++ LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/phyetat0_get_mod.F90	(revision 4482)
@@ -0,0 +1,1 @@
+link ../phylmd/phyetat0_get_mod.F90
Index: LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/phyetat0_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/phyetat0_mod.F90	(revision 4482)
+++ LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/phyetat0_mod.F90	(revision 4482)
@@ -0,0 +1,661 @@
+! $Id: phyetat0.F90 3890 2021-05-05 15:15:06Z jyg $
+
+MODULE phyetat0_mod
+
+  PRIVATE
+  PUBLIC :: phyetat0
+
+CONTAINS
+
+SUBROUTINE phyetat0 (fichnom, clesphy0, tabcntr0)
+
+  USE dimphy, only: klon, zmasq, klev
+  USE iophy, ONLY : init_iophy_new
+  USE ocean_cpl_mod,    ONLY : ocean_cpl_init
+  USE fonte_neige_mod,  ONLY : fonte_neige_init
+  USE pbl_surface_mod,  ONLY : pbl_surface_init
+#ifdef ISO
+  USE fonte_neige_mod,  ONLY : fonte_neige_init_iso
+  USE pbl_surface_mod,  ONLY : pbl_surface_init_iso
+#endif
+  USE phyetat0_get_mod, ONLY : phyetat0_get, phyetat0_srf
+  USE surface_data,     ONLY : type_ocean, version_ocean
+  USE phys_state_var_mod, ONLY : ancien_ok, clwcon, detr_therm, phys_tstep, &
+       qsol, fevap, z0m, z0h, agesno, &
+       du_gwd_rando, du_gwd_front, entr_therm, f0, fm_therm, &
+       falb_dir, falb_dif, prw_ancien, prlw_ancien, prsw_ancien, &
+       ftsol, pbl_tke, pctsrf, q_ancien, ql_ancien, qs_ancien, rneb_ancien, radpas, radsol, rain_fall, ratqs, &
+       rnebcon, rugoro, sig1, snow_fall, solaire_etat0, sollw, sollwdown, &
+       solsw, solswfdiff, t_ancien, u_ancien, v_ancien, w01, wake_cstar, wake_deltaq, &
+       wake_deltat, wake_delta_pbl_TKE, delta_tsurf, beta_aridity, wake_fip, wake_pe, &
+       wake_s, wake_dens, awake_dens, cv_gen, zgam, zmax0, zmea, zpic, zsig, &
+#ifdef ISO
+       fxtevap, xtsol, xt_ancien, xtl_ancien, xts_ancien, wake_deltaxt, &
+       xtrain_fall,xtsnow_fall, &
+#endif
+       zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl, u10m, v10m, treedrg, &
+       ale_wake, ale_bl_stat, ds_ns, dt_ns, delta_sst, delta_sal, dter, dser, &
+       dt_ds, ratqs_inter
+!FC
+  USE geometry_mod,     ONLY: longitude_deg, latitude_deg
+  USE iostart,          ONLY: close_startphy, get_field, get_var, open_startphy
+  USE infotrac_phy,     ONLY: nqtot, nbtr, type_trac, tracers
+  USE readTracFiles_mod,ONLY: maxlen, new2oldH2O
+  USE traclmdz_mod,     ONLY: traclmdz_from_restart
+  USE carbon_cycle_mod, ONLY: carbon_cycle_init, carbon_cycle_cpl, carbon_cycle_tr, carbon_cycle_rad, co2_send, RCO2_glo
+  USE indice_sol_mod,   ONLY: nbsrf, is_ter, epsfra, is_lic, is_oce, is_sic
+  USE ocean_slab_mod,   ONLY: nslay, tslab, seaice, tice, ocean_slab_init
+  USE time_phylmdz_mod, ONLY: init_iteration, pdtphys, itau_phy
+#ifdef CPP_XIOS
+  USE wxios, ONLY: missing_val
+#else
+  use netcdf, only: missing_val => nf90_fill_real
+#endif
+  use config_ocean_skin_m, only: activate_ocean_skin
+#ifdef ISO
+  USE infotrac_phy, ONLY: niso
+  USE isotopes_routines_mod, ONLY: phyisoetat0
+  USE isotopes_mod, ONLY: iso_eau
+#ifdef ISOVERIF
+  USE isotopes_verif_mod, ONLY: iso_verif_egalite_vect2D,iso_verif_egalite
+#endif
+#endif
+
+  IMPLICIT none
+  !======================================================================
+  ! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
+  ! Objet: Lecture de l'etat initial pour la physique
+  !======================================================================
+  include "dimsoil.h"
+  include "clesphys.h"
+  include "alpale.h"
+  include "compbl.h"
+  include "YOMCST.h"
+  !======================================================================
+  CHARACTER*(*) fichnom
+
+  ! les variables globales lues dans le fichier restart
+
+  REAL tsoil(klon, nsoilmx, nbsrf)
+  REAL qsurf(klon, nbsrf)
+  REAL snow(klon, nbsrf)
+  real fder(klon)
+  REAL run_off_lic_0(klon)
+  REAL fractint(klon)
+  REAL trs(klon, nbtr)
+  REAL zts(klon)
+  ! pour drag arbres FC
+  REAL drg_ter(klon,klev)
+
+  CHARACTER*6 ocean_in
+  LOGICAL ok_veget_in
+
+  INTEGER        longcles
+  PARAMETER    ( longcles = 20 )
+  REAL clesphy0( longcles )
+
+  REAL xmin, xmax
+
+  INTEGER nid, nvarid
+  INTEGER ierr, i, nsrf, isoil , k
+  INTEGER length
+  PARAMETER (length=100)
+  INTEGER it, iq, isw
+  REAL tab_cntrl(length), tabcntr0(length)
+  CHARACTER*7 str7
+  CHARACTER*2 str2
+  LOGICAL :: found
+  REAL :: lon_startphy(klon), lat_startphy(klon)
+  CHARACTER(LEN=maxlen) :: tname, t(2)
+
+#ifdef ISO
+  REAL xtsnow(niso,klon, nbsrf)
+  REAL xtrun_off_lic_0(niso,klon)
+  REAL Rland_ice(niso,klon)
+#endif
+  ! FH1D
+  !     real iolat(jjm+1)
+  !real iolat(jjm+1-1/(iim*jjm))
+
+  ! Ouvrir le fichier contenant l'etat initial:
+
+  CALL open_startphy(fichnom)
+
+  ! Lecture des parametres de controle:
+
+  CALL get_var("controle", tab_cntrl)
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  ! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
+  ! Les constantes de la physiques sont lues dans la physique seulement.
+  ! Les egalites du type
+  !             tab_cntrl( 5 )=clesphy0(1)
+  ! sont remplacees par
+  !             clesphy0(1)=tab_cntrl( 5 )
+  ! On inverse aussi la logique.
+  ! On remplit les tab_cntrl avec les parametres lus dans les .def
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  DO i = 1, length
+     tabcntr0( i ) = tab_cntrl( i )
+  ENDDO
+
+  tab_cntrl(1)=pdtphys
+  tab_cntrl(2)=radpas
+
+  ! co2_ppm : value from the previous time step
+
+  ! co2_ppm0 : initial value of atmospheric CO2 (from create_etat0_limit.e .def)
+  co2_ppm0 = 284.32
+  ! when no initial value is available e.g., from a restart 
+  ! this variable must be set  in a .def file which will then be 
+  ! used by the conf_phys_m.F90 routine.
+  ! co2_ppm0 = 284.32 (illustrative example on how to set the variable in .def
+  ! file, for a pre-industrial CO2 concentration value)
+
+  IF (carbon_cycle_tr .OR. carbon_cycle_cpl) THEN
+     co2_ppm = tab_cntrl(3)
+     RCO2    = co2_ppm * 1.0e-06 * RMCO2 / RMD
+     IF (tab_cntrl(17) > 0. .AND. carbon_cycle_rad) THEN
+           RCO2_glo = tab_cntrl(17)
+       ELSE
+           RCO2_glo    = co2_ppm0 * 1.0e-06 * RMCO2 / RMD
+     ENDIF
+     ! ELSE : keep value from .def
+  ENDIF
+
+  solaire_etat0      = tab_cntrl(4)
+  tab_cntrl(5)=iflag_con
+  tab_cntrl(6)=nbapp_rad
+
+  IF (iflag_cycle_diurne.GE.1) tab_cntrl( 7) = iflag_cycle_diurne
+  IF (soil_model) tab_cntrl( 8) =1.
+  IF (new_oliq) tab_cntrl( 9) =1.
+  IF (ok_orodr) tab_cntrl(10) =1.
+  IF (ok_orolf) tab_cntrl(11) =1.
+  IF (ok_limitvrai) tab_cntrl(12) =1.
+
+  itau_phy = tab_cntrl(15)
+
+  clesphy0(1)=tab_cntrl( 5 )
+  clesphy0(2)=tab_cntrl( 6 )
+  clesphy0(3)=tab_cntrl( 7 )
+  clesphy0(4)=tab_cntrl( 8 )
+  clesphy0(5)=tab_cntrl( 9 )
+  clesphy0(6)=tab_cntrl( 10 )
+  clesphy0(7)=tab_cntrl( 11 )
+  clesphy0(8)=tab_cntrl( 12 )
+  clesphy0(9)=tab_cntrl( 17 )
+
+  ! set time iteration
+   CALL init_iteration(itau_phy)
+
+  ! read latitudes and make a sanity check (because already known from dyn)
+  CALL get_field("latitude",lat_startphy)
+  DO i=1,klon
+    IF (ABS(lat_startphy(i)-latitude_deg(i))>=1) THEN
+      WRITE(*,*) "phyetat0: Error! Latitude discrepancy wrt startphy file:",&
+                 " i=",i," lat_startphy(i)=",lat_startphy(i),&
+                 " latitude_deg(i)=",latitude_deg(i)
+      ! This is presumably serious enough to abort run
+      CALL abort_physic("phyetat0","discrepancy in latitudes!",1)
+    ENDIF
+    IF (ABS(lat_startphy(i)-latitude_deg(i))>=0.0001) THEN
+      WRITE(*,*) "phyetat0: Warning! Latitude discrepancy wrt startphy file:",&
+                 " i=",i," lat_startphy(i)=",lat_startphy(i),&
+                 " latitude_deg(i)=",latitude_deg(i)
+    ENDIF
+  ENDDO
+
+  ! read longitudes and make a sanity check (because already known from dyn)
+  CALL get_field("longitude",lon_startphy)
+  DO i=1,klon
+    IF (ABS(lon_startphy(i)-longitude_deg(i))>=1) THEN
+      IF (ABS(360-ABS(lon_startphy(i)-longitude_deg(i)))>=1) THEN
+        WRITE(*,*) "phyetat0: Error! Longitude discrepancy wrt startphy file:",&
+                   " i=",i," lon_startphy(i)=",lon_startphy(i),&
+                   " longitude_deg(i)=",longitude_deg(i)
+        ! This is presumably serious enough to abort run
+        CALL abort_physic("phyetat0","discrepancy in longitudes!",1)
+      ENDIF
+    ENDIF
+    IF (ABS(lon_startphy(i)-longitude_deg(i))>=1) THEN
+      IF (ABS(360-ABS(lon_startphy(i)-longitude_deg(i))) > 0.0001) THEN
+        WRITE(*,*) "phyetat0: Warning! Longitude discrepancy wrt startphy file:",&
+                   " i=",i," lon_startphy(i)=",lon_startphy(i),&
+                   " longitude_deg(i)=",longitude_deg(i)
+      ENDIF 
+    ENDIF
+  ENDDO
+
+  ! Lecture du masque terre mer
+
+  CALL get_field("masque", zmasq, found)
+  IF (.NOT. found) THEN
+     PRINT*, 'phyetat0: Le champ <masque> est absent'
+     PRINT *, 'fichier startphy non compatible avec phyetat0'
+  ENDIF
+
+  ! Lecture des fractions pour chaque sous-surface
+
+  ! initialisation des sous-surfaces
+
+  pctsrf = 0.
+
+  ! fraction de terre
+
+  CALL get_field("FTER", pctsrf(:, is_ter), found)
+  IF (.NOT. found) PRINT*, 'phyetat0: Le champ <FTER> est absent'
+
+  ! fraction de glace de terre
+
+  CALL get_field("FLIC", pctsrf(:, is_lic), found)
+  IF (.NOT. found) PRINT*, 'phyetat0: Le champ <FLIC> est absent'
+
+  ! fraction d'ocean
+
+  CALL get_field("FOCE", pctsrf(:, is_oce), found)
+  IF (.NOT. found) PRINT*, 'phyetat0: Le champ <FOCE> est absent'
+
+  ! fraction glace de mer
+
+  CALL get_field("FSIC", pctsrf(:, is_sic), found)
+  IF (.NOT. found) PRINT*, 'phyetat0: Le champ <FSIC> est absent'
+
+  !  Verification de l'adequation entre le masque et les sous-surfaces
+
+  fractint( 1 : klon) = pctsrf(1 : klon, is_ter)  &
+       + pctsrf(1 : klon, is_lic)
+  DO i = 1 , klon
+     IF ( abs(fractint(i) - zmasq(i) ) .GT. EPSFRA ) THEN
+        WRITE(*, *) 'phyetat0: attention fraction terre pas ',  &
+             'coherente ', i, zmasq(i), pctsrf(i, is_ter) &
+             , pctsrf(i, is_lic)
+        WRITE(*, *) 'Je force la coherence zmasq=fractint'
+        zmasq(i) = fractint(i)
+     ENDIF
+  ENDDO
+  fractint (1 : klon) =  pctsrf(1 : klon, is_oce)  &
+       + pctsrf(1 : klon, is_sic)
+  DO i = 1 , klon
+     IF ( abs( fractint(i) - (1. - zmasq(i))) .GT. EPSFRA ) THEN
+        WRITE(*, *) 'phyetat0 attention fraction ocean pas ',  &
+             'coherente ', i, zmasq(i) , pctsrf(i, is_oce) &
+             , pctsrf(i, is_sic)
+        WRITE(*, *) 'Je force la coherence zmasq=1.-fractint'
+        zmasq(i) = 1. - fractint(i)
+     ENDIF
+  ENDDO
+
+!===================================================================
+! Lecture des temperatures du sol:
+!===================================================================
+
+  found=phyetat0_get(ftsol(:,1),"TS","Surface temperature",283.)
+  IF (found) THEN
+     DO nsrf=2,nbsrf
+        ftsol(:,nsrf)=ftsol(:,1)
+     ENDDO
+  ELSE
+     found=phyetat0_srf(ftsol,"TS","Surface temperature",283.)
+  ENDIF
+
+!===================================================================
+  ! Lecture des albedo difus et direct
+!===================================================================
+
+  DO nsrf = 1, nbsrf
+     DO isw=1, nsw
+        IF (isw.GT.99) THEN
+           PRINT*, "Trop de bandes SW"
+           call abort_physic("phyetat0", "", 1)
+        ENDIF
+        WRITE(str2, '(i2.2)') isw
+        found=phyetat0_srf(falb_dir(:, isw,:),"A_dir_SW"//str2//"srf","Direct Albedo",0.2)
+        found=phyetat0_srf(falb_dif(:, isw,:),"A_dif_SW"//str2//"srf","Direct Albedo",0.2)
+     ENDDO
+  ENDDO
+
+  found=phyetat0_srf(u10m,"U10M","u a 10m",0.)
+  found=phyetat0_srf(v10m,"V10M","v a 10m",0.)
+
+!===================================================================
+  ! Lecture des temperatures du sol profond:
+!===================================================================
+
+   DO isoil=1, nsoilmx
+        IF (isoil.GT.99) THEN
+           PRINT*, "Trop de couches "
+           call abort_physic("phyetat0", "", 1)
+        ENDIF
+        WRITE(str2,'(i2.2)') isoil
+        found=phyetat0_srf(tsoil(:, isoil,:),"Tsoil"//str2//"srf","Temp soil",0.)
+        IF (.NOT. found) THEN
+           PRINT*, "phyetat0: Le champ <Tsoil"//str7//"> est absent"
+           PRINT*, "          Il prend donc la valeur de surface"
+           tsoil(:, isoil, :)=ftsol(:, :)
+        ENDIF
+   ENDDO
+
+!=======================================================================
+! Lecture precipitation/evaporation
+!=======================================================================
+
+  found=phyetat0_srf(qsurf,"QS","Near surface hmidity",0.)
+  found=phyetat0_get(qsol,"QSOL","Surface hmidity / bucket",0.)
+  found=phyetat0_srf(snow,"SNOW","Surface snow",0.)
+  found=phyetat0_srf(fevap,"EVAP","evaporation",0.)
+  found=phyetat0_get(snow_fall,"snow_f","snow fall",0.)
+  found=phyetat0_get(rain_fall,"rain_f","rain fall",0.)
+
+!=======================================================================
+! Radiation
+!=======================================================================
+
+  found=phyetat0_get(solsw,"solsw","net SW radiation surf",0.)
+  found=phyetat0_get(solswfdiff,"solswfdiff","fraction of SW radiation surf that is diffuse",1.)
+  found=phyetat0_get(sollw,"sollw","net LW radiation surf",0.)
+  found=phyetat0_get(sollwdown,"sollwdown","down LW radiation surf",0.)
+  IF (.NOT. found) THEN
+     sollwdown(:) = 0. ;  zts(:)=0.
+     DO nsrf=1,nbsrf
+        zts(:)=zts(:)+ftsol(:,nsrf)*pctsrf(:,nsrf)
+     ENDDO
+     sollwdown(:)=sollw(:)+RSIGMA*zts(:)**4
+  ENDIF
+
+  found=phyetat0_get(radsol,"RADS","Solar radiation",0.)
+  found=phyetat0_get(fder,"fder","Flux derivative",0.) 
+
+
+  ! Lecture de la longueur de rugosite 
+  found=phyetat0_srf(z0m,"RUG","Z0m ancien",0.001)
+  IF (found) THEN
+     z0h(:,1:nbsrf)=z0m(:,1:nbsrf)
+  ELSE
+     found=phyetat0_srf(z0m,"Z0m","Roughness length, momentum ",0.001)
+     found=phyetat0_srf(z0h,"Z0h","Roughness length, enthalpy ",0.001)
+  ENDIF
+!FC
+  IF (ifl_pbltree>0) THEN
+!CALL get_field("FTER", pctsrf(:, is_ter), found)
+    treedrg(:,1:klev,1:nbsrf)= 0.0
+    CALL get_field("treedrg_ter", drg_ter(:,:), found)
+!  found=phyetat0_srf(treedrg,"treedrg","drag from vegetation" , 0.)
+    !lecture du profile de freinage des arbres
+    IF (.not. found ) THEN
+      treedrg(:,1:klev,1:nbsrf)= 0.0
+    ELSE
+      treedrg(:,1:klev,is_ter)= drg_ter(:,:)
+!     found=phyetat0_get(treedrg,"treedrg","freinage arbres",0.)
+    ENDIF
+  ELSE
+    ! initialize treedrg(), because it will be written in restartphy.nc
+    treedrg(:,:,:) = 0.0 
+  ENDIF
+
+  ! Lecture de l'age de la neige:
+  found=phyetat0_srf(agesno,"AGESNO","SNOW AGE",0.001)
+
+  ancien_ok=.true.
+  ancien_ok=ancien_ok.AND.phyetat0_get(t_ancien,"TANCIEN","TANCIEN",0.)
+  ancien_ok=ancien_ok.AND.phyetat0_get(q_ancien,"QANCIEN","QANCIEN",0.)
+  ancien_ok=ancien_ok.AND.phyetat0_get(ql_ancien,"QLANCIEN","QLANCIEN",0.)
+  ancien_ok=ancien_ok.AND.phyetat0_get(qs_ancien,"QSANCIEN","QSANCIEN",0.)
+  ancien_ok=ancien_ok.AND.phyetat0_get(rneb_ancien,"RNEBANCIEN","RNEBANCIEN",0.)
+  ancien_ok=ancien_ok.AND.phyetat0_get(u_ancien,"UANCIEN","UANCIEN",0.)
+  ancien_ok=ancien_ok.AND.phyetat0_get(v_ancien,"VANCIEN","VANCIEN",0.)
+  ancien_ok=ancien_ok.AND.phyetat0_get(prw_ancien,"PRWANCIEN","PRWANCIEN",0.)
+  ancien_ok=ancien_ok.AND.phyetat0_get(prlw_ancien,"PRLWANCIEN","PRLWANCIEN",0.)
+  ancien_ok=ancien_ok.AND.phyetat0_get(prsw_ancien,"PRSWANCIEN","PRSWANCIEN",0.)
+
+  ! Ehouarn: addtional tests to check if t_ancien, q_ancien contain
+  !          dummy values (as is the case when generated by ce0l,
+  !          or by iniaqua)
+  IF ( (maxval(q_ancien).EQ.minval(q_ancien))       .OR. &
+       (maxval(ql_ancien).EQ.minval(ql_ancien))     .OR. &
+       (maxval(qs_ancien).EQ.minval(qs_ancien))     .OR. &
+       (maxval(rneb_ancien).EQ.minval(rneb_ancien)) .OR. &
+       (maxval(prw_ancien).EQ.minval(prw_ancien))   .OR. &
+       (maxval(prlw_ancien).EQ.minval(prlw_ancien)) .OR. &
+       (maxval(prsw_ancien).EQ.minval(prsw_ancien)) .OR. &
+       (maxval(t_ancien).EQ.minval(t_ancien)) ) THEN
+    ancien_ok=.false.
+  ENDIF
+
+  found=phyetat0_get(clwcon,"CLWCON","CLWCON",0.)
+  found=phyetat0_get(rnebcon,"RNEBCON","RNEBCON",0.)
+  found=phyetat0_get(ratqs,"RATQS","RATQS",0.)
+
+  found=phyetat0_get(run_off_lic_0,"RUNOFFLIC0","RUNOFFLIC0",0.)
+
+!==================================
+!  TKE
+!==================================
+!
+  IF (iflag_pbl>1) then
+     found=phyetat0_srf(pbl_tke,"TKE","Turb. Kinetic. Energ. ",1.e-8)
+  ENDIF
+
+  IF (iflag_pbl>1 .AND. iflag_wake>=1  .AND. iflag_pbl_split >=1 ) then
+    found=phyetat0_srf(wake_delta_pbl_tke,"DELTATKE","Del TKE wk/env",0.)
+!!    found=phyetat0_srf(delta_tsurf,"DELTA_TSURF","Delta Ts wk/env ",0.)
+    found=phyetat0_srf(delta_tsurf,"DELTATS","Delta Ts wk/env ",0.)
+!!    found=phyetat0_srf(beta_aridity,"BETA_S","Aridity factor ",1.)
+    found=phyetat0_srf(beta_aridity,"BETAS","Aridity factor ",1.)
+  ENDIF   !(iflag_pbl>1 .AND. iflag_wake>=1 .AND. iflag_pbl_split >=1 )
+
+!==================================
+!  thermiques, poches, convection
+!==================================
+
+! Emanuel
+  found=phyetat0_get(sig1,"sig1","sig1",0.)
+  found=phyetat0_get(w01,"w01","w01",0.)
+
+! Wake
+  found=phyetat0_get(wake_deltat,"WAKE_DELTAT","Delta T wake/env",0.)
+  found=phyetat0_get(wake_deltaq,"WAKE_DELTAQ","Delta hum. wake/env",0.)
+  found=phyetat0_get(wake_s,"WAKE_S","Wake frac. area",0.)
+!jyg<
+!  Set wake_dens to -1000. when there is no restart so that the actual
+!  initialization is made in calwake.
+!!  found=phyetat0_get(1,wake_dens,"WAKE_DENS","Wake num. /unit area",0.)
+  found=phyetat0_get(wake_dens,"WAKE_DENS","Wake num. /unit area",-1000.)
+  found=phyetat0_get(awake_dens,"AWAKE_DENS","Active Wake num. /unit area",0.)
+  found=phyetat0_get(cv_gen,"CV_GEN","CB birth rate",0.)
+!>jyg
+  found=phyetat0_get(wake_cstar,"WAKE_CSTAR","WAKE_CSTAR",0.)
+  found=phyetat0_get(wake_pe,"WAKE_PE","WAKE_PE",0.)
+  found=phyetat0_get(wake_fip,"WAKE_FIP","WAKE_FIP",0.)
+
+! Thermiques
+  found=phyetat0_get(zmax0,"ZMAX0","ZMAX0",40.)
+  found=phyetat0_get(f0,"F0","F0",1.e-5)
+  found=phyetat0_get(fm_therm,"FM_THERM","Thermals mass flux",0.)
+  found=phyetat0_get(entr_therm,"ENTR_THERM","Thermals Entrain.",0.)
+  found=phyetat0_get(detr_therm,"DETR_THERM","Thermals Detrain.",0.)
+
+! ALE/ALP
+  found=phyetat0_get(ale_bl,"ALE_BL","ALE BL",0.)
+  found=phyetat0_get(ale_bl_trig,"ALE_BL_TRIG","ALE BL_TRIG",0.)
+  found=phyetat0_get(alp_bl,"ALP_BL","ALP BL",0.)
+  found=phyetat0_get(ale_wake,"ALE_WAKE","ALE_WAKE",0.)
+  found=phyetat0_get(ale_bl_stat,"ALE_BL_STAT","ALE_BL_STAT",0.)
+
+! fisrtilp/Clouds 0.002 could be ratqsbas. But can stay like this as well
+  found=phyetat0_get(ratqs_inter,"RATQS_INTER","Relative width of the lsc sugrid scale water",0.002)
+
+!===========================================
+  ! Read and send field trs to traclmdz
+!===========================================
+
+!--OB now this is for co2i - ThL: and therefore also for inco
+  IF (ANY(type_trac == ['co2i','inco'])) THEN
+     IF (carbon_cycle_cpl) THEN
+        ALLOCATE(co2_send(klon), stat=ierr)
+        IF (ierr /= 0) CALL abort_physic('phyetat0', 'pb allocation co2_send', 1)
+        found=phyetat0_get(co2_send,"co2_send","co2 send",co2_ppm0)
+     ENDIF
+  ELSE IF (type_trac == 'lmdz') THEN
+     it = 0
+     DO iq = 1, nqtot
+        IF(.NOT.(tracers(iq)%isAdvected .AND. tracers(iq)%isInPhysics)) CYCLE
+        it = it+1
+        tname = tracers(iq)%name
+        t(1) = 'trs_'//TRIM(tname); t(2) = 'trs_'//TRIM(new2oldH2O(tname))
+        found = phyetat0_get(trs(:,it), t(:), "Surf trac"//TRIM(tname), 0.)
+     END DO
+     CALL traclmdz_from_restart(trs)
+  ENDIF
+
+#ifdef ISO
+        ! initialise les isotopes       
+        write(*,*) 'phyetat0 1069' 
+         CALL phyisoetat0 (snow,run_off_lic_0, &
+     &           xtsnow,xtrun_off_lic_0, &
+     &           Rland_ice)
+#ifdef ISOVERIF 
+      write(*,*) 'phyetat0 1074'
+      if (iso_eau.gt.0) then
+      call iso_verif_egalite_vect2D(  &
+     &           xtsnow,snow, &
+     &           'phyetat0 1101a',niso,klon,nbsrf)
+        do i=1,klon  
+              call iso_verif_egalite(Rland_ice(iso_eau,i),1.0, &
+     &         'phyetat0 1101b')
+         enddo
+      endif
+      write(*,*) 'phyetat0 1102'
+#endif
+#endif
+
+!===========================================
+!  ondes de gravite / relief
+!===========================================
+
+!  ondes de gravite non orographiques
+  IF (ok_gwd_rando) found = &
+       phyetat0_get(du_gwd_rando,"du_gwd_rando","du_gwd_rando",0.)
+  IF (.NOT. ok_hines .AND. ok_gwd_rando) found &
+       = phyetat0_get(du_gwd_front,"du_gwd_front","du_gwd_front",0.)
+
+!  prise en compte du relief sous-maille
+  found=phyetat0_get(zmea,"ZMEA","sub grid orography",0.)
+  found=phyetat0_get(zstd,"ZSTD","sub grid orography",0.)
+  found=phyetat0_get(zsig,"ZSIG","sub grid orography",0.)
+  found=phyetat0_get(zgam,"ZGAM","sub grid orography",0.)
+  found=phyetat0_get(zthe,"ZTHE","sub grid orography",0.)
+  found=phyetat0_get(zpic,"ZPIC","sub grid orography",0.)
+  found=phyetat0_get(zval,"ZVAL","sub grid orography",0.)
+  found=phyetat0_get(zmea,"ZMEA","sub grid orography",0.)
+  found=phyetat0_get(rugoro,"RUGSREL","sub grid orography",0.)
+
+!===========================================
+! Initialize ocean
+!===========================================
+
+  IF ( type_ocean == 'slab' ) THEN
+      CALL ocean_slab_init(phys_tstep, pctsrf)
+      IF (nslay.EQ.1) THEN
+        found=phyetat0_get(tslab,["tslab01","tslab  "],"tslab",0.)
+      ELSE
+          DO i=1,nslay
+            WRITE(str2,'(i2.2)') i
+            found=phyetat0_get(tslab(:,i),"tslab"//str2,"tslab",0.)  
+          ENDDO
+      ENDIF
+      IF (.NOT. found) THEN 
+          PRINT*, "phyetat0: Le champ <tslab> est absent"
+          PRINT*, "Initialisation a tsol_oce"
+          DO i=1,nslay
+              tslab(:,i)=MAX(ftsol(:,is_oce),271.35)
+          ENDDO
+      ENDIF 
+
+      ! Sea ice variables
+      IF (version_ocean == 'sicINT') THEN
+          found=phyetat0_get(tice,"slab_tice","slab_tice",0.)
+          IF (.NOT. found) THEN 
+              PRINT*, "phyetat0: Le champ <tice> est absent"
+              PRINT*, "Initialisation a tsol_sic"
+                  tice(:)=ftsol(:,is_sic)
+          ENDIF 
+          found=phyetat0_get(seaice,"seaice","seaice",0.)
+          IF (.NOT. found) THEN
+              PRINT*, "phyetat0: Le champ <seaice> est absent"
+              PRINT*, "Initialisation a 0/1m suivant fraction glace"
+              seaice(:)=0.
+              WHERE (pctsrf(:,is_sic).GT.EPSFRA)
+                  seaice=917.
+              ENDWHERE
+          ENDIF
+      ENDIF !sea ice INT
+  ENDIF ! Slab        
+
+  if (activate_ocean_skin >= 1) then
+     if (activate_ocean_skin == 2 .and. type_ocean == 'couple') then
+        found = phyetat0_get(delta_sal, "delta_sal", &
+             "ocean-air interface salinity minus bulk salinity", 0.)
+        found = phyetat0_get(delta_sst, "delta_SST", &
+             "ocean-air interface temperature minus bulk SST", 0.)
+        found = phyetat0_get(dter, "dter", &
+             "ocean-air interface temperature minus subskin temperature", 0.)
+        found = phyetat0_get(dser, "dser", &
+             "ocean-air interface salinity minus subskin salinity", 0.)
+        found = phyetat0_get(dt_ds, "dt_ds", "(tks / tkt) * dTer", 0.)
+
+        where (pctsrf(:, is_oce) == 0.)
+           delta_sst = missing_val
+           delta_sal = missing_val
+           dter = missing_val
+           dser = missing_val
+           dt_ds = missing_val
+        end where
+     end if
+     
+     found = phyetat0_get(ds_ns, "dS_ns", "delta salinity near surface", 0.)
+     found = phyetat0_get(dt_ns, "dT_ns", "delta temperature near surface", &
+          0.)
+
+     where (pctsrf(:, is_oce) == 0.)
+        ds_ns = missing_val
+        dt_ns = missing_val
+        delta_sst = missing_val
+        delta_sal = missing_val
+     end where
+  end if
+
+  ! on ferme le fichier
+  CALL close_startphy
+
+  ! Initialize module pbl_surface_mod 
+
+! 994 est trop grand quand OpenMP est activé : on sort du tableau
+!#ifdef ISOVERIF
+!        write(*,*) 'phyetat0 572: snow(994,:)=',snow(994,2)
+!        write(*,*) 'xtsnow(:,994,2)=',xtsnow(:,994,2)
+!#endif
+
+  CALL pbl_surface_init(fder, snow, qsurf, tsoil)
+#ifdef ISO
+  CALL pbl_surface_init_iso(xtsnow,Rland_ice)
+#endif
+
+  ! Initialize module ocean_cpl_mod for the case of coupled ocean
+  IF ( type_ocean == 'couple' ) THEN
+     CALL ocean_cpl_init(phys_tstep, longitude_deg, latitude_deg)
+  ENDIF
+
+!  CALL init_iophy_new(latitude_deg, longitude_deg)
+
+  ! Initilialize module fonte_neige_mod      
+  CALL fonte_neige_init(run_off_lic_0)
+#ifdef ISO
+   CALL fonte_neige_init_iso(xtrun_off_lic_0)
+#endif
+
+END SUBROUTINE phyetat0
+
+END MODULE phyetat0_mod
+
Index: LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/phyredem.F90
===================================================================
--- LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/phyredem.F90	(revision 4175)
+++ LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/phyredem.F90	(revision 4482)
@@ -30,5 +30,5 @@
                                 du_gwd_rando, du_gwd_front, u10m, v10m, &
                                 treedrg, solswfdiff, delta_sal, ds_ns, dt_ns, &
-                                delta_sst, ratqs_inter
+                                delta_sst, ratqs_inter, dter, dser, dt_ds
 #ifdef ISO
   USE phys_state_var_mod, ONLY: xtsol, fxtevap,xtrain_fall, xtsnow_fall,     &
@@ -39,5 +39,5 @@
   USE iostart, ONLY: open_restartphy, close_restartphy, enddef_restartphy, put_field, put_var
   USE traclmdz_mod, ONLY : traclmdz_to_restart
-  USE infotrac_phy, ONLY: types_trac, nqtot, tracers, nbtr, niso
+  USE infotrac_phy, ONLY: type_trac, nqtot, tracers, nbtr, niso
 #ifdef ISO
 #ifdef ISOVERIF
@@ -45,5 +45,5 @@
 #endif
 #endif
-  USE carbon_cycle_mod, ONLY : carbon_cycle_cpl, co2_send
+  USE carbon_cycle_mod, ONLY : carbon_cycle_cpl, co2_send, carbon_cycle_rad, RCO2_glo
   USE indice_sol_mod, ONLY: nbsrf, is_oce, is_sic, is_ter, is_lic, epsfra
   USE surface_data, ONLY: type_ocean, version_ocean
@@ -130,5 +130,10 @@
 
   ! co2_ppm0 : initial value of atmospheric CO2
-  tab_cntrl(16) = co2_ppm0
+  ! tab_cntrl(16) = co2_ppm0
+
+  !  PC -- initial value of RCO2 for the radiation scheme
+  !  tab_cntrl(17) = co2_ppm * 1.0e-06 * RMCO2 / RMD
+  IF (carbon_cycle_rad) tab_cntrl(17) = RCO2_glo
+  !PRINT*, "PC : phyredem RCO2_glo =",RCO2_glo
 
   DO pass=1,2   ! pass=1 netcdf definition ; pass=2 netcdf write
@@ -171,13 +176,7 @@
     CALL put_field(pass,"FSIC", "fraction glace mer", pctsrf(:, is_sic))
 
-    IF(nbsrf>99) THEN
-      PRINT*, "Trop de sous-mailles";  CALL abort_physic("phyredem", "", 1)
-    END IF
-    IF(nsoilmx>99) THEN
-      PRINT*, "Trop de sous-surfaces"; CALL abort_physic("phyredem", "", 1)
-    END IF
-    IF(nsw>99) THEN
-      PRINT*, "Trop de bandes"; CALL abort_physic("phyredem", "", 1)
-    END IF
+    IF(nbsrf  >99) CALL abort_physic("phyredem", "Trop de sous-mailles", 1)
+    IF(nsoilmx>99) CALL abort_physic("phyredem", "Trop de sous-mailles", 1)
+    IF(nsw    >99) CALL abort_physic("phyredem", "Trop de bandes", 1)
 
 !    Surface variables 
@@ -345,5 +344,5 @@
 
 
-    IF (ANY(types_trac == 'co2i') .OR. ANY(types_trac == 'inco')) THEN
+    IF (ANY(type_trac == ['co2i','inco'])) THEN
        IF (carbon_cycle_cpl) THEN
           IF (.NOT. ALLOCATED(co2_send)) THEN
@@ -356,5 +355,5 @@
 
     ! trs from traclmdz_mod
-    ELSE IF (ANY(types_trac == 'lmdz')) THEN
+    ELSE IF (type_trac == 'lmdz') THEN
        CALL traclmdz_to_restart(trs)
        it = 0
@@ -394,4 +393,11 @@
           CALL put_field(pass, "delta_SST", &
                "ocean-air interface temperature minus bulk SST", delta_sst)
+          CALL put_field(pass, "dter", &
+               "ocean-air interface temperature minus subskin temperature", &
+               dter)
+          CALL put_field(pass, "dser", &
+               "ocean-air interface salinity minus subskin salinity", dser)
+          CALL put_field(pass, "dt_ds", &
+               "(tks / tkt) * dTer", dt_ds)
        end if
        
Index: LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/phys_output_var_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/phys_output_var_mod.F90	(revision 4175)
+++ LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/phys_output_var_mod.F90	(revision 4482)
@@ -135,25 +135,17 @@
   ! Ocean-atmosphere interface, subskin ocean and near-surface ocean:
   
-  REAL, ALLOCATABLE, SAVE:: dter(:)
-  ! Temperature variation in the diffusive microlayer, that is
-  ! ocean-air interface temperature minus subskin temperature. In K.
-      
-  REAL, SAVE, ALLOCATABLE:: dser(:)
-  ! Temperature variation in the diffusive microlayer, that is
-  ! subskin temperature minus ocean-air interface temperature. In K.
-
-  REAL, SAVE, ALLOCATABLE:: tkt(:)
+  REAL, SAVE, ALLOCATABLE:: tkt(:) ! (klon)
   ! épaisseur (m) de la couche de diffusion thermique (microlayer)
   ! cool skin thickness
 
-  REAL, SAVE, ALLOCATABLE:: tks(:)
+  REAL, SAVE, ALLOCATABLE:: tks(:) ! (klon)
   ! épaisseur (m) de la couche de diffusion de masse (microlayer)
   
-  REAL, SAVE, ALLOCATABLE:: taur(:) ! momentum flux due to rain, in Pa
-
-  REAL, SAVE, ALLOCATABLE:: sss(:)
+  REAL, SAVE, ALLOCATABLE:: taur(:) ! (klon) momentum flux due to rain, in Pa
+
+  REAL, SAVE, ALLOCATABLE:: sss(:) ! (klon)
   ! bulk salinity of the surface layer of the ocean, in ppt
   
-  !$OMP THREADPRIVATE(dter, dser, tkt, tks, taur, sss)
+  !$OMP THREADPRIVATE(tkt, tks, taur, sss)
 
 CONTAINS
@@ -216,6 +208,6 @@
     IF (ok_gwd_rando) allocate(zustr_gwd_rando(klon), zvstr_gwd_rando(klon))
 
-    if (activate_ocean_skin >= 1) allocate(dter(klon), dser(klon), tkt(klon), &
-         tks(klon), taur(klon), sss(klon))
+    if (activate_ocean_skin >= 1) allocate(tkt(klon), tks(klon), taur(klon), &
+         sss(klon))
 
   END SUBROUTINE phys_output_var_init
Index: LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/physiq_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/physiq_mod.F90	(revision 4175)
+++ LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/physiq_mod.F90	(revision 4482)
@@ -39,7 +39,7 @@
     USE ioipsl_getin_p_mod, ONLY : getin_p
     USE indice_sol_mod
-    USE infotrac_phy, ONLY: nqtot, nbtr, nqo, tracers, types_trac, nqCO2
+    USE infotrac_phy, ONLY: nqtot, nbtr, nqo, tracers, type_trac
     USE readTracFiles_mod, ONLY: addPhase
-    USE strings_mod,  ONLY: strIdx, strStack, int2str
+    USE strings_mod,  ONLY: strIdx
     USE iophy
     USE limit_read_mod, ONLY : init_limit_read
@@ -53,7 +53,8 @@
     USE pbl_surface_mod, ONLY : pbl_surface
     USE phyaqua_mod, only: zenang_an
+    USE phyetat0_mod, only: phyetat0
     USE phystokenc_mod, ONLY: offline, phystokenc
     USE phys_cal_mod, only: year_len, mth_len, days_elapsed, jh_1jan, &
-         year_cur, mth_cur,jD_cur, jH_cur, jD_ref, day_cur, hour
+         year_cur, mth_cur,jD_cur, jH_cur, jD_ref, day_cur, hour, calend
 !!  USE phys_local_var_mod, ONLY : a long list of variables
 !!              ==> see below, after "CPP Keys" section
@@ -68,14 +69,19 @@
     USE regr_horiz_time_climoz_m, ONLY: regr_horiz_time_climoz
     USE regr_pr_time_av_m, only: regr_pr_time_av
-    USE surface_data,     ONLY : type_ocean, ok_veget, landice_opt
-    USE time_phylmdz_mod, only: annee_ref, current_time, day_ini, day_ref, &
-          day_step_phy, itau_phy, pdtphys, raz_date, start_time, update_time
+    USE surface_data,     ONLY : type_ocean, ok_veget
+    USE time_phylmdz_mod, only: current_time, itau_phy, pdtphys, raz_date, update_time
     USE tracinca_mod, ONLY: config_inca
     USE tropopause_m,     ONLY: dyn_tropopause
     USE ice_sursat_mod,  ONLY: flight_init, airplane
     USE vampir
-    USE VERTICAL_LAYERS_MOD, ONLY: aps,bps, ap, bp
     USE write_field_phy
+#ifdef CPP_XIOS
+    USE wxios, ONLY: g_ctx, wxios_set_context
+#endif
+    USE atke_turbulence_ini_mod, ONLY : atke_ini
+    USE lscp_ini_mod, ONLY : lscp_ini
     USE lscp_mod, ONLY : lscp
+    USE wake_ini_mod, ONLY : wake_ini
+    USE yamada_ini_mod, ONLY : yamada_ini
     USE thermcell_ini_mod, ONLY : thermcell_ini
 
@@ -96,8 +102,16 @@
 
 
+#ifdef INCA
+    USE geometry_mod,      ONLY: longitude, latitude, boundslon, boundslat, ind_cell_glo
+    USE time_phylmdz_mod,  ONLY: ndays
+    USE infotrac_phy,      ONLY: nqCO2
+#endif
 #ifdef REPROBUS
-    USE CHEM_REP, ONLY : Init_chem_rep_xjour, &
-         d_q_rep,d_ql_rep,d_qi_rep,ptrop,ttrop, &
-         ztrop, gravit,itroprep, Z1,Z2,fac,B
+    USE chem_rep, ONLY: Init_chem_rep_xjour, d_q_rep, d_ql_rep, d_qi_rep, &
+                        ptrop, ttrop, ztrop, gravit, itroprep, Z1, Z2, fac, B
+#endif
+#if defined INCA || defined REPROBUS
+    USE time_phylmdz_mod,    ONLY: annee_ref, day_ini, day_ref, start_time
+    USE vertical_layers_mod, ONLY: aps, bps, ap, bp
 #endif
 
@@ -105,5 +119,5 @@
 #ifdef CPP_RRTM
     USE YOERAD, ONLY : NRADLP
-    USE YOESW, ONLY : RSUN
+!    USE YOESW, ONLY : RSUN
 #endif
 
@@ -115,6 +129,7 @@
 
 #ifdef CPP_XIOS
-    USE xios, ONLY: xios_update_calendar, xios_context_finalize, &
-            xios_get_field_attr, xios_field_is_active
+    USE xios, ONLY: xios_update_calendar, xios_context_finalize
+    USE xios, ONLY: xios_get_field_attr, xios_field_is_active, xios_context
+    USE xios, ONLY: xios_set_current_context
     USE wxios, ONLY: missing_val, missing_val_omp
 #endif
@@ -179,6 +194,6 @@
        d_t_ajsb,d_q_ajsb, &
        d_t_ajs,d_q_ajs,d_u_ajs,d_v_ajs, &
-       d_t_ajs_w,d_q_ajs_w, &
-       d_t_ajs_x,d_q_ajs_x, &
+!       d_t_ajs_w,d_q_ajs_w, &
+!       d_t_ajs_x,d_q_ajs_x, &
        !
        d_t_eva,d_q_eva,d_ql_eva,d_qi_eva, &
@@ -193,5 +208,5 @@
        d_ts, &
        !
-       d_t_oli,d_u_oli,d_v_oli, &
+!       d_t_oli,d_u_oli,d_v_oli, &
        d_t_oro,d_u_oro,d_v_oro, &
        d_t_oro_gw,d_u_oro_gw,d_v_oro_gw, &
@@ -312,4 +327,8 @@
        sij, &
        !
+       rneblsvol, &
+       zqsatl, zqsats, &
+       qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, &
+       Tcontr, qcontr, qcontr2, fcontrN, fcontrP, &
        cldemi,  &
        cldfra, cldtau, fiwc,  &
@@ -365,5 +384,4 @@
 #endif
        !
-
 
     IMPLICIT NONE
@@ -578,5 +596,4 @@
     !
     !
-    INTEGER debug
     INTEGER n
     !ym      INTEGER npoints
@@ -635,5 +652,5 @@
     ! Upmost level reached by deep convection and related variable (jyg)
     !
-    INTEGER izero
+!    INTEGER izero
     INTEGER k_upper_cv
     !------------------------------------------------------------------
@@ -813,5 +830,4 @@
     REAL beta_prec_fisrt(klon,klev) ! taux de conv de l'eau cond (fisrt)
     ! RomP <<<
-    REAL          :: calday
 
     !IM cf FH pour Tiedtke 080604
@@ -869,5 +885,4 @@
     !C      EXTERNAL o3cm      ! initialiser l'ozone
     EXTERNAL orbite    ! calculer l'orbite terrestre
-    EXTERNAL phyetat0  ! lire l'etat initial de la physique
     EXTERNAL phyredem  ! ecrire l'etat de redemarrage de la physique
     EXTERNAL suphel    ! initialiser certaines constantes
@@ -934,5 +949,5 @@
     real zqsat(klon,klev)
     !
-    INTEGER i, k, iq, j, nsrf, ll, l, itr
+    INTEGER i, k, iq, nsrf, l, itr
 #ifdef ISO
     real zxt_apres(ntraciso,klon)
@@ -1056,4 +1071,5 @@
 
     REAL picefra(klon,klev)
+    REAL zrel_oro(klon)
     !IM cf. AM 081204 END
     !
@@ -1265,4 +1281,5 @@
 
 #ifdef INCA
+    REAL :: calday, zxsnow_dummy(klon)
     ! set de variables utilisees pour l'initialisation des valeurs provenant de INCA
     REAL, DIMENSION(klon,klev,naero_grp,nbands) :: init_tauinca
@@ -1313,5 +1330,9 @@
     phys_tstep=NINT(pdtphys)
 #ifdef CPP_XIOS
-    IF (.NOT. debut .AND. is_omp_master) CALL xios_update_calendar(itap+1)
+! switch to XIOS LMDZ physics context 
+    IF (.NOT. debut .AND. is_omp_master) THEN
+       CALL wxios_set_context()
+       CALL xios_update_calendar(itap+1)
+    ENDIF
 #endif
 
@@ -1403,4 +1424,16 @@
           WRITE (lunout, *) ' ok_ice_sursat=y requires 4 H2O tracers ', &
                '(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)
+       ENDIF
+
+       IF (ok_plane_h2o.AND..NOT.ok_ice_sursat) THEN
+          WRITE (lunout, *) ' ok_plane_h2o=y requires ok_ice_sursat=y '
+          abort_message='see above'
+          CALL abort_physic(modname,abort_message,1)
+       ENDIF
+
+       IF (ok_plane_contrail.AND..NOT.ok_ice_sursat) THEN
+          WRITE (lunout, *) ' ok_plane_contrail=y requires ok_ice_sursat=y '
           abort_message='see above'
           CALL abort_physic(modname,abort_message,1)
@@ -1539,5 +1572,5 @@
        tau_overturning_th(:)=0.
 
-       IF (ANY(types_trac == 'inca') .OR. ANY(types_trac == 'inco')) THEN
+       IF (ANY(type_trac == ['inca','inco'])) THEN
           ! jg : initialisation jusqu'au ces variables sont dans restart
           ccm(:,:,:) = 0.
@@ -1847,13 +1880,20 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ! Nouvelle initialisation pour le rayonnement RRTM
-       !
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
        CALL iniradia(klon,klev,paprs(1,1:klev+1))
 
-
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+       CALL wake_ini(rg,rd,rv,prt_level)
+       CALL yamada_ini(klon,lunout,prt_level)
+       CALL atke_ini(prt_level, lunout, RG, RD, RPI, RCPD)
        CALL thermcell_ini(iflag_thermals,prt_level,tau_thermals,lunout, &
    &    RG,RD,RCPD,RKAPPA,RLVTT,RETV)
+       IF (ok_new_lscp) then
+           CALL lscp_ini(pdtphys,ok_ice_sursat)
+       endif
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
        !
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -2151,5 +2191,5 @@
        !c         ENDDO
        !
-       IF (ANY(types_trac == 'inca') .OR. ANY(types_trac == 'inco')) THEN  ! ModThL
+       IF (ANY(type_trac == ['inca','inco'])) THEN ! ModThL
 #ifdef INCA
           CALL VTe(VTphysiq)
@@ -2158,31 +2198,64 @@
           WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
 
-          CALL chemini(  &
-               rg, &
-               ra, &
-               cell_area, &
-               latitude_deg, &
-               longitude_deg, &
-               presnivs, &
-               calday, &
-               klon, &
-               nqtot, &
-               nqo+nqCO2, &
-               pdtphys, &
-               annee_ref, &
-               year_cur, &
-               day_ref,  &
-               day_ini, &
-               start_time, &
-               itau_phy, &
-               date0, &
-               io_lon, &
-               io_lat, &
-               chemistry_couple, &
-               init_source, &
-               init_tauinca, &
-               init_pizinca, &
-               init_cginca, &
-               init_ccminca)
+          call init_const_lmdz( &
+          ndays, nbsrf, is_oce,is_sic, is_ter,is_lic, calend, &
+          config_inca)
+
+          CALL init_inca_geometry( & 
+               longitude, latitude, &
+               boundslon, boundslat, &
+               cell_area, ind_cell_glo) 
+
+          if (grid_type==unstructured) THEN 
+             CALL chemini(  pplay, &
+                  nbp_lon, nbp_lat, &
+                  latitude_deg, &
+                  longitude_deg, &
+                  presnivs, &
+                  calday, &
+                  klon, &
+                  nqtot, &
+                  nqo+nqCO2, &
+                  pdtphys, &
+                  annee_ref, &
+                  year_cur, &
+                  day_ref,  &
+                  day_ini, &
+                  start_time, &
+                  itau_phy, &
+                  date0, &
+                  chemistry_couple, &
+                  init_source, &
+                  init_tauinca, &
+                  init_pizinca, &
+                  init_cginca, &
+                  init_ccminca)
+          ELSE
+             CALL chemini(  pplay, &
+                  nbp_lon, nbp_lat, &
+                  latitude_deg, &
+                  longitude_deg, &
+                  presnivs, &
+                  calday, &
+                  klon, &
+                  nqtot, &
+                  nqo+nqCO2, &
+                  pdtphys, &
+                  annee_ref, &
+                  year_cur, &
+                  day_ref,  &
+                  day_ini, &
+                  start_time, &
+                  itau_phy, &
+                  date0, &
+                  chemistry_couple, &
+                  init_source, &
+                  init_tauinca, &
+                  init_pizinca, &
+                  init_cginca, &
+                  init_ccminca, &
+                  io_lon, &
+                  io_lat)
+          ENDIF
 
 
@@ -2200,5 +2273,5 @@
        ENDIF
        !
-       IF (ANY(types_trac == 'repr')) THEN
+       IF (type_trac == 'repr') THEN
 #ifdef REPROBUS
           CALL chemini_rep(  &
@@ -2286,4 +2359,5 @@
 
 
+
     ENDIF
     !
@@ -2316,5 +2390,5 @@
 
     ! Update time and other variables in Reprobus
-    IF (ANY(types_trac == 'repr')) THEN
+    IF (type_trac == 'repr') THEN
 #ifdef REPROBUS
        CALL Init_chem_rep_xjour(jD_cur-jD_ref+day_ref)
@@ -2390,9 +2464,11 @@
           ql_seri(i,k) = qx(i,k,iliq)
           !CR: ATTENTION, on rajoute la variable glace
-          IF (nqo.eq.2) THEN
+          IF (nqo.EQ.2) THEN             !--vapour and liquid only
              qs_seri(i,k) = 0.
-          ELSE IF (nqo.eq.3) THEN
+             rneb_seri(i,k) = 0.
+          ELSE IF (nqo.EQ.3) THEN        !--vapour, liquid and ice
              qs_seri(i,k) = qx(i,k,isol)
-          ELSE IF (nqo.eq.4) THEN
+             rneb_seri(i,k) = 0.
+          ELSE IF (nqo.EQ.4) THEN        !--vapour, liquid, ice and rneb
              qs_seri(i,k) = qx(i,k,isol)
              rneb_seri(i,k) = qx(i,k,irneb)
@@ -2606,7 +2682,7 @@
              DO i = 1, klon
               if (q_seri(i,k).gt.3e-3) then
-              call iso_verif_positif(deltaD(xt_seri(iso_eau,i,k) &
+              call iso_verif_positif(deltaD(xt_seri(iso_HDO,i,k) &
      &           /q_seri(i,k))+400.0,'physiq 2045a')
-              call iso_verif_positif(deltaD(xt_ancien(iso_eau,i,k) &
+              call iso_verif_positif(deltaD(xt_ancien(iso_HDO,i,k) &
      &            /q_ancien(i,k))+400.0,'physiq 2045b')
 !              call iso_verif_egalite_choix(d_xt_dyn(iso_hdo,i,k) &
@@ -2658,6 +2734,6 @@
        ! !! RomP >>>   td dyn traceur
        IF (nqtot > nqo) d_tr_dyn(:,:,:)= 0.0
+       ! !! RomP <<<
        d_rneb_dyn(:,:)=0.0
-       ! !! RomP <<<
        ancien_ok = .TRUE.
     ENDIF
@@ -3526,5 +3602,5 @@
            endif ! if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then
            if ((iso_O17.gt.0).and.(iso_O18.gt.0)) then
-             if ((q_seri(i,k).gt.ridicule).and.(l.lt.nlevmaxO17)) then
+             if ((q_seri(i,k).gt.ridicule).and.(k.lt.nlevmaxO17)) then
               call iso_verif_aberrant_o17(xt_seri(iso_o17,i,k) &
      &           /q_seri(i,k),xt_seri(iso_o18,i,k) &
@@ -3582,5 +3658,5 @@
           !
           !>jyg
-          IF (ANY(types_trac == 'repr')) THEN
+          IF (type_trac == 'repr') THEN
              nbtr_tmp=ntra
           ELSE
@@ -4608,15 +4684,22 @@
     IF (ok_new_lscp) THEN
 
-    CALL lscp(phys_tstep,missing_val,paprs,pplay, &
+    !--mise à jour de flight_m et flight_h2o dans leur module
+    IF (ok_plane_h2o .OR. ok_plane_contrail) THEN
+      CALL airplane(debut,pphis,pplay,paprs,t_seri)
+    ENDIF
+
+    CALL lscp(klon,klev,phys_tstep,missing_val,paprs,pplay, &
          t_seri, q_seri,ptconv,ratqs, &
-         d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneb_seri, & 
+         d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, rneb_seri, &
          cldliq, picefra, rain_lsc, snow_lsc, &
-         pfrac_impa, pfrac_nucl, pfrac_1nucl, &
          frac_impa, frac_nucl, beta_prec_fisrt, &
          prfl, psfl, rhcl,  &
          zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
-         iflag_ice_thermo, ok_ice_sursat)
+         iflag_ice_thermo, ok_ice_sursat, zqsatl, zqsats, &
+         qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, &
+         Tcontr, qcontr, qcontr2, fcontrN, fcontrP )
 
     ELSE
+
     CALL fisrtilp(phys_tstep,paprs,pplay, &
          t_seri, q_seri,ptconv,ratqs, &
@@ -4690,5 +4773,5 @@
        ENDDO
     ENDDO
-    IF (nqo==3) THEN
+    IF (nqo >= 3) THEN
     DO k = 1, klev
        DO i = 1, klon
@@ -5112,5 +5195,5 @@
     ENDDO
 
-    IF (ANY(types_trac == 'inca') .OR. ANY(types_trac == 'inco')) THEN ! ModThL
+    IF (ANY(type_trac == ['inca','inco'])) THEN ! ModThL
 #ifdef INCA
        CALL VTe(VTphysiq)
@@ -5168,5 +5251,5 @@
 #endif
     ENDIF !type_trac = inca or inco
-    IF (ANY(types_trac == 'repr')) THEN
+    IF (type_trac == 'repr') THEN
 #ifdef REPROBUS
     !CALL chemtime_rep(itap+itau_phy-1, date0, dtime, itap)
@@ -5497,7 +5580,5 @@
           ! 
           !--interactive CO2 in ppm from carbon cycle
-          IF (carbon_cycle_rad.AND..NOT.debut) THEN
-            RCO2=RCO2_glo
-          ENDIF
+          IF (carbon_cycle_rad) RCO2=RCO2_glo
           !
           IF (prt_level .GE.10) THEN
@@ -5555,16 +5636,4 @@
  
 #ifndef CPP_XIOS
-          !--OB 30/05/2016 modified 21/10/2016
-          !--here we return swaero_diag and dryaod_diag to FALSE
-          !--and histdef will switch it back to TRUE if necessary
-          !--this is necessary to get the right swaero at first step
-          !--but only in the case of no XIOS as XIOS is covered elsewhere 
-          IF (debut) swaerofree_diag = .FALSE.
-          IF (debut) swaero_diag = .FALSE.
-          IF (debut) dryaod_diag = .FALSE.
-          !--IM 15/09/2017 here we return ok_4xCO2atm to FALSE
-          !--as for swaero_diag, see above
-          IF (debut) ok_4xCO2atm = .FALSE.
-
           !
           !IM 2eme calcul radiatif pour le cas perturbe ou au moins un
@@ -5734,6 +5803,11 @@
        DO i=1,klon
           itest(i)=0
-          !        IF ((zstd(i).gt.10.0)) THEN
-          IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
+          zrel_oro(i)=zstd(i)/(max(zsig(i),1.E-8)*sqrt(cell_area(i)))
+          !zrel_oro: relative mountain height wrt relief explained by mean slope
+          ! -> condition on zrel_oro can deactivate the drag on tilted planar terrains
+          !    such as ice sheets (work by V. Wiener)
+          ! zpmm_orodr_t and zstd_orodr_t are activation thresholds set by F. Lott to 
+          ! earn computation time but they are not physical.
+          IF (((zpic(i)-zmea(i)).GT.zpmm_orodr_t).AND.(zstd(i).GT.zstd_orodr_t).AND.(zrel_oro(i).LE.zrel_oro_t)) THEN
              itest(i)=1
              igwd=igwd+1
@@ -5788,5 +5862,9 @@
        DO i=1,klon
           itest(i)=0
-          IF ((zpic(i)-zmea(i)).GT.100.) THEN
+          !zrel_oro: relative mountain height wrt relief explained by mean slope
+          ! -> condition on zrel_oro can deactivate the lifting on tilted planar terrains
+          !    such as ice sheets (work by V. Wiener)
+          zrel_oro(i)=zstd(i)/(max(zsig(i),1.E-8)*sqrt(cell_area(i)))
+          IF (((zpic(i)-zmea(i)).GT.zpmm_orolf_t).AND.(zrel_oro(i).LE.zrel_oro_t)) THEN
              itest(i)=1
              igwd=igwd+1
@@ -6004,5 +6082,5 @@
 ! car on peut s'attendre a ce que les petites echelles produisent aussi de la TKE
 ! Mais attention, cela ne va pas dans le sens de la conservation de l'energie! 
-          IF (zstd(i).GT.1.0) THEN
+          IF ((zstd(i).GT.1.0).AND.(zrel_oro(i).LE.zrel_oro_t)) THEN
              itest(i)=1
              igwd=igwd+1
@@ -6016,5 +6094,5 @@
        DO i=1,klon
           itest(i)=0
-        IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
+        IF (((zpic(i)-zmea(i)).GT.zpmm_orodr_t).AND.(zstd(i).GT.zstd_orodr_t).AND.(zrel_oro(i).LE.zrel_oro_t)) THEN
              itest(i)=1
              igwd=igwd+1
@@ -6264,5 +6342,5 @@
     !
 
-    IF (ANY(types_trac=='repr')) THEN
+    IF (type_trac=='repr') THEN
 !MM pas d'impact, car on recupere q_seri,tr_seri,t_seri via phys_local_var_mod
 !MM                               dans Reprobus
@@ -6275,5 +6353,5 @@
     ELSE
        sh_in(:,:) = qx(:,:,ivap)
-       IF (nqo .EQ. 3) THEN 
+       IF (nqo >= 3) THEN 
           ch_in(:,:) = qx(:,:,iliq) + qx(:,:,isol)
        ELSE 
@@ -6362,17 +6440,16 @@
     ! Calculer le transport de l'eau et de l'energie (diagnostique)
     !
-    CALL transp (paprs,zxtsol, &
-         t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, zphi, &
-         ve, vq, ue, uq, vwat, uwat)
+    CALL transp (paprs,zxtsol, t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, zphi, &
+                 ue, ve, uq, vq, uwat, vwat)
     !
     !IM global posePB BEG
     IF(1.EQ.0) THEN
        !
-       CALL transp_lay (paprs,zxtsol, &
-            t_seri, q_seri, u_seri, v_seri, zphi, &
+       CALL transp_lay (paprs,zxtsol, t_seri, q_seri, u_seri, v_seri, zphi, &
             ve_lay, vq_lay, ue_lay, uq_lay)
        !
     ENDIF !(1.EQ.0) THEN
     !IM global posePB END
+    !
     ! Accumuler les variables a stocker dans les fichiers histoire:
     !
@@ -6385,5 +6462,5 @@
     d_t_ec(:,:)=0.
     forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
-    CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap),qx(:,:,iliq),qx(:,:,isol), &
+    CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx,ivap,iliq,isol, &
          u_seri,v_seri,t_seri,q_seri,ql_seri,qs_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
          zmasse,exner,d_t_ec)
@@ -6430,5 +6507,5 @@
 #endif
     !
-    IF (ANY(types_trac == 'inca') .OR. ANY(types_trac == 'inco')) THEN
+    IF (ANY(type_trac == ['inca','inco'])) THEN
 #ifdef INCA
        CALL VTe(VTphysiq)
@@ -6447,5 +6524,5 @@
             pphis, &
             zx_rh, &
-            aps, bps, ap, bp)
+            aps, bps, ap, bp, lafin)
 
        CALL VTe(VTinca)
@@ -6454,4 +6531,9 @@
     ENDIF
 
+    IF (type_trac == 'repr') THEN
+#ifdef REPROBUS
+        CALL coord_hyb_rep(paprs, pplay, aps, bps, ap, bp, cell_area)
+#endif
+    ENDIF
 
     !
@@ -6477,9 +6559,9 @@
           d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / phys_tstep
           !CR: on ajoute le contenu en glace
-          IF (nqo.ge.3) THEN
+          IF (nqo >= 3) THEN
              d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / phys_tstep
           ENDIF
           !--ice_sursat: nqo=4, on ajoute rneb
-          IF (nqo.eq.4) THEN
+          IF (nqo == 4) THEN
              d_qx(i,k,irneb) = ( rneb_seri(i,k) - qx(i,k,irneb) ) / phys_tstep
           ENDIF
@@ -6503,11 +6585,5 @@
       enddo ! DO k = 1, klev
     enddo !do ixt=1,ntraciso
-!#ifdef ISOVERIF
-!        write(*,*) 'physiq 6120: d_qx(1,1,:)=',d_qx(1,1,:)
-!        write(*,*) 'qx(1,1,:)=',qx(1,1,:)
-!        write(*,*) 'xt_seri(:,1,1)=',xt_seri(:,1,1)
-!#endif
-#endif
-! #ifdef ISO
+#endif
     ! DC: All iterations are cycled if nqtot==nqo, so no nqtot>nqo condition required
     itr = 0
@@ -6569,4 +6645,5 @@
     ql_ancien(:,:) = ql_seri(:,:)
     qs_ancien(:,:) = qs_seri(:,:)
+    rneb_ancien(:,:) = rneb_seri(:,:)
 #ifdef ISO
     xt_ancien(:,:,:)=xt_seri(:,:,:)
@@ -6828,4 +6905,13 @@
 #ifdef CPP_XIOS
        IF (is_omp_master) CALL xios_context_finalize
+
+#ifdef INCA
+       if (type_trac == 'inca') then 
+          IF (is_omp_master .and. grid_type==unstructured) THEN 
+             CALL finalize_inca
+          ENDIF
+       endif
+#endif
+
 #endif
        WRITE(lunout,*) ' physiq fin, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
Index: LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/surf_land_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/surf_land_mod.F90	(revision 4175)
+++ LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/surf_land_mod.F90	(revision 4482)
@@ -30,7 +30,5 @@
     USE dimphy
     USE surface_data, ONLY    : ok_veget
-! >> PC
     USE carbon_cycle_mod 
-! << PC
 
     ! See comments in each module surf_land_orchidee_xxx for compatiblity with ORCHIDEE
@@ -51,5 +49,11 @@
     USE surf_land_orchidee_nounstruct_mod
 #else
+#if ORCHIDEE_NOLIC
+    ! Compilation with cpp key ORCHIDEE_NOLIC
+    USE surf_land_orchidee_nolic_mod
+#else
+    ! Default version#else
     USE surf_land_orchidee_mod
+#endif
 #endif
 #endif
@@ -67,8 +71,6 @@
 #endif
 #endif
-
-! >> PC
+    
     USE print_control_mod, ONLY: lunout
-! << PC
 
     INCLUDE "dimsoil.h"
Index: LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/surf_land_orchidee_nolic_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/surf_land_orchidee_nolic_mod.F90	(revision 4482)
+++ LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/surf_land_orchidee_nolic_mod.F90	(revision 4482)
@@ -0,0 +1,1 @@
+link ../phylmd/surf_land_orchidee_nolic_mod.F90
Index: LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/surf_landice_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/surf_landice_mod.F90	(revision 4175)
+++ LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/surf_landice_mod.F90	(revision 4482)
@@ -144,4 +144,5 @@
     INTEGER                  :: i,j,nt
     REAL, DIMENSION(klon)    :: fqfonte,ffonte
+    REAL, DIMENSION(klon)    :: run_off_lic_frac
 #ifdef ISO       
       real, parameter :: t_coup = 273.15
@@ -484,8 +485,15 @@
 ! Send run-off on land-ice to coupler if coupled ocean.
 ! run_off_lic has been calculated in fonte_neige or surf_inlandsis
-!
-!****************************************************************************************
-    IF (type_ocean=='couple') THEN
-       CALL cpl_send_landice_fields(itime, knon, knindex, run_off_lic)
+! If landice_opt>=2, corresponding call is done from surf_land_orchidee
+!****************************************************************************************
+    IF (type_ocean=='couple' .AND. landice_opt .LT. 2) THEN
+       ! Compress fraction where run_off_lic is active (here all pctsrf(is_lic))
+       run_off_lic_frac(:)=0.0
+       DO j = 1, knon
+          i = knindex(j)
+          run_off_lic_frac(j) = pctsrf(i,is_lic)
+       ENDDO
+       
+       CALL cpl_send_landice_fields(itime, knon, knindex, run_off_lic, run_off_lic_frac)
     ENDIF
 
Index: LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/surf_ocean_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/surf_ocean_mod.F90	(revision 4175)
+++ LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/surf_ocean_mod.F90	(revision 4482)
@@ -20,6 +20,6 @@
        z0m, z0h, SFRWL, alb_dir_new, alb_dif_new, evap, fluxsens, fluxlat, & 
        tsurf_new, dflux_s, dflux_l, lmt_bils, &
-       flux_u1, flux_v1, delta_sst, delta_sal, ds_ns, dt_ns, dter, dser, tkt, tks, &
-       taur, sss &
+       flux_u1, flux_v1, delta_sst, delta_sal, ds_ns, dt_ns, dter, dser, &
+       dt_ds, tkt, tks, taur, sss &
 #ifdef ISO
         &       ,xtprecip_rain, xtprecip_snow,xtspechum,Roce, &
@@ -114,4 +114,16 @@
     ! minus foundation temperature. (Can be negative.) In K.
 
+    REAL, intent(inout):: dter(:) ! (knon)
+    ! Temperature variation in the diffusive microlayer, that is
+    ! ocean-air interface temperature minus subskin temperature. In
+    ! K.
+
+    REAL, intent(inout):: dser(:) ! (knon)
+    ! Salinity variation in the diffusive microlayer, that is
+    ! ocean-air interface salinity minus subskin salinity. In ppt.
+
+    real, intent(inout):: dt_ds(:) ! (knon)
+    ! (tks / tkt) * dTer, in K
+
     ! Output variables
     !**************************************************************************
@@ -129,13 +141,4 @@
     REAL, DIMENSION(klon), INTENT(OUT)       :: flux_u1, flux_v1
 
-    REAL, intent(out):: dter(:) ! (knon)
-    ! Temperature variation in the diffusive microlayer, that is
-    ! ocean-air interface temperature minus subskin temperature. In
-    ! K.
-
-    REAL, intent(out):: dser(:) ! (knon)
-    ! Salinity variation in the diffusive microlayer, that is
-    ! ocean-air interface salinity minus subskin salinity. In ppt.
-
     REAL, intent(out):: tkt(:) ! (knon)
     ! épaisseur (m) de la couche de diffusion thermique (microlayer)
@@ -152,4 +155,5 @@
     ! defined for subscripts 1:knon, but we have to declare it with
     ! size klon because of the coupling machinery.)
+
 #ifdef ISO
     REAL, DIMENSION(ntraciso,klon), INTENT(out)        :: xtevap ! isotopes in surface evaporation flux 
@@ -172,5 +176,5 @@
     real s_int(knon) ! ocean-air interface salinity, in ppt
 
-    !******************************************************************************
+    !**************************************************************************
 
 #ifdef ISO
@@ -213,7 +217,5 @@
     ENDIF
 
-
     rhoa = PS(:KNON) / (Rd * temp_air(:knon) * (1. + retv * spechum(:knon)))
-
     !******************************************************************************
     ! Switch according to type of ocean (couple, slab or forced)
@@ -232,5 +234,5 @@
             qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
             tsurf_new, dflux_s, dflux_l, sens_prec_liq, sss, delta_sal, rhoa, &
-            delta_sst)
+            delta_sst, dTer, dSer, dt_ds)
 
     CASE('slab')
@@ -376,7 +378,11 @@
        delta_sst = t_int - tsurf_new(:knon)
        delta_sal = s_int - sss(:knon)
-       if (activate_ocean_skin >= 2) tsurf_new(:knon) = t_int
+
+       if (activate_ocean_skin == 2) then
+          tsurf_new(:knon) = t_int
+          if (type_ocean == 'couple') dt_ds = (tks / tkt) * dter
+       end if
     end if
-
+    
   END SUBROUTINE surf_ocean
   !****************************************************************************
Index: LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/thermcell_down.F90
===================================================================
--- LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/thermcell_down.F90	(revision 4482)
+++ LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/thermcell_down.F90	(revision 4482)
@@ -0,0 +1,1 @@
+link ../phylmd/thermcell_down.F90
Index: LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/wake.F90
===================================================================
--- LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/wake.F90	(revision 4175)
+++ LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/wake.F90	(revision 4482)
@@ -2550,5 +2550,5 @@
     d_deltaqw, sigmaw, d_sigmaw, alpha)
   ! ------------------------------------------------------
-  ! Dtermination du coefficient alpha tel que les tendances
+  ! D\'etermination du coefficient alpha tel que les tendances
   ! corriges alpha*d_G, pour toutes les grandeurs G, correspondent
   ! a une humidite positive dans la zone (x) et dans la zone (w).
Index: LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/wake_ini_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/wake_ini_mod.F90	(revision 4482)
+++ LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/wake_ini_mod.F90	(revision 4482)
@@ -0,0 +1,1 @@
+link ../phylmd/wake_ini_mod.F90
Index: LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/yamada_ini_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/yamada_ini_mod.F90	(revision 4482)
+++ LMDZ6/branches/LMDZ_ECRad/libf/phylmdiso/yamada_ini_mod.F90	(revision 4482)
@@ -0,0 +1,1 @@
+link ../phylmd/yamada_ini_mod.F90
