Index: dynamico_lmdz/simple_physics/phyparam/Makefile
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/Makefile	(revision 4227)
+++ dynamico_lmdz/simple_physics/phyparam/Makefile	(revision 4228)
@@ -30,10 +30,10 @@
 obj/convection.o       : logging.o
 obj/surface.o          : logging.o
-obj/astronomy.o        : logging.o planet.o
+obj/astronomy.o        : logging.o
 obj/solar.o            : logging.o astronomy.o
 obj/radiative_lw.o     : logging.o phys_const.o
 obj/radiative_sw.o     : logging.o phys_const.o
-obj/turbulence.o       : logging.o phys_const.o planet.o
-MAIN                   : callkeys.o phys_const.o astronomy.o turbulence.o surface.o convection.o
+obj/turbulence.o       : logging.o phys_const.o
+MAIN                   : callkeys.o phys_const.o astronomy.o turbulence.o surface.o convection.o planet.o
 obj/iniphyparam_mod.o  : MAIN read_param_mod.o
 obj/phyparam_mod.o     : MAIN solar.o radiative_sw.o radiative_lw.o comgeomfi.o writefield_mod.o
Index: dynamico_lmdz/simple_physics/phyparam/dynphy_lonlat/physiq_mod.F90
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/dynphy_lonlat/physiq_mod.F90	(revision 4227)
+++ dynamico_lmdz/simple_physics/phyparam/dynphy_lonlat/physiq_mod.F90	(revision 4228)
@@ -82,4 +82,6 @@
          d_u,d_v,d_t,d_ps) 
 
+    d_qx(:,:,:)=0.
+
     IF(lafin) THEN
        call iotd_fin
Index: dynamico_lmdz/simple_physics/phyparam/physics/astronomy.F90
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/physics/astronomy.F90	(revision 4227)
+++ dynamico_lmdz/simple_physics/phyparam/physics/astronomy.F90	(revision 4228)
@@ -14,5 +14,4 @@
  
   SUBROUTINE solarlong(pday,psollong)
-    USE planet
     REAL, INTENT(IN) :: pday               ! jour de l'annee (le jour 0 correspondant a l'equinoxe)
     REAL, INTENT(OUT) :: psollong          ! solar longitude
Index: dynamico_lmdz/simple_physics/phyparam/physics/iniphyparam_mod.F90
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/physics/iniphyparam_mod.F90	(revision 4227)
+++ dynamico_lmdz/simple_physics/phyparam/physics/iniphyparam_mod.F90	(revision 4228)
@@ -1,3 +1,4 @@
 MODULE iniphyparam_mod
+#include "use_logging.h"
   IMPLICIT NONE
   PRIVATE
@@ -27,7 +28,4 @@
     REAL, INTENT(IN)  :: ptimestep, prad, pg, pr, pcpp, punjours
     
-    print*,'INIPHYPARAM'
-    print*,'Avant les getpar '
-
     CALL read_param('unjours', 86400., unjours,'unjours')
     CALL read_param('planet_rad',prad,planet_rad,'planet_rad')
@@ -72,12 +70,12 @@
     CALL check_mismatch('cpp', pcpp, cpp)
     
-    print*,'Activation de la physique:'
-    print*,' R=',r
-    print*,' Rayonnement ',callrad
-    print*,' Diffusion verticale turbulente ', calldifv
-    print*,' Ajustement convectif ',calladj
-    print*,' Sol ',callsoil
-    print*,' Cycle diurne ',diurnal
-    
+    WRITELOG(*,*) 'Activation de la physique:'
+    WRITELOG(*,*) ' R=',r
+    WRITELOG(*,*) ' Rayonnement ',callrad
+    WRITELOG(*,*) ' Diffusion verticale turbulente ', calldifv
+    WRITELOG(*,*) ' Ajustement convectif ',calladj
+    WRITELOG(*,*) ' Sol ',callsoil
+    WRITELOG(*,*) ' Cycle diurne ',diurnal
+
     !   choice of the frequency of the computation of radiations
     IF(diurnal) THEN
@@ -87,6 +85,6 @@
     ENDIF
     iradia=1
-    PRINT*,'unjours',punjours
-    PRINT*,'The radiative transfer is computed each ', &
+    WRITELOG(*,*) 'unjours',punjours
+    WRITELOG(*,*) 'The radiative transfer is computed each ', &
          &   iradia,' physical time-step or each ', &
          &   iradia*ptimestep,' seconds'
@@ -94,4 +92,5 @@
     dtphys=ptimestep
 
+    LOG_INFO('iniphyparama')
   END SUBROUTINE iniphyparam
 
Index: dynamico_lmdz/simple_physics/phyparam/physics/phyparam_mod.F90
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/physics/phyparam_mod.F90	(revision 4227)
+++ dynamico_lmdz/simple_physics/phyparam/physics/phyparam_mod.F90	(revision 4228)
@@ -1,3 +1,4 @@
 MODULE phyparam_mod
+#include "use_logging.h"
   USE callkeys
   USE comgeomfi
@@ -6,17 +7,17 @@
   SAVE
 
+  REAL :: nan ! initialized to NaN with adequate compiler options
+
   REAL, PARAMETER :: pi=2*ASIN(1.), solarcst=1370., stephan=5.67e-08, &
        ps_rad=1.e5, height_scale=10000., ref_temp=285., &
-       capcal_nosoil=1e5
-    
-  REAL, ALLOCATABLE :: tsurf(:),tsoil(:,:),rnatur(:), &
-       capcal(:),fluxgrd(:), &
-       dtrad(:,:),fluxrad(:), &
-       q2(:,:),q2l(:,:), &
-       albedo(:),emissiv(:),z0(:),inertie(:)
-  !$OMP THREADPRIVATE( tsurf,tsoil,rnatur)
-  !$OMP THREADPRIVATE( capcal,fluxgrd,dtrad,fluxrad)
-  !$OMP THREADPRIVATE( q2,q2l)
-  !$OMP THREADPRIVATE( albedo,emissiv,z0,inertie)          
+       capcal_nosoil=1e5, tsoil_init=300.
+    
+  ! internal state, written to / read from disk at checkpoint / restart
+  REAL, ALLOCATABLE :: tsurf(:), tsoil(:,:), capcal(:), fluxgrd(:)
+  !$OMP THREADPRIVATE(tsurf, tsoil, capcal, fluxgrd)
+
+  ! precomputed arrays
+  REAL, ALLOCATABLE :: rnatur(:), albedo(:),emissiv(:), z0(:), inertie(:)
+  !$OMP THREADPRIVATE( rnatur, albedo, emissiv, z0, inertie)          
 
   INTEGER :: icount
@@ -74,5 +75,5 @@
     
     !    Local variables :
-    REAL, DIMENSION(ngrid) :: mu0,fract
+    REAL, DIMENSION(ngrid) :: mu0
     INTEGER :: j,l,ig,nlevel,igout
     !
@@ -88,10 +89,10 @@
     REAL zdum3(ngrid,nlayer)
     REAL zdtlw(ngrid,nlayer),zdtsw(ngrid,nlayer)
-    REAL zfluxsw(ngrid),zfluxlw(ngrid)
-    
-    print*,'OK DANS PHYPARAM'
-    print*,'latitude0',ngrid,lati(1:2),lati(ngrid-1:ngrid)
-    print*,'nlayer',nlayer
-    
+    REAL zfluxsw(ngrid),zfluxlw(ngrid), fluxrad(ngrid)
+    
+    WRITELOG(*,*) 'latitude0', ngrid, lati(1:2), lati(ngrid-1:ngrid)
+    WRITELOG(*,*) 'nlayer',nlayer
+    LOG_DBG('phyparam')
+
     IF (ngrid.NE.ngridmax) THEN
        PRINT*,'STOP in inifis'
@@ -113,18 +114,7 @@
        !         zday_last=rjourvrai
        zday_last=zday-ptimestep/unjours
-       CALL alloc_phyparam(ngrid, nlayer, igout)
-
-       !     print*,'OK PHYPARAM 1 '
-       IF(callsoil) THEN
-          CALL soil(ngrid,nsoilmx,firstcall,inertie, &
-               &          ptimestep,tsurf,tsoil,capcal,fluxgrd)
-          ! NB : this call to soil also performs some calculations, see surface.F90
-       ELSE
-          PRINT*,'WARNING!!! Thermal conduction in the soil turned off'
-          DO ig=1,ngrid
-             capcal(ig)  = capcal_nosoil
-             fluxgrd(ig) = 0.
-          ENDDO
-       ENDIF
+       CALL alloc(ngrid, nlayer)
+       CALL precompute
+       CALL coldstart(ngrid, ptimestep)
     ENDIF
 
@@ -172,10 +162,10 @@
     
     !-----------------------------------------------------------------------
-    !    2. Calcul of the radiative tendencies :
+    !    2. Compute radiative tendencies :
     !    ---------------------------------------
     
     IF(callrad) CALL radiative_tendencies(ngrid, igout, nlayer, &
-         gmtime, ptimestep*float(iradia), zday, pplev, pplay, pt, &
-         pdt, zdtlw, zfluxlw, zdtsw, zfluxsw, mu0)
+            gmtime, ptimestep*float(iradia), zday, pplev, pplay, pt, &
+            pdt, zdtlw, zfluxlw, zdtsw, zfluxsw, fluxrad, mu0)
 
     !-----------------------------------------------------------------------
@@ -184,5 +174,5 @@
     !
     IF(calldifv) THEN
-       
+
        DO ig=1,ngrid
           zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
@@ -198,10 +188,10 @@
        enddo
        
-       CALL vdif(ngrid,nlayer,zday,                 &
-       &        ptimestep,capcal,z0,                &
-       &        pplay,pplev,zzlay,zzlev,            &
-       &        pu,pv,zh,tsurf,emissiv,             &
-       &        zdum1,zdum2,zdum3,zflubid,          &
-       &        zdufr,zdvfr,zdhfr,zdtsrfr,q2,q2l,   &
+       CALL vdif(ngrid,nlayer,zday,        &
+       &        ptimestep,capcal,z0,       &
+       &        pplay,pplev,zzlay,zzlev,   &
+       &        pu,pv,zh,tsurf,emissiv,    &
+       &        zdum1,zdum2,zdum3,zflubid, &
+       &        zdufr,zdvfr,zdhfr,zdtsrfr, &
        &        lverbose)
 
@@ -275,7 +265,8 @@
        &          ptimestep,tsurf,tsoil,capcal,fluxgrd)
        IF(lverbose) THEN
-          PRINT*,'Surface Heat capacity,conduction Flux, Ts, dTs, dt'
-          PRINT*,capcal(igout),fluxgrd(igout),tsurf(igout), &
-               &          zdtsrf(igout),ptimestep
+          WRITELOG(*,*) 'Surface Heat capacity,conduction Flux, Ts, dTs, dt'
+          WRITELOG(*,*) capcal(igout), fluxgrd(igout), tsurf(igout), &
+               &        zdtsrf(igout), ptimestep
+          LOG_DBG('phyparam')
        ENDIF
     ENDIF
@@ -284,8 +275,12 @@
     !   sorties:
     !   --------
-    
-    print*,'zday, zday_last ',zday,zday_last,icount
+
+    WRITELOG(*,*) 'zday, zday_last ',zday,zday_last,icount
+    LOG_DBG('phyparam')
+
     if(abs(zday-zday_last-period_sort)<=ptimestep/unjours/10.) then
-       print*,'zday, zday_last SORTIE ',zday,zday_last
+       WRITELOG(*,*) 'zday, zday_last SORTIE ', zday, zday_last
+       LOG_INFO('phyparam')
+
        zday_last=zday
        !  Ecriture/extension de la coordonnee temps
@@ -314,5 +309,4 @@
        call writefield('mu0','mu0',' ',mu0)
        call writefield('alb','alb',' ',albedo)
-       call writefield('fract','fract',' ',fract)
        call writefield('ps','Surface pressure','Pa',pplev(:,1))
        call writefield('slp','Sea level pressure','Pa',zpmer)
@@ -326,5 +320,5 @@
   SUBROUTINE radiative_tendencies(ngrid, igout, nlayer, &
        gmtime, zdtime, zday, pplev, pplay, pt, &
-       pdt, zdtlw, zfluxlw, zdtsw, zfluxsw, mu0)
+       pdt, zdtlw, zfluxlw, zdtsw, zfluxsw, fluxrad, mu0)
     USE planet
     USE phys_const,   ONLY : planet_rad, unjours
@@ -339,7 +333,12 @@
          &                 pt(ngrid, nlayer+1)
     REAL, INTENT(INOUT) :: pdt(ngrid,nlayer)
-    REAL, INTENT(OUT) :: zdtlw(ngrid,nlayer), zfluxlw(ngrid), &
-         zdtsw(ngrid,nlayer), zfluxsw(ngrid), mu0(ngrid)
-    REAL, DIMENSION(ngrid) :: fract
+    REAL, INTENT(OUT)   :: zdtlw(ngrid,nlayer), & ! long-wave temperature tendency
+         &                 zfluxlw(ngrid),      & ! long-wave flux at surface
+         &                 zdtsw(ngrid,nlayer), & ! short-wave temperature tendency
+         &                 zfluxsw(ngrid),      & ! short-wave flux at surface
+         &                 fluxrad(ngrid),      & ! surface flux
+                           mu0(ngrid)             ! ??
+    ! local variables
+    REAL :: fract(ngrid), dtrad(ngrid, nlayer)
     REAL :: zls, zinsol, tsurf2, ztim1,ztim2,ztim3, dist_sol, declin
     REAL :: zplanck(ngrid)
@@ -362,15 +361,16 @@
        ELSE
           CALL zenang(ngrid,zls,gmtime,zdtime,lati,long,mu0,fract)
-          print*,'ZENANG '
        ENDIF
        
        IF(lverbose) THEN
-          PRINT*,'day, declin, sinlon,coslon,sinlat,coslat'
-          PRINT*,zday, declin,                       &
-               &            sinlon(igout),coslon(igout),  &
-               &            sinlat(igout),coslat(igout)
+          WRITELOG(*,*) 'day, declin, sinlon,coslon,sinlat,coslat'
+          WRITELOG(*,*) zday, declin,                       &
+               &        sinlon(igout),coslon(igout),  &
+               &        sinlat(igout),coslat(igout)
+          LOG_DBG('radiative_tendencies')
        ENDIF
     ELSE
-       print*,'declin,ngrid,planet_rad',declin,ngrid,planet_rad
+       WRITELOG(*,*) 'declin,ngrid,planet_rad', declin, ngrid, planet_rad
+       LOG_DBG('radiative_tendencies')
        CALL mucorr(ngrid,declin,lati,mu0,fract,height_scale,planet_rad)
     ENDIF
@@ -414,63 +414,67 @@
     
     IF(lverbose) THEN
-       PRINT*,'Diagnostics for radiation'
-       PRINT*,'albedo, emissiv, mu0,fract,Frad,Planck'
-       PRINT*,albedo(igout),emissiv(igout),mu0(igout), &
-            &           fract(igout),                       &
-            &           fluxrad(igout),zplanck(igout)
-       PRINT*,'Tlay Play Plev dT/dt SW dT/dt LW (K/day)'
-       PRINT*,'unjours',unjours
+       WRITELOG(*,*) 'Diagnostics for radiation'
+       WRITELOG(*,*) 'albedo, emissiv, mu0,fract,Frad,Planck'
+       WRITELOG(*,*) albedo(igout),emissiv(igout),mu0(igout), &
+            &        fract(igout),                       &
+            &        fluxrad(igout),zplanck(igout)
+       WRITELOG(*,*) 'Tlay Play Plev dT/dt SW dT/dt LW (K/day)'
+       WRITELOG(*,*) 'unjours',unjours
        DO l=1,nlayer
-          WRITE(*,'(3f15.5,2e15.2)') pt(igout,l),   &
+          WRITELOG(*,'(3f15.5,2e15.2)') pt(igout,l),     &
                &         pplay(igout,l),pplev(igout,l),  &
                &         zdtsw(igout,l),zdtlw(igout,l)
        ENDDO
+       LOG_DBG('radiative_tendencies')
     ENDIF
     
   END SUBROUTINE radiative_tendencies
 
-  SUBROUTINE alloc_phyparam(ngrid, nlayer, igout)
+  SUBROUTINE alloc(ngrid, nlayer)
+    USE astronomy, ONLY : iniorbit
+    USE surface, ONLY : zc,zd
+    INTEGER, INTENT(IN) :: ngrid, nlayer
+    LOGICAL, PARAMETER :: firstcall=.TRUE.
+    ! allocate arrays for internal state
+    ALLOCATE(tsurf(ngrid))
+    ALLOCATE(tsoil(ngrid,nsoilmx))
+    ! we could avoid the arrays below with a different implementation of surface / radiation / turbulence coupling
+    ALLOCATE(capcal(ngrid),fluxgrd(ngrid))
+    ALLOCATE(zc(ngrid,nsoilmx),zd(ngrid,nsoilmx))
+    ! allocate precomputed arrays
+    ALLOCATE(rnatur(ngrid))
+    ALLOCATE(albedo(ngrid),emissiv(ngrid),z0(ngrid),inertie(ngrid))
+    CALL iniorbit
+  END SUBROUTINE alloc
+
+  SUBROUTINE precompute
     USE surface
-    USE astronomy, ONLY : iniorbit
-    INTEGER, INTENT(IN) :: ngrid, nlayer, igout
-    LOGICAL, PARAMETER :: firstcall=.TRUE.
-
-    print*,'AKk',ngrid,nsoilmx
-    allocate(tsurf(ngrid))
-    print*,'AKa'
-    allocate (tsoil(ngrid,nsoilmx))
-    print*,'AKb'
-    allocate (rnatur(ngrid))
-    print*,'AK2'
-    allocate(capcal(ngrid),fluxgrd(ngrid))
-    print*,'AK3'
-    allocate(dtrad(ngrid,nlayer),fluxrad(ngrid))
-    print*,'AK4'
-    allocate(q2(ngrid,nlayer+1),q2l(ngrid,nlayer+1))
-    print*,'AK5'
-    allocate(albedo(ngrid),emissiv(ngrid),z0(ngrid),inertie(ngrid))
-    print*,'AK6'
-    
-    
-    PRINT*,'FIRSTCALL  '
-    
-    rnatur=1.
-    emissiv(:)=(1.-rnatur(:))*emi_mer+rnatur(:)*emi_ter
-    inertie(:)=(1.-rnatur(:))*I_mer+rnatur(:)*I_ter
-    albedo(:)=(1.-rnatur(:))*alb_mer+rnatur(:)*alb_ter
-    z0(:)=(1.-rnatur(:))*Cd_mer+rnatur(:)*Cd_ter
-    q2=1.e-10
-    q2l=1.e-10
-    tsurf=300.
-    tsoil=300.
-    print*,tsoil(igout,nsoilmx/2+2)
-    print*,'TS ',tsurf(igout),tsoil(igout,5)
-    CALL iniorbit
-    
-    if (.not.callrad) fluxrad(:)=0.
-    
+    ! precompute time-independent arrays
+    rnatur(:)  = 1.
+    emissiv(:) = (1.-rnatur(:))*emi_mer+rnatur(:)*emi_ter
+    inertie(:) = (1.-rnatur(:))*I_mer+rnatur(:)*I_ter
+    albedo(:)  = (1.-rnatur(:))*alb_mer+rnatur(:)*alb_ter
+    z0(:)      = (1.-rnatur(:))*Cd_mer+rnatur(:)*Cd_ter
+  END SUBROUTINE precompute
+
+  SUBROUTINE coldstart(ngrid, ptimestep)
+    ! create internal state to start a run without a restart file
+    USE surface, ONLY : soil
+    INTEGER, INTENT(IN) :: ngrid
+    REAL, INTENT(iN)    :: ptimestep
+    tsurf(:)   = tsoil_init
+    tsoil(:,:) = tsoil_init
     icount=0
-    
-  END SUBROUTINE alloc_phyparam
+    IF(callsoil) THEN
+       ! initializes zc, zd, capcal, fluxgrd
+       CALL soil(ngrid,nsoilmx,.TRUE.,inertie, &
+            &          ptimestep,tsurf,tsoil,capcal,fluxgrd)
+    ELSE
+       WRITELOG(*,*) 'WARNING!!! Thermal conduction in the soil turned off'
+       LOG_WARN('coldstart')
+       capcal(:)  = capcal_nosoil
+       fluxgrd(:) = 0.
+    ENDIF
+  END SUBROUTINE coldstart
 
 END MODULE phyparam_mod
Index: dynamico_lmdz/simple_physics/phyparam/physics/radiative_sw.F90
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/physics/radiative_sw.F90	(revision 4227)
+++ dynamico_lmdz/simple_physics/phyparam/physics/radiative_sw.F90	(revision 4228)
@@ -53,15 +53,11 @@
     !   ----------
     !
-    INTEGER ngrid,nlayer
-    REAL albedo(ngrid),coefvis
-    REAL pmu(ngrid),pfract(ngrid)
-    REAL plevel(ngrid,nlayer+1),ps_rad
-    REAL psolarf0
-    REAL fsrfvis(ngrid),dtsw(ngrid,nlayer)
-    LOGICAL lwrite,ldiurn
-    !
-    !   variables locales:
-    !   ------------------
-    !
+    INTEGER, INTENT(IN) :: ngrid,nlayer
+    LOGICAL, INTENT(IN) :: lwrite,ldiurn
+    REAL, INTENT(IN)    :: albedo(ngrid), coefvis
+    REAL, INTENT(IN)    :: pmu(ngrid), pfract(ngrid)
+    REAL, INTENT(IN)    :: plevel(ngrid,nlayer+1), ps_rad
+    REAL, INTENT(IN)    :: psolarf0
+    REAL, INTENT(OUT)   :: fsrfvis(ngrid),dtsw(ngrid,nlayer)
     
     REAL zalb(ngrid),zmu(ngrid),zfract(ngrid)
@@ -80,7 +76,6 @@
     !   -------------------
     
-    
     nlevel=nlayer+1
-    
+
     !-----------------------------------------------------------------------
     !   Definitions des tableaux locaux pour les points ensoleilles:
Index: dynamico_lmdz/simple_physics/phyparam/physics/surface.F90
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/physics/surface.F90	(revision 4227)
+++ dynamico_lmdz/simple_physics/phyparam/physics/surface.F90	(revision 4228)
@@ -13,11 +13,14 @@
        &           alb_mer,alb_ter,emi_mer,emi_ter
 
-  !   local saved variables:                                              
-  !   ----------------------                                              
+  ! precomputed variables
   REAL :: lambda 
-  REAL,ALLOCATABLE :: dz1(:),dz2(:),zc(:,:),zd(:,:) 
-  !$OMP THREADPRIVATE(dz1,dz2,zc,zd,lambda)                               
-
-  PUBLIC :: soil
+  REAL,ALLOCATABLE :: dz1(:),dz2(:)
+  !$OMP THREADPRIVATE(dz1,dz2,zc,lambda)
+
+  ! internal state variables needed for checkpoint/restart
+  REAL, ALLOCATABLE :: zc(:,:),zd(:,:) 
+  !$OMP THREADPRIVATE(zc,zd)
+
+  PUBLIC :: soil, zc, zd
 
 CONTAINS
@@ -40,5 +43,4 @@
 
     ALLOCATE(dz1(nsoil),dz2(nsoil)) 
-    ALLOCATE(zc(ngrid,nsoil),zd(ngrid,nsoil)) 
     
     min_period=20000. 
@@ -59,4 +61,5 @@
     ENDDO
     lambda=fz(.5)*dz1(1) 
+
     WRITELOG(*,*) 'full layers, intermediate layers (secoonds)' 
     DO jk=1,nsoil 
@@ -67,6 +70,6 @@
             &        fz(rk)*fz(rk)*pi                                          
     ENDDO
-    
     LOG_INFO('init_soil')
+
   END SUBROUTINE init_soil
   
@@ -134,5 +137,4 @@
     
     INTEGER ig,jk 
-    REAL za(ngrid),zb(ngrid) 
     REAL zdz2(nsoil),z1(ngrid) 
         
@@ -154,5 +156,5 @@
        DO jk=1,nsoil-1 
           DO ig=1,ngrid 
-             ptsoil(ig,jk+1)=zc(ig,jk)+zd(ig,jk)*ptsoil(ig,jk) 
+             ptsoil(ig,jk+1)=zc(ig,jk)+zd(ig,jk)*ptsoil(ig,jk)
           ENDDO
        ENDDO
Index: dynamico_lmdz/simple_physics/phyparam/physics/turbulence.F90
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/physics/turbulence.F90	(revision 4227)
+++ dynamico_lmdz/simple_physics/phyparam/physics/turbulence.F90	(revision 4228)
@@ -101,5 +101,4 @@
        ptimestep,pg,pzlev,pzlay,pz0,pu,pv,ph,pkv,pkh)
     ! FIXME : pkh := pkv
-    USE planet
     INTEGER, INTENT(IN) :: ngrid,nlay
     REAL, INTENT(IN) ::  ptimestep, pg, &
@@ -160,5 +159,5 @@
        pu,pv,ph,ptsrf,pemis, &
        pdufi,pdvfi,pdhfi,pfluxsrf, &
-       pdudif,pdvdif,pdhdif,pdtsrf,pq2,pq2l, &
+       pdudif,pdvdif,pdhdif,pdtsrf, &
        lwrite)
     USE phys_const
@@ -195,5 +194,4 @@
     REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay),pdhdif(ngrid,nlay)
     REAL pdtsrf(ngrid),pcapcal(ngrid),pz0(ngrid)
-    REAL pq2(ngrid,nlay+1),pq2l(ngrid,nlay+1)
     LOGICAL lwrite
     !
