Index: dynamico_lmdz/simple_physics/phyparam/physics/surface.F90
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/physics/surface.F90	(revision 4242)
+++ 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
+   
