Index: dynamico_lmdz/simple_physics/config/DYNAMICO/TEST/file_def_physics.xml
===================================================================
--- dynamico_lmdz/simple_physics/config/DYNAMICO/TEST/file_def_physics.xml	(revision 4244)
+++ dynamico_lmdz/simple_physics/config/DYNAMICO/TEST/file_def_physics.xml	(revision 4245)
@@ -10,12 +10,7 @@
       <field field_ref="phyparam_swtop" />
       <field field_ref="phyparam_lwsurf" />
+      <field field_ref="ps"/>
 
-      <field field_ref="phyparam_u" />
-      <field field_ref="phyparam_v" />
       <field field_ref="phyparam_temp" />
-      <field field_ref="phyparam_du" />
-      <field field_ref="phyparam_dv" />
-      <field field_ref="phyparam_dtsw" />
-      <field field_ref="phyparam_dtlw" />
       <field field_ref="phyparam_swflux_down" />
       <field field_ref="phyparam_swflux_up" />
@@ -23,5 +18,14 @@
       <field field_ref="phyparam_lwflux_up" />
 
-    </field_group>
+<!--
+      <field field_ref="phyparam_u" />
+      <field field_ref="phyparam_v" />
+      <field field_ref="phyparam_du" />
+      <field field_ref="phyparam_dv" />
+      <field field_ref="phyparam_dtsw" />
+      <field field_ref="phyparam_dtlw" />
+-->
+
+    </field_group> 
 
     <variable name="model" type="string" > dynamico </variable>
Index: dynamico_lmdz/simple_physics/config/DYNAMICO/TEST/job_CICLAD.sh
===================================================================
--- dynamico_lmdz/simple_physics/config/DYNAMICO/TEST/job_CICLAD.sh	(revision 4244)
+++ dynamico_lmdz/simple_physics/config/DYNAMICO/TEST/job_CICLAD.sh	(revision 4245)
@@ -4,5 +4,5 @@
 #PBS -n
 #PBS -l nodes=1:ppn=60
-#PBS -l walltime=01:00:00
+#PBS -l walltime=04:00:00
 #PBS -l mem=120gb
 #PBS -l vmem=120gb
Index: dynamico_lmdz/simple_physics/config/DYNAMICO/TEST/run.def
===================================================================
--- dynamico_lmdz/simple_physics/config/DYNAMICO/TEST/run.def	(revision 4244)
+++ dynamico_lmdz/simple_physics/config/DYNAMICO/TEST/run.def	(revision 4245)
@@ -42,5 +42,6 @@
 
 #---------------- Run ---------------
-run_length=2592000
+run_length=8640000
+#run_length=2592000
 # 21600s = 6h
 write_period=21600
Index: dynamico_lmdz/simple_physics/phyparam/physics/iniphyparam_mod.F90
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/physics/iniphyparam_mod.F90	(revision 4244)
+++ dynamico_lmdz/simple_physics/phyparam/physics/iniphyparam_mod.F90	(revision 4245)
@@ -71,4 +71,6 @@
 
   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
@@ -76,5 +78,4 @@
 
     CALL read_params(ptimestep)
-
     CALL check_mismatch('unjours', punjours, unjours)
     CALL check_mismatch('rad', prad, planet_rad)
@@ -98,4 +99,6 @@
 
     LOG_INFO('iniphyparam')
+
+    CALL init_soil(nsoilmx)
   END SUBROUTINE iniphyparam
 
Index: dynamico_lmdz/simple_physics/phyparam/physics/phyparam_mod.F90
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/physics/phyparam_mod.F90	(revision 4244)
+++ dynamico_lmdz/simple_physics/phyparam/physics/phyparam_mod.F90	(revision 4245)
@@ -3,5 +3,4 @@
   USE callkeys
   USE comgeomfi
-  USE surface
   IMPLICIT NONE
   PRIVATE
@@ -25,4 +24,8 @@
        &            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 turbulence, ONLY : vdif
     USE convection, ONLY : convadj
@@ -66,19 +69,27 @@
 
     !    Local variables :
-    REAL, DIMENSION(ngrid) :: mu0
-    INTEGER :: j,l,ig,nlevel,igout
-    LOGICAL :: lwrite
-    !
-    REAL :: zday, zdtime
-    REAL zh(ngrid,nlayer),z1,z2
-    REAL zzlev(ngrid,nlayer+1),zzlay(ngrid,nlayer)
-    REAL zdvfr(ngrid,nlayer),zdufr(ngrid,nlayer)
-    REAL zdhfr(ngrid,nlayer),zdtsrf(ngrid),zdtsrfr(ngrid)
-    REAL zflubid(ngrid),zpmer(ngrid)
-    REAL zpopsk(ngrid,nlayer)
+    REAL :: zh(ngrid,nlayer),      & ! potential temperature
+         &  zpopsk(ngrid,nlayer),  & ! Exner function
+         &  zzlev(ngrid,nlayer+1), & ! altitude of interfaces
+         &  zzlay(ngrid,nlayer),   & ! altitude of full levels
+         &  fluxrad(ngrid),        & ! radiative flux at surface
+         &  zc(ngrid, nsoilmx),    & ! LU coefficients for soil implicit solve
+         &  zd(ngrid, nsoilmx),    &
+         &  fluxgrd(ngrid),        & ! heat flux from deep soil
+         &  capcal(ngrid)          & ! effective heat capacity of soil
+         &  zdufr(ngrid,nlayer),   & ! partial tendencies for zonal velocity,
+         &  zdvfr(ngrid,nlayer),   & !   meridional velocity,
+         &  zdhfr(ngrid,nlayer),   & !   potential temperature,
+         &  zdtsrfr(ngrid),        & !   surface temperature
+         &  zdtsrf(ngrid),         & ! total tendency of surface temperature
+         &  zflubid(ngrid),        & ! radiative + deep soil fluxes
+         &  zpmer(ngrid),          & ! sea-level pressure
     REAL zdum1(ngrid,nlayer)
     REAL zdum2(ngrid,nlayer)
     REAL zdum3(ngrid,nlayer)
-    REAL fluxrad(ngrid)
+
+    INTEGER :: j,l,ig,nlevel,igout
+    LOGICAL :: lwrite
+    REAL    :: zday, zdtime, z1, z2
 
     WRITELOG(*,*) 'latitude0', ngrid, lati(1:2), lati(ngrid-1:ngrid)
@@ -88,5 +99,5 @@
     IF (ngrid.NE.ngridmax) THEN
        PRINT*,'STOP in inifis'
-       PRINT*,'Probleme de dimenesions :'
+       PRINT*,'Probleme de dimensions :'
        PRINT*,'ngrid     = ',ngrid
        PRINT*,'ngridmax   = ',ngridmax
@@ -159,4 +170,29 @@
     ENDDO
 
+    !-------------------------------------------------------------
+    !  soil temperatures : 1st half of implicit time integration
+    !  forward sweep from deep ground to surface
+    !  yields LU coefficients zc,zd and capcal, fluxgrd
+    !   ----------------------------------------------------------
+
+    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)
+    ELSE
+       capcal(:)  = capcal_nosoil
+       fluxgrd(:) = 0.
+    ENDIF
+
+    IF(lverbose) THEN
+       WRITELOG(*,*) 'Surface Heat capacity, conduction Flux, Ts'
+       WRITELOG(*,*) capcal(igout), fluxgrd(igout), tsurf(igout)
+       LOG_DBG('phyparam')
+    ENDIF
+
     !-----------------------------------------------------------------------
     !    2. Compute radiative tendencies :
@@ -169,5 +205,8 @@
     !-----------------------------------------------------------------------
     !    3. Vertical diffusion (turbulent mixing):
-    !    -----------------------------------------
+    ! Kz is computed then vertical diffusion is integrated in time implicitly
+    ! using a linear relationship between surface heat flux and air temperature
+    ! in lowest level (Robin-type BC)
+    !    -------------------------------------------------------------------
     !
     IF(calldifv) THEN
@@ -212,4 +251,26 @@
        ENDDO
     ENDIF
+
+    !-------------------------------------------------------------
+    !   soil temperatures : 2nd half of implicit time integration
+    !   using updated tsurf as input
+    !   ----------------------------------------------------------
+
+    DO ig=1,ngrid
+       tsurf(ig)=tsurf(ig)+ptimestep*zdtsrf(ig)
+    ENDDO
+
+    WRITE(55,'(2e15.5)') zday,tsurf(ngrid/2+1)
+
+    IF (callsoil) THEN
+       CALL soil_backward(ngrid,nsoilmx, zc,zd, tsurf,tsoil)
+       IF(lverbose) THEN
+          WRITELOG(*,*) 'Surface Ts, dTs, dt'
+          WRITELOG(*,*) tsurf(igout), zdtsrf(igout), ptimestep
+          LOG_DBG('phyparam')
+       ENDIF
+    END IF
+
+
     !
     !-----------------------------------------------------------------------
@@ -243,29 +304,4 @@
        ENDDO
 
-    ENDIF
-
-    !-----------------------------------------------------------------------
-    !   On ajoute les tendances physiques a la temperature du sol:
-    !   ---------------------------------------------------------------
-
-    DO ig=1,ngrid
-       tsurf(ig)=tsurf(ig)+ptimestep*zdtsrf(ig)
-    ENDDO
-
-    WRITE(55,'(2e15.5)') zday,tsurf(ngrid/2+1)
-
-    !-----------------------------------------------------------------------
-    !   soil temperatures:
-    !   --------------------
-
-    IF (callsoil) THEN
-       CALL soil(ngrid,nsoilmx,.false.,inertie, &
-            &          ptimestep,tsurf,tsoil,capcal,fluxgrd)
-       IF(lverbose) THEN
-          WRITELOG(*,*) 'Surface Heat capacity,conduction Flux, Ts, dTs, dt'
-          WRITELOG(*,*) capcal(igout), fluxgrd(igout), tsurf(igout), &
-               &        zdtsrf(igout), ptimestep
-          LOG_DBG('phyparam')
-       ENDIF
     ENDIF
 
@@ -309,14 +345,16 @@
     !$cython wrapper def alloc(ngrid, nlayer) : phy.phyparam_alloc(ngrid, nlayer)
     USE astronomy, ONLY : iniorbit
+    USE surface
     INTEGER, INTENT(IN), VALUE :: ngrid, nlayer
+    ! allocate precomputed arrays
+    ALLOCATE(rnatur(ngrid), albedo(ngrid), emissiv(ngrid))
+    ALLOCATE(z0(ngrid),inertie(ngrid))
     ! 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), albedo(ngrid), emissiv(ngrid))
-    ALLOCATE(z0(ngrid),inertie(ngrid))
+    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
@@ -326,4 +364,5 @@
     !$cython wrapper def precompute() : phy.phyparam_precompute()
     ! precompute time-independent arrays
+    USE surface
     rnatur(:)  = 1.
     inertie(:) = (1.-rnatur(:))*I_mer+rnatur(:)*I_ter
@@ -337,4 +376,5 @@
     !$cython wrapper def coldstart (ngrid, timestep): phy.phyparam_coldstart(ngrid, timestep)
     ! create internal state to start a run without a restart file
+    USE surface
     INTEGER, INTENT(IN), VALUE :: ngrid
     REAL, INTENT(IN),    VALUE :: ptimestep
@@ -343,7 +383,14 @@
     icount=0
     IF(callsoil) THEN
-       ! initializes zc, zd, capcal, fluxgrd
-       CALL soil(ngrid,nsoilmx,.TRUE.,inertie, &
-            &          ptimestep,tsurf,tsoil,capcal,fluxgrd)
+       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
        WRITELOG(*,*) 'WARNING!!! Thermal conduction in the soil turned off'
Index: dynamico_lmdz/simple_physics/phyparam/physics/radiative_sw.F90
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/physics/radiative_sw.F90	(revision 4244)
+++ dynamico_lmdz/simple_physics/phyparam/physics/radiative_sw.F90	(revision 4245)
@@ -69,5 +69,5 @@
     REAL :: buf1(ngrid), buf2(ngrid, nlayer+1) ! buffers for output
     ! fluxes are non-zero only on those points where the sun shines (mu0>0)
-    ! We compute only on those ncount points and gather them to vectorize 
+    ! We compute only on those ncount points and gather them to vectorize
     INTEGER :: ncount, index(ngrid)
     ! In the work arrays below, ngrid should be ncount but ncount is not known yet
Index: dynamico_lmdz/simple_physics/phyparam/physics/surface.F90
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/physics/surface.F90	(revision 4244)
+++ dynamico_lmdz/simple_physics/phyparam/physics/surface.F90	(revision 4245)
@@ -16,15 +16,18 @@
   REAL :: lambda
   REAL, ALLOCATABLE :: dz1(:),dz2(:)
-  !$OMP THREADPRIVATE(dz1,dz2,zc)
+  !$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 :: zc(:,:),zd(:,:)
-  !$OMP THREADPRIVATE(zc,zd)
-  REAL, ALLOCATABLE :: tsurf(:), tsoil(:,:), capcal(:), fluxgrd(:)
-  !$OMP THREADPRIVATE(tsurf, tsoil, capcal, fluxgrd)
-
-  PUBLIC :: soil, zc, zd, &
+  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
@@ -32,6 +35,6 @@
 CONTAINS
 
-  SUBROUTINE init_soil(ngrid,nsoil)
-    INTEGER, INTENT(IN) :: ngrid, nsoil
+  SUBROUTINE init_soil(nsoil)
+    INTEGER, INTENT(IN) :: nsoil
     REAL :: min_period,dalph_soil, rk,fz1,rk1,rk2
     INTEGER :: jk
@@ -42,5 +45,5 @@
     !   --------------------------------------------------------
 
-    WRITELOG(*,*) 'nsoil,ngrid,firstcall=',nsoil,ngrid, .TRUE.
+    WRITELOG(*,*) 'nsoil,firstcall=',nsoil, .TRUE.
 
     ALLOCATE(dz1(nsoil),dz2(nsoil))
@@ -82,4 +85,97 @@
 
   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,          &
@@ -149,66 +245,76 @@
 
     IF (firstcall) THEN
-       CALL init_soil(ngrid, nsoil)
+       ! 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 soil temperatures using the Cgrd and Dgrd
-       !  coefficient computed at the previous time-step:
-       !  -----------------------------------------------
-
-       !    surface temperature
+       !   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
-          ptsoil(ig,1)=(lambda*zc(ig,1)+ptsrf(ig))/                   &
-               &      (lambda*(1.-zd(ig,1))+1.)
-       ENDDO
-
-       !   other temperatures
-       DO jk=1,nsoil-1
+          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
-             ptsoil(ig,jk+1)=zc(ig,jk)+zd(ig,jk)*ptsoil(ig,jk)
+             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
-
-    ENDIF
-
-    !-----------------------------------------------------------------------
-    !   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
+       
+       !-----------------------------------------------------------------------
+       !   computation of the surface diffusive flux from ground and
+       !   calorific capacity of the ground:
+       !   ---------------------------------
+       
        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
-
+          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
+     
+   END MODULE surface
+   
