Index: /LMDZ6/trunk/libf/dyn3d/conf_gcm.F90
===================================================================
--- /LMDZ6/trunk/libf/dyn3d/conf_gcm.F90	(revision 4357)
+++ /LMDZ6/trunk/libf/dyn3d/conf_gcm.F90	(revision 4358)
@@ -599,14 +599,4 @@
      type_trac = 'lmdz'
      CALL getin('type_trac',type_trac)
-
-     !Config  Key  = config_inca
-     !Config  Desc = Choix de configuration de INCA
-     !Config  Def  = none
-     !Config  Help = Choix de configuration de INCA :
-     !Config         'none' = sans INCA
-     !Config         'chem' = INCA avec calcul de chemie
-     !Config         'aero' = INCA avec calcul des aerosols 
-     config_inca = 'none'
-     CALL getin('config_inca',config_inca)
 
      !Config  Key  = ok_dynzon 
@@ -672,5 +662,4 @@
      write(lunout,*)' offline = ', offline
      write(lunout,*)' type_trac = ', type_trac
-     write(lunout,*)' config_inca = ', config_inca
      write(lunout,*)' ok_dynzon = ', ok_dynzon 
      write(lunout,*)' ok_dyn_ins = ', ok_dyn_ins 
@@ -795,14 +784,4 @@
      type_trac = 'lmdz'
      CALL getin('type_trac',type_trac)
-
-     !Config  Key  = config_inca
-     !Config  Desc = Choix de configuration de INCA
-     !Config  Def  = none
-     !Config  Help = Choix de configuration de INCA :
-     !Config         'none' = sans INCA
-     !Config         'chem' = INCA avec calcul de chemie
-     !Config         'aero' = INCA avec calcul des aerosols 
-     config_inca = 'none'
-     CALL getin('config_inca',config_inca)
 
      !Config  Key  = ok_dynzon 
@@ -912,5 +891,4 @@
      write(lunout,*)' offline = ', offline
      write(lunout,*)' type_trac = ', type_trac
-     write(lunout,*)' config_inca = ', config_inca
      write(lunout,*)' ok_dynzon = ', ok_dynzon
      write(lunout,*)' ok_dyn_ins = ', ok_dyn_ins 
Index: /LMDZ6/trunk/libf/dyn3d_common/control_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/dyn3d_common/control_mod.F90	(revision 4357)
+++ /LMDZ6/trunk/libf/dyn3d_common/control_mod.F90	(revision 4358)
@@ -29,5 +29,4 @@
   INTEGER,SAVE :: ip_ebil_dyn
   LOGICAL,SAVE :: offline
-  CHARACTER(len=4),SAVE :: config_inca
   CHARACTER(len=10),SAVE :: planet_type ! planet type ('earth','mars',...)
   LOGICAL,SAVE :: output_grads_dyn ! output dynamics diagnostics in
Index: /LMDZ6/trunk/libf/dyn3d_common/infotrac.F90
===================================================================
--- /LMDZ6/trunk/libf/dyn3d_common/infotrac.F90	(revision 4357)
+++ /LMDZ6/trunk/libf/dyn3d_common/infotrac.F90	(revision 4358)
@@ -67,17 +67,17 @@
 !  | phase       | Phases list ("g"as / "l"iquid / "s"olid)             | /           | [g][l][s]              |
 !  | component   | Name(s) of the merged/cumulated section(s)           | /           | coma-separated names   |
-!  | iadv        | Advection scheme number                              | iadv        | 1-20,30 exc. 3-9,15,19 |
 !  | iGeneration | Generation (>=1)                                     | /           |                        |
-!  | isAdvected  | advected tracers flag (.TRUE. if iadv >= 0)          | /           | nqtrue  .TRUE. values  |
-!  | isInPhysics | tracers not extracted from the main table in physics | /           | nqtottr .TRUE. values  |
 !  | iqParent    | Index of the parent tracer                           | iqpere      | 1:nqtot                |
 !  | iqDescen    | Indexes of the childs       (all generations)        | iqfils      | 1:nqtot                |
 !  | nqDescen    | Number of the descendants   (all generations)        | nqdesc      | 1:nqtot                |
 !  | nqChildren  | Number of childs            (1st generation only)    | nqfils      | 1:nqtot                |
+!  | keys        | key/val pairs accessible with "getKey" routine       | /           |                        |
+!  | iadv        | Advection scheme number                              | iadv        | 1,2,10-20(exc.15,19),30|
+!  | isAdvected  | advected tracers flag (.TRUE. if iadv >= 0)          | /           | nqtrue  .TRUE. values  |
+!  | isInPhysics | tracers not extracted from the main table in physics | /           | nqtottr .TRUE. values  |
 !  | iso_iGroup  | Isotopes group index in isotopes(:)                  | /           | 1:nbIso                |
 !  | iso_iName   | Isotope  name  index in isotopes(iso_iGroup)%trac(:) | iso_indnum  | 1:niso                 |
 !  | iso_iZone   | Isotope  zone  index in isotopes(iso_iGroup)%zone(:) | zone_num    | 1:nzone                |
 !  | iso_iPhas   | Isotope  phase index in isotopes(iso_iGroup)%phas(:) | phase_num   | 1:nphas                |
-!  | keys        | key/val pairs accessible with "getKey" routine       | /           |                        |
 !  +-------------+------------------------------------------------------+-------------+------------------------+
 !
@@ -114,5 +114,5 @@
 
 SUBROUTINE init_infotrac
-   USE control_mod, ONLY: planet_type, config_inca
+   USE control_mod, ONLY: planet_type
 #ifdef REPROBUS
    USE CHEM_REP,    ONLY: Init_chem_rep_trac
@@ -182,5 +182,5 @@
       msg1 = 'For type_trac = "'//TRIM(types_trac(it))//'":'
       SELECT CASE(types_trac(it))
-         CASE('inca'); CALL msg(TRIM(msg1)//' coupling with INCA chemistry model, config_inca='//config_inca, modname)
+         CASE('inca'); CALL msg(TRIM(msg1)//' coupling with INCA chemistry model',        modname)
          CASE('inco'); CALL msg(TRIM(msg1)//' coupling jointly with INCA and CO2 cycle',  modname)
          CASE('repr'); CALL msg(TRIM(msg1)//' coupling with REPROBUS chemistry model',    modname)
@@ -191,8 +191,4 @@
       END SELECT
 
-      !--- COHERENCE TEST BETWEEN "type_trac" AND "config_inca"
-      IF(ANY(['inca', 'inco'] == types_trac(it)) .AND. ALL(['aero', 'aeNP', 'chem'] /= config_inca)) &
-         CALL abort_gcm(modname, 'Incoherence between type_trac and config_inca. Please modify "run.def"',1)
-
       !--- COHERENCE TEST BETWEEN "type_trac" AND PREPROCESSING KEYS
       SELECT CASE(types_trac(it))
@@ -214,8 +210,4 @@
    END DO
    !---------------------------------------------------------------------------------------------------------------------------
-
-   !--- DISABLE "config_inca" OPTION FOR A RUN WITHOUT "INCA" IF IT DIFFERS FROM "none"
-   IF(fmsg('Setting config_inca="none" as you do not couple with INCA model', &
-         modname, ALL(types_trac /= 'inca') .AND. ALL(types_trac /= 'inco') .AND. config_inca /= 'none')) config_inca = 'none'
 
    nqCO2 = COUNT( [ANY(types_trac == 'inco') .OR. (ANY(types_trac == 'co2i') .AND. ANY(types_trac == 'inca'))] )
Index: /LMDZ6/trunk/libf/dyn3dmem/conf_gcm.F90
===================================================================
--- /LMDZ6/trunk/libf/dyn3dmem/conf_gcm.F90	(revision 4357)
+++ /LMDZ6/trunk/libf/dyn3dmem/conf_gcm.F90	(revision 4358)
@@ -644,14 +644,4 @@
      type_trac = 'lmdz'
      CALL getin('type_trac',type_trac)
-
-     !Config  Key  = config_inca
-     !Config  Desc = Choix de configuration de INCA
-     !Config  Def  = none
-     !Config  Help = Choix de configuration de INCA :
-     !Config         'none' = sans INCA
-     !Config         'chem' = INCA avec calcul de chemie
-     !Config         'aero' = INCA avec calcul des aerosols 
-     config_inca = 'none'
-     CALL getin('config_inca',config_inca)
 
      !Config  Key  = ok_dynzon 
@@ -725,5 +715,4 @@
      write(lunout,*)' offline = ', offline
      write(lunout,*)' type_trac = ', type_trac
-     write(lunout,*)' config_inca = ', config_inca
      write(lunout,*)' ok_dynzon = ', ok_dynzon 
      write(lunout,*)' ok_dyn_ins = ', ok_dyn_ins 
@@ -855,14 +844,4 @@
      type_trac = 'lmdz'
      CALL getin('type_trac',type_trac)
-
-     !Config  Key  = config_inca
-     !Config  Desc = Choix de configuration de INCA
-     !Config  Def  = none
-     !Config  Help = Choix de configuration de INCA :
-     !Config         'none' = sans INCA
-     !Config         'chem' = INCA avec calcul de chemie
-     !Config         'aero' = INCA avec calcul des aerosols 
-     config_inca = 'none'
-     CALL getin('config_inca',config_inca)
 
      !Config  Key  = ok_dynzon 
@@ -1006,5 +985,4 @@
      write(lunout,*)' offline = ', offline
      write(lunout,*)' type_trac = ', type_trac
-     write(lunout,*)' config_inca = ', config_inca
      write(lunout,*)' ok_dynzon = ', ok_dynzon
      write(lunout,*)' ok_dyn_ins = ', ok_dyn_ins 
Index: /LMDZ6/trunk/libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90	(revision 4357)
+++ /LMDZ6/trunk/libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90	(revision 4358)
@@ -30,5 +30,5 @@
   USE mod_phys_lmdz_omp_data, ONLY: klon_omp
 #endif
-  USE control_mod, ONLY: dayref,anneeref,day_step,nday,offline, iphysiq, config_inca
+  USE control_mod, ONLY: dayref,anneeref,day_step,nday,offline, iphysiq
   USE inifis_mod, ONLY: inifis
   USE time_phylmdz_mod, ONLY: init_time
Index: /LMDZ6/trunk/libf/misc/readTracFiles_mod.f90
===================================================================
--- /LMDZ6/trunk/libf/misc/readTracFiles_mod.f90	(revision 4357)
+++ /LMDZ6/trunk/libf/misc/readTracFiles_mod.f90	(revision 4358)
@@ -54,17 +54,17 @@
     CHARACTER(LEN=maxlen) :: phase       = 'g'             !--- Phase ('g'as / 'l'iquid / 's'olid)
     CHARACTER(LEN=maxlen) :: component   = ''              !--- Coma-separated list of components (Ex: lmdz,inca)
-    INTEGER               :: iadv        = 10              !--- Advection scheme used
     INTEGER               :: iGeneration = -1              !--- Generation number (>=0)
-    LOGICAL               :: isAdvected  = .FALSE.         !--- "true" tracers: iadv > 0.   COUNT(isAdvected )=nqtrue
-    LOGICAL               :: isInPhysics = .TRUE.          !--- "true" tracers: in tr_seri. COUNT(isInPhysics)=nqtottr
     INTEGER               :: iqParent    = 0               !--- Parent index
     INTEGER,  ALLOCATABLE :: iqDescen(:)                   !--- Descendants index (in growing generation order)
     INTEGER               :: nqDescen    = 0               !--- Number of descendants (all generations)
     INTEGER               :: nqChildren  = 0               !--- Number of children  (first generation)
+    TYPE(keys_type)       :: keys                          !--- <key>=<val> pairs vector
+    INTEGER               :: iadv        = 10              !--- Advection scheme used
+    LOGICAL               :: isAdvected  = .FALSE.         !--- "true" tracers: iadv > 0.   COUNT(isAdvected )=nqtrue
+    LOGICAL               :: isInPhysics = .TRUE.          !--- "true" tracers: in tr_seri. COUNT(isInPhysics)=nqtottr
     INTEGER               :: iso_iGroup  = 0               !--- Isotopes group index in isotopes(:)
     INTEGER               :: iso_iName   = 0               !--- Isotope  name  index in isotopes(iso_iGroup)%trac(:)
     INTEGER               :: iso_iZone   = 0               !--- Isotope  zone  index in isotopes(iso_iGroup)%zone(:)
     INTEGER               :: iso_iPhase  = 0               !--- Isotope  phase index in isotopes(iso_iGroup)%phase
-    TYPE(keys_type)       :: keys                          !--- <key>=<val> pairs vector
   END TYPE trac_type
 !------------------------------------------------------------------------------------------------------------------------------
@@ -1895,6 +1895,9 @@
   ip = getiPhase(newName)                                                      !--- Phase index
   IF(ip /= 0 .AND. PRESENT(iPhase)) iPhase = ip                                !--- Returning phase index
-  ix = strIdx(newH2OIso, newName)                                              !--- Index in the known H2O isotopes list
-  IF(ix /= 0) oldName = 'H2O'//'_'//TRIM(oldH2OIso(ix))                        !=== WATER ISOTOPE WITHOUT PHASE
+  ix = strIdx(newH2OIso, strHead(newName, '_'))                                !--- Index in the known H2O isotopes list
+  IF(ix /= 0) THEN
+    oldName = 'H2O'//'_'//TRIM(oldH2OIso(ix))                                  !=== WATER ISOTOPE WITHOUT PHASE
+    IF(newH2OIso(ix)/=newName) oldName=TRIM(oldName)//'_'//strTail(newName,'_')!=== WATER ISOTOPIC TAGGING TRACER WITHOUT PHASE
+  END IF
   IF(ix /= 0 .OR. ip == 0)           RETURN
   oldName = 'H2O'//old_phases(ip:ip)
Index: /LMDZ6/trunk/libf/misc/strings_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/misc/strings_mod.F90	(revision 4357)
+++ /LMDZ6/trunk/libf/misc/strings_mod.F90	(revision 4358)
@@ -505,5 +505,5 @@
 !=== Count the number of elements separated by "delimiter" in list "rawList". =================================================
 !==============================================================================================================================
-LOGICAL FUNCTION strCount_11(rawList, delimiter, nb, lSc) RESULT(out)
+LOGICAL FUNCTION strCount_11(rawList, delimiter, nb, lSc) RESULT(lerr)
   CHARACTER(LEN=*),  INTENT(IN)  :: rawList
   CHARACTER(LEN=*),  INTENT(IN)  :: delimiter
@@ -513,8 +513,8 @@
   LOGICAL :: ll
   ll = .FALSE.; IF(PRESENT(lSc)) ll = lSc
-  out = strCount_1m(rawList, [delimiter], nb, ll)
+  lerr = strCount_1m(rawList, [delimiter], nb, ll)
 END FUNCTION strCount_11
 !==============================================================================================================================
-LOGICAL FUNCTION strCount_m1(rawList, delimiter, nb, lSc) RESULT(out)
+LOGICAL FUNCTION strCount_m1(rawList, delimiter, nb, lSc) RESULT(lerr)
   CHARACTER(LEN=*),     INTENT(IN)  :: rawList(:)
   CHARACTER(LEN=*),     INTENT(IN)  :: delimiter
@@ -525,8 +525,8 @@
   INTEGER :: id
   ll  = .FALSE.; IF(PRESENT(lSc)) ll = lSc .AND. INDEX('+-', delimiter) /= 0
-  out = .TRUE.
+  lerr = .TRUE.
   ALLOCATE(nb(SIZE(rawList)))
   DO id = 1, SIZE(rawList)
-    out = out .AND. strCount_1m(rawList(id), [delimiter], nb(id), ll)
+    lerr = lerr .AND. strCount_1m(rawList(id), [delimiter], nb(id), ll)
   END DO
 END FUNCTION strCount_m1
Index: /LMDZ6/trunk/libf/phylmd/infotrac_phy.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/infotrac_phy.F90	(revision 4357)
+++ /LMDZ6/trunk/libf/phylmd/infotrac_phy.F90	(revision 4358)
@@ -76,9 +76,9 @@
 !  | nqDescen    | Number of the descendants   (all generations)        | nqdesc      | 1:nqtot                |
 !  | nqChildren  | Number of childs            (1st generation only)    | nqfils      | 1:nqtot                |
+!  | keys        | key/val pairs accessible with "getKey" routine       | /           |                        |
 !  | iso_iGroup  | Isotopes group index in isotopes(:)                  | /           | 1:nbIso                |
 !  | iso_iName   | Isotope  name  index in isotopes(iso_iGroup)%trac(:) | iso_indnum  | 1:niso                 |
 !  | iso_iZone   | Isotope  zone  index in isotopes(iso_iGroup)%zone(:) | zone_num    | 1:nzone                |
 !  | iso_iPhas   | Isotope  phase index in isotopes(iso_iGroup)%phas(:) | phase_num   | 1:nphas                |
-!  | keys        | key/val pairs accessible with "getKey" routine       | /           |                        |
 !  +-------------+------------------------------------------------------+-------------+------------------------+
 !
@@ -123,8 +123,7 @@
 
 SUBROUTINE init_infotrac_phy
-   USE control_mod, ONLY: planet_type, config_inca
    USE ioipsl_getin_p_mod, ONLY: getin_p
 #ifdef REPROBUS
-   USE CHEM_REP,    ONLY: Init_chem_rep_trac
+   USE CHEM_REP, ONLY: Init_chem_rep_trac
 #endif
    IMPLICIT NONE
@@ -200,5 +199,5 @@
       msg1 = 'For type_trac = "'//TRIM(types_trac(it))//'":'
       SELECT CASE(types_trac(it))
-         CASE('inca'); CALL msg(TRIM(msg1)//' coupling with INCA chemistry model, config_inca='//config_inca, modname)
+         CASE('inca'); CALL msg(TRIM(msg1)//' coupling with INCA chemistry model',        modname)
          CASE('inco'); CALL msg(TRIM(msg1)//' coupling jointly with INCA and CO2 cycle',  modname)
          CASE('repr'); CALL msg(TRIM(msg1)//' coupling with REPROBUS chemistry model',    modname)
@@ -209,8 +208,4 @@
       END SELECT
 
-      !--- COHERENCE TEST BETWEEN "type_trac" AND "config_inca"
-      IF(ANY(['inca', 'inco'] == types_trac(it)) .AND. ALL(['aero', 'aeNP', 'chem'] /= config_inca)) &
-         CALL abort_gcm(modname, 'Incoherence between type_trac and config_inca. Please modify "run.def"',1)
-
       !--- COHERENCE TEST BETWEEN "type_trac" AND PREPROCESSING KEYS
       SELECT CASE(types_trac(it))
@@ -236,8 +231,4 @@
    END IF
 !##############################################################################################################################
-
-   !--- DISABLE "config_inca" OPTION FOR A RUN WITHOUT "INCA" IF IT DIFFERS FROM "none"
-   IF(fmsg('Setting config_inca="none" as you do not couple with INCA model', &
-         modname, ALL(types_trac /= 'inca') .AND. ALL(types_trac /= 'inco') .AND. config_inca /= 'none')) config_inca = 'none'
 
    nqCO2 = COUNT( [ANY(types_trac == 'inco') .OR. (ANY(types_trac == 'co2i') .AND. ANY(types_trac == 'inca'))] )
Index: DZ6/trunk/libf/phylmd/phyetat0.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phyetat0.F90	(revision 4357)
+++ 	(revision )
@@ -1,667 +1,0 @@
-! $Id$
-
-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
-  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, &
-       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, type_trac, types_trac, tracers
-  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
-
-  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)
-
-  ! 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(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,rneb_ancien,"RNEBANCIEN","RNEBANCIEN",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(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(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)
-        found=phyetat0_get(1,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
-        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
-
-
-!===========================================
-!  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 
-
-  CALL pbl_surface_init(fder, snow, qsurf, tsoil)
-
-  ! 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)
-
-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/trunk/libf/phylmd/phyetat0_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phyetat0_mod.F90	(revision 4358)
+++ /LMDZ6/trunk/libf/phylmd/phyetat0_mod.F90	(revision 4358)
@@ -0,0 +1,729 @@
+! $Id$
+
+MODULE phyetat0_mod
+
+  PRIVATE
+  PUBLIC :: phyetat0, phyetat0_get, phyetat0_srf
+
+  INTERFACE phyetat0_get
+    MODULE PROCEDURE phyetat0_get10, phyetat0_get20, phyetat0_get11, phyetat0_get21
+  END INTERFACE phyetat0_get
+  INTERFACE phyetat0_srf
+    MODULE PROCEDURE phyetat0_srf20, phyetat0_srf30, phyetat0_srf21, phyetat0_srf31
+  END INTERFACE phyetat0_srf
+
+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
+  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, &
+       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, type_trac, types_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
+
+  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)
+
+  ! 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(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(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
+
+
+!===========================================
+!  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.)
+     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 
+
+  CALL pbl_surface_init(fder, snow, qsurf, tsoil)
+
+  ! 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)
+
+END SUBROUTINE phyetat0
+
+!==============================================================================
+LOGICAL FUNCTION phyetat0_get10(field, name, descr, default) RESULT(lFound)
+! Read a field. Check whether reading succeded and use default value if not.
+  IMPLICIT NONE
+  REAL,             INTENT(INOUT) :: field(:) ! klon
+  CHARACTER(LEN=*), INTENT(IN)    :: name
+  CHARACTER(LEN=*), INTENT(IN)    :: descr
+  REAL,             INTENT(IN)    :: default
+!------------------------------------------------------------------------------
+  REAL :: fld(SIZE(field),1)
+  lFound = phyetat0_get21(fld, [name], descr, default); field = fld(:,1)
+END FUNCTION phyetat0_get10
+!==============================================================================
+LOGICAL FUNCTION phyetat0_get20(field, name, descr, default) RESULT(lFound)
+! Same as phyetat0_get11, field on multiple levels.
+  IMPLICIT NONE
+  REAL,             INTENT(INOUT) :: field(:,:) ! klon, nlev
+  CHARACTER(LEN=*), INTENT(IN)    :: name
+  CHARACTER(LEN=*), INTENT(IN)    :: descr
+  REAL,             INTENT(IN)    :: default
+!-----------------------------------------------------------------------------
+  lFound = phyetat0_get21(field, [name], descr, default)
+END FUNCTION phyetat0_get20
+!==============================================================================
+LOGICAL FUNCTION phyetat0_get11(field, name, descr, default) RESULT(lFound)
+! Same as phyetat0_get11, multiple names.
+  IMPLICIT NONE
+  REAL,             INTENT(INOUT) :: field(:) ! klon
+  CHARACTER(LEN=*), INTENT(IN)    :: name(:)
+  CHARACTER(LEN=*), INTENT(IN)    :: descr
+  REAL,             INTENT(IN)    :: default
+!-----------------------------------------------------------------------------
+  REAL :: fld(SIZE(field),1)
+  lFound = phyetat0_get21(fld, name, descr, default); field = fld(:,1)
+END FUNCTION phyetat0_get11
+!==============================================================================
+LOGICAL FUNCTION phyetat0_get21(field, name, descr, default, tname) RESULT(lFound)
+! Same as phyetat0_get11, field on multiple levels, multiple names.
+  USE iostart,           ONLY: get_field
+  USE print_control_mod, ONLY: lunout
+  IMPLICIT NONE
+  REAL,             INTENT(INOUT) :: field(:,:) ! klon, nlev
+  CHARACTER(LEN=*), INTENT(IN)    :: name(:)
+  CHARACTER(LEN=*), INTENT(IN)    :: descr
+  REAL,             INTENT(IN)    :: default
+  CHARACTER(LEN=*), OPTIONAL, INTENT(OUT) :: tname
+!-----------------------------------------------------------------------------
+  CHARACTER(LEN=LEN(name)) :: tnam
+  INTEGER :: i
+  DO i = 1, SIZE(name)
+    CALL get_field(TRIM(name(i)), field, lFound)
+    IF(lFound) EXIT
+    WRITE(lunout,*) "phyetat0: Missing field <",TRIM(name(i)),"> "
+  END DO
+  IF(.NOT.lFound) THEN
+    WRITE(lunout,*) "Slightly distorted start ; continuing."
+    field(:,:) = default
+    tnam = name(1)
+  ELSE
+    tnam = name(i)
+  END IF
+  WRITE(lunout,'(2(a,ES14.7))') 'phyetat0: '//TRIM(tnam)//' ('//TRIM(descr)//') min/max=', &
+    MINval(field),' ',MAXval(field)
+  IF(PRESENT(tname)) tname = tnam
+END FUNCTION phyetat0_get21
+!==============================================================================
+LOGICAL FUNCTION phyetat0_srf20(field, name, descr, default) RESULT(lFound)
+! Read a field per sub-surface.
+! Check whether reading succeded and use default value if not.
+  IMPLICIT NONE
+  REAL,             INTENT(INOUT) :: field(:,:)
+  CHARACTER(LEN=*), INTENT(IN)    :: name
+  CHARACTER(LEN=*), INTENT(IN)    :: descr
+  REAL,             INTENT(IN)    :: default
+!-----------------------------------------------------------------------------
+  REAL :: fld(SIZE(field,1),1,SIZE(field,2))
+  lFound = phyetat0_srf31(fld, [name], descr, default); field = fld(:,1,:)
+END FUNCTION phyetat0_srf20
+
+!==============================================================================
+LOGICAL FUNCTION phyetat0_srf30(field, name, descr, default) RESULT(lFound)
+! Same as phyetat0_sfr11, multiple names tested one after the other.
+  IMPLICIT NONE
+  REAL,             INTENT(INOUT) :: field(:,:,:)
+  CHARACTER(LEN=*), INTENT(IN)    :: name
+  CHARACTER(LEN=*), INTENT(IN)    :: descr
+  REAL,             INTENT(IN)    :: default
+!-----------------------------------------------------------------------------
+  lFound = phyetat0_srf31(field, [name], descr, default)
+END FUNCTION phyetat0_srf30
+
+!==============================================================================
+LOGICAL FUNCTION phyetat0_srf21(field, name, descr, default) RESULT(lFound)
+! Same as phyetat0_sfr11, field on multiple levels.
+  IMPLICIT NONE
+  REAL,             INTENT(INOUT) :: field(:,:)
+  CHARACTER(LEN=*), INTENT(IN)    :: name(:)
+  CHARACTER(LEN=*), INTENT(IN)    :: descr
+  REAL,             INTENT(IN)    :: default
+!-----------------------------------------------------------------------------
+  REAL :: fld(SIZE(field,1),1,SIZE(field,2))
+  lFound = phyetat0_srf31(fld, name, descr, default); field = fld(:,1,:)
+END FUNCTION phyetat0_srf21
+
+!==============================================================================
+LOGICAL FUNCTION phyetat0_srf31(field, name, descr, default) RESULT(lFound)
+! Same as phyetat0_sfr11, field on multiple levels, multiple names tested one after the other.
+  USE iostart,           ONLY: get_field
+  USE print_control_mod, ONLY: lunout
+  USE strings_mod,       ONLY: int2str, maxlen
+  IMPLICIT NONE
+  REAL,             INTENT(INOUT) :: field(:,:,:)
+  CHARACTER(LEN=*), INTENT(IN)    :: name(:)
+  CHARACTER(LEN=*), INTENT(IN)    :: descr
+  REAL,             INTENT(IN)    :: default
+!-----------------------------------------------------------------------------
+  INTEGER :: nsrf, i
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: nam(:)
+  CHARACTER(LEN=maxlen) :: tname, des
+  IF(SIZE(field,3)>99) CALL abort_physic("phyetat0", "Too much sub-cells", 1)
+  DO nsrf = 1, SIZE(field,3)
+    nam = [(TRIM(name(i))//TRIM(int2str(nsrf,2)), i=1, SIZE(name))]
+    des = TRIM(descr)//" srf:"//int2str(nsrf,2)
+    lFound = phyetat0_get21(field(:,:,nsrf), nam, TRIM(des), default, tname)
+  END DO
+  WRITE(lunout,'(2(a,ES14.7))') 'phyetat0: '//TRIM(tname)//' ('//TRIM(descr)//') min/max=', &
+    MINval(field),' ',MAXval(field)
+END FUNCTION phyetat0_srf31
+
+END MODULE phyetat0_mod
+
Index: /LMDZ6/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4357)
+++ /LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4358)
@@ -54,4 +54,5 @@
     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, &
@@ -802,5 +803,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
Index: /LMDZ6/trunk/libf/phylmd/tracinca_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/tracinca_mod.F90	(revision 4357)
+++ /LMDZ6/trunk/libf/phylmd/tracinca_mod.F90	(revision 4358)
@@ -12,4 +12,5 @@
                      ! config_inca='chem' => INCA with chemistry
                      ! config_inca='aero' => INCA with aerosols
+                     ! config_inca='aeNP' => INCA with aerosols NP (?)
 CONTAINS
 
@@ -17,6 +18,5 @@
     ! This subroutine initialize some control varaibles. 
 
-    USE infotrac_phy, ONLY: nbtr
-    USE ioipsl_getin_p_mod, ONLY: getin_p
+    USE infotrac_phy, ONLY: nbtr, types_trac
     IMPLICIT NONE
     
@@ -25,8 +25,17 @@
     LOGICAL,INTENT(OUT) :: lessivage
     
-    
     ! Initialization
     lessivage  =.FALSE.
     aerosol(:) = .FALSE.
+
+    !--- COHERENCE TEST BETWEEN "type_trac" AND "config_inca"
+    IF((ANY(types_trac == 'inca') .OR. ANY(types_trac == 'inco')) .AND. ALL(config_inca /= ['aero', 'aeNP', 'chem']))&
+       CALL abort_gcm('tracinca_init', 'INCA enabled, but unknown config_inca = "'//TRIM(config_inca)//'".'          &
+                             //'Please modify "run.def"', 1)
+
+    !--- PROBLEM IF "config_inca" DIFFERS FROM "none" AND INCA HAS NOT BEEN ACTIVATED
+    IF(ALL(types_trac /= 'inca') .AND. ALL(types_trac /= 'inco')  .AND.     config_inca /= 'none') &
+       CALL abort_gcm('tracinca_init', 'INCA disabled, but config_inca = "'//TRIM(config_inca)//'" should be "none".'&
+                             //'Please modify "run.def"', 1)
 
   END SUBROUTINE tracinca_init
Index: /LMDZ6/trunk/libf/phylmdiso/isotopes_routines_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/isotopes_routines_mod.F90	(revision 4357)
+++ /LMDZ6/trunk/libf/phylmdiso/isotopes_routines_mod.F90	(revision 4358)
@@ -16508,261 +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_mod,      ONLY: phyetat0_get, phyetat0_srf
+   USE readTracFiles_mod, ONLY: new2oldH2O
+   USE strings_mod,       ONLY: strIdx, strHead, 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
- USE readTracFiles_mod, ONLY: newH2Oiso, oldH2Oiso
- USE strings_mod, ONLY: strIdx, strHead, strTail
-
-#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, oldIso
-        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
+   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)
-      k = strIdx(newH2Oiso, strHead(outiso, '_'))
-      oldIso = outiso; IF(k /= 0) oldIso = oldH2Oiso(k)
-      IF(INDEX(outiso, '_') /= 0) THEN
-        outiso = TRIM(outiso)//TRIM(strTail(outiso, '_'))
-        oldIso = TRIM(oldIso)//TRIM(strTail(outiso, '_'))
+      oldIso = strTail(new2oldH2O(outiso), '_', lFirst=.TRUE.)
+      ! on lit seulement si ixt<=niso ou si on initialise les traceurs d'après fichier:
+#ifdef ISOTRAC
+      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(xtrain_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,  "XTASNCIEN", "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
-           
-      ! on lit seulement si ixt<=niso ou si on initialise les traceurs d'après
-      ! fichier:
+      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
-      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.AND.k/=0) found=phyetat0_srf(1,iso_tmp_lonsrf,"XTSNOW"//TRIM(oldIso),"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.)
-      if (.NOT.found.AND.k/=0) found=phyetat0_srf(1,iso_tmp_lonsrf,"XTEVAP"//TRIM(oldIso),"evaporation",0.)
-      fxtevap(ixt,:,:)=iso_tmp_lonsrf(:,:)
-
-      found=phyetat0_get(1,iso_tmp,"xtrain_f"//TRIM(outiso),"xrain fall",0.)
-      if (.NOT.found.AND.k/=0) found=phyetat0_get(1,iso_tmp,"xtrain_f"//TRIM(oldIso),"xrain fall",0.)
-      xtrain_fall(ixt,:)=iso_tmp(:)
-
-      found=phyetat0_get(1,iso_tmp,"xtsnow_f"//TRIM(outiso),"snow fall",0.)
-      if (.NOT.found.AND.k/=0) found=phyetat0_get(1,iso_tmp,"xtsnow_f"//TRIM(oldIso),"snow fall",0.)
-      xtsnow_fall(ixt,:)=iso_tmp(:)
-
-      found=phyetat0_get(klev,iso_tmp_lonlev,"XTANCIEN"//TRIM(outiso),"QANCIEN",0.)
-      if (.NOT.found.AND.k/=0) found=phyetat0_get(1,iso_tmp_lonlev,"XTANCIEN"//TRIM(oldIso),"QANCIEN",0.)
-      xt_ancien(ixt,:,:)=iso_tmp_lonlev(:,:)
-
-      found=phyetat0_get(klev,iso_tmp_lonlev,"XTLANCIEN"//TRIM(outiso),"QLANCIEN",0.)
-      if (.NOT.found.AND.k/=0) found=phyetat0_get(1,iso_tmp_lonlev,"XTLANCIEN"//TRIM(oldIso),"QLANCIEN",0.)
-      xtl_ancien(ixt,:,:)=iso_tmp_lonlev(:,:)
-
-      found=phyetat0_get(klev,iso_tmp_lonlev,"XTSANCIEN"//TRIM(outiso),"QSANCIEN",0.)
-      if (.NOT.found.AND.k/=0) found=phyetat0_get(1,iso_tmp_lonlev,"XTSANCIEN"//TRIM(oldIso),"QSANCIEN",0.)
-      xts_ancien(ixt,:,:)=iso_tmp_lonlev(:,:)
-
-      found=phyetat0_get(1,iso_tmp,"XTRUNOFFLIC0"//TRIM(outiso),"RUNOFFLIC0",0.)  
-      if (.NOT.found.AND.k/=0) found=phyetat0_get(1,iso_tmp,"XTRUNOFFLIC0"//TRIM(oldIso),"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.)  
-      if (.NOT.found.AND.k/=0) found=phyetat0_get(1,iso_tmp_lonlev,"WAKE_DELTAXT"//TRIM(oldIso),"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.)  
-        if (.NOT.found.AND.k/=0) found=phyetat0_get(1,iso_tmp,"XTSOL"//TRIM(oldIso),"Surface hmidity / bucket",0.)
-        xtsol(ixt,:)=iso_tmp(:)
-
-        found=phyetat0_get(1,iso_tmp,"Rland_ice"//TRIM(outiso),"R land ice",0.)
-        if (.NOT.found.AND.k/=0) found=phyetat0_get(1,iso_tmp,"Rland_ice"//TRIM(oldIso),"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
+      END IF ! IF(ixt > niso .OR. initialisation_isotrac == 0))
+#endif
+
+   END DO
 
 #ifdef ISOTRAC
-     endif !if ((ixt.le.niso).or.(initialisation_isotrac.eq.0)) then
-#endif
-
-  enddo !do ixt=1,ntraciso
-
-#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:nsrf) = fxtevap(iiso,i,1:nsrf)
+               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       
@@ -16770,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: DZ6/trunk/libf/phylmdiso/phyetat0.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/phyetat0.F90	(revision 4357)
+++ 	(revision )
@@ -1,716 +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_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,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
-
-  ! 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(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_ppm0)
-     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/trunk/libf/phylmdiso/phyetat0_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/phyetat0_mod.F90	(revision 4358)
+++ /LMDZ6/trunk/libf/phylmdiso/phyetat0_mod.F90	(revision 4358)
@@ -0,0 +1,781 @@
+! $Id: phyetat0.F90 3890 2021-05-05 15:15:06Z jyg $
+
+MODULE phyetat0_mod
+
+  PRIVATE
+  PUBLIC :: phyetat0, phyetat0_get, phyetat0_srf
+
+  INTERFACE phyetat0_get
+    MODULE PROCEDURE phyetat0_get10, phyetat0_get20, phyetat0_get11, phyetat0_get21
+  END INTERFACE phyetat0_get
+  INTERFACE phyetat0_srf
+    MODULE PROCEDURE phyetat0_srf20, phyetat0_srf30, phyetat0_srf21, phyetat0_srf31
+  END INTERFACE phyetat0_srf
+
+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 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, 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, types_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(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(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.)
+     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 
+
+#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
+
+!==============================================================================
+LOGICAL FUNCTION phyetat0_get10(field, name, descr, default) RESULT(lFound)
+! Read a field. Check whether reading succeded and use default value if not.
+  IMPLICIT NONE
+  REAL,             INTENT(INOUT) :: field(:) ! klon
+  CHARACTER(LEN=*), INTENT(IN)    :: name
+  CHARACTER(LEN=*), INTENT(IN)    :: descr
+  REAL,             INTENT(IN)    :: default
+!------------------------------------------------------------------------------
+  REAL :: fld(SIZE(field),1)
+  lFound = phyetat0_get21(fld, [name], descr, default); field = fld(:,1)
+END FUNCTION phyetat0_get10
+!==============================================================================
+LOGICAL FUNCTION phyetat0_get20(field, name, descr, default) RESULT(lFound)
+! Same as phyetat0_get11, field on multiple levels.
+  IMPLICIT NONE
+  REAL,             INTENT(INOUT) :: field(:,:) ! klon, nlev
+  CHARACTER(LEN=*), INTENT(IN)    :: name
+  CHARACTER(LEN=*), INTENT(IN)    :: descr
+  REAL,             INTENT(IN)    :: default
+!-----------------------------------------------------------------------------
+  lFound = phyetat0_get21(field, [name], descr, default)
+END FUNCTION phyetat0_get20
+!==============================================================================
+LOGICAL FUNCTION phyetat0_get11(field, name, descr, default) RESULT(lFound)
+! Same as phyetat0_get11, multiple names.
+  IMPLICIT NONE
+  REAL,             INTENT(INOUT) :: field(:) ! klon
+  CHARACTER(LEN=*), INTENT(IN)    :: name(:)
+  CHARACTER(LEN=*), INTENT(IN)    :: descr
+  REAL,             INTENT(IN)    :: default
+!-----------------------------------------------------------------------------
+  REAL :: fld(SIZE(field),1)
+  lFound = phyetat0_get21(fld, name, descr, default); field = fld(:,1)
+END FUNCTION phyetat0_get11
+!==============================================================================
+LOGICAL FUNCTION phyetat0_get21(field, name, descr, default, tname) RESULT(lFound)
+! Same as phyetat0_get11, field on multiple levels, multiple names.
+  USE iostart,           ONLY: get_field
+  USE print_control_mod, ONLY: lunout
+  IMPLICIT NONE
+  REAL,             INTENT(INOUT) :: field(:,:) ! klon, nlev
+  CHARACTER(LEN=*), INTENT(IN)    :: name(:)
+  CHARACTER(LEN=*), INTENT(IN)    :: descr
+  REAL,             INTENT(IN)    :: default
+  CHARACTER(LEN=*), OPTIONAL, INTENT(OUT) :: tname
+!-----------------------------------------------------------------------------
+  CHARACTER(LEN=LEN(name)) :: tnam
+  INTEGER :: i
+  DO i = 1, SIZE(name)
+    CALL get_field(TRIM(name(i)), field, lFound)
+    IF(lFound) EXIT
+    WRITE(lunout,*) "phyetat0: Missing field <",TRIM(name(i)),"> "
+  END DO
+  IF(.NOT.lFound) THEN
+    WRITE(lunout,*) "Slightly distorted start ; continuing."
+    field(:,:) = default
+    tnam = name(1)
+  ELSE
+    tnam = name(i)
+  END IF
+  WRITE(lunout,'(2(a,ES14.7))') 'phyetat0: '//TRIM(tnam)//' ('//TRIM(descr)//') min/max=', &
+    MINval(field),' ',MAXval(field)
+  IF(PRESENT(tname)) tname = tnam
+END FUNCTION phyetat0_get21
+!==============================================================================
+LOGICAL FUNCTION phyetat0_srf20(field, name, descr, default) RESULT(lFound)
+! Read a field per sub-surface.
+! Check whether reading succeded and use default value if not.
+  IMPLICIT NONE
+  REAL,             INTENT(INOUT) :: field(:,:)
+  CHARACTER(LEN=*), INTENT(IN)    :: name
+  CHARACTER(LEN=*), INTENT(IN)    :: descr
+  REAL,             INTENT(IN)    :: default
+!-----------------------------------------------------------------------------
+  REAL :: fld(SIZE(field,1),1,SIZE(field,2))
+  lFound = phyetat0_srf31(fld, [name], descr, default); field = fld(:,1,:)
+END FUNCTION phyetat0_srf20
+
+!==============================================================================
+LOGICAL FUNCTION phyetat0_srf30(field, name, descr, default) RESULT(lFound)
+! Same as phyetat0_sfr11, multiple names tested one after the other.
+  IMPLICIT NONE
+  REAL,             INTENT(INOUT) :: field(:,:,:)
+  CHARACTER(LEN=*), INTENT(IN)    :: name
+  CHARACTER(LEN=*), INTENT(IN)    :: descr
+  REAL,             INTENT(IN)    :: default
+!-----------------------------------------------------------------------------
+  lFound = phyetat0_srf31(field, [name], descr, default)
+END FUNCTION phyetat0_srf30
+
+!==============================================================================
+LOGICAL FUNCTION phyetat0_srf21(field, name, descr, default) RESULT(lFound)
+! Same as phyetat0_sfr11, field on multiple levels.
+  IMPLICIT NONE
+  REAL,             INTENT(INOUT) :: field(:,:)
+  CHARACTER(LEN=*), INTENT(IN)    :: name(:)
+  CHARACTER(LEN=*), INTENT(IN)    :: descr
+  REAL,             INTENT(IN)    :: default
+!-----------------------------------------------------------------------------
+  REAL :: fld(SIZE(field,1),1,SIZE(field,2))
+  lFound = phyetat0_srf31(fld, name, descr, default); field = fld(:,1,:)
+END FUNCTION phyetat0_srf21
+
+!==============================================================================
+LOGICAL FUNCTION phyetat0_srf31(field, name, descr, default) RESULT(lFound)
+! Same as phyetat0_sfr11, field on multiple levels, multiple names tested one after the other.
+  USE iostart,           ONLY: get_field
+  USE print_control_mod, ONLY: lunout
+  USE strings_mod,       ONLY: int2str, maxlen
+  IMPLICIT NONE
+  REAL,             INTENT(INOUT) :: field(:,:,:)
+  CHARACTER(LEN=*), INTENT(IN)    :: name(:)
+  CHARACTER(LEN=*), INTENT(IN)    :: descr
+  REAL,             INTENT(IN)    :: default
+!-----------------------------------------------------------------------------
+  INTEGER :: nsrf, i
+  CHARACTER(LEN=maxlen), ALLOCATABLE :: nam(:)
+  CHARACTER(LEN=maxlen) :: tname, des
+  IF(SIZE(field,3)>99) CALL abort_physic("phyetat0", "Too much sub-cells", 1)
+  DO nsrf = 1, SIZE(field,3)
+    nam = [(TRIM(name(i))//TRIM(int2str(nsrf,2)), i=1, SIZE(name))]
+    des = TRIM(descr)//" srf:"//int2str(nsrf,2)
+    lFound = phyetat0_get21(field(:,:,nsrf), nam, TRIM(des), default, tname)
+  END DO
+  WRITE(lunout,'(2(a,ES14.7))') 'phyetat0: '//TRIM(tname)//' ('//TRIM(descr)//') min/max=', &
+    MINval(field),' ',MAXval(field)
+END FUNCTION phyetat0_srf31
+
+END MODULE phyetat0_mod
+
Index: /LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 4357)
+++ /LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 4358)
@@ -53,4 +53,5 @@
     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, &
@@ -869,5 +870,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
