Index: dynamico_lmdz/simple_physics/config/DYNAMICO/TEST/run.def
===================================================================
--- dynamico_lmdz/simple_physics/config/DYNAMICO/TEST/run.def	(revision 4245)
+++ dynamico_lmdz/simple_physics/config/DYNAMICO/TEST/run.def	(revision 4246)
@@ -41,4 +41,7 @@
 diurnal=.TRUE.
 
+planet_rad=6.37122e6
+mugaz=28.97
+
 #---------------- Run ---------------
 run_length=8640000
Index: dynamico_lmdz/simple_physics/config/python/test_python.py
===================================================================
--- dynamico_lmdz/simple_physics/config/python/test_python.py	(revision 4246)
+++ dynamico_lmdz/simple_physics/config/python/test_python.py	(revision 4246)
@@ -0,0 +1,5 @@
+import phyparam as phys 
+
+phys.setup()
+phys.alloc(30,100)
+phys.coldstart(30)
Index: dynamico_lmdz/simple_physics/config/python/test_python.sh
===================================================================
--- dynamico_lmdz/simple_physics/config/python/test_python.sh	(revision 4246)
+++ dynamico_lmdz/simple_physics/config/python/test_python.sh	(revision 4246)
@@ -0,0 +1,14 @@
+#!/bin/env bash
+
+# modules for camelot :
+#     module load python/3.6-anaconda50gnu/7.2.0
+
+cd ../../phyparam
+make py
+cd -
+
+rm -f lib phyparam*.so
+ln -s ../../phyparam/python/phyparam*.so .
+ln -s ../../phyparam/lib .
+
+LD_LIBRARY_PATH=lib python test_python.py
Index: dynamico_lmdz/simple_physics/phyparam/Makefile
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/Makefile	(revision 4245)
+++ dynamico_lmdz/simple_physics/phyparam/Makefile	(revision 4246)
@@ -19,5 +19,4 @@
 	bash/extract_cython.sh
 	cd python ; rm -rf build *.c ; python setup.py build_ext --inplace
-	cd python ; python -c 'import phyparam as phys ; phys.setup(100.); phys.alloc(30,100) ; phys.coldstart(30, 100.)'
 obj :
 	@rm -rf obj include lib xml tmp *~ */*~
Index: dynamico_lmdz/simple_physics/phyparam/physics/iniphyparam_mod.F90
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/physics/iniphyparam_mod.F90	(revision 4245)
+++ dynamico_lmdz/simple_physics/phyparam/physics/iniphyparam_mod.F90	(revision 4246)
@@ -10,7 +10,7 @@
 CONTAINS
 
-  SUBROUTINE read_params(ptimestep) BIND(C, name='phyparam_setup')
-    !$cython header void phyparam_setup(double);
-    !$cython wrapper def setup(timestep) : phy.phyparam_setup(timestep)
+  SUBROUTINE read_params() BIND(C, name='phyparam_setup')
+    !$cython header void phyparam_setup();
+    !$cython wrapper def setup() : phy.phyparam_setup()
     USE read_param_mod
     USE phys_const, ONLY : planet_rad,g,r,cpp,rcp,dtphys,unjours,mugaz
@@ -18,7 +18,6 @@
     USE planet, ONLY : coefir, coefvis
     USE turbulence, ONLY : lmixmin, emin_turb
-    USE surface
+    USE soil_mod
     USE callkeys
-    REAL, INTENT(IN), VALUE :: ptimestep
 
     CALL read_param('planet_rad',6.4e6 ,planet_rad,'planet_rad')
@@ -59,4 +58,14 @@
     CALL read_param('period_sort', 1., period_sort, 'period sorties en jour')
 
+  END SUBROUTINE read_params
+
+  SUBROUTINE iniphyparam(ptimestep, punjours, prad, pg, pr, pcpp)
+    USE comgeomfi, ONLY : nsoilmx
+    USE soil_mod, ONLY : init_soil
+    USE phys_const, ONLY : planet_rad,g,r,cpp,rcp,dtphys,unjours
+    USE callkeys
+    REAL, INTENT(IN)  :: ptimestep, punjours, prad, pg, pr, pcpp
+
+    CALL read_params
     !   choice of the frequency of the computation of radiations
     IF(diurnal) THEN
@@ -66,21 +75,11 @@
     ENDIF
     iradia=1
+    dtphys=ptimestep
 
-    dtphys=ptimestep
-  END SUBROUTINE read_params
-
-  SUBROUTINE iniphyparam(ptimestep, punjours, prad, pg, pr, pcpp)
-    USE comgeomfi, ONLY : nsoilmx
-    USE surface, ONLY : init_soil
-    USE phys_const, ONLY : planet_rad,g,r,cpp,rcp,dtphys,unjours
-    USE callkeys
-    REAL, INTENT(IN)  :: ptimestep, punjours, prad, pg, pr, pcpp
-
-    CALL read_params(ptimestep)
-    CALL check_mismatch('unjours', punjours, unjours)
-    CALL check_mismatch('rad', prad, planet_rad)
-    CALL check_mismatch('g', pg, g)
-    CALL check_mismatch('R', pr, r)
-    CALL check_mismatch('cpp', pcpp, cpp)
+    CALL check_mismatch('day lenght (s)', punjours, unjours)
+    CALL check_mismatch('planetary radius (km)', prad/1000., planet_rad/1000.)
+    CALL check_mismatch('gravity', pg, g)
+    CALL check_mismatch('specific R', pr, r)
+    CALL check_mismatch('specific heat capacity', pcpp, cpp)
     LOG_WARN('iniphyparam')
 
Index: dynamico_lmdz/simple_physics/phyparam/physics/phyparam_mod.F90
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/physics/phyparam_mod.F90	(revision 4245)
+++ dynamico_lmdz/simple_physics/phyparam/physics/phyparam_mod.F90	(revision 4246)
@@ -24,8 +24,7 @@
        &            pdu,pdv,pdt,pdpsrf) BIND(C, name='phyparam_phyparam')
     USE phys_const, ONLY : g, rcp, r, unjours
-    USE surface, ONLY : soil_forward, soil_backward
-    USE surface, ONLY : z0, inertie, emissiv, albedo  ! precomputed
-    USE surface, ONLY : tsurf, tsoil ! state variables
-!    USE surface, ONLY : capcal, fluxgrd, zc, zd ! should be temporaries
+    USE soil_mod, ONLY : soil_forward, soil_backward
+    USE soil_mod, ONLY : z0, inertie, emissiv, albedo  ! precomputed
+    USE soil_mod, ONLY : tsurf, tsoil ! state variables
     USE turbulence, ONLY : vdif
     USE convection, ONLY : convadj
@@ -77,5 +76,5 @@
          &  zd(ngrid, nsoilmx),    &
          &  fluxgrd(ngrid),        & ! heat flux from deep soil
-         &  capcal(ngrid)          & ! effective heat capacity of soil
+         &  capcal(ngrid),         & ! effective heat capacity of soil
          &  zdufr(ngrid,nlayer),   & ! partial tendencies for zonal velocity,
          &  zdvfr(ngrid,nlayer),   & !   meridional velocity,
@@ -84,5 +83,5 @@
          &  zdtsrf(ngrid),         & ! total tendency of surface temperature
          &  zflubid(ngrid),        & ! radiative + deep soil fluxes
-         &  zpmer(ngrid),          & ! sea-level pressure
+         &  zpmer(ngrid)            ! sea-level pressure
     REAL zdum1(ngrid,nlayer)
     REAL zdum2(ngrid,nlayer)
@@ -118,5 +117,5 @@
        CALL alloc(ngrid, nlayer)
        CALL precompute
-       CALL coldstart(ngrid, ptimestep)
+       CALL coldstart(ngrid)
     ENDIF
 
@@ -178,14 +177,14 @@
     IF (callsoil) THEN
        CALL soil_forward(ngrid,nsoilmx, ptimestep, inertie, tsurf, tsoil, &
-         &            zc, zd, capcal, fluxgrd)
-
-!       CALL soil_new(ngrid,nsoilmx,ptimestep,inertie, &
-!            tsurf, tsoil, capcal,fluxgrd)
-!       CALL soil(ngrid,nsoilmx,.false.,inertie, &
-!            &          ptimestep,tsurf,tsoil,capcal,fluxgrd)
+            &            zc, zd, capcal, fluxgrd)
+
+       !       CALL soil_new(ngrid,nsoilmx,ptimestep,inertie, &
+       !            tsurf, tsoil, capcal,fluxgrd)
+       !       CALL soil(ngrid,nsoilmx,.false.,inertie, &
+       !            &          ptimestep,tsurf,tsoil,capcal,fluxgrd)
     ELSE
        capcal(:)  = capcal_nosoil
        fluxgrd(:) = 0.
-    ENDIF
+    END IF
 
     IF(lverbose) THEN
@@ -193,5 +192,5 @@
        WRITELOG(*,*) capcal(igout), fluxgrd(igout), tsurf(igout)
        LOG_DBG('phyparam')
-    ENDIF
+    END IF
 
     !-----------------------------------------------------------------------
@@ -345,5 +344,5 @@
     !$cython wrapper def alloc(ngrid, nlayer) : phy.phyparam_alloc(ngrid, nlayer)
     USE astronomy, ONLY : iniorbit
-    USE surface
+    USE soil_mod
     INTEGER, INTENT(IN), VALUE :: ngrid, nlayer
     ! allocate precomputed arrays
@@ -353,8 +352,4 @@
     ALLOCATE(tsurf(ngrid))
     ALLOCATE(tsoil(ngrid,nsoilmx))
-    IF(.FALSE.) THEN ! arrays below are now local temporaries in phyparam
-       ALLOCATE(capcal(ngrid),fluxgrd(ngrid))
-       ALLOCATE(zc(ngrid,nsoilmx),zd(ngrid,nsoilmx))
-    END IF
     CALL iniorbit
   END SUBROUTINE alloc
@@ -364,5 +359,5 @@
     !$cython wrapper def precompute() : phy.phyparam_precompute()
     ! precompute time-independent arrays
-    USE surface
+    USE soil_mod
     rnatur(:)  = 1.
     inertie(:) = (1.-rnatur(:))*I_mer+rnatur(:)*I_ter
@@ -372,30 +367,16 @@
   END SUBROUTINE precompute
 
-  SUBROUTINE coldstart(ngrid, ptimestep) BIND(C, name='phyparam_coldstart')
-    !$cython header void phyparam_coldstart(int, double);
-    !$cython wrapper def coldstart (ngrid, timestep): phy.phyparam_coldstart(ngrid, timestep)
+  SUBROUTINE coldstart(ngrid) BIND(C, name='phyparam_coldstart')
+    !$cython header void phyparam_coldstart(int);
+    !$cython wrapper def coldstart (ngrid): phy.phyparam_coldstart(ngrid)
     ! create internal state to start a run without a restart file
-    USE surface
+    USE soil_mod, ONLY : tsurf, tsoil
     INTEGER, INTENT(IN), VALUE :: ngrid
-    REAL, INTENT(IN),    VALUE :: ptimestep
     tsurf(:)   = tsoil_init
     tsoil(:,:) = tsoil_init
     icount=0
-    IF(callsoil) THEN
-       IF(.FALSE.) THEN
-          ! init_soil is now called by iniphyparam
-          ! initializes zc, zd, capcal, fluxgrd
-          CALL soil(ngrid,nsoilmx,.TRUE.,inertie, &
-               &          ptimestep,tsurf,tsoil,capcal,fluxgrd)
-       END IF
-       IF(.FALSE.) THEN ! soil_forward is now called by phyparam
-          CALL soil_forward(ngrid, nsoilmx, ptimestep, inertie, tsurf, tsoil, &
-               &            zc, zd, capcal, fluxgrd)
-       END IF
-    ELSE
+    IF(.NOT. callsoil) THEN
        WRITELOG(*,*) 'WARNING!!! Thermal conduction in the soil turned off'
        LOG_WARN('coldstart')
-       capcal(:)  = capcal_nosoil
-       fluxgrd(:) = 0.
     ENDIF
   END SUBROUTINE coldstart
Index: dynamico_lmdz/simple_physics/phyparam/physics/radiative_mod.F90
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/physics/radiative_mod.F90	(revision 4245)
+++ dynamico_lmdz/simple_physics/phyparam/physics/radiative_mod.F90	(revision 4246)
@@ -21,5 +21,5 @@
     USE astronomy,    ONLY : orbite, solarlong
     USE solar,        ONLY : solang, zenang, mucorr
-    USE surface,      ONLY : albedo, emissiv, tsurf
+    USE soil_mod,     ONLY : albedo, emissiv, tsurf
     USE radiative_sw, ONLY : sw
     USE radiative_lw, ONLY : lw
Index: dynamico_lmdz/simple_physics/phyparam/physics/soil_mod.F90
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/physics/soil_mod.F90	(revision 4246)
+++ dynamico_lmdz/simple_physics/phyparam/physics/soil_mod.F90	(revision 4246)
@@ -0,0 +1,190 @@
+MODULE soil_mod
+
+#include "use_logging.h"
+
+  IMPLICIT NONE
+  PRIVATE
+  SAVE
+
+  REAL, PARAMETER :: pi=2.*ASIN(1.)
+
+  ! common variables
+  REAL, PUBLIC ::  I_mer,I_ter,Cd_mer,Cd_ter, &
+       &           alb_mer,alb_ter,emi_mer,emi_ter
+
+  ! precomputed variables
+  REAL :: lambda
+  REAL, ALLOCATABLE :: dz1(:),dz2(:)
+  !$OMP THREADPRIVATE(dz1,dz2)
+  REAL, ALLOCATABLE :: rnatur(:), albedo(:),emissiv(:), z0(:), inertie(:)
+  !$OMP THREADPRIVATE( rnatur, albedo, emissiv, z0, inertie)
+
+  ! internal state, written to / read from disk at checkpoint / restart
+  REAL, ALLOCATABLE :: tsurf(:), tsoil(:,:)
+  !$OMP THREADPRIVATE(tsurf, tsoil)
+
+  PUBLIC :: init_soil, soil_forward, soil_backward, &
+       rnatur, albedo, emissiv, z0, inertie, &
+       tsurf, tsoil
+
+CONTAINS
+
+  SUBROUTINE init_soil(nsoil)
+    INTEGER, INTENT(IN) :: nsoil
+    REAL :: min_period,dalph_soil, rk,fz1,rk1,rk2
+    INTEGER :: jk
+
+    !-----------------------------------------------------------------------
+    !   ground levels
+    !   grnd=z/l where l is the skin depth of the diurnal cycle:
+    !   --------------------------------------------------------
+
+    WRITELOG(*,*) 'nsoil,firstcall=',nsoil, .TRUE.
+
+    ALLOCATE(dz1(nsoil),dz2(nsoil))
+
+    min_period=20000.
+    dalph_soil=2.
+
+    !   la premiere couche represente un dixieme de cycle diurne
+    fz1=sqrt(min_period/pi)
+
+    DO jk=1,nsoil
+       rk1=jk
+       rk2=jk-1
+       dz2(jk)=fz(rk1)-fz(rk2)
+    ENDDO
+    DO jk=1,nsoil-1
+       rk1=jk+.5
+       rk2=jk-.5
+       dz1(jk)=1./(fz(rk1)-fz(rk2))
+    ENDDO
+    lambda=fz(.5)*dz1(1)
+
+    WRITELOG(*,*) 'full layers, intermediate layers (secoonds)'
+    DO jk=1,nsoil
+       rk=jk
+       rk1=jk+.5
+       rk2=jk-.5
+       WRITELOG(*,*) fz(rk1)*fz(rk2)*pi,        &
+            &        fz(rk)*fz(rk)*pi
+    ENDDO
+    LOG_INFO('init_soil')
+
+  CONTAINS
+
+    FUNCTION fz(rk) RESULT(val)
+      REAL :: val, rk
+      val = fz1*(dalph_soil**rk-1.)/(dalph_soil-1.)
+    END FUNCTION fz
+
+  END SUBROUTINE init_soil
+
+  !=======================================================================
+  !
+  !   Auteur:  Frederic Hourdin     30/01/92
+  !   -------
+  !
+  !   objet:  computation of : the soil temperature evolution
+  !   ------                   the surfacic heat capacity "Capcal"
+  !                            the surface conduction flux pcapcal
+  !
+  !
+  !   Method: implicit time integration
+  !   -------
+  !   Consecutive ground temperatures are related by:
+  !           T(k+1) = C(k) + D(k)*T(k)  (1)
+  !   the coefficients C and D are computed at the t-dt time-step.
+  !   Routine structure:
+  !   1)new temperatures are computed  using (1)
+  !   2)C and D coefficients are computed from the new temperature
+  !     profile for the t+dt time-step
+  !   3)the coefficients A and B are computed where the diffusive
+  !     fluxes at the t+dt time-step is given by
+  !            Fdiff = A + B Ts(t+dt)
+  !     or     Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt
+  !            with F0 = A + B (Ts(t))
+  !                 Capcal = B*dt
+  !
+
+  PURE SUBROUTINE soil_backward(ngrid,nsoil, zc,zd, ptsrf,ptsoil)
+    INTEGER, INTENT(IN) :: ngrid, nsoil         ! number of columns, of soil layers
+    REAL, INTENT(IN)    :: zc(ngrid, nsoil), zd(ngrid, nsoil) ! LU factorization
+    REAL, INTENT(IN)    :: ptsrf(ngrid)         ! new surface temperature
+    REAL, INTENT(INOUT) :: ptsoil(ngrid,nsoil)  ! soil temperature
+    INTEGER :: ig, jk
+
+    !-----------------------------------------------------------------------
+    !  Computation of the soil temperatures using the Cgrd and Dgrd
+    !  coefficient computed during the forward sweep
+    !  -----------------------------------------------
+
+    !  surface temperature => temperature in first soil layer
+    DO ig=1,ngrid
+       ptsoil(ig,1)=(lambda*zc(ig,1)+ptsrf(ig))/                   &
+            &      (lambda*(1.-zd(ig,1))+1.)
+    ENDDO
+
+    !   other temperatures
+    DO jk=1,nsoil-1
+       DO ig=1,ngrid
+          ptsoil(ig,jk+1)=zc(ig,jk)+zd(ig,jk)*ptsoil(ig,jk)
+       ENDDO
+    ENDDO
+  END SUBROUTINE Soil_backward
+
+  PURE SUBROUTINE soil_forward(ngrid, nsoil, ptimestep, ptherm_i, ptsrf, ptsoil, &
+       &                       zc, zd, pcapcal, pfluxgrd)
+    INTEGER, INTENT(IN) :: ngrid, nsoil         ! number of columns, of soil layers
+    REAL, INTENT(IN)    :: ptimestep,         & ! time step
+         &                 ptherm_i(ngrid),   & ! thermal inertia ??
+         &                 ptsrf(ngrid),      & ! surface temperature before heat conduction
+         &                 ptsoil(ngrid, nsoil) ! soil temperature before heat conduction
+    REAL, INTENT(OUT)   :: zc(ngrid,nsoil),   &
+         &                 zd(ngrid, nsoil),  & ! LU factorization for backward sweep
+         &                 pcapcal(ngrid),    & ! effective calorific capacity
+         &                 pfluxgrd(ngrid)      ! conductive heat flux at the ground
+    REAL :: z1, zdz2(ngrid)
+    INTEGER :: jk, ig
+
+    !-----------------------------------------------------------------------
+    !   Computation of the Cgrd and Dgrd coefficients the backward sweep :
+    !   ---------------------------------------------------------------
+
+    DO jk=1,nsoil
+       zdz2(jk)=dz2(jk)/ptimestep
+    ENDDO
+
+    DO ig=1,ngrid
+       z1=zdz2(nsoil)+dz1(nsoil-1)
+       zc(ig,nsoil-1)=zdz2(nsoil)*ptsoil(ig,nsoil)/z1
+       zd(ig,nsoil-1)=dz1(nsoil-1)/z1
+    ENDDO
+
+    DO jk=nsoil-1,2,-1
+       DO ig=1,ngrid
+          z1=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk)))
+          zc(ig,jk-1)=                                                &
+               &      (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk))*z1
+          zd(ig,jk-1)=dz1(jk-1)*z1
+       ENDDO
+    ENDDO
+
+    !-----------------------------------------------------------------------
+    !   computation of the surface diffusive flux from ground and
+    !   calorific capacity of the ground:
+    !   ---------------------------------
+
+    DO ig=1,ngrid
+       pfluxgrd(ig)=ptherm_i(ig)*dz1(1)*                              &
+            &   (zc(ig,1)+(zd(ig,1)-1.)*ptsoil(ig,1))
+       z1=lambda*(1.-zd(ig,1))+1.
+       pcapcal(ig)=ptherm_i(ig)*                                      &
+            &   ptimestep*(zdz2(1)+(1.-zd(ig,1))*dz1(1))/z1
+       pfluxgrd(ig)=pfluxgrd(ig)                                      &
+            &   +pcapcal(ig)*(ptsoil(ig,1)*z1-lambda*zc(ig,1)-ptsrf(ig))   &
+            &   /ptimestep
+    ENDDO
+  END SUBROUTINE soil_forward
+
+END MODULE soil_mod
Index: dynamico_lmdz/simple_physics/phyparam/physics/surface.F90
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/physics/surface.F90	(revision 4245)
+++ 	(revision )
@@ -1,320 +1,0 @@
-MODULE surface
-
-#include "use_logging.h"
-
-  IMPLICIT NONE
-  PRIVATE
-  SAVE
-
-  REAL, PARAMETER :: pi=2.*ASIN(1.)
-
-  ! common variables
-  REAL, PUBLIC ::  I_mer,I_ter,Cd_mer,Cd_ter, &
-       &           alb_mer,alb_ter,emi_mer,emi_ter
-
-  ! precomputed variables
-  REAL :: lambda
-  REAL, ALLOCATABLE :: dz1(:),dz2(:)
-  !$OMP THREADPRIVATE(dz1,dz2)
-  REAL, ALLOCATABLE :: rnatur(:), albedo(:),emissiv(:), z0(:), inertie(:)
-  !$OMP THREADPRIVATE( rnatur, albedo, emissiv, z0, inertie)
-
-  ! internal state, written to / read from disk at checkpoint / restart
-  REAL, ALLOCATABLE :: tsurf(:), tsoil(:,:)
-  !$OMP THREADPRIVATE(tsurf, tsoil)
-  ! variables below should be temporary arrays, not persistent
-  REAL, ALLOCATABLE :: zc(:,:),zd(:,:), capcal(:), fluxgrd(:)
-  !$OMP THREADPRIVATE(zc,zd, capcal, fluxgrd)
-
-  PUBLIC :: init_soil, &
-       soil, soil_new, soil_forward, soil_backward, &
-       zc, zd, &
-       rnatur, albedo, emissiv, z0, inertie, &
-       tsurf, tsoil, capcal, fluxgrd
-
-CONTAINS
-
-  SUBROUTINE init_soil(nsoil)
-    INTEGER, INTENT(IN) :: nsoil
-    REAL :: min_period,dalph_soil, rk,fz1,rk1,rk2
-    INTEGER :: jk
-
-    !-----------------------------------------------------------------------
-    !   ground levels
-    !   grnd=z/l where l is the skin depth of the diurnal cycle:
-    !   --------------------------------------------------------
-
-    WRITELOG(*,*) 'nsoil,firstcall=',nsoil, .TRUE.
-
-    ALLOCATE(dz1(nsoil),dz2(nsoil))
-
-    min_period=20000.
-    dalph_soil=2.
-
-    !   la premiere couche represente un dixieme de cycle diurne
-    fz1=sqrt(min_period/pi)
-
-    DO jk=1,nsoil
-       rk1=jk
-       rk2=jk-1
-       dz2(jk)=fz(rk1)-fz(rk2)
-    ENDDO
-    DO jk=1,nsoil-1
-       rk1=jk+.5
-       rk2=jk-.5
-       dz1(jk)=1./(fz(rk1)-fz(rk2))
-    ENDDO
-    lambda=fz(.5)*dz1(1)
-
-    WRITELOG(*,*) 'full layers, intermediate layers (secoonds)'
-    DO jk=1,nsoil
-       rk=jk
-       rk1=jk+.5
-       rk2=jk-.5
-       WRITELOG(*,*) fz(rk1)*fz(rk2)*pi,        &
-            &        fz(rk)*fz(rk)*pi
-    ENDDO
-    LOG_INFO('init_soil')
-
-  CONTAINS
-
-    FUNCTION fz(rk) RESULT(val)
-      REAL :: val, rk
-      val = fz1*(dalph_soil**rk-1.)/(dalph_soil-1.)
-    END FUNCTION fz
-
-  END SUBROUTINE init_soil
-
-  PURE SUBROUTINE soil_backward(ngrid,nsoil, zc,zd, ptsrf,ptsoil)
-    INTEGER, INTENT(IN) :: ngrid, nsoil         ! number of columns, of soil layers
-    REAL, INTENT(IN)    :: zc(ngrid, nsoil), zd(ngrid, nsoil) ! LU factorization
-    REAL, INTENT(IN)    :: ptsrf(ngrid)         ! new surface temperature
-    REAL, INTENT(INOUT) :: ptsoil(ngrid,nsoil)  ! soil temperature
-    INTEGER :: ig, jk
-
-    !-----------------------------------------------------------------------
-    !  Computation of the soil temperatures using the Cgrd and Dgrd
-    !  coefficient computed during the forward sweep
-    !  -----------------------------------------------
-    
-    !  surface temperature => temperature in first soil layer
-    DO ig=1,ngrid
-       ptsoil(ig,1)=(lambda*zc(ig,1)+ptsrf(ig))/                   &
-            &      (lambda*(1.-zd(ig,1))+1.)
-    ENDDO
-    
-    !   other temperatures
-    DO jk=1,nsoil-1
-       DO ig=1,ngrid
-          ptsoil(ig,jk+1)=zc(ig,jk)+zd(ig,jk)*ptsoil(ig,jk)
-       ENDDO
-    ENDDO
-  END SUBROUTINE Soil_backward
-
-  PURE SUBROUTINE soil_forward(ngrid, nsoil, ptimestep, ptherm_i, ptsrf, ptsoil, &
-       &                       zc, zd, pcapcal, pfluxgrd)
-    INTEGER, INTENT(IN) :: ngrid, nsoil         ! number of columns, of soil layers
-    REAL, INTENT(IN)    :: ptimestep,         & ! time step
-         &                 ptherm_i(ngrid),   & ! thermal inertia ??
-         &                 ptsrf(ngrid),      & ! surface temperature before heat conduction
-         &                 ptsoil(ngrid, nsoil) ! soil temperature before heat conduction
-    REAL, INTENT(OUT)   :: zc(ngrid,nsoil),   &
-         &                 zd(ngrid, nsoil),  & ! LU factorization for backward sweep
-         &                 pcapcal(ngrid),    & ! effective calorific capacity
-         &                 pfluxgrd(ngrid)      ! conductive heat flux at the ground
-    REAL :: z1, zdz2(ngrid)
-    INTEGER :: jk, ig
-
-    !-----------------------------------------------------------------------
-    !   Computation of the Cgrd and Dgrd coefficients the backward sweep :
-    !   ---------------------------------------------------------------
-    
-    DO jk=1,nsoil
-       zdz2(jk)=dz2(jk)/ptimestep
-    ENDDO
-    
-    DO ig=1,ngrid
-       z1=zdz2(nsoil)+dz1(nsoil-1)
-       zc(ig,nsoil-1)=zdz2(nsoil)*ptsoil(ig,nsoil)/z1
-       zd(ig,nsoil-1)=dz1(nsoil-1)/z1
-    ENDDO
-    
-    DO jk=nsoil-1,2,-1
-       DO ig=1,ngrid
-          z1=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk)))
-          zc(ig,jk-1)=                                                &
-               &      (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk))*z1
-          zd(ig,jk-1)=dz1(jk-1)*z1
-       ENDDO
-    ENDDO
-    
-    !-----------------------------------------------------------------------
-    !   computation of the surface diffusive flux from ground and
-    !   calorific capacity of the ground:
-    !   ---------------------------------
-    
-    DO ig=1,ngrid
-       pfluxgrd(ig)=ptherm_i(ig)*dz1(1)*                              &
-            &   (zc(ig,1)+(zd(ig,1)-1.)*ptsoil(ig,1))
-       z1=lambda*(1.-zd(ig,1))+1.
-       pcapcal(ig)=ptherm_i(ig)*                                      &
-            &   ptimestep*(zdz2(1)+(1.-zd(ig,1))*dz1(1))/z1
-       pfluxgrd(ig)=pfluxgrd(ig)                                      &
-            &   +pcapcal(ig)*(ptsoil(ig,1)*z1-lambda*zc(ig,1)-ptsrf(ig))   &
-            &   /ptimestep
-    ENDDO
-  END SUBROUTINE soil_forward
-
-  SUBROUTINE soil_new(ngrid,nsoil,ptimestep,ptherm_i, ptsrf,ptsoil, pcapcal,pfluxgrd)
-    INTEGER, INTENT(IN) :: ngrid, nsoil         ! number of columns, of soil layers
-    REAL, INTENT(IN)    :: ptimestep,         & ! time step
-         &                 ptherm_i(ngrid)      ! thermal inertia ??
-    REAL, INTENT(INOUT) :: ptsrf(ngrid),      & ! surface temperature
-         &                 ptsoil(ngrid,nsoil)  ! soil temperature
-    REAL, INTENT(OUT)   :: pcapcal(ngrid),    & ! effective calorific capacity
-         &                 pfluxgrd(ngrid)      ! conductive heat flux at the ground
-    CALL soil_backward(ngrid,nsoil, zc,zd, ptsrf,ptsoil)
-    CALL soil_forward(ngrid, nsoil, ptimestep, ptherm_i, ptsrf, ptsoil, &
-            &            zc, zd, pcapcal, pfluxgrd)
-  END SUBROUTINE soil_new
-
-  SUBROUTINE soil(ngrid,nsoil,firstcall,ptherm_i,          &
-       &          ptimestep,ptsrf,ptsoil,                  &
-       &          pcapcal,pfluxgrd)
-
-    !=======================================================================
-    !
-    !   Auteur:  Frederic Hourdin     30/01/92
-    !   -------
-    !
-    !   objet:  computation of : the soil temperature evolution
-    !   ------                   the surfacic heat capacity "Capcal"
-    !                            the surface conduction flux pcapcal
-    !
-    !
-    !   Method: implicit time integration
-    !   -------
-    !   Consecutive ground temperatures are related by:
-    !           T(k+1) = C(k) + D(k)*T(k)  (1)
-    !   the coefficients C and D are computed at the t-dt time-step.
-    !   Routine structure:
-    !   1)new temperatures are computed  using (1)
-    !   2)C and D coefficients are computed from the new temperature
-    !     profile for the t+dt time-step
-    !   3)the coefficients A and B are computed where the diffusive
-    !     fluxes at the t+dt time-step is given by
-    !            Fdiff = A + B Ts(t+dt)
-    !     or     Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt
-    !            with F0 = A + B (Ts(t))
-    !                 Capcal = B*dt
-    !
-    !   Interface:
-    !   ----------
-    !
-    !   Arguments:
-    !   ----------
-    !   ngrid               number of grid-points
-    !   ptimestep              physical timestep (s)
-    !   pto(ngrid,nsoil)     temperature at time-step t (K)
-    !   ptn(ngrid,nsoil)     temperature at time step t+dt (K)
-    !   pcapcal(ngrid)      specific heat (W*m-2*s*K-1)
-    !   pfluxgrd(ngrid)      surface diffusive flux from ground (Wm-2)
-    !
-    !=======================================================================
-    !   declarations:
-    !   -------------
-
-
-    !-----------------------------------------------------------------------
-    !  arguments
-    !  ---------
-
-    INTEGER ngrid,nsoil
-    REAL ptimestep
-    REAL ptsrf(ngrid),ptsoil(ngrid,nsoil),ptherm_i(ngrid)
-    REAL pcapcal(ngrid),pfluxgrd(ngrid)
-    LOGICAL firstcall
-
-
-    !-----------------------------------------------------------------------
-    !  local arrays
-    !  ------------
-
-    INTEGER ig,jk
-    REAL zdz2(nsoil),z1(ngrid)
-
-    IF (firstcall) THEN
-       ! init_soil is now called by iniphyparam
-       ! CALL init_soil(ngrid, nsoil)
-    ELSE
-       IF(.FALSE.) THEN
-          !-----------------------------------------------------------------------
-          !   Computation of the soil temperatures using the Cgrd and Dgrd
-          !  coefficient computed at the previous time-step:
-          !  -----------------------------------------------
-          
-          !    surface temperature
-          DO ig=1,ngrid
-             ptsoil(ig,1)=(lambda*zc(ig,1)+ptsrf(ig))/                   &
-                  &      (lambda*(1.-zd(ig,1))+1.)
-          ENDDO
-          
-          !   other temperatures
-          DO jk=1,nsoil-1
-             DO ig=1,ngrid
-                ptsoil(ig,jk+1)=zc(ig,jk)+zd(ig,jk)*ptsoil(ig,jk)
-             ENDDO
-          ENDDO
-       ELSE
-          CALL soil_backward(ngrid,nsoil, zc,zd, ptsrf,ptsoil)
-       END IF
-       
-    ENDIF
-
-    IF(.FALSE.) THEN
-       !-----------------------------------------------------------------------
-       !   Computation of the Cgrd and Dgrd coefficient for the next step:
-       !   ---------------------------------------------------------------
-       
-       DO jk=1,nsoil
-          zdz2(jk)=dz2(jk)/ptimestep
-       ENDDO
-       
-       DO ig=1,ngrid
-          z1(ig)=zdz2(nsoil)+dz1(nsoil-1)
-          zc(ig,nsoil-1)=zdz2(nsoil)*ptsoil(ig,nsoil)/z1(ig)
-          zd(ig,nsoil-1)=dz1(nsoil-1)/z1(ig)
-       ENDDO
-       
-       DO jk=nsoil-1,2,-1
-          DO ig=1,ngrid
-             z1(ig)=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk)))
-             zc(ig,jk-1)=                                                &
-                  &      (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk))*z1(ig)
-             zd(ig,jk-1)=dz1(jk-1)*z1(ig)
-          ENDDO
-       ENDDO
-       
-       !-----------------------------------------------------------------------
-       !   computation of the surface diffusive flux from ground and
-       !   calorific capacity of the ground:
-       !   ---------------------------------
-       
-       DO ig=1,ngrid
-          pfluxgrd(ig)=ptherm_i(ig)*dz1(1)*                              &
-               &   (zc(ig,1)+(zd(ig,1)-1.)*ptsoil(ig,1))
-          z1(ig)=lambda*(1.-zd(ig,1))+1.
-          pcapcal(ig)=ptherm_i(ig)*                                      &
-               &   ptimestep*(zdz2(1)+(1.-zd(ig,1))*dz1(1))/z1(ig)
-          pfluxgrd(ig)=pfluxgrd(ig)                                      &
-               &   +pcapcal(ig)*(ptsoil(ig,1)*z1(ig)-lambda*zc(ig,1)-ptsrf(ig))   &
-               &   /ptimestep
-       ENDDO
-    ELSE
-       CALL soil_forward(ngrid, nsoil, ptimestep, ptherm_i, ptsrf, ptsoil, &
-            &            zc, zd, pcapcal, pfluxgrd)
-    END IF
-  END SUBROUTINE soil
-     
-   END MODULE surface
-   
Index: dynamico_lmdz/simple_physics/phyparam/python/setup.py
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/python/setup.py	(revision 4245)
+++ dynamico_lmdz/simple_physics/phyparam/python/setup.py	(revision 4246)
@@ -6,4 +6,4 @@
        ext_modules = cythonize([ Extension("phyparam",["phyparam.pyx"], 
                                            libraries=['phyparam'],
-                                           library_dirs=['.'])] )
+                                           library_dirs=['./lib'])] )
        )
